The Failure Mechanism of Methane Hydrate-Bearing Specimen Based on Energy Analysis Using Discrete Element Method

Abstract

Studying the failure mechanism of methane hydrate specimens (MHSs) is of great significance to the exploitation of methane hydrate. Most previous studies have focused on the macro or micromechanical response of MHS under different conditions. However, there are a few studies that have investigated the mechanical response mechanism of MHS based on energy evolution. Therefore, in this study, a numerical model of the methane hydrate-bearing sediments was constructed in the particle flow code (PFC) environment. Then, the numerical model was validated using the conducted laboratory tests; and a series of numerical tests were conducted under different methane hydrate saturation conditions, and the obtained results were analyzed. These results qualitatively describe the main mechanical properties of the methane hydrate-bearing sediments from the viewpoint of energy evolution. The simulation results indicated that during the shear test, the bond breaks at first. Then, the soil particles (sediments) start to roll and rarely slid before shear strength arrives at the highest value. Around the highest shear strength value, more soil particles begin to roll until they occlude with each other. Strain softening is induced by the combined action of the breakage of the hydrate bond and the slipping of soil particles. The higher the hydrate saturation is, the more obvious the strain softening is. Considering that a good agreement was observed between the numerical simulation results and the laboratory test results, it can be concluded that the numerical simulation approach can complement the existing experimental techniques, and also can further clarify the deformation and failure mechanism of various methane hydrate-bearing sediments. The results obtained from the present study will contribute to a better understanding of the mechanical behavior of the gas hydrate-bearing sediments during hydrate dissociation and gas exploitation.

1. Introduction

Methane hydrate, in which gas molecules occupy the cage of water molecules [1,2,3], is one of the most important new energy resources in the world and is being closely considered by more and more countries. Previous studies have shown that in the presence of hydrate, the mechanical properties of sediment could vary while using different test techniques such as the triaxial shear test, direct shear test, bending test, etc. [4,5,6,7,8]. Nevertheless, the strength of the methane hydrate-bearing sediments will decrease after the dissociation of the hydrate and some catastrophic geohazard, such as a landslide [9,10,11], and wellbore instability [12,13,14] will occur [15]. Moreover, considering that methane is a greenhouse gas and its release to the atmosphere while exploiting gas hydrate reservoirs would intensify global warming, environmental impacts are also a very important issue that should be taken care of appropriately. Hence, it is essential to study the deformation mechanism of the hydrate-bearing sediments to avoid the occurrence of such hazards during the commercial production process in the future.

Many researchers have conducted a series of experiments on the mechanical properties of the methane hydrate-bearing specimen (MHS). Winter et al. [16] conducted acoustic and triaxial shear tests using natural specimens drilled from the Mackenzie Delta and synthetic specimens made in the laboratory. The results showed that the velocity of the compressional wave and shear strength is increased due to the presence of methane hydrate. Masui et al. [17] conducted a series of triaxial tests of synthetic methane hydrate specimens in which methane hydrate was generated in an ice–sand or water–sand mixture and studied shear strength, scant elastic modulus, cohesive force, strain softening, and hydrate saturation (Sh). Miyazaki [18] studied the deformation mechanism of the MHS and found that the volume of the MHS was compressed at first and then expanded.

According to the above discussion, the failure mechanism of the MHS has been studied using laboratory tests at the macro level. The macro level indicates studying the mechanical properties of natural gas hydrate specimens by observing their overall deformation and failure. However, making artificial specimens or obtaining drill samples is complicated because of the nature of the hydrate-bearing sediments. Many factors, such as the dissociation of the hydrate, cannot be avoided in the laboratory, which can influence the accuracy of the experiments. Meanwhile, it is almost impossible to make two identical specimens with the same methane hydrate saturation in the laboratory. Therefore, characterizing the mechanical properties of the MHS is not a straightforward job using laboratory testing. With this in mind, the importance of numerical simulation methods is highlighted.

