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 [912]. 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 [1315]. 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.

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/cm3, 7.4 MPa, 0.5 GPa, 1.5 MPa, 45°, and 0.22, respectively.

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

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.

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

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

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.

The statistical results show that when rock burst occurs, the energy level of mine tremor is generally 105∼106 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 aswhere 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:

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

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.

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.

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.

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.

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.

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)

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.


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