Advances in Civil Engineering
Volume 2021 (2021), Article ID 6697507, 13 pages
https://doi.org/10.1155/2021/6697507

Research Article

Dynamic Behavior of Coal Pillar under Different Load Percentage by Numerical Simulation

Guo-Wei Dong, Hai-Yang Liu, and Guang-an Zhu

School of Energy Engineering, Xi’an University of Science and Technology, Xi’an, Shaanxi 710054, China

Correspondence should be addressed to Guo-Wei Dong; leng285@163.com

Received 15 December 2020; Revised 17 March 2021; Accepted 13 April 2021; Published 27 April 2021

Academic Editor: Hualei Zhang

Copyright © 2021 Guo-Wei Dong et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The stability and dynamic response of coal pillar is of great importance in underground coal mining. In this paper, a series of uniaxial compressive experiments were first carried out to investigate the mechanical properties of coal. Subsequently, a statistical constitutive damage model for coal was proposed and applied to the numerical simulation. The proposed strain damage softening function showed almost the same goodness-of-fit on the experimental curve. According to this investigation, a numerical model FLAC3Dwas created to investigate the dynamic behavior of the coal pillar under different load percentage (LP). Modelling suggests that the incident and transmitted wave stress evolution observes similar rule and its process can be divided into three stages, namely, static preload, dynamic disturbance, and stabilized stages. The effects of dynamic disturbance intensity are also studied at 10 MPa, 20 MPa, and 30 MPa of peak stress, respectively. The results indicate that under the same load percentage, the peak incident and transmitted wave stress increase with the increase of dynamic disturbance intensity. On the contrary, the attenuation decreases. It is also observed that the failure zone interior the coal can be predicted by the wave propagation.

1. Introduction

Rock bursts, which are sudden dynamic disasters induced by deformation and fracture of the coal-rock mass, pose a serious threat to the production and safety of underground coal mining throughout the word [1, 2]. Pillar burst is one of the typical types of rock burst called strain burst, which is caused by violent fracturing of intact rock and, in many cases, occurs in the vicinity of mining openings or pillars [3, 4].

Practice and studies have shown that both the pillar and hard roof will accumulate a large amount of energy, a dynamic response is easily induced in the surrounding coal and rock mass, which will be served under their coupling effect. In the mining areas of Shaanxi, Inner Mongolia and Xinjiang in West China, many rock burst accidents with huge energies have occurred under the conditions of wide coal pillar and hard roof [5, 6].

Coal pillars play a key role in underground coal mines, which are left to provide support to the overlying strata and control surface subsidence, such as in the bord-and-pillar mining method and the strip mining method. The coal pillar strength has been a focus of research for many years. Traditionally, the coal pillar design is considered in terms of the relationship between width () and height (h) of the pillar. Based on databases of actual coal mine pillars, many empirical and theoretical formulas related to the width-to-height ratio were developed by Salamon [7] and Galvin [8]. Numerous authors have tried to quantify coal pillar strength using analytical approaches, numerical models, and field measurements [9–12]. However, the above research studies seldom take into account the material properties in their estimations, so their applications have limitations.

Numerical simulation provides us with an opportunity to include a good number of relevant factors in the analysis so that the results may be more realistic. Furthermore, the effects of single factors can be easily studied numerically, which is helpful to identify the most important factors for use in practice. Previous studies of this kind have proven the validity of the method [13–15]. Jaiswal and Shrivastva [13] proposed an approach for estimation of strain-softening constitutive behavior of coal mass through calibration of a numerical model with field cases, considering various permutation and combination of strain-softening parameters. Wang et al. [14] presented a numerical investigation on the dynamic mechanical state of a coal pillar, and the results show that the intact core of the coal pillar plays a significant role in its loading capacity; that is, the wider the coal pillar, the greater the volume of the elastic core. Poulsen et al. [15] estimated the pillar strength in both dry and saturated conditions by numerical methods and found that the pillar strength with saturation is approximated by weighting the pillar lithological component.

However, the determination of dynamic behavior of coal pillar under different LP is less widely reported in the literature. Therefore, the dynamic stress evolution is a topic of this study. In contrast to laboratory experiments and the empirical method, numerical simulation is capable of capturing the stress state and damage evolution in coal mass and is therefore helpful in clarifying the associated coal failure mechanism under static combined with dynamic condition [16].

