Coal is a typical sedimentary deposit principally composed of carbon and hydrocarbon that is easily combustible. So far, coal is still the most plentiful fossil fuel in the world. The extraction of coal via underground mining (opposite to surface mining) currently accounts for an obviously larger share of global coal production than open pit. For the past several decades, China has dominated the coal mining and steadily maintained its position as the largest coal-production country, and accounted for approximately half of the global coal production in 2021 (Garside 2022). Due to the exhaustion and depletion of shallow coal in east China, more and more mega mines with thick coal seams (Wang et al. 2021) are constructed and shifted to the central and western parts of China and the coal production in these regions now can account for over 70% of total domestic coal production (Chi et al. 2019). However, the precipitation volume is unexpectedly low in these regions, e.g., in provinces of Shaanxi, Inner Mongolia, and Ningxia, where are almost arid and semiarid climates. The precipitation volume of these regions merely accounts for approximately 4% of China’s total precipitation volume (Chi et al. 2022). The low precipitation volume results in an ultra-fragile natural and ecological environment, the coal mining in these areas can easily incur permanent damage to natural circumstances and exacerbate the environmental vulnerability, such as desertification, surface water loss, vegetation degradation. To tackle such a serious issue, an engineering conception and framework: the coal mine-based underground water reservoir (CMUWR) was proposed by Gu (2015) to utilize the coal pillars (here refer to the pillars which are intentionally left intact to protect against roof collapse, see the coal pillars indicated in Fig. 1a) in the gob as the dams with the connection of reinforced concrete artificial dam (Yao et al. 2020a, b), a schematic diagram for CMUWR is shown in Fig. 1a. It is noticed that the dam composed of coal pillars is subject to the load of overlying strata, which can be either static or dynamic: e.g., strata weighting, seismic-events, drilling, and mining operations in adjunct mining stopes. With development of seismic technology, the monitoring of seismicity at the in-situ scale is critical to improving mining efficiency and operation safety (Cai et al. 2020). Unravelling the fundamental pattern of seismic waves is a prerequisite for data analysis and numerical model construction. A real seismic wave in the form of vibration velocity vs. time is shown in Fig. 1d. Two areas (A and B) are zoomed-in as shown in Fig. 1e, f, where the distinct rates (slopes) can be noticed at the rising and falling phases. This represents diverse loading and unloading rates and this kind of cyclic stress is termed “Differential Cyclic Loading (DCL)” in our previous research (Song et al. 2022c, b). Although the coal’s response under compressive cyclic loading has a long record of research (Li et al. 2019b; Zhang et al. 2020; Duan et al. 2021), the attempts to understand the fatigue behaviors by incorporating the diverse loading and unloading rates still have massive unresolved issues and remains unclarified. The coal dams of CMUWR suffer from oscillation and dynamic load due to the effect both natural or human-induced activities: e.g., earthquakes (Brown and Hudson 1974), seismic events (Dang et al. 2018), hydraulic drilling (Zhuang et al. 2019; Yu et al. 2021), shearer cutting (Wang et al. 2015), underground blasting (Yu et al. 2021). The previous studies from Duan et al. (2021) and Nikolenko et al. (2021) have revealed that the compressive cyclic loading can obviously deteriorate the strength and modulus of coal samples, a low failure stress is always observed in cyclic loading. Whether the DCL can exacerbate or mitigate the fatigue damage of coal is still not sufficiently revealed.

a An overview of coal mine-based underground water reservoir (CMUWR), after Zhang et al. (2021a), the mechanical model of the coal dam can be thought of as a coal sample exposed to compressive cyclic loading due to the imposed vertical load from strata. b Passive seismic velocity tomography presented by with the aid of seismic technology, from Cai et al. (2016). c A layout of the sensor station to monitor the seismicity and record the seismic wave, after Cao et al. (2016). d A real specific seismic wave. e A zoomed-in view of section A in d, a rapid loading and slow unloading pattern is noted. d A zoomed-in view of section B in d, a slow loading and rapid unloading pattern is noted

