Damage evolution of surrounding sandstone rock under charging–discharging cyclic loading in the natural gas storage of abandoned mines based on the discrete element method


Abstract

Gas storage in abandoned mines is one way to reuse waste space resources. The surrounding rock of gas storage reservoirs in underground roadways undergoes damage and deformation under the cyclic loading of gas charging and discharging, which can pose a risk to the safety of the reservoirs. This study establishes a true triaxial numerical model of rock mass with the discrete element method (DEM) and explores the crack evolution of surrounding rock of underground gas storage during cyclic loading and unloading. Also, a damage evolution model in numerical analysis considering residual deformation is developed to explain the experimental results. As was revealed, cyclic loading and unloading resulted in fatigue damage in the specimen and caused strength deterioration of the specimen. During the loading process, the uniformly distributed force chains of the rock mass redistributed, evolving gradually to mostly transverse force chains. This contributed to the appearance of blank areas in the force chains when through cracks appear. The ratio of tensile cracks to shear cracks gradually decreases and finally stabilizes at 7:1. The damage evolution model considering residual strain can be mutually verified with the numerical simulation results. Based on the DEM model, it was found that there was a certain threshold of confining pressure. When the confining pressure exceeded 30 MPa, the deformation to ductility of sandstone samples began to accelerate, with a greater residual strength. This study provides a theoretical basis for analyzing the long-term mechanical behavior of surrounding rock of gas storage in abandoned mines.

Highlights


  • This study established a true triaxial numerical model of rock mass with the discrete element method and explored the crack evolution of surrounding rock of underground gas storage during charging–discharging cycles.

  • Cyclic loading and unloading result in fatigue damage development in the specimen and strength deterioration of the specimen.

  • The damage evolution model considering residual strain can be mutually verified with the numerical simulation results.

  • When the confining pressure exceeds a threshold (30 MPa), deformation to ductility of sandstone samples begins to accelerate, with greater residual strength.

  • This study provided a theoretical basis for analyzing the long-term mechanical behavior of surrounding rock of gas storage in abandoned mines.


1 INTRODUCTION

The Paris Climate Agreement requires governments to take measures to reduce the extraction and use of fossil fuels (Fan, Liu, et al., 2020). China has implemented a policy for reduction of coal production capacity, eliminating outdated-capacity and depleted mines (Pu et al., 2021). China has invested considerably in the construction and maintenance of these mines. Accordingly, direct collapse of these mines is not only a waste of resources but also a loss to the national economy (Hu, 2019). In addition, China, the United Kingdom, and Brazil still face challenges in terms of gas and hydrogen storage facilities in terms of high construction costs and long construction periods (AbuAisha & Billiotte, 2021; Balland et al., 2018; Firme et al., 2019; Gao et al., 2022; Vandeginste et al., 2023). Mining production has high requirements of support strength and air-tightness of roadways. Therefore, when a mine is abandoned, a large amount of roadway space can be converted into gas and hydrogen storage spaces after essential upgradation and on meeting the standards and requirements of support strength and air-tightness. Specifically, the stability of surrounding rock under the action of charging and discharging of gas and hydrogen storage plays a crucial role in whether the abandoned mine can store gas effectively and efficiently.

The stability of surrounding rock is the critical factor to ensure the normal operation of a gas storage facility in abandoned mines (Hu, 2019; Wei, 2021; Wu et al., 2019; Zhang et al., 2021). However, in the long-term operation of underground gas storage, the gas charging–discharging cycle will cause fatigue damage to the surrounding rock of gas storage (Grgic et al., 2022; Yuan & Yang, 2021). The magnitude of the upper limit stress and amplitude of cyclic loading are the main parameters that need to be studied to assess the remaining service life of surrounding rock and gas storage (Chen et al., 2007; Ge et al., 2003; Xi et al., 2006; Xu et al., 2006; Zhang et al., 2006). As for the study of damage criteria of surrounding rock, many scholars investigated the lateral deformation of rock as the basis of analysis of macroscopic damage criteria (Bagde & Petros, 2005; Ge & Lu, 1992; Jiang et al., 2004; Ren et al., 2005; Xu et al., 2005). Recently, many scholars have utilized electron microscope scanning and CT real-time scanning to study rock damage under static loading and cyclic loading (Evans & Wong, 1992; Lu & Zhang, 1990; Wu et al., 1998; Yang & Zhang, 1998; Zhao & Huang, 1992). However, scanning methods are only applicable to the damage condition of a specific rock mass with standard dimensions and loading conditions. More studies need to be conducted to reveal the general damage evolution under cyclic loading and unloading.

