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

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.


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 serve 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 strataand control surface subsidence,such as in thebord and pillar mining method, the strip mining method, etc.The coal pillar strength has been a focus of research for many years. Traditionally, the coal pillars design is considered in terms of the relationship between width (w) and height (h) of the pillar.Based on databases of actual coal mine pillars, many Empirical and theoreticalformulas related to the width-to-height ratio were developed by Salamon [7] and Galvin [8].Numerous authors have tried toquantify coal pillar strength using analytical approaches, numerical models and field measurements [9][10][11][12].However, the above researchesseldom 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, whichis helpful to identify the most important factors for use in practice.Previous studies of this kind have proven the validity of the method [13][14][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, which 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, the results show that the intact core of the coal pillar plays a significant role in its loading capacity, i.e., 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 LD is less widely reported in the literature. Therefore, the dynamic stress evolution is a topic of this study. In contrast to laboratory experiments and 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 was 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 LD will be investigated.

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

Set up of constitutive model
The coalspecimen 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 micro-unit 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 micro-unit (such as strength or Young's modulus); the scale parameter F is related to the average of the micro-unit 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 micro-units to N, the number of all the micro-units. 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 Eq. (3) into Eq. (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 micro-unit, the damage variable of all the directions is D, and obtain [19]   where Cn is the damage proportional coefficient, ranges from 0 to 1, which reflects the residual strength of coal.
Substituting Eq.(5) into Eq.(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 micro-units will reduce with increasing plastic strain, wherec0 and 0  represent the initial cohesion and friction angle of the micro-units, respectively.

3.2Discussion of model parameters and their physical meaning
According to Eq. (8), the strain damage softening factor CDis defined as follows: Thestrain damage softening factor CDreflects the damage softening degree of the material, and the physical meanings of the parameters are as follows: (1) when the parametersF and m are constant, as shown in Fig. 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 Fig. 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 post peak 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 brittlenessand rock burst tendency of coal and rock mass.
(3) when the parameters m and Cn are constant, as shown in Fig. 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.3Physical experiment verification of the model
It can be seen from Eq. (7) that the model parameters are determined by the following calculation methods:taking the parameters E, F, mand Cn as variables and adopting the optimization method, the objective function is established as follows wherenis the sample number of experimental data;σiandεiare the stress and strain of the ith experimental data.
Coal samples are collected from Litang, Xuzhuang, Longgu, Zhangshuanglou coal mines, and Fig. 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.

Boundary conditions of numerical model
In order to obtain the condition of uniaxial compression, reasonable boundary conditions should be set around

Determination of relevant material parameters
In this study, formula (6)

Dynamic disturbance
As regards the combined static and dynamic loading conditions (as defined in Fig. 5.), the dynamic incident stress pulse pd(t), is applied at the top surface of a model domain after the static boundary stressesis 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 Fig. 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 LD.
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 10 5~ 10 6 J. According to the Table 2, the sine wave with the peak value of 30MPa 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.5Calibration of local damping coefficient in the dynamic analysis
According to the manual for FLAC 3D , 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 FLAC 3D , the local damping coefficient L  can be expressed as: whereD is the critical damping coefficient.  Fig.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: y=-0.0036x 2 -0.5316x+13.316 (12) While the ratio of PPV at point A to that at point B decreases,and the resulting polynomial fit can be expressed as follows: y= 0.028x 2 +0.1355x+2.5296 (13) It can be seen from Fig. 7  To verify the calibrated damping coefficient, an additional run was performed, with the PPV and error valuesat different distances depicted in Fig. 8. It can be seen that the PPVs recorded 2 m and 12 m from the center of the equivalent cavity agree well with those calculated form theoretical calculation, 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.

Comparison of experiment and numerical simulation
The comparison of the stress-strain curves obtained by the numerical model and that of the laboratory test is shown in Fig. 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 can't effectively simulate the initial stage of the curve, a large number of experiments show that this stage is called the crack compression stage.
The numerical simulation doesn't take into account thenaturally existing cracksand 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 Fig. 9).Then, it can be seen that the numerical simulation stress-strain curve of the model matches the laboratory test data very well. Fig. 9.Stress vs strain experimental and numerical curves for coal specimens.

Effect of load percentage
If the dynamic disturbance intensity is constant (10MPa), the evolution of the incident and transmitted wave stress with different LD is shown in Fig. 11.In the entire process (take Fig. 11.(d) as an example), the stress wave evolution for various LD 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, stabilized stage. .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 Fig. 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,theduration of the static preload stage decreases.
The relationship betweenwave attenuation andLD is shown inFig. 12. Statistic show that for models with different LD, the curves present their periodic changes. When LD≤30%,for an increase in LD, 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 LD=40%~70%, the attenuation coefficient is almost constant (about 50%). When LD≥80%, compression-induced cracks form again and parts of micro-units start to damage.
Therefore, the increase in LD makes the concentration more obvious and leads to further attenuation. The stress wave evolution recorded at the monitoring points for a constant LD (40%) is presented in Fig. 8.In general, the attenuation of stress wave in coal and rock mass shows a power relation with the propagation distance, i.e., the longer the distance is, themore evident the attenuation isat a certain distance [26]. However, due to various LD, the wave propagation is no longer uniform. By comparing the recorded data, we can see thatthe recordedpeak vertical stressof 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 micro-unit 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.

Effect of dynamic disturbance intensity
The change in vertical stress for different dynamic disturbance intensity is shown in Fig. 14. 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 varies greatly. To better understand these relationships, the detailed data are plotted in Fig. 15.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.28MPa to 3.81MPa. However, the attenuation of the incident wave stress varies slightly, which trends to be zero when pdm≥30MPa.

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 experimental result were generated by FLAC 3D software. The following main conclusions can be drawn: (1) 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.
(2) When LD =0~30%,the attenuation coefficient decreases almost linearly and sharply. While LD=40%~70%, the attenuation coefficient is almost constant. If LD≥80%, compression-induced cracks form again and parts of micro-units start to damage, which leads to further attenuation.
(3) With increasing dynamic disturbance intensity, the peak incident and transmitted wave stress increase

Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.