Except for in-situ monitoring and lab test, numerical characterization of coal’s response under DCL is also of fundamental relevance to evaluate the stability of CMUWR especially when the coal dams are frequently subject to dynamic disturbances. The particle based numerical modelling can deliver a more intricate and comprehensive analysis on damage mechanisms and fracturing of porous materials from crack initiation, coalescence, and final rupture (Hazzard and Young 2002b; Potyondy and Cundall 2004; Hoek and Martin 2014). This can facilitate the calibration of model not only from the view of its constitutive relations but also the associated mechanical responses, e.g., AE behaviors and crack patterns. The numerical characterization of the time-stress-dependent damage evolution incurred by constant, monotonic or cyclic loading via particle model has been studied by a plethora of academics (Chang et al. 2002; Hossain et al. 2007; Potyondy 2007; Li et al. 2019a; Song et al. 2019a, 2022a). The recent investigations on numerical cyclic loading are highlighted as follows: Liu et al. (2018) used the particle model numerically to reproduce the results of cyclic-Brazilian test based the flattened Brazilian disks, the development of microcracks is consistent with the observed experimental result, and a macrocrack lies vertically along the central line. The damage in numerical model is based on the degradation of elastic modulus and bond strength with ongoing of cycles. Li et al. (2022a, b) numerically investigated the behavior of mono/hybrid pile foundation under horizontal cyclic stress, the results indicate that the rise of pile stiffness and application of hybrid pile foundation could enhance the secant stiffness and restrict cumulative lateral displacement. Bai et al. (2022a) conducted a particle model-based numerical investigation to reveal the effect of cyclic normal stress on frictional behaviors of jointed rocks, the phase shift incurred by dynamic load was interpreted and was found related to the changes of the local slipping and frozen contacts. Our former work (Song et al. 2019a, b, 2022a, 2023a) developed various numerical models to replicate the mechanical and micro-seismic behaviors of concrete and coal in constant and variable amplitude cyclic loading. The damage is characterized via the stress-induced bond property degradation. These abovementioned articles numerically examined the applicability of PFC in simulating mechanical responses of different materials under various loads. However, these studies are based on identical loading/unloading rate. Up to now, the numerical simulation on DCL tests is still rarely reported due to the lack of robust lab data concerning such a novel loading regime.

This work is motivated by the formerly conducted DCL tests on coal and strives to more in-depth understand its underlying damage mechanisms induced by distinct loading and unloading rates (Song et al. 2022d). An explicit modelling is conducted based on 3D particle model to simulate the mechanical behaviors. Then, from the microscale, a specific attention is paid to how to numerically replicate the stress–strain curves, micro-fracturing, acoustic emission behaviors observed in lab test. The study can facilitate to deepen the understanding of the mechanism how the DCL impacts the coal’s mechanical responses and can provide practical guidance to the stability assessment of coal pillars in underground coal mines and infrastructures.

In this section, we briefly present the major results of previous laboratory tests (Song et al. 2022d). The compressive DCL tests were conducted based on an advanced closed-loop servo-controlled hydraulic machine and the geometry of cylindrical anthracite coal sample is [50 mm diameter × 100 mm length]. The characterization of section morphology of a coal sample used in lab testing is presented in Fig. 2 in three scales (Song et al. 2022d). Figure 2c shows that the intersectant cracks are densely distributed and the size of coal matrix is close to the order of millimeters. If we characterize the coal matrix and fracture with disc in a 2D view (Fig. 2d), it has a similar contact with the 2D particle model. The particles intersected by cracks are subject to restricted movements in normal and shear directions. The micro-structural characteristics indicate that the coal is suitable to be modelled via discrete method due to densely distributed cracks and matrixes, where normal contact and frictional behaviors can be well characterized via assigned normal and shear stiffness and frictional parameters. To avoid repeating published content, we only present the relevant data and compare with the numerical results in this study. The calibrations in numerical model with the experimental data are mainly from two parts: monotonic loading (see C16, C18 in Fig. 3a) and DCL tests (see Regime 1: C5, C11, C19; Regime 2:C6, C22, C21 in Fig. 3a). The physical properties for these samples are presented in Table 1. Figure 3 shows the distribution of microcracks and microvoids via thin section of anthracite. The monotonic tests provide a reasonable reference UCS, and the stiffness parameters are recommended to be separately calibrated with each dataset of cyclic loaded sample. The AE characteristics (AE count, b-value, magnitudes of AE events), stress–strain relations, fracture evolution are considered and calibrated with both monotonic and DCL tests.