Aiming to reveal the surrounding rock damage evolution under the action of charging–discharging cycles in abandoned mine gas storage, this study explored the damage evolution process of rock mass in the process of cyclic loading and unloading and developed a numerical model by the discrete element method (DEM) of the software PFC3D. The findings of this article can provide theoretical support for the reuse of the underground space of abandoned mines and promote environmentally friendly utilization of mining resources.

2 ESTABLISHMENT OF DEM MODELS OF SANDSTONE SAMPLES UNDER TRIAXIAL LOADS

Since cracks in rock mass mostly start at the paraxial boundary of quartz particles (Hu, 2019; Vandeginste et al., 2023), it is reasonable to assume that the failure of rock mass materials is related to the internal micro-structure of mineral particles (Alkan et al., 2007; Li et al., 2002; Zhang et al., 2022). The response of stone samples under cyclic triaxial loading was investigated by Alkan et al. (2007) based on combined acoustic emission and triaxial compression tests. Combined with the analysis of microstructure, Zhang et al. (2022) discussed the empirical fatigue model of salt rock under cyclic loading.

DEM is a numerical simulation method that focuses on the microstructure of rock materials and has advantages such as small displacement limitations and convenient dynamic monitoring. DEM has been widely accepted and used in research on rock crack development under static and cyclic loads (Fan, Xie, et al., 2020; He & Wang, 2018; Meng et al., 2017, 2021; Yuan et al., 2021).

In order to study the crack damage development and stability of surrounding rock under the cyclic action of gas charging and discharging in abandoned mines microscopically, a DEM model of cyclic loading and unloading test was established. To simulate the sandstone sample, the size of the model designed was 50 mm × 50 mm × 100 mm (length × width × height). This numerical model is shown in Figure 1. The contacts in the numerical model involve parallel bonding, with a total of 11 713 particles and 34 325 contacts.

Details are in the caption following the image
Numerical model of sandstone in discrete element method simulation.

From Figure 2, it is clear that the displacement of the specimen is basically symmetrical about the midplane before macroscopic failure occurs, with larger displacement on both sides (top and bottom). After the macroscopic fracture surface appears in the sample, the trend of displacement of the spherical element along the crack surface is significant.

Details are in the caption following the image
Diagram of displacement and failure form.

For the crack development process, there are no cracks in the numerical model at the initial state, and the sample is intact. When the sample is loaded to about 70% of the peak strength, randomly distributed cracks arise within the sample, but the distribution of cracks inside the sample could not be observed during indoor testing. When loaded to the peak strength, a macroscopic fracture surface can be seen on the surface of the sample, and the surface failure form can also be seen in the simulation results. As the load reaches the peak, the cracks continue to expand along the macroscopic fracture surface, and a significant failure surface can be seen in the simulation diagram.

By conducting the uniaxial compression test, the parameters in the DEM model were determined for sandstone samples. The microscopic parameters used in the numerical model of the cyclic loading and unloading test and the static loading test are shown in Table 1. As shown in Figure 2, the DEM model based on the parameters shown in Table 1 can provide simulation results that are consistent with the experimental uniaxial compression test results.

Table 1. Simulated microscopic parameters of sandstone samples.

image

