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.
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.
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.
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 (
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.
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
As can be seen from Figure 3, with the increase of
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
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.
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.
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 (
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.
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
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.
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.
Confining pressure, σ3 (MPa) | E0 (GPa) | b | k | m | n |
|
---|---|---|---|---|---|---|
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
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
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.
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
5 CONCLUSIONS
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
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.
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.