This paper presents a study on the coal pillar strength and dynamic behavior of the coal pillar under different load percentage. The arrangement is listed as follows: in Section 2, a series of experimental tests were carried out, and in Section 3, a statistical constitutive damage model for coal was derived. In Sections 4 and 5, the preparation of this numerical study and the comparison between experimental and numerical results will be firstly conducted. Subsequently, dynamic behavior of the coal pillar under different LP will be investigated.

2. Prior Work

Before starting to develop the numerical models, a series of tasks and tests had to be carried out to characterize the coal mass. Experimental samples were taken from coal seam No.5 of the 2130 Coal Mine in Xinjiang Province, China. At the beginning of the experiment, we cut the samples into cuboid specimens. All the specimens were about 50 mm long, 50 mm wide, and 100 mm high. Details of the specimen properties of the tests carried out are shown in Table 1.

Table 1: Results of the uniaxial compressive strength test in coal.

Sample no.
Density/(g/cm³)
Failure load/kN
Compression strength/MPa
Young’s modulus/GPa
Tensile strength/MPa

1
1.37
20.25
7.57
0.62
0.41

2
1.42
22.12
8.42
0.55

3
1.38
16.59
6.22
0.35

Average value
1.39
19.65
7.40
0.5

The overall stress-strain curves obtained from experiments are shown in Figure 1. According to experimental studies, the average density, uniaxial compressive strength, elastic modulus, cohesion, friction angle, and Poisson’s ratio of the coal are 1.39 g/cm³, 7.4 MPa, 0.5 GPa, 1.5 MPa, 45°, and 0.22, respectively.

Figure 1: Full stress-strain curves of the coal specimens.

3. Constitutive Model Analysis

3.1. Setup of Constitutive Model

The coal specimen has obvious heterogeneity, and many defects exist in the interior. There is a big difference among the mechanical properties of the defects, and they are distributed randomly. The mechanical properties of these microunits are assumed to conform to a given Weibull distribution as defined by the following probability density function [17]:

where is the mechanical parameter of the microunit (such as strength or Young’s modulus); the scale parameter is related to the average of the microunit parameters, and the parameter m defines the shape of the distribution function.

According to the damage mechanics theory [18], the following constitutive relation can be obtained:

Assume D, the statistical damage variable, as the ratio of c, the number of damaged microunits to N, the number of all the microunits. The damage variable D, which reflects the degree of internal damaged of the coal, ranges from 0 to 1. Thus, the damage variable of the coal can be expressed as

After substituting equations (3) into (2), we obtain the following equation:

However, in the process of coal mass under compression, the effective area of normal and shear stress is the same after the destruction of the microunit, the damage variable of all the directions is D, and we obtain [19]

where Cn is the damage proportional coefficient, ranges from 0 to 1, which reflects the residual strength of coal:

Substituting equations (5) into (4), we obtain the statistical model for damage in coal under uniaxial compression:

For the softening constitutive model after the peak strength, the load-bearing capacity of the microunits will reduce with increasing plastic strain,

where c₀ and represent the initial cohesion and friction angle of the microunits, respectively.

3.2. Discussion of Model Parameters and Their Physical Meaning

According to equation (8), the strain damage softening factor CD is defined as follows:

The strain damage softening factor CD reflects the damage softening degree of the material, and the physical meanings of the parameters are as follows:

(1)When the parameters F and m are constant, as shown in Figure 2(a), with increase of Cn, CD decreases, i.e., the higher the Cn, the higher the final softening degree of the material is. Therefore, the parameter Cn reflects the residual strength characteristics of coal and rock mass.(2)When the parameters F and Cn are constant, as shown in Figure 2(b), the higher the m, the steeper the strain damage softening coefficient curve is, which indicates the higher the softening rate, the higher the brittleness. From the perspective of rock burst tendency, the faster the postpeak softening rate is, the shorter the time required for dynamic failure of coal and rock mass, that is, the stronger the rock burst tendency is. Therefore, the parameter Cn reflects the brittleness and rock burst tendency of coal and rock mass.(3)When the parameters m and Cn are constant, as shown in Figure 2(c), with the increase of F, the decreasing section of strain damage softening coefficient curves shifts to the right part, i.e., the higher the m, the more plastic strain needs to enter the softening stage. Therefore, the parameter F reflects the sensitivity of the material begins to softening.