In order to compare and analyze the cyclic loading and unloading model and obtain the peak value of the compression-bearing capacity of sandstone samples, conventional triaxial experiments with static loading were first conducted on the numerical model. The conventional triaxial loading process of the model is as follows: through the equilibrium of multiple time steps, the overlap between the particles reaches the pre-set reasonable values to stabilize the force and displacement.

The abandoned mine-type underground gas storage in Denver, the United States, can hold about 73.6 × 106 m3 of natural gas under a storage pressure of 17 MPa (Brown, 1978); the storage of the abandoned mine-type underground gas storage in Peronnes, Belgium, is 500 × 106 m3, and the allowable maximum gas storage pressure is 20 MPa (Sa, 1981). Accordingly, the confining pressure of the numerical model ( 𝜎 2 = 𝜎 3 ) is controlled by using the servo mechanism and set to 10, 20, 30, and 40 MPa. For example, if the elevation of the top of the model is 800 m, and the surrounding rock is sandstone, the upper boundary load would be 800 m × 2500 kg/m3 × 9.81 m/s2 ≈ 20 MPa, according to a depth of 800 m, density of surrounding rock of 2500 kg/m3, and acceleration of gravity of 9.81 m/s2. The confining pressures 10, 20, 30, and 40 MPa correspond to rock depths of 400, 800, 1200, and 1600 m, respectively, which are consistent with the mining depth of coal mines investigated in this study. After the σ3 is stabilized, the axial displacement is realized by controlling the displacement of wall units until a convergence condition is reached.

The cyclic loading process of the model is as follows: the initial loading is the same as that of the conventional triaxial test; after the σ3 is stabilized by the servo mechanism, the model is loaded to 70% of its peak strength and then unloaded until the axial stress is equal to the confining pressure. The loading and unloading process is then repeated cyclically. The strain of each loading is 100.1% of the previous strain level.

3 DEM SIMULATION RESULTS OF TRIAXIAL TESTS

3.1 Deformation of sandstone samples under conventional triaxial loading and cyclic loading

The stress–strain results of sandstone samples of the DEM model under conventional triaxial loading and cyclic loading are shown in Figure 3.

Details are in the caption following the image
Stress–strain curves of sandstone samples under conventional triaxial loading and cyclic loading. (a) σ 3 = 10 MPa, (b) σ 3 = 20 MPa, (c) σ 3 = 30 MPa, (d) σ 3 = 40 MPa.

As shown in Figure 3, by plotting the stress–strain curves of sandstone samples under conventional triaxial loading and cyclic loading, the mechanical responses of sandstone samples under the two loading modes can be observed and compared. According to Figure 3, when the 𝜎 3 is 10, 20, 30, and 40 MPa, there is a certain deviation between the stress–strain curve under cyclic loading and that of the conventional triaxial stress–strain curve. With the same strain, 𝜎 1 𝜎 3 under cyclic loading is always smaller than that under the conventional triaxial condition. This phenomenon confirms that the cyclic loading process causes certain fatigue damage to sandstone samples, which reduces the compressive strength. Bagde and Petros (2005) observed the impact of fatigue damage on the stiffness deterioration of rock samples under cyclic triaxial loading.

As can be seen from Figure 3, with the increase of 𝜎 3 from 10 to 40 MPa, the compressive strength of sandstone samples 𝜎 1 𝜎 3 in the conventional triaxial test is 91, 129, 167, and 210 MPa, showing a gradual upward trend. This shows that with the increase of σ3, the rock samples are subjected to greater confinement, and accordingly display a gradual increase in compressive strength. Considering the similarity between the conventional triaxial test results and the cyclic loading test results, the bearing capacity of rock samples under cyclic loading also increase with the increase of σ3.

