Numerical Simulation Study on the Mechanics and Pore Characteristics of Tectonically Deformed Coal under Multi-Level and Multi-Cycle Loading and Unloading Conditions

Abstract

Horizontal well cavern completion and stress release is considered a potential technique for efficient development of coalbed methane in tectonically deformed coal (TDC). Pulsating loading and unloading is a key technique for the controlled expansion of caverns and broader stress release within the reservoir. However, current understanding of the mechanical characteristics and pore network structure evolution of TDC under cyclic loading and unloading conditions is still limited. This paper employs numerical simulation methods to study the mechanical behavior and damage characteristics of TDC under cyclic loading and unloading. After obtaining a set of micromechanical parameters reflecting the behavior of TDC samples under triaxial compression in high-stress states, the effects of different stress gradients and cyclic amplitudes on the stress–strain curve, porosity changes, and crack propagation in TDC samples were analyzed. The study results indicate that under various cyclic loading and unloading conditions, the mechanical response characteristics of TDC samples are broadly similar, primarily divided into compression, slow expansion, and accelerated expansion phases. Under low unloading level conditions, the volume expansion of TDC samples is minimal. Also, at the same unloading level, the strain increment decreases with an increasing number of cycles. Correspondingly, under these conditions, the porosity and microcrack expansion in TDC are less than in high-stress gradient scenarios. Under the same unloading level but different amplitudes, the volume expansion rate at 50% unloading amplitude is higher than at 1 MPa unloading amplitude for TDC, with an increased number of crack expansions. Therefore, under cyclic loading conditions, the sensitivity of crack propagation within TDC samples to amplitude is greater than that to unloading level. Under actual pulsating excitation conditions, a low-amplitude, low-stress gradient pulsation method should be used to maintain the stability of horizontal well caverns, and gradually increase the cyclic amplitude to achieve the efficient extraction of coalbed methane in TDC reservoirs. The findings of this study can serve as an important reference for optimizing process parameters in cyclic pulsating stress release engineering for TDC.

1. Introduction

China has abundant reserves of coal and associated coalbed methane, with significant development potential [1,2]. However, many coal seams in China have been affected by one or more phases of tectonic movements during their formation, resulting in a wide distribution of tectonically deformed coal reservoirs [3,4]. Additionally, the gradual depletion of high-quality shallow coal reserves [5] hampers the development of coalbed methane resources in China [6]. Therefore, investigating the efficient surface development of coalbed methane from tectonically deformed coal reservoirs is crucial for coalbed methane extraction, mine gas management, and methane emission reduction.

It is well known that tectonically deformed coal, especially typical varieties, is characterized by low-strength, extremely low permeability, and difficulties in dewatering [5]. Despite extensive development projects in tectonically deformed coal reservoirs across regions such as Anhui and Guizhou, China, and achieving varied production increases, challenges such as unstable production, rapid decline rates, and low extraction efficiency still prevail [7,8]. Statistics indicate that direct fracturing remains the main technique for enhancing production in tectonically deformed coal reservoirs, yet due to their poor mechanical strength, significant achievements are difficult to obtain [9,10]. Numerous scholars have proposed indirect fracturing techniques targeting the roof [11], interburden [12], and inherent structure coal seams [13] in tectonically deformed coal reservoirs, achieving moderate production enhancements. Yuan et al. have established a technique for protective layer mining and depressurization to enhance permeability and gas extraction in tectonically deformed coal reservoirs [14,15]. Following this, Sang et al. have suggested stimulating extensive stress release in these reservoirs through cavity-induced collapse, inspired by cave completion technology [3]. Engineering practices have shown that stress release techniques in tectonically deformed coal reservoirs significantly increase single well production, reaching up to 2000 m³/d [16]. This demonstrates the feasibility and practical value of horizontal well caving completion and pulsating stimulation for inducing stress release.

The technology primarily involves “U-shaped well positioning, large-diameter drilling in horizontal wells, stress release and reservoir stimulation under horizontal well pulsating pumping, high coal dust content fluid lifting, efficient recovery of output, and fluid circulation with pumping”. Specifically, the pulsating pumping segment facilitates the stabilization of cavities in tectonically deformed coal reservoirs and controlled stress release and collapse [17]. Figure 1 presents the schematic diagram of the technology. During the pulsating stimulation, tectonically deformed coal reservoirs endure cumulative cyclic loading, leading to potential coal rock fatigue and deformation, which may result in instability and failure [18]. Additionally, research by Niu et al. suggests that increased moisture content can soften coal rock [19], alter its macroscopic and microscopic surface morphology, and ultimately affect pore distribution [20]. Therefore, investigating the micromechanical damage patterns and the evolution of pore-fracture network structures under cyclic loading and unloading, as well as analyzing porosity expansion under cyclic pulsating loading conditions, is crucial for further application of this technology and the evaluation of its enhanced permeability effects.