(a)

(b)

(c)

Figure 2: Variation of the damage strain-softening coefficient CD with the parameters of Cn, m, and F changing (a) F = 0.01, m = 3, (b) F = 0.01, Cn = 0.9, and (c) m = 7, Cn = 0.9.

3.3. Physical Experiment Verification of the Model

It can be seen from equation (7) that the model parameters are determined by the following calculation methods: taking the parameters E, F, m, and Cn as variables, and adopting the optimization method, the objective function is established as follows:

where n is the sample number of experimental data; σi and εi are the stress and strain of the i th experimental data.

Coal samples are collected from Litang, Xuzhuang, Longgu, and Zhangshuanglou coal mines, and Figure 3 shows the comparison results of the above model calculation, literature calculation [20, 21], and experimental results. It can be seen from the figure that the calculation model used in this paper is more reasonable than that in the literature, which can better fit the uniaxial stress-strain curve of coal and rock mass.

(a)

(b)

(c)

(d)

Figure 3: Comparison between experimental curves and theoretical curves: (a) coal sample 1, (b) coal sample 2, (c) coal sample 3, and (d) coal sample 4.

4. Numerical Simulation

Having obtained the strength properties of the coal, we shall now describe the working method followed in analyzing dynamic behavior of a coal pillar under different LP.

4.1. Numerical Model and Failure Criterion

FLAC3D (Fast Lagrangian Analysis of Continua in 3 Dimensions) [22] was employed to create a model of the coal pillar with dimensions of 50 mm × 50 mm × 100 mm, and the element size of pillar is kept constant as 2.5 mm × 2.5 mm × 2.5 mm, as shown in Figure 4. The chosen failure criterion is the strain-softening model in FLAC3D, in which the coal pillar was modelled as a nonlinear strain-softening material with cohesion and friction angle softening as functions of the plastic strain (refer to equation (7)).

Figure 4: A sketch of FLAC3D mesh and the stress-strain curve of the numerical model.

4.2. Boundary Conditions of Numerical Model

In order to obtain the condition of uniaxial compression, reasonable boundary conditions should be set around the numerical model and the displacement of four vertical symmetry planes of the model is restricted in the normal direction. A constant velocity is applied to the top and bottom of the model in the z-direction to generate a vertical loading on this model. In this study, the magnitude of constant velocity has been set as 2.5 × 10⁻⁷ m/s.

4.3. Determination of Relevant Material Parameters

In this study, formula (6) was used to calibrate the material parameter for the model. The final calibration results give a coal uniaxial compression strength (UCS) of 7.4 MPa, and tensile strength is set as 0.41 MPa. Finally, after repeatedly verification and commission, the shape parameter m, the scale parameter F, and the damage proportional coefficient Cn are set as 30, 0.8, and 0.002, respectively.

4.4. Dynamic Disturbance

As regards the combined static and dynamic loading conditions (as defined in Figure 5), the dynamic incident stress pulse pd (t) is applied at the top surface of a model domain after the static boundary stresses are applied. Under the dynamic loading condition, the bottom surface of the domain that has been restrained during the static loading is free from stress. The half-sine stress pulse (as defined in Figure 5) is used because this incident waveform is similar to the half-sine shape. For this waveform, two parameters (including amplitude (pdm) and duration (tdm)) are involved in the numerical simulation in order to analyze dynamic behavior of the coal pillar under different LP.

Figure 5: The waveform of dynamic stress pd (t) with amplitude of pdm and duration of tdm.

According to the seismology theory, the P-wave of seismic source can be simplified as the superposition of sine waves in three directions perpendicular to each other in space. Therefore, the waveform used in the numerical simulation can be replaced by sine wave. At the same time, based on the statistics of experimental and field observation data, different mine tremors energy corresponds to different maximum amplitude of stress wave, and the specific corresponding relationship is shown in Table 2.

Table 2: Relationship between mine tremor energy and peak dynamic load.

No.
Energy/J
Wave velocity/(m/s)
Density/(kg/m³)
PPV/(m/s)
Peak dynamic load/MPa

1
10³∼10⁴
4500
2500
1.00
10