As shown in Figure 3a, when σ3 is 10 MPa, the fatigue damage is small in the first few cycles after the peak under cyclic loading, but with an increase of the number of cycles, the strength decreases rapidly and the sample presents the characteristics of brittle failure. As σ3 increases from 10 MPa to 30 MPa, the values of 𝜎 1 𝜎 3 for both conventional triaxial and cyclic loading do not show a sudden decrease with the increase of strain, but stabilize after a decrease by 20%–30%. This shows that the failure mode of sandstone samples gradually changes from brittle failure to ductile failure with the increase of σ3, and the difference of 𝜎 1 𝜎 3 between each cycle decreases with the increase of loading times. This phenomenon was put forth by Bagde and Petros (2005) and Meng et al. (2021) based on the results of the experimental tests and DEM simulations that they carried out.

3.2 Microscopic damage characteristics of sandstone under cyclic loading

3.2.1 Deformation of the sandstone force chain under cyclic loading

Since the force chains are distributed along the axial direction of particles in the DEM model, the axial contact between particles can be represented by force chains. Also, the wide lines are applied to represent the strong contact forces (Nie et al., 2022).

In order to further study the microscopic mechanical mechanism of damage propagation within sandstone samples under cyclic loading and unloading, the internal force chain of the samples was analyzed using DEM, with σ3 being 10 MPa. The results are shown in Figure 4.

Details are in the caption following the image
3-D Distribution of the force chain under cyclic loading and unloading. (a) initial equilibrium phase, (b) the first cyclic loading to the peak value σ 1σ 3, (c) the first cyclic unloading to σ 1σ 3 = 0, (d) the fourth cyclic unloading to σ 1σ 3 = 0, (e) the seventh cyclic unloading to σ 1σ 3 = 0, (f) the eleventh cyclic unloading to σ 1σ 3 = 0, (g) comparison of force chains and cracks after the fourth cyclic unloading to σ 1σ 3 = 0, (h) comparison of force chains and cracks after the seventh cyclic unloading to σ 1σ 3 = 0.

According to Figure 4a, in the initial equilibrium phase, due to the internal equilibrium of the sample and no micro-damage, the force chains are distributed relatively uniformly. The force chain diagram after the first unloading is less dense than that in the initial phase (Figure 4c). After each cyclic loading, the number of transverse force chains decreases, resulting in a relatively sparse distribution of force chains (Figure 4c–f).

As can be seen from Figure 4b, at the beginning of loading, there are mainly tensile force chains in the transverse direction and compressive force chains in the longitudinal direction. Consequently, after loading, the transverse force chains are more susceptible to damage than the longitudinal force chains. This indicates indirectly that at this time, the vertical stress of the first principal stress is much larger than the lateral stress, and the load is transmitted vertically, but the cracks in the sample are not completely through at this time, so the force chains are relatively uniform. When the first cycle is unloaded to σ1 − σ3 = 0, as shown in Figure 4c, the sample is in the state of σ1 = σ2 = σ3, the force chain is redistributed, and the transverse force chains are restored to a certain extent, but not fully to the pre-loading state (Figure 4a). Therefore, after the first loading–unloading cycle, due to the transverse tension process of the internal particle group, the sample accumulates corresponding fatigue damage, resulting in the decline of the transverse tensile strength between the particles.

Figure 4d shows the distribution of the force chains of the sample after unloading in the fourth loading cycle. Compared with the first loading cycle (Figure 4c), the transverse force chains in the sample appear to be locally concentrated. At the same time, it can be seen through further analysis of Figure 4g that concentrated damage occurs in the sample, and some adjacent microscopic cracks are connected with each other, resulting in the appearance of macroscopic cracks. Figure 4g shows that the sample can resist further crack propagation by increasing the bonding force between spherical elements at the corresponding position.

After the seventh cyclic loading of the sample (Figure 4e,h), the damage of the rock mass further increases, with more microscopic cracks interconnecting and macroscopic cracks appearing. Meanwhile, the internal damage is aggravated. At this time, the bonding force at some cracks is not enough to resist crack propagation and accordingly bond failure occurs, and the bonding force returns to zero. This is manifested as a blank area of the transverse force chain in the force chain diagram.