a An overview of section morphology of a coal sample in lab test (Song et al. 2022d). b A zoomed-in view with the plotting scale of 1 mm. c A zoomed-in view with the plotting scale of 200 μm, the cracks are indicted with the red solid lines. d A characterization of coal matrix with particle models and the order of the size is indicated

a An overview of the cylindrical coal samples (Song et al. 2022d), the data of samples used in this study is indicated. b A graph of the micro structure of coal thin section, the microvoids and microcracks are densely distributed

The numerical simulation is based on PFC^3D, the micro contact model used in this work is the linear parallel bonded model (LPBM). The readers want to understand fundamentals of LPBM can visit the issued PFC official manual (Itasca Consulting Group Inc. 2014) or the article from Potyondy and Cundall (2004).

3.1 Monotonic loading and calibration

The coal samples were subject to wetting–drying (WD) treatments preceding cyclic testing. We choose six samples exposed to two DCL regimes (Regime 1: C5, C11, C19; Regime 2: C6, C22, C21) to calibrate and replicate the mechanical responses. Prior to the DCL simulation, two monotonically loaded samples (C16, C18) were selected to calibrate the quasi-static monotonic strength. As documented in our previous work (Song et al. 2023b), the model with heterogeneous microparameters exhibits a more realistic AE behavior as well as a reasonable development process of microcracks. In the calibration of monotonic loading, we compare the results obtained in homogeneous and heterogeneous models. Due to the lack of AE monitoring in cyclic test (Song et al. 2022d), the calibration of AE parameters is conducted based on a dataset from Zhang et al. (2020), where the coal sample holds a similar monotonic strength with the samples in our work (Song et al. 2022d). The AE data were captured via 8 nano sensors. Only AE events detected by at least four sensors are used for source location. As recommended in Song et al. (2023b), a Weibull distribution with the shape parameter m = 2.5 is implemented in both bond tensile strength (τ_ten) and normal and shear parallel bond stiffness (k_n’ and k_s’,the ratio of k_s’/ k_n’ is fixed to 2). The homogeneous model with identical (τ_ten) and (k_n’ and k_s’) was established as a control group. To assign Weibull distribution based-micro parameters, we used a method introduced by Liu et al. (2004) and Wang et al. (2019) and was also used in our former work (Song et al. 2023b). First of all, the parameter in the homogeneous model, here we call it Q₀, should be determined with calibration of lab data. Then the parameter under Weibull distribution is termed Q_i. The Q₀ and Q_i conform to such a relation: Q_i = Q₀*x, the x is called the random variable. The cumulative distribution function (CDF) and probability density function (PDF) for x are expressed in Eqs. (1) and (2), where x_u is threshold and is recommended to be 0 (Vales et al. 1999; Liu et al. 2004; Song et al. 2023a). m is shape parameter, and x₀ is an expectation related scale parameter (Wang et al. 2019). The Eqs. (1) and (2) are turned to Eqs. (3) and (4) when x_u = 0. The mean E(x) and variance Var(x) are given in Eqs. (5) and (6). The generation of x is realized via the solution of inverse function of which is controlled by m and x₀ (Wang et al. 2019). F(x) varies between [0–1], therefore can be denoted by a random number u between [0–1] The random number can be obtained from a universal random number generator. The x and x₀ can determined from Eqs. (7) and (8):