2
10⁴∼10⁵
4500
2500
2.90
30

3
10⁵∼10⁶
4500
2500
5.60
60

4
10⁶∼10⁷
4500
2500
9.40
100

The statistical results show that when rock burst occurs, the energy level of mine tremor is generally 10⁵∼10⁶ J. According to Table 2, the sine wave with the peak value of 30 MPa can be selected as the dynamic load source. Therefore, in the numerical simulation, the vibration frequency of longitudinal stress wave is set as 50 Hz (the corresponding vibration period is 0.02 s), the peak stress is 30 MPa, and the monitoring time is 1.2 s.

4.5. Calibration of Local Damping Coefficient in the Dynamic Analysis

According to the manual for FLAC3D, there were two aspects that the user should consider for a dynamic analysis: (1) boundary conditions and (2) mechanical damping. Here, the quiet boundary, developed by Lysmer and Kuhlemeyer [23], was chosen to reduce the reflection of propagating waves. The model is almost entirely effective at absorbing body waves approaching the boundary at angles of incidence greater than 30° [22]. Damping is due, in part, to energy loss as a result of internal friction in the intact material or slippage along interfaces. Regarding the attenuation of elastic waves, local damping was considered in the dynamic analysis. In the numerical modelling of FLAC3D, the local damping coefficient can be expressed as

where D is the critical damping coefficient.

The local damping coefficient can be calibrated in the numerical model with respect to calculated and measured PPVs [24, 25]. To determine the local damping parameters through this numerical analysis, a three-dimension model, including coal, roof, and floor, was developed for the calculation. The model had the following dimensions: 40 m × 66 m × 40 m (length × width × height), as shown in Figure 6. In this section, the numerical model was run five times under conditions commensurate with different critical damping coefficients, i.e., the critical damping coefficient was set to 1%, 2%, 3%, 4%, and 5%, respectively. Afterwards, the PPVs at these monitoring points were recorded, and the points A and B are set 2 m and 12 m from the numerical model’s centre, respectively. Based on the theoretical and numerical PPVs, the damping coefficient can be obtained and the calculated results are depicted in Figure 6. It can be seen that, with the increase of the critical damping coefficient, the PPV at point A increases: the resulting polynomial fit can be expressed as follows:

(a)

(b)

(c)

Figure 6: Numerical simulation for dynamic analysis of local damping. (a) Numerical model geometry. (b) Sectional panel of the calibrated model (z = 0). (c) Sectional panel of the calibrated model (y = 0).

While the ratio of PPV at point A to that at point B decreases, the resulting polynomial fit can be expressed as follows:

It can be seen from Figure 7 that the ratio of the PPV at point A to that at point B was 3.694: substituting y = 3.694 into equation (13) gives x = 4.468, which means that the critical damping coefficient in this dynamic analysis is 4.468%.

Figure 7: Graphs for calibration with respect to dynamic stress and local damping coefﬁcient.

To verify the calibrated damping coefficient, an additional run was performed, with the PPV and error values at different distances depicted in Figure 8. It can be seen that the PPVs recorded 2 m and 12 m from the centre of the equivalent cavity agree well with those calculated form theoretical calculation, and the error can be limited to less than 12.7%. The calibrated model thus represents real effects on the coal and rock mass, and the damping coefficient and modelling scheme are adopted in the following simulations.

Figure 8: Error analysis of PPV with increasing distance.

5. Numerical Simulation Results and Discussion

5.1. Comparison of Experiment and Numerical Simulation

The comparison of the stress-strain curves obtained by the numerical model and that of the laboratory test are shown in Figure 4. Readers may find that the proposed numerical curve is not well fitted to the experimental curves of stress-strain, especially in the initial stage. It is because that the numerical calculation cannot effectively simulate the initial stage of the curve, and a large number of experiments show that this stage is called the crack compression stage. The numerical simulation does not take into account the naturally existing cracks and the mechanical properties of the model. Therefore, eliminating the influence of the initial stage, i.e., translating the simulation curve, the “numerical-correction” curve can be obtained (red dotted line in Figure 9). Then, it can be seen that the numerical simulation stress-strain curve of the model matches the laboratory test data very well.

Figure 9: Stress vs. strain experimental and numerical curves for coal specimens.