By the time of the 11th cyclic loading of the sample (Figure 4f), the main crack of the sample is completely through and the rupture surface is fully developed. There is a continuous absence of transverse force chains in the force chain diagram, with the missing force chain regions representing the damaged areas that are connected from top to bottom.

3.2.2 Microscopic fracture propagation of sandstone under cyclic loading

Since the force chains shown in Figure 3 are distributed along the axial direction of particles in the DEM model, the axial contact between particles can be represented by force chains. Also, the wide lines are applied to represent the strong contact forces. Additional statements are added in the manuscript to discuss this aspect of the model (Nie et al., 2022). To further investigate the crack propagation in sandstone samples under cyclic loading, the change of the number of micro-cracks in the samples with axial strain during loading and unloading should be statistically counted and observed. The σ3 of 10 MPa can be taken as an example. In the process of cyclic loading, the relationship between the number of micro-cracks and axial strain is shown in Figure 5, and the relationship between the number of micro-cracks and loading timestep is shown in Figure 6. The envelope line of the loading–unloading curve in Figure 5 denotes the connecting line among the peak strain points.

Details are in the caption following the image
Curves of the total micro-crack number to strain.
Details are in the caption following the image
Curves of the micro-crack number to timestep.

According to Figure 6, the total number of micro-cracks is proportional to the loading times, and the total number of micro-cracks continuously changes in a stepwise manner during the whole loading process. The number of micro-cracks increases exponentially during the first cycle of pressurization, which is consistent with the crack propagation observed in Figure 4a–c. As also shown in Figure 6, before the sample reaches the residual strength stage, the growth rate of micro-cracks increases first and then decreases with the increase of cyclic loading times, and the growth rate of all micro-cracks reaches the maximum ( > 800 ) during the fifth to eighth loading. Starting from the ninth loading cycle, the growth rate of micro-cracks gradually becomes stabilized ( < 500 ) after the sample reaches the residual strength.

Figure 7 shows the distribution of micro-cracks generated in sandstone samples during cyclic loading and the propagation process observed and recorded in this study. Figure 7a–f correspond to the distribution of tensile and shear cracks in the specimen during the 1st, 3rd, 5th, 6th, 8th, and 11th loading cycles after unloading to σ1 − σ3 = 0, respectively. In these figures, the blue thin disk denotes the tensile crack and the red thin disk represents the shear crack.

Details are in the caption following the image
Micro-crack evolution of sandstone sample under cyclic loading and unloading. (a) the first cyclic unloading to σ 1σ 3 = 0, (b) the third cyclic unloading to σ 1σ 3 = 0, (c) the fifth cyclic unloading to σ 1σ 3 = 0, (d) the sixth cyclic unloading to σ 1σ 3 = 0, (e) the eighth cyclic unloading to σ 1σ 3 = 0, (f) the eleventh cyclic unloading to σ 1σ 3 = 0.

It can also be seen from Figure 6 that the cracks generated in the loading process are mainly tensile cracks, accounting for more than 80% of the total cracks. At the same time, according to the analysis of Figure 7, at the initial stage of crack initiation in the first loading cycle, with the gradual increase of tensile cracks, some cracks start to connect, leading to the generation of shear cracks, and the ratio of tensile to shear cracks reaches 40∶1. With the emergence of shear bands, the proportion of shear cracks gradually increases, and the ratio of tensile–shear cracks is 10∶1 when the cyclic shear bands appear at the fifth to eighth loading cycles. Toward the end of the 9th–11th loading cycles, corresponding to the residual strength stage of the sample, the ratio stabilizes at 7∶1.

Therefore, in the first four loading cycles, the growth number of micro-cracks in each stage is relatively smaller than that in the preceding stage. Therefore, during the fifth to eighth simulated loading cycles, due to the formation of shear bands, micro-cracks appear and expand rapidly, and the number of micro-cracks shows an obvious upward trend. Figure 7 shows the change process in detail. When the loading process enters the ninth loading cycle, the generation of micro-cracks gradually decreases and enters a gentle stage corresponding to the residual strength stage.