Figure 4 details the calibration of monotonic strength. Figure 4a, b respectively show the statistical distributions of τ_ten and k_n’ under Weibull distribution (m = 2.5) in numerical model via Python. In Fig. 4c, it is seen that the UCS for C16 and C18 are respectively close to 10 and 16 MPa. To guarantee an identical sample strength, a comprise is made on strength diversity and the reference UCS is set to 16 MPa, the higher strength allows for a more reasonable range for the variation of bond radius multiplier in the following cyclic loading. The homogeneous and heterogeneous models both show a similar elastic modulus and UCS with respect to the results of C18 in Fig. 4c, which validates the rationality of parameters. The corresponding discrete (AE counts recorded in the fixed time steps) and cumulative AE counts in two models are plotted in Fig. 4d, where the first AE event occurs at approximately 5% UCS in heterogeneous model whereas the value is much higher in homogeneous model and almost reaches 60% UCS. Due to the lack of AE monitoring in samples C16 and C18, we adopt an AE dataset (H4) from the same series of coal (Zhang et al. 2020) which has a similar UCS. The stress–strain curves and the AE results are shown in Fig. 4e, f. The first located AE event occurs around 3% UCS which is close to the level in heterogeneous numerical model (5%). The pattern of cumulative AE counts for H4 is highly consistent with the heterogeneous model in Fig. 4c. Figure 4g, h respectively present locations of AE sources in heterogeneous model and lab test. Three levels (small, medium, and large) of AE events are classified according to the magnitudes of AE energy. The amount of AE energy is proportional to the released energy at failure of the bond (Song et al. 2019b, 2023b). The evolutions of three levels of cumulative AE events are respectively indicated by black, red and blue solid lines. The general pattern and proportions of three levels of cumulative AE events in heterogeneous model highly coincide with lab data which further validates the reasonability of heterogeneous model.

Illustration of Weibull distribution in a tensile strength of bond and b parallel bond stiffness in heterogeneous model, the shape parameter of m is set to 2.5. c Stress–strain curves of C18, C16, homogeneous and heterogeneous models in calibration of uniaxial compressive loading. d AE events and cumulative AE events in homogeneous and heterogeneous models. e Stress–strain curves of H4, where the coal samples is from the same series (Zhang et al. 2020). f AE events and cumulative AE events in lab test of H4. g Evolution of three levels of AE events in heterogeneous model. h Evolution of three levels of AE events in lab test of H4

We run 5 sets of numerical models with different random seeds and the stress–strain curves are shown in Fig. 5a. The comparisons between numerical models (18.02 ± 1.59 MPa for UCS; 1.125 ± 0.182 GPa for secant modulus) and experimental results (17.87 MPa for UCS; 1.13 GPa for secant modulus) in terms of secant moduli and UCS are close to each other. The scatters in numerical models indicate that the change of random seeds can cause the moderate variation of the UCS and secant modulus. The random features can directly influence the evolutions of the bond radius multiplier in the following cyclic simulation. Usually, a higher strength is associated with a smaller bond radius multiplier at failure and vice versa. The uncertainties incurred by particle fabric should be considered in simulation, which is also another heterogeneous characteristic in PFC model. In LPBM, a part of stored energy will be released at failure of a bond which shares a similar physical mechanism with a realistic AE event (Hazzard and Young 2002a; Bai et al. 2022b). Some parameters (e.g., stored bond energy (Song et al. 2023b), moment tensor (Bai et al. 2022b)) were used to characterize the magnitudes of AE energy. The Gutenberg-Richter relationship (Gutenberg and Richter 1944) has been used as a basic characterization of the size distribution of earthquake events and rock failure. It is mathematically expressed as below:

where N is the number of earthquakes with magnitude larger than M, a and b are fitting parameters. b describes the size distribution scaling, namely, b-value. Due to the nature of M is related to emitted energy. The M in PFC can be characterized by a scaled released energy at bond failure. M in lab tests is normalized via a scaled absolute energy, which is consistent with M calculated in numerical simulation. The scaling factors for normalization in PFC and lab tests are respectively, 1 × 10³ and 1 × 10⁻⁷. The sample sizes for calculation in PFC and lab test are 400 and 40 respectively. 4 time points are selected to calculate the b-value in numerical model and lab test. Due to different scaling factors, the for b-value may exhibit significant disparity. However, the declining of b-value with the rise of stress in numerical model is consistent with lab data ahead of failure. The details of calculation of b-value of the four timepoints are presented in Fig. 5d, e. The calibrated micro-parameters subject to monotonic loading are plotted in Table 2.

Comparisons of b-value in heterogeneous model and H4 data based on 4 timepoints (1–4): a Stress–strain curves for heterogeneous model under 5 random seeds and lab results for H4. b b-value of the four timepoints in heterogeneous model, the results is based on the case of seed 1. c b-value of the four timepoints in lab test of H4. d The calculation of b-value in heterogeneous model, the magnitude is normalized by a factor of 1e3 based on the AE energy captured in the code. e The calculation of b-value in lab test of H4, the magnitude is normalized by a factor of 1e-7 based on the absolute AE energy captured in lab test

3.2 DCL test and calibration