As a well-known numerical simulation method, the discrete element method (DEM) has been widely applied to many fields of rock and soil engineering. DEM offers a new way to study the desired geotechnical problems. In recent years, several research works have studied the mechanical properties of the MHS using the DEM. Brugada et al. [19] conducted a series of DEM simulations of triaxial compression tests for the pore-filling methane hydrate using PFC3D. The influence of methane hydrate saturation on the stress–strain relationship, the volumetric response, and the friction and dilation angle was studied. Jung et al. [20] simulated the MHS by considering two different kinds of hydrate particle distribution, namely small randomly distributed boned grains and clusters. They examined the relationship between porosity, hydrate saturation, stiffness, strength, and dilative tendency. Jiang et al. [21] proposed a micro bond model to study the mechanical response of bonding-type methane hydrate and conducted a biaxial test using PFC2D. Yu et al. [22] investigated the mechanical response of MHS, considering the particle shape and shear wave velocity using PFC3D. These studies prove the capability of the DEM as well as PFC2D/PFC3D for examining hydrate-bearing sediments.

There are three major microscopic models to describe hydrate distribution in the pore space among the coarse-sized soil grains [23], which are called hydrate morphology and include (i) pore-filling type, where hydrates nucleate on sediment grain boundaries and grow freely into pore spaces without bridging two or more particles together, (ii) load-bearing type, where hydrates bridge neighboring grains and contribute mechanical stability to the granular skeleton by becoming part of the load-bearing framework, and finally (iii) cementation type, in which the soil skeleton is cemented due to the existence of hydrate nucleated at intergranular contacts. The aforementioned studies of Brugada et al. [19] and Jung et al. [20] only considered pore-filling type hydrates. Jiang et al. [21] only considered cementation-type hydrates. However, the distribution of methane hydrate inside the sediments or soil particles is not always in a single type. In the meanwhile, there should be rolling resistance between soil grains. If these factors are considered during simulations, the conducted tests and results will be much closer to the real hydrate-bearing sediments. Moreover, considering that the movement of the soil particles can more clearly describe the failure mechanism of MHS at the meso level (analyzing the mechanical properties of natural gas hydrate-bearing sediments through the relationship between particles), the present work established the discrete element numerical model considering different distributions and saturations of methane hydrate inside the sediment particles. While the existing literature has addressed various aspects of the macro- or micro-mechanical response of MHS under different conditions, a detailed analysis of the mechanical response based on the alteration and evolution of the energy inside the hydrate dissociation system has not been performed yet. If the mechanical response is considered the effect, the energy evolution would be considered the corresponding cause, which highlights the importance of its accurate understanding to correctly interpret the resultant response of the gas hydrate reservoirs during exploitation. While the several techniques for the exploitation of gas hydrates, such as depressurization, thermal stimulation, chemical injection, or gas replacement, differ in their practical principles, the energy evolution is the common concept among all of them, which further highlights the contribution of the present study to the relevant engineering field. Therefore, in the current research, a numerical model of the methane hydrate-bearing sediments was established in the particle flow code (PFC) environment. Then, the particle movements and energy evolution were analyzed and its fundamental patterns and governing rules were obtained. Thus, the presented method and the obtained results will help to understand the failure mechanism of the hydrate-bearing sediments more clearly.

2. Numerical MHS Models

2.1. Simulation Mechanism of Particle Flow Code (PFC)