5.2. Dynamic Behavior

The dynamic module was used to investigate the response features of the coal pillar under dynamic disturbance. In this section, 6 monitoring points were set up in the interior of the model (x = 0.025 m, y = 0.025 m, and z = 0-0.1 m), as shown in Figure 10.

Figure 10: Cross-section of the model, in which 6 typical points where the stresses will be retrieved.

5.2.1. Effect of Load Percentage

If the dynamic disturbance intensity is constant (10 MPa), the evolution of the incident and transmitted wave stress with different LP is shown in Figure 11. In the entire process (take Figure 11(d) as an example), the stress wave evolution for various LP follows a similar rule and the time history curve of stress wave variation can be divided into three stages: namely, static preload stage, dynamic disturbance stage, and stabilized stage.

(1)Static preload stage: at the beginning, the stresses on the two ends of the model are different until it reaches the uniaxial compression strength (about 7 MPa). The reason is that before the dynamic disturbance, the simulation model is subject to the static preload through the boundary, which is already compressed, and thus, the stress gauge interior the model measured the preload.(2)Dynamic disturbance stage: this stage lasts for ∼0.2 s, which is the main dynamic disturbance period. The wave stress reached a peak value at 0.1 s at this stage.(3)Stabilized stage: when the calculation time exceeds 0.2 s, the wave stress is maintained at a low level and the models reach equilibrium.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 11: Waveform of the incident and transmission for the model with different LP. From Figures 11(a)–11(j), the LP is from 10% to 100%, respectively. (a) C, (b) LP = 20%, (c) LP = 30%, (d) LP = 40%, (e) LP = 50%, (f) LP = 60%, (g) LP = 70%, (h) LP = 80%, (i) LP = 90%, and (j) LP = 100%.

Another interesting phenomenon found in from Figure 11 is that the duration of the first stage gets shorter and shorter. With an increase in load percentage, the static preload in the coal increases. Consequently, the duration of the static preload stage decreases.

The relationship between wave attenuation and LP is shown in Figure 12. Statistics show that for models with different LP, the curves present their periodic changes. When LP ≤ 30%, for an increase in LP, the attenuation coefficient decreases almost linearly. The reason for this phenomenon is that the cracks in the coal models close; thus, the wave attenuation decreases sharply. When LP = 40%–70%, the attenuation coefficient is almost constant (about 50%). When LP ≥ 80%, compression-induced cracks form again and parts of microunits start to damage. Therefore, the increase in LP makes the concentration more obvious and leads to further attenuation.

Figure 12: Attenuation coefficient versus time curve for different load percentage.

The stress wave evolution recorded at the monitoring points for a constant LP (40%) is presented in Figure 8. In general, the attenuation of stress wave in coal and rock mass shows a power relation with the propagation distance; that is, the longer the distance is, the more evident the attenuation is at a certain distance [26]. However, due to various LP, the wave propagation is no longer uniform. By comparing the recorded data, we can see that the recorded peak vertical stress of monitoring points C and D are larger than the ones of points B and E. This means dynamic disturbance is obvious as the monitoring points gets close to the two ends of the model, where the load applied to the microunit is much higher and leads to varying degrees of damage. The results indicate that the failure zone interior the coal can be predicted by the wave propagation.

5.2.2. Effect of Dynamic Disturbance Intensity

The change in vertical stress for different dynamic disturbance intensity is shown in Figure 13. It can be seen that the dynamic disturbance intensity has a large influence on the evolution of stress. At the specific value of load percentage (40%), with the increase of the amplitude of dynamic disturbance intensity, the recorded stresses vary greatly. To better understand these relationships, the detailed data are plotted in Figure 14. With increasing dynamic disturbance intensity, the peak incident wave stress increases from 2.33 MPa to 4.02 MPa and the peak transmitted wave stress increases from 1.28 MPa to 3.81 MPa. However, the attenuation of the incident wave stress varies slightly, which trends to be zero when pdm ≥ 30 MPa. (Figure 15)

(a)

(b)

(c)

Figure 13: Vertical stress versus time curves for different dynamic disturbance intensity (LD = 0.4): (a) 10 MPa; (b) 20 MPa; (c) 30 MPa.

Figure 14: Effect of dynamic disturbance intensity on the peak value of incident and transmitted wave stress.