3.2.1 Stress propagation

The frequency in DCL test is 0.1 Hz, the low stress rate makes the DCL tests fall in the range of quasi-static loading (Zhang and Zhao 2014), where the inertial forces are negligible in this case according to Cerfontaine and Collin (2018). The checking of stress propagation in numerical model is essential to evaluate the accuracy of loading state. The stress-controlled loading regime is similar to the one used in our former study concerning concrete fatigue (Song et al. 2019b, 2023b), where the vertical stress is imposed via the clump composed of densely packed (with overlap) pebbles set on the top surface of assembly and the rotation and horizontal movements are completely restricted. The vertical stress monitoring is conducted at five parts: namely, the input signal, contact stress at two ends of sample, stress in three measuring spheres, the positions of three measuring spheres are shown in Fig. 6c. Figure 6a, b present the stress paths of two modes in previous DCL tests: (1) rapid loading and slow unloading (2) slow loading and rapid unloading (Song et al. 2022d). In Fig. 6c, d, the monitored vertical stresses for C5 and C6 are shown, almost no inertial forces are generated during cyclic loading, which means that the vertical stress recorded at different parts are consistent with the pattern of input stress wave, and meets the requirement of stress propagation under quasi-static loading. This substantiates the accuracy of stress-controlled loading regime constructed in PFC.

a Stress paths for the samples (C5, C11, C19) exposed to DCL of rapid loading and slow unloading, after Song et al. (2022d). b Stress paths for the samples (C6, C22, C21) exposed to DCL of slow loading and rapid unloading, after Song et al. (2022d). c Stress monitoring in simulation in the case (C5) of rapid loading and slow unloading. d Stress monitoring in simulation in the case (C11) of slow loading and rapid unloading, the stresses at five parts (two sample ends, input signal, measuring spheres 1, 2, 3) are monitored in simulation

3.2.2 Evolution of bond radius multiplier and axial strain

The micro-scale damage of particle assembly characterized by the deterioration of bond size, termed “parallel-bonded stress corrosion (PSC) model”, was first proposed by Potyondy (2007) to simulate a stress-time-dependent deformation due to the corrosion reaction of rocks in presence of water. Since then, the following research (Song et al. 2019a, 2022a; Shen et al. 2022; Xu et al. 2022) implemented the time-dependent variation of bond radius multiplier to numerically replicate mechanical responses of rocks under fixed/variable load amplitudes or environmental actions. The original PSC model idealizes the rock as a cemented granular material and then removes the connection based on a stress-level dependent rate. The corrosion rate is expressed as shown in Eq. 10 (Potyondy 2007), where σ_c is the bond strength, σ_a is the crack-initiation stress, and β₁ and β₂ are two constants to determine the corrosion rate. The original PSC model assumes that the rate of bond size reduction remains constant when the load level is fixed and exceeds the σ_a. Our previous work on fatigue behavior of concrete indicates that the variation of bond size can be in various manners, the multiplier can be even larger than 1 to mimic an enhancement of elasticity. The reduction in bond diameter is specified by the scaled ratio β, the diameter is reduced to βD, both the effective stiffness and the bond force were reduced (Potyondy 2007). The strain rate of the assembly is determined by bond radius, stress level and crack number. The accurate evolution of bond radius multiplier should be determined with calibration of stress–strain relations.