In recent years, numerous scholars have conducted research on the physical experiments and numerical simulations of cyclic loading on coal rock. In terms of physical experiments, cyclic loading and unloading typically involve sinusoidal and triangular wave patterns [21]. There are three loading modes: constant initial stress, and progressive pressurization or depressurization. Some scholars have studied the mechanical damage and permeability of samples under progressive cyclic loading and unloading paths through experiments. Regarding mechanical damage, Hou et al. believe that microcrack formation is the primary cause of sample instability during cyclic loading [22], while Yang Yang et al. have shown that irreversible strain increases with continued cyclic unloading [23]. In terms of micro-damage mechanisms, physical experiments often analyze through online permeability tests, acoustic emission signals, and energy evolution. Nasseri et al. utilized acoustic emission technology to gather information on rock damage processes and predict the types of failure [24]. Gao et al. described the mechanism of microcrack formation in tectonically deformed coal from an energy evolution perspective. However, this method does not allow for direct observation of crack evolution and stress distribution [18].

In numerical simulations, the discrete element method (DEM) has been widely used in cyclic loading and unloading experiments [25,26,27]. Xu et al. observed that with continuous cycling, samples exhibited slight expansion, and their mechanical properties were somewhat enhanced after multiple cycles [28]. Lv et al. noted that multistage cyclic stress paths of different waveforms significantly affect the deformation of samples [29]. Cao et al. conducted a comparative study on the destruction and mechanical behavior of transversely isotropic rocks through indoor experiments and DEM simulations, finding a high concordance between DEM results and experimental outcomes [30]. The application of these numerical simulations provides new insights for in-depth studies on the evolution of fractures and changes in mechanical properties in rocks under progressive cyclic loading and unloading.

Despite the diversity in research methods and topics, several issues remain: (1) Few cases achieve the real-time online monitoring of stress–strain and pore-fracture network evolution under cyclic loading and unloading; (2) Most studies focus on the loading rate of cyclic loading, with less consideration for conditions arising from different stress levels and amplitudes; (3) Studies on the effects of cyclic loading on porosity under different operational conditions are rare.

In this study, using typical Huainan tectonically deformed coal reservoir samples, we calibrated micromechanical parameters of particle flow to match macroscopic mechanical properties, established a numerical model of coal samples under multi-level cyclic loading and unloading paths, and analyzed characteristics such as force chains, stress–strain, porosity, and crack evolution. The effects of different cyclic loading conditions on the mechanical properties and pore-fracture network structure of tectonically deformed coal samples were discussed (see Figure 2 for the research process). The experimental results confirm that cyclic loading conditions are beneficial for in situ coalbed methane development in tectonically deformed coal reservoirs, providing valuable insights into the mechanical behavior, damage and deformation patterns, and evolution of pore-fracture structures during horizontal well caving and depressurization extraction processes.

2. Parameter Sensitivity Analysis and Calibration of Micromechanical Parameters

As illustrated in Figure 2, to closely replicate the cyclic loading and unloading physical experiments and ensure the accuracy of discrete element method (DEM) simulation results, it is necessary to calibrate the micromechanical parameters before conducting the cyclic loading and unloading numerical simulation. This calibration aims to align the numerical model more closely with the physical and mechanical properties of the actual samples. To further reduce the workload of parameter calibration and achieve rapid calibration, it is advisable to conduct a parameter sensitivity analysis under high confining pressure conditions. This analysis identifies the impact of each micromechanical parameter on macroscopic mechanical parameters and establishes corresponding mathematical relationships. Based on this, combined with values measured in physical experiments, perform calculations to determine the approximate range of micromechanical parameters. Afterward, make appropriate adjustments to finally calibrate the specific micromechanical parameters, serving as the foundation for subsequent experimental designs.

2.1. Introduction to the Discrete Element Method

Particle Flow Code (PFC), rooted in the discrete element method, decomposes geotechnical materials into aggregates of abstract particle units. By attributing micromechanical parameters to these particles, PFC accurately simulates macroscopic mechanical characteristics of materials, encompassing the spatiotemporal development of cracks and rock mechanics, from a microscale standpoint. PFC is extensively applied in a range of studies. In standard PFC2D simulations, aligning micromechanical parameters with mechanical properties typically requires a trial-and-error approach [31]. This process includes performing indoor triaxial compression tests on mylonite coal samples, accompanied by a sensitivity analysis of parameters to determine the relationship between micromechanical and macroscopic mechanical properties. Fine-tuning these parameters based on triaxial test results subsequently enables the alignment of simulated micromechanical parameters with actual experimental data. The calibrated parameters serve as a foundation for subsequent experiments.

2.2. Parameter Sensitivity Analysis

This section details extensive numerical simulation experiments conducted through a trial-and-error approach. It examines the influence of micromechanical parameters within a discrete element model on macroscopic parameters. The study concentrates the parallel bond contact model, assessing the sensitivity of five micromechanical parameters related to coals deformation parameters (Elastic modulus, Poisson’s ratio) and strength parameters (Peak strength, Cohesion, Internal friction angle). The study involves two key processes. First, it includes fitting significant macro-/micro-parameters. Second, it focuses on creating multivariate function fittings. These functions are specifically for micromechanical parameters that have a strong influence on macroscopic parameters. In the numerical simulations, a confining pressure of 4 MPa was applied, setting seven values for each micromechanical parameter. The parameters for these experiments are detailed in

Table 1. Experimental design for micromechanical sensitivity.