In Figure 7a, after the first cycle of the loading, there are few micro-cracks in the sample. After the third loading cycle (Figure 7b), the number of cracks increases, and some adjacent cracks are connected with each other. After the fifth loading cycle (as shown in Figure 7c), cracks in the sample are gradually connected, displaying a trend of two shear through-going cracks. The tensile cracks are relatively evenly distributed in the sample, while the shear cracks tend to gather only in the shear bands. At the sixth cycle (Figure 7d), a new shear band appears in the upper right corner of the sample. In the eighth cycle (Figure 7e), the fracture zone increases significantly and new shear zones also penetrate the main fracture. The upper right part of the whole sample shows severe damage and macroscopic through-fracture. After that, the sample enters the stage of residual strength. At the 11th loading cycle (Figure 7f), the cracks develop fully, and the new cracks mainly occur on the macroscopic fracture surface. Overall, during the cyclic loading and unloading process, more tensile cracks are generated, and they are more widely distributed. Shear cracks tend to gather within macroscopic shear fractures. Meng et al. (2017) also observed the concentration of cracks within macroscopic fractures by conducting DEM simulations.

4 EVOLUTION MODEL OF ROCK MASS DAMAGE UNDER CYCLIC LOADING AND UNLOADING CONDITIONS

In order to further validate the numerical simulation results, a damage calculation model considering residual strain is needed to further discuss the damage evolution of sandstone samples under different σ3 and analyze the numerical simulation results more comprehensively.

4.1 Damage evolution model of sandstone

Assuming that the damage rate conforms to a Weibull distribution, the damage evolution rate can be calculated using the following equation (Cao et al., 2003):
𝑣 d ( 𝜀 1 ) = 𝑚 𝑛 ( 𝜀 1 𝑛 ) 𝑚 1 e ( 𝜀 1 𝑛 ) 𝑚 , (1)
where 𝑣 d denotes the rock damage evolution rate; 𝜀 1 is the maximum principal strain; and m and n are constants, which are related to σ 3 and lithology.
If the strain value is ε t, the cumulative damage value can be expressed as
𝐷 𝑡 = 0 𝜀 𝑡 𝑣 d ( 𝜀 1 ) d 𝜀 1 , (2)
where D t denotes the damage variable at time t and ε t is the axial strain value at time t.
Solving Equation ( 2), damage variable D can be calculated from axial strain as follows:
𝐷 = 1 exp [ ( 𝜀 1 𝑛 ) 𝑚 ] . (3)
The dynamic elastic modulus E of the rock can be calculated using the equation
𝐸 = 𝐸 0 ( 𝑏 𝜀 1 1 + 𝑏 𝜀 1 ) , (4)
where 𝐸 0 represents the ideal elastic modulus of rock mass, which can be obtained by fitting the stress–strain curve at the linear elastic stage, and b is a constant.
Considering the damage evolution of a rock mass, the dynamic elastic modulus E can be expressed as
𝐸 = 𝐸 0 ( 𝑏 𝜀 1 1 + 𝑏 𝜀 1 ) [ ( 1 𝐷 ) + 𝑘 𝐷 ] , (5)
where k is the parameter, and 0 ≤  k ≤ 1, and the value of k is 0 under uniaxial loading, indicating that there is no residual elastic modulus after failure; in triaxial loading, the k value increases with the increase of σ 3. When a certain threshold is reached, the elastic modulus basically does not change.
The total strain of rock can be expressed by elastic strain and plastic strain:
𝜀 1 = 𝜀 1 e + 𝜀 1 p , (6)
where 𝜀 1 e represents the elastic strain of rock in the axial direction and 𝜀 1 p is the plastic strain of rock in the axial direction.
As can be seen from the cyclic loading and unloading curve in Figure 2, residual elastic strain still exists in the sandstone sample under the triaxial loading condition. As the plastic strain is caused by damage, therefore,
𝜀 1 p = ( 𝜀 1 𝜀 10 ) 𝐷 , (7)
where 𝜀 10 is the residual elastic strain of rock, which is affected by σ 3 and lithology.
By combining Equations ( 5) and ( 7) with elastic Hook's law, the following equation can be obtained:
𝜎 1 = 𝐸 0 ( 𝑏 𝜀 1 1 + 𝑏 𝜀 1 ) [ ( 1 𝐷 ) + kD ] [ 𝜀 1 ( 1 𝐷 ) + 𝜀 10 𝐷 ] . (8)