To reach a balance between computation cost and model complexity, the radii of particles (1.8–2.3 mm) in DCL simulation are slightly larger than the values in monotonic simulation (1.25–2 mm). The strength-related micro parameters in DCL simulation is consistent with the calibrated parameters in monotonic loading. Figure 7 shows the path of bond radius multiplier after calibration with stress–strain relations of the six selected samples. The calibrated numerical stress–strain curves are indicated by orange dash lines and shown in Fig. 8 together with its corresponding lab test results. Here we want to clarify a phase shift is imposed to each sample to better compare the stress–strain curves during cyclic loading, the magnitudes of phase shift are indicted at the top left corner of Fig. 8. In numerical simulation, a user-defined phase shift is usually adopted to offset the cycle gaps when only partial cycles were used to calibrate simulation or to indicate other parameters: e.g., the strain rate and secant modulus (Song et al. 2019b, 2023b). As known from Fig. 8, the stress–strain curves exposed to rapid loading and slow unloading display a more evident “loose” characteristic which is a manifest indication of larger energy dissipation and a higher growth rate of axial strain. This well correlates with the low fatigue life occurred in this regime and implies that the violent energy dissipation and drastic strain growth can occur much earlier in the first regime than the other one. The evolution path of bond radius multiplier in Fig. 7 between two loading regimes are completely distinct. For the first regime (C5, C11, C19), the radius multiplier remains constant or is monotonically reduced with ongoing of cycles. This coincides with the principle of PSC model, which means that the corrosion gradually deteriorates with rise of stress and increase of cycles. The secant moduli are indicated via colorful bubble diagram (Song et al. 2022d), a high consistency is noted at the moment between drastic drop of moduli and bond multiplier reduction for the samples exposed to rapid loading and slow unloading. For the second loading regime (slow loading and rapid unloading), it is noted that the multiplier first exhibits a step rise with elevation of stress, then the multiplier is reduced and the failure is reached. It can be interpreted from the view of two competing mechanisms: stiffness hardening and softening. Normally, the growth rate of axial strain is increased when the stress level is elevated and the bond radius multiplier remains unchanged. As known from the stress–strain results shown in Fig. 8, the growth rate of the axial strain doesn’t increase significantly after the elevation of stress level, which means that the bond radius multiplier should be raised to offset or restrain the effect of stress elevation. For Figs. 8d, e, two phases respectively correlating stiffness hardening and softening are clear, the peak of the secant modulus is the connection point of the two phases, the variation of radius multiplier varies distinctively before and after the point. For the samples in the second regimes, the failure all occurred in monotonic loading phase after 62 cycles. Therefore, a linear reduction of bond radius multiplier is implemented in this phase to characterize an intensive development of damage prior to failure. Figure 9 presents the comparisons of zeroed axial strains (values of the axial strain of the first cycle is zeroed) at maximum (P) and minimum (R) stresses between lab tests and numerical simulation. The results clearly indicate the applied path of bond radius multiplier in numerical simulation delivers an excellent match with lab testing results and the accuracy of the model is validated.

Evolution path of applied bond radius multiplier versus timestep and the evolutions of secant loading moduli (colorful bubble diagram) and stress levels are indicated a C5 b C11 c C19 d C6 e C22 e C21

Comparisons of stress–strain curves between lab test and numerical simulation, a phase shift is added between lab test and simulation results to offset the strain gap, the phase shift is indicated at the left top corner of each figure a C5 b C11 c C12 d C6 e C22 f C21, the lab data are from Song et al. (2022d)

Comparisons of the zeroed axial peak and residual strains between lab test and numerical simulation, for the convenience of comparison, the axial strains were zeroed from the first cycle, the lab data are from Song et al. (2022d)

3.2.3 AE behaviors in DCL tests

3.2.3.1 AE energy evolution

Numerical AE system can reduce the labor intensity, especially when AE source location is used. A 3D real-time AE system is implemented in PFC^3D during both monotonic and DCL tests. For DCL tests, Fig. 10 plots the evolution of AE parameters from beginning to failure based on the samples exposed to DCL. Seven representative moments (a-g) are indicated in Fig. 10 which are used thereafter to show the spatial distribution of AE event magnitude and to calculate the numerical b-value. The AE evolutions for the samples under the first DCL regime (rapid loading and slow unloading) obviously exhibit several “active” phases particularly when the cyclic stress is raised or bond radius multiplier is sharply reduced. For the results in the second regime (slow loading and rapid unloading), the bond radius multiplier is raised since the second cyclic loading stage (CLS), the enhanced multiplier evidently boost the active growth of AE energy in the phase with lower stress level. A transient spike of AE energy is noticed in Fig. 10d, this occurs without reduction of bond radius multiplier and is possibly induced by the elevation of localized contact force in the heterogeneous model.

Evolution of AE energy in numerical model a C5 b C11 c C19 d C6 e C22 f C21

3.2.3.2 Spatial distribution of AE events magnitude