Figure 15: Curves of vertical stress varying with time of monitoring points (LP = 40%).

6. Conclusions

The focus of this study is to investigate dynamic behavior of the coal pillar under different load percentage. Based on laboratory coal strength data, a statistical constitutive damage model for coal under uniaxial compression was established. Numerical models considering the experimental result were generated by FLAC3D software. The following main conclusions can be drawn.

The heterogeneous strain damage softening constitutive model is established, and the numerical test of heterogeneous coal and rock material is carried out. It is revealed that the heterogeneity is the source of coal and rock failure precursor. The observed constitutive model can reflect the elasticity, plasticity, brittleness, and residual strength of coal and rock materials, which fits well with the whole process of stress-strain curve of coal and rock in physical experiment.

For the entire process, the stress wave evolution for various load percentage follows a similar rule and the time history curve of stress wave variation can be divided into three stages: namely, static preload, dynamic disturbance, and stabilized stage.

When LP = 0–30%, the attenuation coefficient decreases almost linearly and sharply. While LP = 40%–70%, the attenuation coefficient is almost constant. If LP ≥ 80%, compression-induced cracks form again and parts of microunits start to damage, which leads to further attenuation.

The dynamic disturbance intensity has a significant influence on the dynamic behavior. With increasing dynamic disturbance intensity, the peak incident and transmitted wave stress increase whilst the attenuation decreases.

Data Availability

All the data supporting the research in this paper are shown in the figures and tables.

Conflicts of Interest

The authors declare that there are no conflicts of interests regarding the publication of this paper.

Authors’ Contributions

Dong, G. and Zhu, G. conceived and designed the experiments. Liu, H. performed the experiments and numerical simulation. Dong, G. and Zhu, G. analyzed the data. Dong, G. and Zhu, G. wrote the paper.

Acknowledgments

This work was supported and financed by the Natural Science Foundation of Shaanxi Province (grant nos. 2020JM-519 and 2019JQ-487), China Postdoctoral Science Foundation (grant no. 2018M643692), and National Natural Science Foundation of China (grant nos. 51904235 and 51974231).

References

L. M. Dou, Z. L. MU, Z. L. LI, A. Y. Cao, and S. Y. Gango, “Research progress of monitoring, forecasting, and prevention of rockburst in underground coal mining in China,” International Journal of Coal Science and Technology, vol. 1, no. 3, pp. 278–288, 2014.

C. Zhang, I. Canbult, B. Hebblewhite, and C. R. Ward, “Assessing coal burst phenomena in mining and insights into directions for future research,” International Journal of Coal Geology, vol. 179, pp. 28–24, 2017.

S.-Y. Wang, K.-C. Lam, S.-K. Au, C. A. Tang, W. C. Zhu, and T. H. Yang, “Analytical and numerical study on the pillar rockbursts mechanism,” Mechanics and Rock Engineering, vol. 39, no. 3, pp. 445–467, 2006.

G.-A. Zhu, L.-M. Dou, A.-Y. Cao, C.-B. Wang, Z.-G. Liu, and J. Li, “Assessment and analysis of strata movement with special reference torock burst mechanism in island longwall panel,” Journal of Central South University, vol. 24, no. 12, pp. 2951–2960, 2017.

B. Wang, F. X. Jiang, S. T. Zhu et al., “Mechanism and prevention of rock burst induced by segment pillars in the deep mining areas of shaanxi-inner Mongolia adjacent regions,” Journal of Mining & Safety Engineering, vol. 37, no. 3, pp. 505–513, 2020, in Chinese.

T. Wang, S. You, F. Pei et al., “Instability mechanism and control technology of coal pillar bumps under hard roof,” Journal of Mining & Safety Engineering, vol. 34, no. 1, pp. 54–59, 2017.

M. D. G. Salamon, “Stability, instability and design of pillar workings,” International Journal of Rock Mechanics & Mining Sciences, vol. 7, no. 6, pp. 613–631, 1970.

J. M. Galvin, “Considerations associated with the application of the UNSW and other pillar design formulae,” in Proceedings of the 41st US Symposium on Rock Mechanics, pp. 1129–1137, Golden, CO, USA, June 2006.