4.2 Validation of the damage evolution model

In this numerical simulation, the initial monitoring time of strain is set to be the moment when the axial pressure σ1 is equal to the confining pressure σ3 (i.e., σ1 – σ3 = 0). For more rigorous fitting with Equation (8), the confining pressure correction term σ3 should be added to the right-hand side of the equation. Figure 8 shows the stress–strain curve fitting results of the DEM model of sandstone with Equation (8) during cyclic loading and unloading.

Details are in the caption following the image
Fitting curves of stress–strain under cyclic loading–unloading. (a) σ 3 = 10 MPa, (b) σ 3 = 20 MPa, (c) σ 3 = 30 MPa, (d) σ 3 = 40 MPa.

The fitting parameters are shown in Table 2. According to Figure 8, the sandstone damage evolution model (Equation 8) can well fit the stress–strain curve of rock under cyclic loading and unloading; the greater the σ3, the better the fitting results reflect the simulated situation of the DEM model. With the increase of σ3, the residual elastic strain increases and the sample gradually shows a trend of ductility transformation. After the stress peak, the axial stress remains stable as the axial strain continues to increase.

Table 2. Fitting parameters of the damage evolution equation.
Confining pressure, σ3 (MPa) E0 (GPa) b k m n 𝜀 10 (%)
10 9.97 40 0.90 1.17 9 0.50
20 10.73 90 0.88 1.67 6 0.70
30 12.89 60 0.95 1.72 6 0.90
40 14.25 60 0.99 1.84 4 1.28

According to Table 2, the ideal elastic modulus 𝐸 0 of the sandstone sample during the first cycle of cyclic loading is the same as that of conventional triaxial loading, which can be obtained from Figure 3. The 𝐸 0 increases with the increase of σ3. However, the value of b shows no obvious rule and its effect on the results is negligible. The value of k does not vary significantly within the range of σ3 from 10 to 40 MPa. The values of m and n are related to σ3 and the rock properties. The value of m is smaller at low σ3 and becomes stable after the σ3 exceeds 20 MPa. The value of n shows a gradually decreasing trend. The trends of the fitting parameters of this study are consistent with the findings of Cao et al. (2003).

When the σ3 is low (10 MPa), due to the setting of loading convergence conditions, the residual elastic strain of the sample is small (0.50%) and the sample does not show significant residual strength. As the σ3 increases, the residual elastic strain 𝜀 10 gradually increases, indicating that the residual strength of the sample increases gradually (refer to Equation [8] and Figure 8). When the σ3 of the sample is 40 MPa, the sample shows pronounced residual strength characteristics with a larger residual elastic strain 𝜀 10 (1.28%) and its deformation shows ductile characteristics.

4.3 Damage evolution based on fitting parameters

The variation of damage variable D with axial strain is obtained by fitting parameters, as shown in Figure 9.

Details are in the caption following the image
Curves of D to axial strain.