We particularly select seven moments (indicated in Fig. 10) correlating different stress levels and cycle numbers to present the spatial distribution of AE energy magnitudes, see Fig. 11. The size of a cubic in Fig. 11 is proportional to the amount of released energy at bond failure. The results clearly demonstrate that the sample under the loading regime with slow loading and rapid unloading holds a significantly larger magnitude of AE energy shortly ahead of rupture than the samples subject to the other regime. This is attributed to the higher stress that the sample experienced under the regime 2, the higher stress will result in more stored energy at contact of two particles. Therefore, the amount of released energy is increased accordingly. Apart from the effect of stress level, the value of bond radius multiplier is another factor to impact the released energy. For regime 2, the stiffness hardening is observed for all 3 samples at earlier phase, the increased bond radius multiplier (larger than 1) is used instead of remaining constant or gradually reduced. The increase in bond radius multiplier can make the assembly survive under a higher stress via enhanced elastic moduli (Potyondy 2007), the stiffened assembly shows a good performance in energy storage and displays a better elasticity.

Spatial distribution of AE events with indication of the magnitudes of scaled AE energy a C5 b C11 c C19 d C6 e C22 f C21

3.2.3.3 Evolution of numerical b-value

We calculated the numerical b-value based on the released energy of AE events of the last five moments shown in Fig. 11, the scaling factor of the released energy is set to 1×10³ and is consistent with the value used in monotonic loading in Sect. 3.1. As documented by Bai et al. (2022b), the b-value can be computed based on moment tensor of AE events, the events that fall within a specific space window and time window (selected timestep) are considered to affect the size of event magnitudes and need to be carefully calibrated. In our study, the sample size for calculation of b-value is 400 AE events. Figure 12 shows the relations of normalized magnitude and the distribution probability. Figure 13a summarized the interrelations between the b-value and the stress level as well as the bond radius multiplier in a 3D diagram. Figure 13b, c respectively plot the impact of stress level and bond radius multiplier on b-value. The results in Fig. 13b indicate that the stress increase ahead of failure (particularly for regime 2) will result in a reduction in b-value, the samples in regime 1 obviously holds a much higher b-value due to the moments selected are far from final rupture. The relation between the bond radius multiplier and b-value revealed in Fig. 13c shows that the reduction of bond multiplier can also incur a smaller b-value. To sum up, the elevation of stress level and reduction of bond radius multiplier can both intensify the damage of particle assembly from insights of b-value. The quantification of damage in numerical model should simultaneously consider the effect of both stress and bond radius multiplier.

Numerical b-values of the last five moments shown Fig. 10

Distribution of b-values versus stress level and bond radius multiplier a 3D distribution b b-value versus stress level c b-value versus bond radius multiplier. The solid lines in Fig. 13b, c are just plotted for the indication of the general distribution pattern of b-value and are not fitting lines

4.1 Discussion on energy dissipation and damping

As revealed in previous DCL tests, the anthracite samples exposed to rapid loading and slow unloading can highly likely incur a lower fatigue life than the samples subjected to a contrary loading pattern (Song et al. 2022d). In addition, more energy was dissipated when the loading rate is larger than the counterpart of unloading. To examine whether this characteristic is well replicated, we present and compare the evolution of discrete energy pillars under two loading regimes, the calculation of the energy of a cycle is based on the summation of discrete energy pillars. Each energy pillar quantifies the energy (positive value for shortened specimen in loading phase and negative value for lengthened specimen in unloading phase) between two consecutive data points (Song et al. 2022d). The results of 15th cycle of C11 and C22 are shown as the example in Fig. 14. The data of the lab test is shown in Fig. 14a (Song et al. 2022d), it is noticed that a 0.36 s phase shift is seen between the peak of stress and cumulative energy pillars. The phase shift is obvious for C11 in lab test, and the counterpart is marginal in C22 which is exposed to an opposite loading pattern. The results plotted in Fig. 14b is from the numerical simulation in this work, the data for the same cycle (15th) is used to make a comparison. The numerical result for C11 shows a 470 timestep (TS) phase shift at the peak stress, and the magnitude of phase shift for C22 is lower. This coincides with the lab tests documented in Fig. 14. Moreover, the significant phase shift in lab tests between stress and cumulative energy pillars can incur excessive discrete positive energy after reaching the peak stress, which is the root reason for the more dissipated energy observed in the first loading regime at the end of the cycle, see the area III. This specific behavior is well replicated in our simulation, see the area III of Fig. 14b. Even the differences of absolute values of energy dissipation between two regimes in lab tests and numerical simulation are not close to each other, the performed simulation still succeeds in replicating the most dominant characteristics observed in lab tests: e.g., phase shift, excessive discrete positive energy, and difference in energy dissipation between two regimes.