Cundall and Strack (1979) [24] used the discrete element method to analyze a continuous medium divided into several discrete parts such as spherical or block elements. The properties of the continuous medium can be described using the microscopic model defined in discrete elements. Particle flow code (PFC) is one of the essential software for discrete element method simulations. Because of the simple contact calculation, PFC can run calculations more efficiently than other discrete methods. PFC can simulate the break of blocks composed of bonded particles, and there is no limit to extending the model. The two main elements in PFC2D are the disk-shaped particle and wall, both of which are rigid and undeformable. When a simulation is calculated using PFC2D, the displacement, contact, and interaction among particles are simulated by Newton’s second law and force–displacement law. The force among particles is updated using the force–displacement law. The Newtonian law of motion is used to find the position between particles and boundaries and form new contacts. Figure 1 shows the fundamental principle of this phenomenon. In this work, the numerical MHS was simulated by using disks, the rolling resistance model (RRM), and the parallel bond contact model (PBM) [25]. The RRM (Figure 2a) has rolling and sliding resistance friction to resist the particle rotation or slide. The PBM (Figure 2b) can resist such particle movements, since the PBM acts like a beam that resists the bending moment occurring within the bonded area [26]. The bond stiffness will lose its effect when the bond is broken. This study assembled the PBM between the soil particles and hydrate particles, while the RRM was assembled among the soil particles. The default contact model among the particles was set to the RRM when the PBM was broken. A comprehensive review of different numerical simulation methods including DEM can be found in Bakhshi et al. [27].

2.2. Validation of the Numerical Model

The constructed numerical model was first validated reasonably and then adopted to study the mechanical response of the hydrate-bearing sediments. The validation job was performed using the test results of real hydrate-bearing sediments as the reference and comparing the obtained results of the numerical model with the reference values. The details are as follows.

The microscopic mechanical parameters of the particles and the contact model parameters must be set before running numerical simulation tests using PFC. However, these parameters cannot be directly acquired from laboratory tests. The microscopic parameters of particles and contact models must be selected and verified before the numerical simulations. When setting the parameters, many numerical simulation tests were first conducted under test conditions similar to laboratory tests or in situ tests. Then, the numerical simulation results were compared with the results of experimental tests or field tests. The microscopic mechanical parameters were adjusted repeatedly by the “trial and error” method until they met the required simulation conditions. Herein, by comparing the results of the numerical test and indoor mechanical tests and continuously adjusting the micro parameters, a basic match would be achieved between the numerical simulations and indoor real tests. The “trial and error” method used to check the PFC simulation process (version 5.0) is shown in Figure 3 [28]. The procedure for calibrating parameters by the trial and error method included the following steps:

(a)

Determining the initial value of the parameters, and then adjusting $k_n$ and $k_s$ so that $v$ is close to the test value;

(b)

Adjusting $E_c$ to make $E$ close to the test value;

(c)

Adjusting ${\sigma }_b$ to make ${\sigma }_t$ close to the test value;

(d)

Adjusting and $\mu$ to make $tan\phi$ close to the test value.

PFC3D can simulate the experimental triaxial test of specimens. When the specimen is cylindrical, the confining pressure is a constant value. Thus, the shear tests can be simulated by using PFC2D. For generating the numerical specimen more realistically and closer to the natural specimen, a new method of making numerical MHS is proposed. The grain size distribution of sediments is shown in Figure 4 and coincides with the size distribution of Toyoura sand used in experimental tests [6]. The host sample is consolidated at the confining pressure of 0.5 MPa. The porosity of the host sample is detected, and the number of methane hydrate particles is calculated according to the methane hydrate saturation. Then, the radius of methane hydrate particles is minimized 1.5 times and generated in the pore. The minimized methane hydrate particles are distributed in the pore randomly. They are just like the condensation nodule when methane hydrate generates in natural and laboratory environments. The radius of the methane hydrate particle is magnified to the goal number several times. The process in which the radius of methane hydrate particles enlarges generally can realistically simulate the growing process of methane hydrate in the pore space of the sediments. The last step is the consolidation of the mixture specimen under 1 MPa confining pressure.

These numerical samples were rectangular in two dimensions with a height of 5 mm and a width of 2.5 mm. The material densities used in the simulation correspond to those characterized by the Nankai Trough [29]. The grain size distribution of soil corresponds to that described by Masui et al. [6]. The mechanical parameters for the PFC simulations are listed in

Table 1. Mechanical parameters.