As can be seen from Figure 9, the variation curves of damage variables with axial strain in the model can be divided into three stages: damage accumulation stage (axial strain 0–1%), rapid increase stage, and stable stage (D ≥ 0.9). When σ3 is 10 MPa, the damage accumulation stage is relatively short, and the sample quickly enters the rapid increase stage. When the σ3 is 30 MPa , the damage accumulation rate in the rapid increase stage decreases as the σ3 increases, although the magnitude of the decrease varies. There is a certain threshold of σ3. When the σ3 is below this threshold (30 MPa), the decreasing amplitude of the damage accumulation rate gradually decreases with the increase of σ3. However, when the σ3 exceeds the threshold, the decreasing amplitude of the damage accumulation rate of the sample increases synchronously with the σ3. Therefore, when the σ3 exceeds a certain threshold, the rate of deformation transformation into ductility in the sandstone sample accelerates, while showing greater residual strength (Figure 8). In order to obtain the threshold value mentioned above, Hu (2019) also conducted a conventional triaxial experiment and true triaxial experiments. The data of these experiments showed that when the confining pressure was close to 30 MPa, the full stress–strain curves transformed into ductile failure, indicating that 30 MPa was very close to the threshold value. Therefore, with the verification from the experimental results, the σ3 of 30 MPa can be taken as a threshold of deformation of sandstone samples.

5 CONCLUSIONS

By using DEM, this study established a numerical model of surrounding rock of gas storage in abandoned mines and simulated the fatigue damage development and deformation of sandstone rock mass during gas charging–discharging cycles. Meanwhile, the damage evolution model considering residual elastic strain was developed to verify the simulation results. The following conclusions were obtained.
  • 1.

    Cyclic loading and unloading on the sandstone sample caused fatigue damage in the sample, resulting in a decrease in its strength. With the increase of the number of loading–unloading cycles, the proportion of brittle failure increased and the strength decreased rapidly. With the increase of confining pressure σ3, ductile failure occurred gradually. During the loading process, the uniformly distributed internal force chains of the rock mass sample redistributed and transformed into mostly transverse force chains. In the failure stage of the specimen, there were concentrated blank areas in the force chains when the through crack appeared.

  • 2.

    With the increase of gas charging–discharging cycles, the crack growth rate increased first and then decreased until reaching the stage of residual strength, and the crack development gradually became steady. Meanwhile, tensile cracks appeared earlier than shear cracks. In the early stage, the number of tensile cracks was 30–40 times that of shear cracks. With the emergence of shear bands, the proportion of shear cracks gradually increased.

  • 3.

    The damage evolution model considering residual strain fit well with the numerical simulation results, and the fitting effect was better with the increase of σ3. After that, there was a certain threshold for the σ3. When the σ3 was below this threshold (30 MPa), the decrease rate of damage accumulation reduced as σ3 increased. However, when the σ3 exceeded 30 MPa, the decrease rate of damage accumulation increased simultaneously with σ3. Therefore, after reaching a certain threshold for the σ3, the rate of deformation to ductility in the sandstone sample accelerated, while showing larger residual strength.

ACKNOWLEDGMENTS

The authors would like to acknowledge the financial support received from the National Natural Science Foundation of China (No. U22A20598; No. 52104107), the National Key Research and Development Program of China (No. 2023YFC2907300; 2019YFE0118500; No. 2019YFC1904304), and the Natural Science Foundation of Jiangsu Province (No. BK20200634).

    CONFLICT OF INTEREST STATEMENT

    The authors declare no conflict of interest.

    Biographies

    • image

      Zhanguo Ma, PhD, is a professor in Civil Engineering at the China University of Mining and Technology. Prof. Ma specializes in the fields of geotechnical engineering and mining engineering. He teaches geotechnical engineering and rock mechanics to undergraduate students and advanced rock mechanics to postgraduate students.

    • image

      Junyu Sun, PhD, is a postdoctoral research fellow at the China University of Mining and Technology. He obtained his PhD degree from the School of Engineering and Built Environment, Griffith University, Gold Coast, Australia, in 2023. In addition, his current research interest includes underground pavement structural evaluation. He is also a multidisciplinary design engineer specializing in airport pavement design and underground pavement design.