Evolution of energy pillars and phase shift between two DCL regimes, the 15th cycle of C11 and C22 is used as an example a Lab test results for C11 and C22, after Song et al. (2022d) b Numerical simulation results for C11 and C22

The damping in LPBM is belong to the dashpot component which can provide viscous behaviors and transmits a dashpot force. The normal F_n^d and shear F_s^d dashpot forces are generated via normal and shear critical-damping ratios, β_n and β_s, see Eqs. (11–12). Where the \(\dot{{\delta }_{\text{n}}}\) and \(\dot{{\delta }_{\text{s}}}\) are relative normal and shear translational velocities. For compact assemblies, non-zero damping is recommended to use to reach the equilibrium and to conduct quasi-static deformation simulations (Itasca Consulting Group Inc. 2014). The low frequency DCL tests makes the loading lie in the range of quasi-static state. Here we discussed the effect of normal and shear damping parameters on energy dissipation under the two DCL regimes. We assign four levels of damping ratios (both β_n and β_s): 0.1, 0.25, 0.5, 0.75 in the first cycle. Figure 15a indicates the results for rapid loading and slow unloading and vice versa for Fig. 15b. It is observed that with increase of damping ratios: the fluctuations of the discrete energy pillars are obviously mitigated in both loading and unloading phases under the two regimes. This indicates that a highly-damped system can smooth the evolution of axial strain and can result in a less scattered discrete energy pillars. Here we count the phase shift between the peak stress and the peak of cumulative energy, it is seen that a greater phase shift is observed when the larger damping ratio is used in both two loading regimes. The effect of damping on dissipated energy is not clearly observed from the values indicated in the zoomed-in view of area II, it seems that the dissipated energy during the first cycle is independent on the damping ratios for the values used in this work.

Evolution of energy pillars and phase shift exposed to different damping ratios (0.1, 0.25, 0.5, 0.75) a Rapid loading and slow unloading b Slow loading and rapid unloading

4.2 Outlooks

The recent numerical simulations on fatigue behaviors of brittle geomaterials mainly focus on the macroscopic mechanical response, such as the replication of hysteresis loops and failure patterns. The characterization and replication of progressive damage in meso-micro-scales under cyclic loading are still rarely considered and reported. To better present the innovative points and major advantages of the work, a comparison is made in a tabular form between our work and previous studies in literatures, see Table 3. Furthermore, we also discuss the limitations of current model and the outlooks and future research directions are highlighted to provide a research framework for readers, see Table 4.

A Weibull-distribution-based heterogeneous 3D particle model is employed to simulate the mechanical and micro-seismic responses of coal samples exposed to multi-level differential cyclic loading (DCL). A microscale damage model is implemented by considering the time-stress-dependent evolution of radius multiplier, which incorporates the characterization of stiffness hardening and softening. The proposed model strives to replicate the major results obtained in previous DCL test (Song et al. 2022d) and numerically show its AE behaviors after model validation and calibration with the monotonic loading and AE data. The following conclusions can be drawn:

1. (1)

   The Weibull heterogeneous model exhibits more realistic behaviors in terms of AE count as well as the ratio of stress level to UCS when the first located AE count occurs. The shape parameter mis set to 2.5. The employed numerical model highly replicates the axial strain evolution under the DCL regime. The damage is impacted by bond radius multiplier and stress level.

2. (2)

   The numerical b-value is calculated via the released bond energy in the numerical model, the b-value is reduced with stress elevation and reduction of bond radius multiplier shortly preceding the failure. Two different mechanisms: stiffness hardening and softening coexist in regime of slow loading and rapid unloading. Two loading regimes have totally different damage mechanisms, this also results in the distinct evolution of AE energy in numerical models.

3. (3)

   The characteristics of energy dissipation within a cycle under two regimes are simulated, the phase shift and larger energy dissipated in regime 1 is well replicated. The effect of damping on energy dissipation and phase shift between stress and cumulation of energy is discussed. The highly damped system holds a larger phase shift as well as the less scattered energy variation.