Property	Soil	Methane Hydrate
Density (Kg/m3)	2500	320
Particle sizes, D (mm)	0.01–0.4	0.006
Normal stiffness kn (N/m)	1 × 108	1 × 105
Shear stiffness ks (N/m)	1 × 108	1 × 105
Inter-particle friction μ	0.7	0.75

, and the contact parameters are listed in

Table 2. Contact parameters.

Property	Soil-Hydrate	Soil-Soil	Hydrate-Hydrate
Friction μ	0.15	0.5	0.15
Normal stiffness kn (N/m)	1 × 105	3 × 108	1 × 105
Shear stiffness ks (N/m)	1 × 104	3 × 107	1 × 104
Tension strength (N)	3 × 106		3 × 106
Cohesion (N)	5 × 106		5 × 106
Friction angle	10		10
Rolling resistance coefficient (μr)	0.6

. The smooth lateral cylindrical wall was given a standard stiffness of one-tenth of the mean particle stiffness (Knw = 1 × 10⁴ N/m, and Ksw = 0 N/m, µw = 0) to simulate soft confinement. An example numerical specimen is shown in Figure 5 for which the hydrate saturation is 65%. The total number of soil particles (represented by yellow circles) and hydrate particles (represented by red circles) is 17,562.

2.3. Contact Model

In order to consider the rolling of particles, a rolling resistance model suitable for soil–soil contact and fracture is applied between soil and soil. The specific formulation is as follows [30]:

The force–displacement law for the rolling resistance linear model updates the contact force and moment as:

$$\left\{\begin{array}{c}{\mathit{F}}_c={\mathit{F}}^1+{\mathit{F}}^d\\ M_c={\mathit{M}}^r\end{array}\right.$$

(1)

The magnitude of the updated rolling resistance moment is then checked against a threshold limit:

$${\mathit{M}}^r=\left\{\begin{array}{c} {\mathit{M}}^r, \Vert {\mathit{M}}^r\Vert \le M^{*}\\ M^{*}\left({\mathit{M}}^r/\Vert {\mathit{M}}^r\Vert \right), \mathrm{otherwise}\end{array}\right.$$

(2)

where:

$$\left\{\begin{array}{c}{\mathit{M}}^r={\mathit{M}}^r-k_r\Delta {\mathit{\theta }}_b\\ k^r=k_s{\overline{R}}^2\\ M^{*}={\mu }_r\overline{R}{F}_n^l\end{array}\right.$$

(3)

where ${\mathit{F}}_c$ is the total force, ${\mathit{F}}^1$ is the linear force, ${\mathit{F}}^d$ is the dashpot force, ${\mathit{M}}^r$ is the rolling resistance moment, $\Delta {\mathit{\theta }}_b$ is the relative bend–rotation increment, $M^{*}$ is the limiting torque, ${\mu }_r$ is the rolling friction coefficient, $k_r$ is the rolling resistance coefficient, and $F_n^l$ is the normal force of the linear part.

In this study, in order to consider the cementation of hydrate, the parallel bond contact model is applied between hydrate–hydrate contact and soil–hydrate contact, and the specific formula is as follows [30]:

The force–displacement law for the linear parallel bond model updates the contact force and moment:

$$\left\{\begin{array}{c}{\mathit{F}}_c={\mathit{F}}^1+{\mathit{F}}^d+\overline{\mathit{F}}\\ {\mathit{M}}_c=\overline{\mathit{M}}\end{array}\right.$$

(4)

where:

$$\left\{\begin{array}{c}\overline{\mathit{F}}={\overline{\mathit{F}}}_n+{\overline{\mathit{F}}}_s\\ \overline{\mathit{M}}={\overline{\mathit{M}}}_b\\ {\overline{\mathit{F}}}_s={\overline{\mathit{F}}}_s-{\overline{k}}_s\overline{A}\Delta {\mathit{\delta }}_s\\ {\overline{\mathit{M}}}_b={\overline{\mathit{M}}}_b-{\overline{k}}_n\overline{I}\Delta {\mathit{\theta }}_b\end{array}\right.$$

(5)

For the 2D model, in the case of cementation, update the maximum normal and shear stresses at the parallel–bond periphery:

$$\overline{\sigma }=\frac{\overline{F_n}}{A}+\frac{\Vert {\overline{\mathit{M}}}_b\Vert \overline{R}}{\overline{I}}\overline{\beta }$$

(6)

$$\overline{\tau }=\frac{\Vert {\overline{\mathit{F}}}_s\Vert }{\overline{A}}$$

(7)

After cementation failure, shear stresses are:

$${\overline{\tau }}_c=\overline{c}-\frac{{\overline{F}}_n}{\overline{A}}tan\overline{\phi }$$

(8)

where ${\mathit{F}}_c$ is the total force, ${\mathit{F}}^1$ is the linear force, ${\mathit{F}}^d$ is the dashpot force, $\overline{\mathit{F}}$ is the parallel–bond force, ${\overline{\mathit{F}}}_n$ is the parallel–bond normal force, ${\overline{\mathit{F}}}_s$ is the parallel–bond tangential force, $\overline{\mathit{M}}$ is the parallel–bond moment, ${\overline{\mathit{M}}}_b$ is the parallel–bond moment bending component, $\Delta {\mathsf{\theta }}_b$ is the relative bend–rotation increment, $\Delta {\mathsf{\delta }}_s$ is the relative shear–displacement increment, $\overline{A}$ is the cross-sectional area, $\overline{I}$ is the moment of inertia of the parallel bond cross section, ${\overline{k}}_n$ is the normal stiffness, and ${\overline{k}}_s$ is the shear stiffness. $\overline{\beta }$ is the moment–contribution factor, $\overline{R}$ is the bond radius, $\overline{\sigma }$ is the normal stress at the bond periphery, and $\overline{\tau }$ is the shear stress at the bond periphery.

3. Numerical Simulation Results and Analysis

The deviator stress–axial strain curves with four different hydrate saturations (Sh = 25, 40, 50, and 65%) obtained from numerical simulations are shown in Figure 6. The shear strength increased with hydrate saturation from 3.13 MPa (Sh = 25%) to 6.32 MPa (Sh = 65%). To evaluate the simulation results, a comparison with laboratory experimental results on real hydrate-bearing sediment samples from the target area was conducted and an acceptable agreement was observed. That is, in the laboratory tests, the shear strength also increased from 3.18 MPa (Sh = 26.4%) to 7.22 MPa (Sh = 67.8%) as shown in Figure 6 [6]. This observation confirms the reliability of the simulation results.

In Figure 6, the deviation between experimental and simulation results is slightly larger at the lowest saturation of gas hydrate. The reason for this phenomenon lies in the fact that when Sh is at a low level, the bonding interaction of hydrate particles between sand grains is insignificant or very weak, and the MHS characteristics are rather close to those of a sand specimen. In this case, the friction force (rolling friction force and sliding friction force) between sand grains mostly controls the MHS shear strength. Because of the low level of void filling, the MHS represents compaction. However, hydrate particles bond together and cement sand grains in the MHS as Sh increases. The increasing bonding force makes it increasingly difficult to damage the MHS. In this case, the peak strengths are mostly controlled by the bonding force of hydrate particles and partly contribute to the friction force. Therefore, a higher Sh results in a larger bonding force among hydrate particles and larger differences between the peak strengths. In addition, the residual strengths are enhanced with increasing Sh. This enhancement may be attributed to the complicated structure of MHS with a high Sh, which has a tight void space and lower tendency to be destroyed because of the dendritic geometry of hydrates growing into the pore spaces.

The simulated Young’s modulus (E50) generally increased with hydrate saturation, as shown in Figure 7. Here, a linear relationship is observed between Young’s modulus and hydrate saturation. A similar linear relationship between Young’s modulus and hydrate saturation was observed for laboratory test results as well. The regular pattern between the two lines shown in Figure 7 is the same, but the values are slightly different, which does not question the reliability of the numerical model.

Figure 8 represents the relationship between peak stress and hydrate saturation for the numerical simulation sample, and the experimental results, respectively. Comparing the two figures, it is undeniable that the linear relationships between peak stress and hydrate saturation in the two figures are almost the same. Although there is a numerical value disparity, the simulated and experimental test results are qualitatively similar.

Contrasting the stress–strain curve of simulations and experimental tests indicate that the mechanical response of the numerical and real sample are very similar in terms of the following aspects: (1) the strain softening was enhanced with an increase in Sh; (2) the elastic modulus increased with increased Sh; (3) the peak shear strength increased with increasing Sh; and (4) the axial strain corresponding to the peak shear strength was about 2–4% both in simulations and experimental tests. Therefore, the DEM simulation presented in this research can quantitatively explain the laboratory experiment.

It can be observed in Figure 9 that the volume strain shows a dilatation phenomenon at the beginning of the simulation. The greater the hydrate saturation is, the more obvious the dilatation phenomenon is. This observation is very similar to the experimental test results [29]. As is known, the damage to the MHS occurs in the form of compression in the axial direction and enlargement in the lateral direction. Therefore, simulation results can reveal the real damage to a well drilled into the hydrate-bearing sediments layer.

The analysis above makes it very clear that the simulated and experimental results have the same qualitative regularity. Axial stress increases rapidly with the enlargement of axial strain at the beginning of the stress–strain curve. Axial stress arrives at the top value of about 3% of axial strain for different hydrate saturations, and then there is strain softening after the maximum value. The larger the hydrate saturation is, the more obvious the strain softening is. The elastic modulus values of both the simulated sample and experimental tests are a linear relationship with hydrate saturation. The relationship between peak stress and hydrate saturation is linear both in PFC simulation and laboratory experiments. The volume strain shows apparent dilation against axial strain at the beginning of the simulation, and the specimen is compressed in the axial direction. Specimens are not destroyed into fragments in various hydrate saturations and show lateral expansion only. All these observations lead to the conclusion that the PFC hydrate-bearing model is reliable for studying the mechanical and geomechanical properties of the hydrate-bearing layers. Hence, the characteristics of the MHS can be researched qualitatively by using the numerical model proposed in this work.

4. Results and Discussions

Failure Pattern of MHS

In this section, the failure pattern of MHS is analyzed from the viewpoint of energy dissipation. The energy dissipation in the process of hydrate sediment failure is irreversible. Therefore, according to the law of conservation of energy

$$E_w=E_a+E_i$$

(9)

$$E_i=E_k+E_c+E_p+E_r+E_{rs}+E_s$$

(10)

where E_w is the input energy, E_a is the energy dissipation by confining pressure, E_i is energy excluded confining pressure dissipation, E_k is the kinetic energy, E_c is the strain energy, E_p is the bond strain energy, E_r is the rolling strain energy, E_s is the energy dissipated by slip, and E_rs is the energy dissipated by rolling slip.

Figure 10 shows the rotation condition of sand particles in different axial-strain stages under various hydrate saturation levels. The rotation-degree distribution of sand particles is plotted at axial strains of 2%, 4%, 6%, 8%, 10%, and 12%. Few sand particles rotate when the axial strain is 2%. At this stage, the sand particles may be mainly constrained by the friction force. As the axial strain increases, the impact force between two sand particles gradually increases. A sand particle will rotate when the friction force exceeds a critical value. When the axial strain is 4%, sand particles rotate randomly. The rotated sand particles are distributed in a band in a certain direction. A negative correlation exists between the number of rotated sand particles and hydrate saturation. The number of rotated sand particles for Sh = 5% is greater than that for Sh = 45% and 70%, which indicates that the constraint force between sand particles increases with the increase in Sh. Furthermore, the distribution of rotated sand particles is similar when the axial strain is 15%, and the MHS has been destroyed. It illustrates that friction force, which includes rolling friction force and sliding friction force, dominates the residual strength of the MHS.

Figure 11 shows the change in the number of cracks with increasing axial strain for different hydrate saturation levels. Under the same hydrate saturation level, as the axial strain increases, the number of cracks gradually increases; the growth trend tends to be flat, and it is predicted that it will not increase after reaching the maximum value. At the same axial strain, the greater the hydrate saturation, the faster the increase in cracks, and the more cracks that are generated. The loss of bonding between particles constitutes a crack. The higher the saturation of natural gas hydrate, the more the bonding between particles is lost; hence, more cracks will be generated. The number of cracks is positively correlated with hydrate saturation level and axial strain. The difference in the number of cracks between Sh = 70% and Sh = 65% is larger than the difference in the number of cracks between other adjacent saturations in the figure. The results show that Sh = 70% has a great influence on the cracks of hydrate specimens, indicating that when Sh = 70%, the hydrate saturation level has a significant influence on the mechanical strength of hydrate-bearing sediment specimens.

Figure 12 shows how much energy is required throughout the simulation process as the axial strain increases for different hydrate saturation levels. Under the condition of constant hydrate saturation, the greater the axial strain, the more the required energy is. This shows that the greater the axial strain, the more strain energy dissipation and friction energy dissipation occur, and the more energy input is required. At the same axial strain, as the hydrate saturation increases, more energy is required. It shows that the higher the hydrate saturation, the stronger the resistance to deformation, and the less likely the sample is to be damaged.

Figure 13 shows the energy-excluded confining pressure dissipation at different hydrate saturations under different axial strains. The results showed that when the hydrate saturation remains unchanged, by the increase in axial strain, the energy excluded confining pressure dissipation also increases, and there is a proportional relationship. This relationship is more obvious when the hydrate saturation increases. When the axial strain is constant, by the increase in hydrate saturation, the energy excluded confining pressure dissipation also increases. When the axial strain increases, the magnitude of the increase in energy also increases.

Figure 14 shows the variation of the total energy dissipated by confining pressure with increasing axial strain under different hydrate saturations. It can be seen from the figure that the hydrate saturation has little effect on the total energy consumed by the confining pressure. However, the total energy dissipated by confining pressure is significantly affected by the axial strain and increases with the increase in the axial strain.

From Figure 13 and Figure 14, the input energy in Figure 12 is just consumed by the confining pressure and various movements inside the specimen. Through comparison, it is found that the energy consumption of confining pressure is not very different under different hydrate saturations. The difference between the energy consumption of Sh = 5% and Sh = 70% is only about 10 KJ. It can be seen that the hydrate saturation mainly affects the energy consumption inside the specimen. The energy dissipation inside the hydrate sample is mainly composed of strain energy and friction energy dissipation between particles, so hydrate saturation can affect the energy dissipation of the hydrate sample.

Figure 15 shows the relationship between kinetic energy and strain under different hydrate saturations. From this figure, it can be seen that there are suspended hydrate particles in the sediment and they have been moving irregularly.

Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20 show the relationship between different kinds of energy consumption and axial strain in samples under different hydrate saturations. From Figure 16, it can be understood that when the axial strain is around 0–4%, the strain energy increases with the increase in the axial strain; when the axial strain is 4–8%, the strain energy begins to decrease with increasing axial strain; and the strain energy gradually tends to a stable value after the axial strain of 8%. The higher the hydrate saturation, the larger the peak of the strain energy.

In Figure 17, the bond strain energy increases with the increase in the axial strain around 0–4%, and it decreases with the increase in the axial strain after 4%. By analyzing this phenomenon, it can be concluded that before the axial strain is 4%, the input energy is mainly converted into strain energy. After the axial strain is 4%, it enters the strain-softening stage. The input energy is mainly converted into friction energy consumption, and the strain energy is gradually reduced. Here, the higher the hydrate saturation, the more bond strain energy is consumed.

Figure 18 shows the total energy dissipated by slip. It can be seen from the figure that the total energy dissipated by slip increases with the increase in the axial strain, and the hydrate saturation has a significant effect on the total energy dissipated by slip and its increasing trend. When Sh = 5%, Sh = 15%, Sh = 20%, Sh = 40%, Sh = 50%, Sh = 55%, Sh = 60%, and Sh = 70%, the total energy dissipated by slip is less and the growth trend is stable. However, when Sh = 10%, Sh = 25%, Sh = 30%, Sh = 35%, Sh = 45%, and Sh = 65%, the total energy dissipated by slip changes abruptly between 4% and 12% of the axial strain rate and consumes a lot of energy. The total energy dissipated by slip is the most when Sh = 25%. The results show that before the axial strain is 4%, it is in the elastic stage and basically does not produce slip friction energy consumption. After the axial strain is 4%, the particles begin to move and begin to produce slip friction energy consumption.

Figure 19 shows the total energy dissipated by the rolling slip. It can be seen from the figure that starting from the axial strain of 4%, the rolling slip begins to consume energy and has an increasing trend. Hydrate saturation has a great influence on the energy consumed by rolling slips, but there is no fixed law of influence.

It can be seen from Figure 20 that the hydrate saturation has little effect on the rolling strain energy. Starting from an axial strain rate of 4%, the rolling strain energy increases with increasing axial strain rate.

From the above-mentioned five figures (Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20), the total energy consumed by slippage is the main part of the energy consumed inside the specimen and is greatly affected by the hydrate saturation.

5. Conclusions

The present work established the discrete element numerical model considering different distribution types of methane hydrate. This model is credible. It can describe the mechanical properties of a real well in the hydrate-bearing sediment layer. The results of the PFC simulations show similar trends to the laboratory measurement results reported by Masui et al. [6], where the shear strength increases with the increase in the hydrate saturation. The relationships of peak shear stress against axial strain and elastic modulus against axial strain are linear, and they are similar to the relationships observed for the real specimens.

By comparing Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20, it can be understood that in the shear test process, when the axial strain is less than 4%, the strain energy is mainly produced, the sample is in the elastic stage, and the slip or rolling between particles is hard. After the axial strain of 4%, the strain energy began to decrease slowly, and the energy consumed by slip and rolling began to increase, indicating that the interaction between particles provided the strength of the sample during the strain-softening stage.

From the analysis of energy consumption, the failure mechanism of the MHS can be obtained. The bond is broken at first, and particles slide or roll hard when axial strain is less than 1%. In this stage, the specimen is under elastic conditions. After the elastic stage, the great mass of soil particles rolls until shear stress arrives at the biggest numerical value. When hydrate saturation is low, the slipping of soil particles dominates the mechanical performance of the MHS. Moreover, strain softening is not evident while strain hardening is performed slightly. When hydrate saturation is high, strain softening is more apparent, and the breakage of hydrate bonds and slipping of soil particles dominate the mechanical performance. The breakage of the bond causes the strain softening of the hydrate specimen.

In this study, the failure mechanism of the MHS was examined using numerical simulation and experimental test results. The proposed method and findings of the present study will help us to understand the complex failure process of methane hydrate-bearing sediments. Additionally, these findings will help us construct more reliable failure models to better control the stability of the seabed while exploiting methane hydrate. After all, the joint application of the numerical simulation approach introduced in the present study with fractal geometry theory, digital rock technology, and deep learning method is suggested as future extensions to the present study [31,32,33,34,35,36,37]. Meanwhile, a good collection of information and data on the basic geomechanical properties of the gas hydrate-bearing sediments as well as the mechanical behavior of gas hydrate layers, as presented in the literature [38,39,40], offer the availability of developing more accurate numerical models based on the fundamental properties of these reservoirs. Herein, for simulations using numerical specimens, the difference in the intensity of particle movement in the specimen, and the energy consumption inside the specimen are of particular importance [41,42].