L. Yuan, Y. D. Jiang, X. Q. He et al., “Research progress of precise riskaccurate identification and monitoring early warning on typical dynamic disasters in coal mine,” Journal of China Coal Society, vol. 43, no. 2, pp. 306–318, 2018.

B. Robert, D. Kurt, and M. Garry, “An approach to model the strength of coal pillars,” International Journal of Rock Mechanics & Mining Sciences, vol. 89, pp. 165–175, 2016.

G.-F. Wang, S.-Y. Gong, Z.-L. Li, L.-M. Dou, W. Cai, and Y. Mao, “Evolution of stress concentration and energy release before rock bursts: two case studies from Xingan coal mine, Hegang, China,” Rock Mechanics and Rock Engineering, vol. 49, no. 8, pp. 3393–3401, 2016.

S. Sankhaneel and W. Gabriel, “Investigation of longwallheadgate stress distribution with an emphasis on pillar behavior,” International Journal of Rock Mechanics & Mining Sciences, vol. 121, pp. 104–409, 2019.

A. Jaiswal and B. K. Shrivastva, “Numerical simulation of coal pillar strength,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 4, pp. 779–788, 2009.

H. Wang, Y. Jiang, Y. Zhao, J. Zhu, and S. Liu, “Numerical investigation of the dynamic mechanical state of a coal pillar during longwall mining panel extraction,” Rock Mechanics and Rock Engineering, vol. 46, no. 5, pp. 1211–1221, 2013.

B. A. Poulsen, B. Shen, D. J. Williams, C. Huddlestone-Holmes, N. Erarslan, and J. Qin, “Strength reduction on saturation of coal and coal measures rocks with implications for coal pillar strength,” International Journal of Rock Mechanics and Mining Sciences, vol. 71, pp. 41–52, 2014.

S.-L. Wang, S.-P. Hao, Y. Chen, J.-B. Bai, X.-Y. Wang, and Y. Xu, “Numerical investigation of coal pillar failure under simultaneous static and dynamic loading,” International Journal of Rock Mechanics and Mining Sciences, vol. 84, pp. 59–68, 2016.

S. C. Yuan and J. P. Harrison, “A review of the state of the art in modelling progressive mechanical breakdown and associated fluid flow in intact heterogeneous rocks,” International Journal of Rock Mechanics and Mining Sciences, vol. 43, no. 7, pp. 1001–1022, 2006.

Z. X. Li, Damage Mechanics and its Application, Science Press, Beijing, China, 2002.

W. Y. Xu and L. D. Wei, “Study on statistical damage constitutive model of rock,” Chinese Journal of Rock Mechanics and Engineering, vol. 21, no. 6, pp. 787–791, 2002.

Z. Wu and C. J. Zhang, “Investigation of rock damage model and its mechanical behaviour,” Chinese Journal of Rock Mechanics and Engineering, vol. 15, no. 1, pp. 55–61, 1996.

S. Q. Yang, W. Y. Xu, L. D. De et al., “Statistical constitutive model for rock damage under uniaxial compression and its experimental study,” Journal of Hohai University (Natural Sciences), vol. 32, no. 2, pp. 200–203, 2004.

Itasca Consulting Group Inc, FLAC3D Fast Lagrangian Analysis of Continua in 3 Dimensions User’s Guide, Itasca ConsultingGroup Inc, Minneapolis, MN, USA, 2012.

J. Lysmer and R. L. Kuhlemeyer, “Finite dynamic model for infinite media,” Journal of the Engineering Mechanics Division, vol. 95, no. 4, pp. 859–877, 1969.

A. Sainoki and H. S. Mitri, “Dynamic modelling of fault slip induced by stress waves due to stope production blasts,” Rock Mechanics and Rock Engineering, vol. 49, no. 1, pp. 165–181, 2016.

G. A. Zhu, L. M. Dou, C. B. Wang, J. Li, W. Cai, and Z.-W. Ding, “Numerical investigation of the evolution of overlying strata and distribution of static and dynamic loads in a deep island coal panel,” Arabian Journal of Geosciences, vol. 10, no. 24, p. 549, 2017.

B. Liu, C. P. Lu, L. M. Dou et al., “Simulation study on shock wave propagation character in coal and rock,” Journal of China Coal Society, vol. 36, no. 2, pp. 247–253, 2011.