Test Number	EM-01	EM-02	EM-03	EM-04	EM-05	EM-06	EM-07
Emod/GPa	0.2	0.5	2	5	10	20	50
Test Number	KR-01	KR-02	KR-03	KR-04	KR-05	KR-06	KR-07
Kraito	0.5	0.8	1	1.2	1.5	5	10
Test Number	FR-01	FR-02	FR-03	FR-04	FR-05	FR-06	FR-07
Fric	0.01	0.05	0.1	0.25	1	2	5
Test Number	CO-01	CO-02	CO-03	CO-04	CO-05	CO-06	CO-07
PB_Coh/MPa	0.1	0.5	0.8	2.5	5	7	10
Test Number	FA-01	FA-02	FA-03	FA-04	FA-05	FA-06	FA-07
PB_Fa/°	0.1	1	10	20	25	40	50

Notes: Emod represents effective modulus; Kraito represents Normal-to-shear stiffness ratio, PB_Coh represents cohesion; PB_Fa represents friction coefficient.

.

Significant variations were observed in macroscopic mechanical parameters measured under different micromechanical parameter experimental conditions. Due to the inconsistency in the units of these macroscopic mechanical parameters, it is challenging to intuitively demonstrate their mathematical relationships and the impact of micromechanical parameters on them. Therefore, to enable a clearer comparison of the differences and response characteristics of macroscopic mechanical parameters when micromechanical parameters change, we utilized the initialization method for these parameters, essentially normalizing them. The calculation formula is:

$$P_{\mathrm{i}}=\frac{X_{\mathrm{i}}}{X_0}$$

(1)

In the equation, $P_{\mathrm{i}}$ represents the initialization result; $X_{\mathrm{i}}$ represents macro mechanical parameter results; and $X_0$ represents the first non-empty data.

Figure 3 illustrates the response patterns of five macroscopic mechanical parameters to various micromechanical parameters. As indicated by Figure 3, it is evident that macroscopic mechanical parameters are influenced by various micromechanical parameters, each with distinct effects and influenced parameters. Normalized data in Figure 3a show that the stiffness ratio significantly impacts the Poisson’s ratio, following an exponential relationship. Cohesion, internal friction angle, elastic modulus, and peak stress demonstrate characteristics of segmented functions. This variation is particularly noticeable with the stiffness ratio around a value of 1. For instance, cohesion linearly increases with the stiffness ratio up to a value of 1.5. However, the normalized data indicate a relatively minor impact of the stiffness ratio on peak stress, elastic modulus, and internal friction angle, all showing a trend of slight increase before stabilizing.

Figure 3b indicates that although the friction coefficient affects Poisson’s ratio, this effect is relatively minor, with the Poisson’s ratio slightly decreasing as the friction coefficient increases. The friction coefficient has a more pronounced effect on the internal friction angle, which shows an initial increase and then stabilizes at higher values of the friction coefficient. This trend is also mirrored in the peak stress. Additionally, the friction coefficient influences cohesion, demonstrating an exponential growth pattern. There is a logarithmic relationship between the friction coefficient and the elastic modulus. The elastic modulus has a significant impact on peak stress, exhibiting an exponential relationship with high accuracy. It also shows a linear relationship with cohesion but has negligible effects on the internal friction angle and Poisson’s ratio. According to Figure 3c, the internal friction angle has a considerable impact on cohesion, with the elastic modulus and peak stress following in terms of influence, but negligible impact on the Poisson’s ratio. Figure 3d indicates that tangential bond strength mainly influences cohesion and peak stress, with a moderate effect on the internal friction angle, but less so on the elastic modulus and Poisson’s ratio. Figure 3e demonstrates that the parallel bond elastic modulus significantly impacts the elastic modulus, more than peak strength, internal friction angle, cohesion, and Poisson’s ratio, and shows a linear correlation.

Table 2. Single-factor macromolecular parameters fitting results.

Micromechanical Parameters	Macromechanical Parameters	Fitting Formula	Fitting Coefficients (R2)
Effective Modulus	Peak Stress	Y = 19.65X0.1	0.994
	Elastic Modulus	Y = 1.148X − 0.109	0.996
	Internal Friction Angle	Y = −12.06 × e (−X/0.435) + 31.47	0.943
	Force of Cohesion	Y = 2.62 + 0.075X	0.995
Friction Angle	Peak Stress	0.218X + 12.813	0.946
	Elastic Modulus	Y = −1.12 × e (−X/17.31) + 1.57	0.96
	Force of Cohesion	Y = −2.77 × e (−X/12.99) + 3.16	0.996
Friction Coefficient	Peak Stress	Y = −18.41 × e (−X/0.93) + 29.82	0.999
	Internal Friction Angle	Y = −25.53 × e (−X/0.72) + 40.9	0.984
Stiffness Ratio	Poisson Ratio	Y = 0.06X + 0.25	0.979
	Force of Cohesion	Y = 1.08(lnX)/X + 3.22	0.869
Cohesion	Peak Stress	Y = 2.87X + 18.07	0.998
	Force of Cohesion	Y = 1.03X + 3.32	0.941

presents the mathematical relationships between various micromechanical and macroscopic mechanical parameters, elucidating the interrelations depicted in these figures.

Based on the results of the numerical experiments, the analysis of the correlations using multivariate function fitting is as follows:

