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 ${𝜎}_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 $\ge$ 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.

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 ( ) during the fifth to eighth loading. Starting from the ninth loading cycle, the growth rate of micro-cracks gradually becomes stabilized ( ) 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.

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.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.

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.

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 $\ge 30\text{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.