$$\{\begin{array}{ll}E=0.12E_{\mathrm{c}}-0.88{\mathrm{e}}^{(\frac{-{\phi }_c}{17.31})}+0.3& {\mathrm{R}}^2=0.996\\ \phi =-10.97{\mathrm{e}}^{(\frac{-E_{\mathrm{c}}}{0.44})}-29.61{\mathrm{e}}^{(\frac{-\mu }{0.72})}+44.44& {\mathrm{R}}^2=0.9504\\ c=0.05E_{\mathrm{c}}+0.91\frac{\ln (kn/ks)}{kn/ks}+1.04{\tau }_c-3.43{\mathrm{e}}^{(\frac{-{\phi }_c}{12.99})}+2.32& {\mathrm{R}}^2=0.9556\\ V=0.06\frac{kn}{ks}+0.25& {\mathrm{R}}^2=0.979\\ {\sigma }_c=18.67E_{\mathrm{c}}+0.24{\phi }_c+2.87{\tau }_c-17.66{\mathrm{e}}^{(\frac{-\mu }{0.93})}+0.98& {\mathrm{R}}^2=0.9983\end{array}$$

(2)

In the equation, macromechanical parameters are defined as follows: $E$ represents elastic modulus (GPa); $\phi$ represents internal friction angle (°); $c$ represents force of cohesion (MPa); $V$ represents Poisson ratio; ${\sigma }_c$ represents peak stress (MPa); micromechanical parameters: $E_{\mathrm{c}}$ represents effective modulus (GPa); ${\phi }_{\mathrm{c}}$ represents friction angle (°); ${\tau }_c$ represents cohesion (MPa); $\frac{kn}{ks}$ represents stiffness ratio; and $\mu$ represents friction coefficient.

2.3. Calibration of Microscopic Mechanical Parameters and Their Macroscopic Mechanical Characteristics

To facilitate a deeper investigation into the evolution of pore-fracture networks in tectonically deformed coal samples under cyclic loading and unloading, we utilized the widely recognized Parallel Bond (PB) model [32] in discrete element method simulations. Before performing the biaxial cyclic loading and unloading simulations, calibrating the micromechanical parameters using data from indoor [18] was essential. The numerical simulation initially replicated these triaxial tests using the PB model with discrete element cylindrical samples measuring 50 mm by 100 mm, containing 5354 particles with radii ranging from 0.5 to 0.7 mm. The model’s porosity was set to the measured value of 0.14977, with the confining pressure fixed at 4 MPa to align with the indoor tests. Simulations were conducted under various deviatoric stresses (q values of 4 MPa, 8 MPa, 12 MPa). To reduce simulation time, the duration of the numerical simulation was proportionally scaled down relative to the indoor test time. The micro-parameters of the contact model in the DEM simulation are presented in

Table 3. Basic micromechanical parameters of the model.

Microscopic Mechanical Parameters of PFC Model	Value
Minimum Particle Size, dmin (mm)	0.5
Ratio of Maximum to Minimum Particle Size, dmax/dmin	1.4
Particle Density, kg/m3	1610
Effective modulus, Ec (GPa)	0.36
Normal-to-shear stiffness ratio, kn/ks	0.49
Friction coefficient, μ	0.8
Bond Gap	0.25
$\mathrm{Friction}\mathrm{angle},{\phi }_c$ (°)	30
Tensile strength/Cohesion	0.7
$\mathrm{Cohesion},{\tau }_c$ (MPa)	0.16

. The sample model, as shown in Figure 4, features rigid boundary elements, termed “walls”, on its four sides (top, bottom, left, and right). During the experiment, the PFC servo mechanism is employed to maintain the walls at the left and right ends within a reasonable fluctuation range, ensuring constant confining pressure. The top and bottom walls are set at specific velocities to act as loading platforms, enabling load application through displacement control.

To verify the accuracy of simulation results, it is necessary to compare these with laboratory experiment outcomes. Figure 5 displays the stress–strain curves obtained from laboratory physical experiments (dashed lines) alongside those from PFC^2D numerical simulations of conventional triaxial compression tests (solid lines). The data from physical experiments indicate that the tectonically deformed coal samples exhibit significant elastoplastic strain characteristics, with mechanical properties similar to those of soft rock and relatively low peak strains. The numerical simulation results closely match the actual laboratory experiment outcomes. Under a confining pressure of 4 MPa, the peak strength measured in the laboratory was 9.32 MPa, compared to 9.4 MPa in the numerical simulations, with a maximum absolute error of 0.08 MPa.

Table 4. Comparison of numerical simulation experimental results with indoor experimental results.

Macroscopic Mechanical Properties	Indoor Experiment Results	Simulations Result	Maximum Absolute Error
Peak strength	9.32 MPa	9.40 MPa	0.08 MPa
Elastic modulus	281.67 MPa	287.11 MPa	5.44
Poisson’s ratio	0.249	0.316	0.067
Internal friction angle	28.5°	29.2°	0.7
Cohesion	0.34 MPa	0.37 MPa	0.03

presents the results of the mechanical property parameters from both the laboratory physical experiments and numerical simulations.

A comparative analysis demonstrates that the PFC^2D program can realistically simulate the mechanical properties of tectonically deformed coal samples. Building on this, conducting cyclic loading and unloading simulation experiments with PFC^2D to analyze the evolution of macroscopic mechanical parameters, cracks, and porosity under various conditions is feasible.

2.4. Cyclic Loading and Unloading Experiment Design

The study on horizontal well cavity completion for coalbed methane development, several stress application modes were considered, such as pulsating pressurization, depressurization, and cyclic pulsation. Initially, the research outlined four stress paths for cyclic loading and unloading to assess their impact on fracture expansion. This evaluation was crucial in selecting the most advantageous stress path.

Subsequently, experiments were carried out to determine the most effective stress loading path. The cyclic loading and unloading process was divided into two phases, during which the axial pressure maintained constant. Cyclic loading and unloading commenced from a predefined initial stress value and were conducted at a fixed amplitude. If the impact on fracture expansion was minimal after 5–6 cycles, the initial stress value was adjusted for the next cycle. This procedure was repeated at a consistent 50% amplitude until sample failure. To ensure consistency between numerical simulations and physical experiments, the simulation were designed to closely replicate actual conditions (Table 4). This was achieved by setting the stress differential between unloading levels equal to the cyclic amplitude and maintaining other parameters constant. As presented in

Table 5. Experimental parameters for cycling with unloading under various conditions.

Test Number	Initial Pressure	Pressure Gradient	Cycle Amplitude	Frequency
BB-01	16 MPa	1 MPa	8 MPa	50 Hz
BB-02	16 MPa	2 MPa	8 MPa	50 Hz
DF-01	16 MPa	1 MPa	1 MPa	50 Hz
DF-02	16 MPa	2 MPa	2 MPa	50 Hz

, for Pathway 1, with a constant axial pressure of 16 MPa, the upper limit confining pressure is reduced, and cyclic loading and unloading are conducted according to a set pressure gradient of either 1 MPa or 2 MPa; and for Pathway 2, with the axial pressure held constant at 16 MPa, the cyclic amplitude is altered. Instead of using 50% of the upper limit confining pressure as the amplitude, cycling is conducted under conditions based on a pressure gradient for loading and unloading.

3. Results

3.1. Stress–Strain Curve

Across various cyclic loading and unloading modes, specimens demonstrated consistent deformation and failure patterns, as Figure 6 shows. Specifically, under continuous equal-amplitude cyclic unloading of confining pressure, the deformation process comprises three stages: compression, gradual expansion, and rapid expansion. In the compression stage, increased axial stress and decreased confining pressure caused the gaps between particles to compress gradually, leading to continuous volume reduction. At this stage, axial strain exceeded radial strain. In the second and third stages, as the gaps between particles became further compacted, the mechanical strength of the sample underwent a change, becoming more sensitive to increases in deviatoric stress. Consequently, the sample’s deformation shifted from compressive to expansive.

Under varying conditions, the initial cycle’s strain during cyclic loading at a constant level is notably larger than in later cycles. Significantly, during the first stage, the increase in axial strain is more pronounced than in radial strain, indicating rapid compaction of the sample. However, in the second and third stages, the rate of increase in radial strain exceeds that of axial strain. Additionally, the width of the delay curves for both axial and radial strains gradually widens with each cycle.

Coal, particularly tectonically deformed coal, is not a perfectly elastic material. Consequently, the deformation incurred during loading is partly irreversible, resulting in irreversible plastic strain. This phenomenon is manifested as delay loops in the stress–strain curves, which reveal the accumulation of irreversible plastic strain during cyclic loading and unloading. As the number of cycles increases, samples undergo volumetric expansion and eventually fail.

3.2. Fracture Expansion and Porosity Evolution Characteristics

Cyclic loading clearly caused periodic changes in the samples’ fracture characteristics. Experimental results showed that fractures increased with cycle count, especially in early cycle stages, but stabilized after approximately five cycles. Consequently, only the first five cycles of each stage were analyzed in this paper. Figure 7 illustrates that larger amplitudes led to a notable increase in cracks after the first cycle, due to higher deviatoric stress at minimal confining pressure. Conversely, lower amplitude cycles resulted in fewer fractures due to smaller stress fluctuations, causing compaction instead of fracturing. With more cycles, fracture growth approached saturation, and the specimen’s fracture structure adapted, showing stable growth until reaching dynamic equilibrium at the specific confining pressure level. Additionally, the development of fractures varied noticeably at different confining pressures, correlating with each pressure decrease. At higher confining pressures (e.g., 16 MPa to 15 MPa), fractures were denser yet smaller, while at lower pressures (e.g., 12 MPa to 10 MPa), larger but fewer fractures formed and expanded over cycles. The number of fractures increased by 216% from the first to the seventh stage of cyclic loading and unloading, with the fracture growth in each stage fluctuating around 20% ± 2%. These findings indicate continuous volumetric expansion resulted from gradually reduced stress levels in the second stage. Both the stress cycle path and decompression effectively enhanced reservoir stimulation while maintaining sample stability and preventing collapse, resulting in controlled stress release and improved reservoir stimulation effects.

To quantify porosity evolution in the samples, the built-in ‘measurement circle’ program in PFC^2D was utilized to set 400 measurement points within the specimen. This approach allowed for the collection of porosity data both in its original state and at various stages of the cyclic process. The porosity distribution was then visualized with a cloud diagram. As observed in Figure 8, sample porosity increases from the center to the ends. This trend is likely related to the shear dilation effect observed in granular materials under shear load. In the central part, which is less compacted, particles are looser and show a decrease in pore volume under shear load. At the ends, the compaction effect leads to tighter particle compression. Upon load application, particles displace and roll, increasing pore volume and indicating shear dilation. With further cyclic loading and unloading, the porosity distribution pattern remains consistent across different stages, but an overall increase in porosity is noted with the reduction in confining pressure.

4. Discussion

4.1. The Impact of Unloading Rate Differential on Strain and Porosity

4.1.1. Effects on Strain

Given the uniformity of post-failure trends, this study concentrates on the stress–strain curves during the cyclic phase. Figure 9 presents these curves under cyclic amplitudes of 50% confining pressure, contrasting unload differentials of 1 MPa and 2 MPa. The figure illustrates, despite varying unload differentials, the samples’ axial and radial strain curves follow a similar pattern. Significantly, at a confining pressure of 14 MPa, the sample with a larger unload (BB-02) differential shows more axial strain and less radial expansion than the one with a smaller differential (BB-01). This observation suggests that a larger unload differential under cyclic unloading conditions facilitates more pronounced deformation and volumetric expansion in the coal sample. With advancing cycles and decreasing confining pressure, differences in the stress–strain curves become more pronounced, signifying increased deformation.

Axial strain is more affected by cyclic unloading than radial strain across different unload differentials. At the 16 MPa stage, axial strain increases more than radial strain during cyclic unloading. Moreover, unloading across different pressure levels results in a greater strain increase than observed during the entire cyclic loading and unloading process. Regarding the curve shape, the hysteresis loop of the axial strain curve for specimen BB-02 is slightly larger than that for BB-01, particularly in the final cyclic loading and unloading stage. In this stage, the gap between the two curves gradually widens, eventually surpassing the axial strain of BB-01. A similar pattern occurs in radial strain, showing that as pressure differential increases, the sample’s irreversible strain rises, increasing susceptibility to damage.

4.1.2. Effects on Porosity

To address the limitations of using a 2D model in numerical simulations, and to reduce particle count and computational complexity, it is necessary to convert the 2D porosity into an equivalent 3D porosity. This conversion typically relies on an established relationship derived from structures with equal particle sizes. A commonly used conversion formula relates 2D porosity to 3D porosity, ensuring the 2D model’s mechanical properties and pore structures accurately reflect those of a 3D solid. Using this approach is vital for the reliability and accuracy of results from 2D numerical simulations [33].

$${\epsilon }_{2\mathrm{D}}=1-\frac{1}{\xi }\left(1-{\epsilon }_{3\mathrm{D}}\right)$$

(3)

$$\xi =\frac{\sqrt{2}}{\sqrt{\pi \sqrt{3}}}+{\rho }_{\mathrm{d}}\left(\frac{2}{\sqrt{\pi \sqrt{3}}}-\frac{\sqrt{2}}{\sqrt{\pi \sqrt{3}}}\right)$$

(4)

In the equation, ${\epsilon }_{2\mathrm{D}}$ represents two-dimensional porosity, ${\epsilon }_{3\mathrm{D}}$ represents three-dimensional porosity, ξ is the correction factor, and ${\rho }_{\mathrm{d}}$ stands for relative density.

Figure 10 shows a clear trend of increasing porosity axially and radially as stress is released. Axial porosity evolves more markedly than radial porosity. Even under the same confining pressure, the porosity at different stages does not overlap. Initially, high confining pressure restricts radial strain, while axial strain increases. Decreasing confining pressure causes microfractures to form and expand due to increased deviatoric stress. Consequently, both fracture expansion and axial porosity significantly increase.

Within a single cycle, porosity initially decreases with more cycles, showing a narrowing change range.

Table 6. Strain values under different stress paths.

Pressure Levels	Axial Strain (%)				Radial Strain (%)				Volumetric Strain (%)
	BB-01		BB-02		BB-01		BB-02		BB-01		BB-02
	First	Last	First	Last	First	Last	First	Last	First	Last	First	Last
16	0	0.29	0	0.30	0	−0.03	0	0.11	0	0.23	0	0.52
14	0.42	0.45	0.40	0.45	−0.21	−0.2	−0.22	−0.21	0	0.05	−0.03	0.04
12	0.54	0.72	0.45	0.73	−0.58	−0.59	−0.58	−0.60	−0.51	−0.47	−0.71	−0.46
10	1.17	1.4	0.99	1.54	−1.17	−1.31	−1.08	−1.42	−1.18	−1.22	−1.17	−1.31

indicates that volume strain increases with each cycle, suggesting progressive sample compression during cyclic loading and unloading at the same pressure level. The plastic failure of the rock mainly occurs when the unloading level changes. Meanwhile, the proportion of plastic deformation decreases with cyclic loading at the same level, leading to primarily elastic deformation [34] and, consequently, smaller changes in porosity. Since the numerical model treats particles as rigid bodies, unable to create new fractures or expand internal pores, subsequent cycles cannot produce the same voids or microfractures. Thus, under cyclic loading, only particle compaction occurs, gradually reducing porosity over time.

Comparing porosity evolution under various stress differentials reveals that average radial porosity consistently stays below axial porosity. Initially, there is no significant difference in porosity between the two sample types. However, as unloading progresses, BB-02 consistently exhibits higher porosity than BB-01. The increase in porosity at different pressure levels is more pronounced with a larger stress differential. This observation suggests that greater stress differentials lead to enhanced porosity and improved permeability in the samples.

4.2. The Impact of Different Cyclic Amplitudes on Strain and Porosity

4.2.1. Effects on Strain

Figure 11 illustrates the strain behavior of samples under different amplitudes of cyclic loading and unloading within the same stress differential. Although strain curves differ significantly under various cyclic amplitudes, the final strain values at cycle end are quite similar. With decreasing cyclic amplitude, the hysteresis loop area in the stress–strain curve shrinks, indicating a gradual rise in axial strain increment across stress differentials. Consequently, the initial cycles show that axial strain at 1 MPa amplitude is smaller than at 50% amplitude. However, as the axial strain values diverge across stress differentials, the discrepancy in final axial strain under the two amplitudes narrows. During the 12 MPa cyclic loading and unloading stage, the axial strain difference in DF-01 increases from 0.23% to 0.6%.

Radial strain is more affected by varying amplitudes than axial strain. Before the initial stress drops to 10 MPa, the increase in radial strain during equal-amplitude cyclic loading and unloading exceeds that during the 50% amplitude stepwise process. The radial strain increment consistently remains within −0.16%. In terms of volume strain, equal-amplitude cyclic loading and unloading does not lead to volume compression. As indicated in

Table 7. Strain values under different cyclic amplitudes.

Pressure Levels	Axial Strain (%)				Radial Strain (%)				Volumetric Strain (%)
	BB-01		DF-01		BB-01		DF-01		BB-01		DF-01
	First	Last	First	Last	First	Last	First	Last	First	Last	First	Last
16	0	0.29	0	0.03	0	−0.03	0	−0.16	0	0.23	0	−0.29
14	0.42	0.45	0.06	0.17	−0.21	−0.2	−0.32	−0.48	0	0.05	−0.58	−0.79
12	0.54	0.72	0.3	0.53	−0.58	−0.59	−0.64	−0.8	−0.51	−0.47	−0.98	−1.07
10	1.17	1.4	0.79	1.39	−1.17	−1.31	−0.96	−1.27	−1.18	−1.22	−1.13	−1.15

, after the first cycle, the sample expands, with the volume strain reaching −0.29%, and increases to −1.15% by the end of the stepwise cyclic process. In the 50% amplitude stepwise process, initial compression is followed by expansion at each unloading stage. Ultimately, BB-01 exhibits greater expansion than DF-01, as observed in volume strain.

The results show that with varying amplitudes, tectonically deformed coal exhibits more stable macroscopic damage characteristics at smaller amplitudes. With increasing amplitude, the range of volume strain changes significantly, resulting in pronounced expansion effects. This trend increases the coal sample’s susceptibility to damage, ultimately causing instability.

4.2.2. Effects on Porosity

This section compares porosity evolution under various cyclic amplitudes to examine their effects on porosity expansion in tectonically deformed coal across different loading and unloading scenarios. The results demonstrate as shown in Figure 12 that under a 1 MPa amplitude, the increase in porosity of tectonically deformed coal is more pronounced than under the 50% amplitude loading and unloading method. Additionally, repeated cycles at the same stress level do not reduce porosity. Volume strain data also reveal that with equal-amplitude stepwise cyclic loading and unloading, the sample mainly expands across different stress differentials. There is also some expansion within the same stress level, contrasting with the compressive effect observed in the 50% amplitude cyclic loading and unloading at the same stress level.

The effects of equal-amplitude cyclic loading and unloading indicate enhanced stability in sample behavior. Specifically, the DF-01 sample demonstrates a more stable porosity evolution within each cycle, without experiencing compressive deformation. Additionally, the rate of porosity increase in DF-01 is higher than in BB-01 at different stress stages. Notably, at the fourth stress stage (13 MPa), the porosity of DF-01 surpasses that of BB-01 after completing all amplitude cycles. This observation suggests that in practical scenarios, cyclic loading with smaller amplitude and lower pressure differential leads to more pronounced expansion effects in tectonically deformed coal layers.

Previous research indicates that the permeability of brittle materials like primary structure coal tends to decrease with increased cycles of loading and unloading. However, mylonite coal exhibits significant elastic-plastic characteristics. The stress path employed in cyclic loading and unloading helps to mitigate major damage during rapid stress release phases in mylonite coal. This enhances the mechanical strength of tectonically deformed coal and facilitates controlled stress release, contributing to the stability of horizontal wells. High-amplitude cyclic loading and unloading are more effective in quickly increasing porosity and creating interconnected fracture networks, which benefits coalbed methane extraction. Additionally, as the stress differential between cycles increases, the porosity expansion effect also increases, but this leads to greater volumetric expansion and a higher likelihood of sample damage.

4.3. Crack Expansion Process and Force Chain Analysis

Physical experiments demonstrate that during cyclic loading and unloading, deformation tends to rebound and expand towards the direction of unloading. Figure 13 illustrates this crack expansion process. Initially, cracks are randomly distributed across the sample, but as the process continues, they gradually increase in number and form shear bands. Notably, in areas near these shear bands and in the direction of unloading, the number of cracks continues to grow significantly.

At the end of the first cycle, microcracks show a random distribution but soon start to exhibit an X-shaped shear dilation trend. As illustrated in Figure 7, during multiple cycles at the same stress level, the number of microcracks increases slowly and then stabilizes. This suggests that with repeated cycling, the tectonically deformed coal sample becomes progressively more compacted, and the elastic strain begins to surpass the plastic strain. In the stepwise cyclic loading process, as the confining pressure decreases, the previously compacted coal matrix becomes further relaxed. This relaxation not only leads to an increase in strain and volume strain, but also to a rise in the number of microcracks and porosity. The sample eventually reaches an equilibrium after several cycles at the next stress level.

Under cyclic loading and unloading, Figure 13 depicts the evolution of the micromechanical force field within the tectonically deformed coal. This figure is crucial for studying the micromechanical mechanisms of internal damage expansion under stepwise cyclic loading and unloading. It shows the size and distribution of internal contact force chains in the sample throughout the cyclic process. The thickness of these force chains represents the magnitude of the force, with thicker chains indicating larger forces. Initially, the force chains are uniformly distributed, reflecting an intact sample. However, after the first cycle of the first stage of unloading, the stress concentration and an X-shaped shear dilation pattern become evident in the sample. This pattern aligns with the locations of damage, suggesting that as local damage and microcracks develop, stress increasingly concentrates around these areas.

According to the research by Hou et al. [22], microcracks can amplify the applied stress, resulting in stress concentration phenomena. As the stepwise cyclic loading and unloading continues, the force chain distribution cloud diagrams in Figure 13, ranging from stage 1 to stage 7, demonstrate that with a further reduction in confining pressure, the internal contact force in the sample gradually decreases. However, the areas of stress concentration continue to primarily surround the main shear bands.

In actual stress release projects for tectonically deformed coal, the controllable conditions primarily involve managing the cyclic pressure differential between the pulsating pressure and the static hydrostatic pressure of the in situ coal seam, along with the amplitude of each pulsation. Initially, a low-amplitude, high-pressure differential cyclic loading mode is employed to achieve controlled induced expansion and rapid permeability enhancement. Subsequently, as the pulsating pressure approaches the peak deviatoric stress measured in triaxial tests of the coal, a mode of pulsation with low amplitude and pressure differential is utilized. This approach aims to minimize damage to the coal seam, preventing issues such as well collapse, and ensuring sustainable in situ coalbed methane extraction.

All our preliminary results throw light on the nature of coal pore-fracture network evolution under complex stress paths. The research has resulted in a solution for setting engineering parameters during the application of horizontal well cavity completion technology. Meanwhile, the findings of this study have to be seen in light of some limitations. First, the study is limited by the simulation tools available, precluding further simulation of permeability characteristics. Second, the evolution of pore-fracture networks from the end of cyclic loading to failure state was not further investigated. Future work should enhance research methods to better conduct numerical simulations related to fluid migration.

In summary, through cyclic loading and unloading particle flow numerical simulation analysis, it has been effectively confirmed that horizontal well cavity completion and stimulation (HWCCS) can increase porosity and form well-connected pore-fracture networks, facilitating future coalbed methane production.

However, as a promising new technology, much work remains to be done in the future: (1) Dynamic evolution of permeability in tectonically deformed coal samples based on fluid-solid coupling theory; (2) Reservoir characteristic evolution and multiphysics numerical simulation coupling methods after horizontal well completion; (3) Experimental and simulation studies on the in situ fluid migration mechanisms in tectonically deformed coal reservoirs.

5. Summary and Conclusions

(1)

The Particle Flow Code was employed to simulate the mechanical behavior and the evolution of the internal pore-fracture network structure in tectonically deformed coal samples under various cyclic loading and unloading scenarios. The outcomes of these simulations closely align with the results of indoor experimental studies. Furthermore, through extensive trial-and-error data analysis, a mathematical relationship was established between the micromechanical parameters of tectonically deformed coal and its macroscopic mechanical properties.

(2)

During the initial stages of cyclic loading, the internal contact forces within the sample are evenly distributed. However, as cyclic unloading progresses, localized areas begin to develop cracks accompanied by stress concentration. With continued expansion of the sample’s volume, the distribution of cracks in mylonite coal samples assumes an ‘X’ shape, becoming particularly dense near shear bands. This pattern mirrors the shear dilation effect observed in physical experiments. Before the complete failure of the sample, cracks eventually permeate throughout the entire specimen.

(3)

Under cyclic loading and unloading, the number of internal cracks in the sample increases stepwise with each level of stress reduction. In terms of porosity, significant improvements primarily occur between different stress differentials. As the cyclic amplitude increases, more cycles within the same stress level result in a compressive effect, leading to a higher likelihood of porosity decrease within these cycles. Conversely, with smaller amplitude single-stage cycles, porosity remains almost unchanged, and the final expansion effect is significantly greater than that observed with larger amplitude stepwise cycles.

(4)

Large-amplitude stepwise unloading proves more effective in enhancing the fluid network connectivity within tectonically deformed coal reservoirs. Conversely, smaller amplitude stepwise cyclic loading is better suited for maintaining wellbore stability. In practical engineering scenarios, employing a combination of different amplitude stepwise cyclic loading strategies can be highly beneficial. This approach allows for controlled extraction under stress release conditions in horizontal wells of tectonically deformed coal, thereby improving extraction efficiency.