Inverse analysis of dynamic failure characteristics of roadway surrounding rock under rock burst

Rock burst is one of the most serious dynamic disasters in the process of coal resources mining, in order to study its occurrence mechanism and the characteristics of its impact on roadways. Based on the No. 25110 working face of Yuejin Mine in Yi Coal Mine and No. 21032 return air uphill event in Qianqiu Mine, combined with geological conditions, shock location, and source waveform, this paper uses numerical simulation to inverse the dynamic response characteristics of roadway surrounding rock under rock burst. The results show that the plastic strain energy density of No. 25110 working face in Yuejin Mine is mainly distributed around the roadway, and the maximum plastic strain energy density is 3.30 × 107 J at 75 ms, all rock masses around the roadway are damaged, and the maximum damage variable value is 0.44. The plastic strain energy density of No. 21032 return air uphill in Qianqiu mine is mainly distributed around the lower side of the roadway, and the maximum plastic strain energy density is 2.68 × 108 J at 80 ms, all rock masses around the roadway are damaged, and the maximum damage variable value is 1.76, indicating that the rock masses at the roof and two sides of the roadway have been damaged after the impact. Based on the results of numerical simulation, the energy gradient criterion of rock burst tendency was proposed, the degree of energy accumulation of surrounding rock is measured by the change amplitude of energy gradient, and the tendency index of rock burst is quantified, which provides a new method for the prediction of rock burst.

impacts leads to roadway deformation and personnel injury, which seriously threatens the safe and efficient mining of mines. 1,2 At present, there are 133 rock burst mines in China, 3 of which as many as 50 mines have a mining depth of more than one kilometer.
Many scholars at home and abroad have conducted a lot of research on the mechanism and prevention of mine rock burst disasters. By means of acoustic emission and on-site microseismic monitoring in the laboratory, the precursor identification of coal and rock sudden failure process is carried out, and a comprehensive early warning model of rock burst based on microseismic precursor index system is proposed. [4][5][6] The coal rock mechanics and impact characteristics experiments under triaxial loading and unloading path are carried out, the dynamic characteristics of triaxial coal under impact load are studied, 7 the energy evolution process of coal and gas outburst is analyzed, and the relationship index between acoustic emission parameter characteristics and coal and gas outburst precursor information is established. 8,9 Xia et al 10 established the static comprehensive evaluation index and method of impact risk in view of the critical mutation and difficulty in weight quantification in the impact risk evaluation by the comprehensive index method and multifactor coupling method. When studying the mechanism of rock burst and coal and gas outburst, Qian et al considered that the four types of typical rock burst in coal mine are mainly coal seam material instability type, coal seam structure instability type, roof fracture type, and fault slip dislocation type. 11 Li et al 12 used high-precision microseismic monitoring system and roof dynamic monitor to study the microseismic activity characteristics of coal seam floor and the influence of dynamic pressure on roadway deformation and failure. The laboratory used the drop hammer impact test device to carry out the mechanical response test of common pallets and their composite components in coal mines and revealed the impact resistance mechanical properties of bolt surface protection components. 13 Li et al [14][15][16] analyzed the rock burst mechanism and microseismic (MS) and electromagnetic radiation (EMR) responses in coal and rock with structural planes. Based on the particularity of graben structure and the occurrence characteristics of overburden, Wu et al 17 analyzed the characteristics of overburden movement and the stress evolution law of coal and rock mass in the graben structure area, established the corresponding mechanical model, and studied the mechanism of rock burst in the graben structure area. Based on the rock burst experiment system, the fractal characteristics of sandstone rock burst debris under dynamic and static combined load are analyzed. 18 21,22 To sum up, it can be seen that during the excavation of the rock mass, due to the sudden release of the surrounding rock high stress and the combined action of the dynamic disturbance load, the mine dynamic disaster has complex dynamic characteristics. In order to study its occurrence mechanism and the characteristics of its influence on the roadway, based on the No. 25110 working face of Yuejin Mine and No. 21032 return wind uphill event of Qianqiu Mine, combined with geological conditions, shock location, and source waveform, this paper uses numerical simulation to invert the dynamic response characteristics of roadway surrounding rock under the action of rock burst. Based on the analysis of the energy density distribution characteristics of the surrounding rock, an energy gradient criterion for the tendency of rock burst is proposed, in order to provide theoretical guidance for the control of the surrounding rock of the roadway under dynamic load.

CHARACTERISTICS OF THE ROADWAY UNDER ROCK BURST IN YUEJIN MINE
ANSYS/LS-DYNA finite element analysis was used to simulate the dynamic response of the surrounding rock under the influence of ground pressure signals, such as stress state, plastic strain energy density, damage distribution, and other parameters.

| Geological conditions
The depth of No. 25110 working face in Yuejin Coal mine, the Yimei mine area, is about 1041 m, and the coal seam of the working face tends to be 191 m long with a strike length of 1000 m. The upper roadway was excavated along with the air, and the support of "anchor net cable +Ishaped steel ladder-shaped shed" was adopted. The lower roadway was arranged in the solid coal seam, and the roof of the coal seam was excavated to adopt a combined support method of "anchor net cable +36 U oval shed." The original section of the roadway is a straight wall arc arch, with a clear height of 4.0 m and a width of 6.2 m. The coal seam shows an inclination angle of 10° to 15° and a net cross-sectional area of 24 m 2 . The coal seam thickness is 8.4 m-13.2 m, with an average thickness of 11.5 m. The pseudo-top is sandy mudstone with a thickness of 0.2 m, layered structure, easy to fall off, and quartz sandstone is hard in part. The direct top is fragile mudstone with a thickness of about 20 m. The old roof is mainly conglomerate with thin sandstone, which belongs to the weak aquifer containing water-rich abnormal areas. The direct bottom is fragile mudstone with a thickness of about 4 m.
The impact tendency of the coal seam in No. 25110 fully mechanized caving face is strong impact, and the roof is weak impact tendency.

| Calculation model
After excavating the tunnel, the source is the velocitytime history curve collected by the microseismic monitoring system when the impact ground pressure occurs on the spot. The velocity-time curve is shown in Figure 1, adopted from HJC constitutive model. Table 1 shows the HJC constitutive model's parameters.
Holmquist Johnson cook (HJC) 23 model is a constitutive model for concrete materials, which comprehensively considers the effects of strain rate effect, damage evolution effect, confining pressure effect, crushing, and compaction effect. The HJC yield surface equation is used to describe the equivalent strength of its material. The expression is: where * = ∕f c and P * = P∕f c are dimensionless equivalent stress and hydrostatic pressure, respectively; f c is the uniaxial compressive strength of the quasi-static cylinder of the material; P is the actual hydrostatic pressure; ε * is the ratio of the actual strain rate ε to the reference strain rate ε 0 =1.0 s −1 . A, B, N, and S max are the strength parameters of the material model: A is the characteristic viscosity strength coefficient; B and N are the characteristic pressure hardening coefficient and pressure hardening index, respectively; D is the damage variable; C is the strain rate effect coefficient.
The simulated seismic source is located 20 m away from the roadway floor. The roadway section is a straight-wall arc arch with a width of 6.2 m and a net height of 4 m. A Note: ρ 0 is the density; G is shear modulus; K 1 , K 2 , K 3 are pressure constants; μ lock is the volumetric strain during compaction; ε fmin is the minimum plastic strain at fracture; T is uniaxial tensile strength; D 1 and D2 are damage parameters; P Crush is the hydrostatic pressure at the elastic limit; μ Crush is the volume strain at the elastic limit; P lock is the compaction hydrostatic pressure.
T A B L E 1 Material parameters of HJC constitutive model combined support mode of "anchor net cable +36 U oval shed" was adopted as shown in Figure 2. Figure 3 shows the stress nephogram in the Y direction around the roadway at different times. As shown in Figure 3, the stress in the Y direction around the roadway gradually increases with the continuation of the disturbance time. The monitoring points were arranged on the roof, bottom, and two sides of the roadway. Figure 4 presents the curve of stress change for each measuring point with time.

| Result analysis
As shown in Figure 4, the Y-direction stress of the roof measuring point is fluctuating with a rapid increase initially and a decrease subsequently. At 4 ms, the maximum Y-direction stress of the roof reaches 96 MPa, and the roadway returns to steady-state stress with a value of 40 MPa at 75 ms The Y-direction stress changes from tensile stress to compressive stress while the time is gradually increased to reach a stable state with a maximum value of 250 MPa, the Y-direction stress in a stable state at 10-75 ms It can be seen that the bottom rock mass experiences a high-stress state, indicating a more significant impact on the bottom than the top. This paper only analyzes the stress of the left side since the stress changes of the left and the right sides are the same. At 5 ms, the maximum tensile stress of the left side is 60 MPa. However, at 75 ms, the stress value of the left side drops to 0, and the left rock mass is damaged. Figure 5 presents a nephogram of the plastic strain energy density vs. different times around the roadway.
From Figure 5, at 20 ms, the plastic strain energy density is mainly concentrated at the arch toe of the roadway. It is extended to the deep part of the rock mass with a maximum value of 5.34 × 10 6 J. At 40 ms, the plastic strain energy density is distributed at the arch toe and the two sides of the roadway with a maximum value of 2.53 × 10 7 J. At this time, the anchor rod at the arch toe of the roadway was damaged and deformed. It caused to shrink the roadway section. At 60 ms, the plastic strain energy density is mainly distributed at the arch toe and floor of the roadway with the maximum value of 6.99 × 10 7 J. In this case, the anchor rod was deformed at the arch toe of the roadway. At 75 ms, the plastic strain energy density is mainly distributed around the roadway with a maximum value of 1.73 × 10 8 J. At this time, the anchor rod at the arch foot of the roadway was seriously deformed, and the roadway section was shrunk about 0.75 times of the initial section area. Figure 6 presents the maximum value of plastic strain energy density vs. time by applying Figure 5.  As shown in Figure 7, the rock mass around the roadway was damaged after applying the dynamic load. At 20 ms, the roof and floor of the roadway were damaged with a damage variation of 0.098. At 40 ms, the left and right arches of the roadway were damaged, with the maximum damage value of 0.36. At 60 ms, damage was also appeared on the two sides of the roadway floor and arch corners and extended to the outside, with the maximum damage value of 0.63. At 75 ms, all damages occurred around the roadway, and the maximum damage value was 1.42, indicating that the rock masses were broken at the left and right arch feet. Figure 8 shows the curve of damage parameter D and time based on the results from The 21-mining area was identified as medium to medium-strong rock burst zone.

| Calculation model
The simulation analysis was based on applying the seismic sources after excavating the roadway. The microseismic monitoring system collected the impact source for the rate time history curve once the rock burst occurred on The simulated seismic source is located 10 m away from the roadway floor. The roadway section is a straightwall arc arch, with a width of 5 m and a net height of 3 m. The supportive measures have been adopted by "spray anchor net cable +5154 type 36 U steel shed." Table 2 presents the model calculation parameters. Figure 10 presents the nephogram of Y-direction stress around the roadway at different times.

| Analysis results
The results show that the Y-direction stress gradually increases around the roadway. Monitoring points were arranged on the roof, bottom, and two sides of the roadway. The characterizations of the stress variation for each measuring point with time are plotted in Figure 11. As shown in Figure 11, the Y-direction stress of the roof first decreases in 80 ms and then increases. Afterward, the stress in the Y direction decreases and returns to a steady state. At 65 ms, the stress in the Y direction of the roof changes from compression to tension. At 75 ms, the stress in the Y direction of the roof is 20 MPa when the roadway returns to a static state. Also, the Y-direction stress of the floor first decreases, then increases, and finally tends to be stable. At 0 to 4 ms, the floor rock mass bears both compression and tensile stresses. At 80 ms, the maximum stress in the Y direction of the floor is 300 MPa. At 10-75 ms, the stress trend in the Y direction of the floor is unchanged. It can be seen that the floor rock mass is in a high-stress state, indicating the higher impact of the floor roof. The Y-direction stress of the left side measuring point increases with increasing time, while the Y-direction stress of the right side first increases and then decreases to 0, indicating that the right side is damaged before the left side. Figure 12 shows the nephogram of plastic strain energy density change around the roadway at different times.
From the results shown in Figure 12, the maximum values of plastic strain energy vs. time are plotted in Figure 13.
At 20 ms, the plastic strain energy density was mainly concentrated on the left and right sides of the roadway and extended outward at the maximum value of 1.93 × 10 7 J. At this time, the bolts on the two sides of the roadway have been damaged and deformed. At 40 ms, the plastic strain energy density was distributed on the left and right sides at the maximum value of 4.79 × 10 7 J. At 60 ms, the plastic strain energy density was mainly distributed in the left and right arch and both sides of the roadway at the maximum value of 1.68 × 10 8 J. Also, at 80 ms, the plastic strain energy density was mainly distributed around the roadway and the position of the floor. The maximum value of the plastic strain energy density at 80 ms was 1.70 × 10 8 J.
The two sides of the roadway and the arch toe were deformed, and the section was shrunk about 0.5 times the initial section. Figure 14 shows a variation nephogram of the damage variable D at different times.
From the results shown in Figure 14, the damage variable D values with time are plotted in Figure 15.
The results show that the maximum damage of the roadway floor and the left and right sides of the rock mass is 0.1 at 20 ms Besides, the left and right sides of the roadway and the arch toe were damaged at 40 ms The damage value of the left arch was greater than that of the right arch. The maximum value of the damage was 0.22 at 40 ms Also, at 60 ms, the damage has appeared on the left and right sides and arch corners of the roadway floor. The maximum value of the damage was 0.82 at 60 ms Moreover, at 80 ms, the damage occurred at the bottom of the roadway and the arch toe. The maximum value of the damage was 1.16 at 80 ms When D is greater than 1, the rock mass is broken, indicating the rock mass destruction at the bottom of the roadway and the arch during 60 ms-80 ms

GRADIENT OF ROCK BURST
The size of the energy density U d determines the energy accumulation characteristics of a certain area. 24 When the U d in the surrounding area reaches the maximum value, it is the high energy accumulation area. This also shows that the greater the (U d ) max , the higher the degree of energy accumulation in the area, and the more energy may be released. On the other hand, if the distance between a point of the maximum energy density (U d ) max and the free surface of the roadway is smaller, that is, the value of d is smaller, the energy consumed for the failure and instability of the surrounding rock is less, and the sudden failure and instability of the surrounding rock are more likely to occur, that is, the greater the tendency of rock burst. Using the combination of the (U d ) max point position and the roadway free surface distance d, the surrounding rock energy density distribution factor can be determined by the tendency of the roadway surrounding rock to burst.
The value of K is related to the stress state of the surrounding rock, lithology, and the characteristics of the disturbing stress wave. When the value of K exceeds a certain critical value, the surrounding rocks will have a rock burst. The critical value (K*) is the critical energy density distribution factor, which reflects the relationship between the energy accumulation degree of the surrounding rock and the relative distance of the free surface. The following condition (Equation 3) is necessary to occur the rock burst.
The critical energy density distribution factor is only related to the inherent properties, such as lithology and structural characteristics of the surrounding rocks. It must be obtained through experiments, field measurements, and much statistical data of rock bursts. According to the actual situation, the energy density distribution factor K* of the surrounding rock can be obtained by numerical analysis.
The value of energy density (U d ) max reflects a certain extent that the surrounding rock has a higher degree of energy accumulation at this location. However, the energy density U d value of a point can only reflect the energy accumulation characteristics of the point in the surrounding rock mass. However, it cannot explain the rock mass's total energy characteristics and accumulation degree. It can be seen that there are obvious limitations in judging the tendency of rock burst at this location only by the magnitude of the energy density value of a point. Based on the rock burst energy density criterion, the energy gradient criterion is proposed to judge the shocking tendency. The energy gradient is expressed as follows: where V is the volume of a point energy accumulation zone. U and U d are the total energy and the energy density of a point, respectively. (U) max and (U d ) max are the maximum energy and the maximum energy density, respectively. where K is the energy distribution factor, △K are the energy gradient, d is the distance from the center of the roadway, d 1 and (U d ) max are the position and maximum energy of one point from the center of the roadway, and d 2 and U 2 are the other points, respectively. The position and energy from the center of the roadway, △K*, are the critical impact tendency index.
From the expression of the energy gradient, a significant variation in the ∆K means that the energy value of the rock mass in the vicinity of the point drops rapidly. However, the energy accumulated in the vicinity of this point is not necessarily significant. The changing range of the ∆K is relatively small, indicating that the energy value of the rock mass near this point decreases slowly. On the other hand, the rock accumulates a large amount of energy in a certain area, and the energy accumulation area is relatively close to the roadway boundary. In this case, a large amount of energy may be released instantaneously, resulting in a rock burst accident.
This paper proposes an energy gradient criterion for rock burst tendency based on the energy density distribution analysis of the surrounding rocks. Also, the energy accumulation of surrounding rock is measured by determining the changing range of energy gradient. Besides, this paper provides a new method for predicting rock burst and makes the energy gradient criterion more reasonable by quantifying the tendency index of rock burst.

| CONCLUSION
Based on the analysis of rock burst events in the Yuejin Mine and Qianqiu Mine of Yimei mini area, combined with geological conditions, shock location, and source waveform, numerical simulation methods are used to inverse analyze the dynamic response characteristics of the surrounding rock of the roadway and proposed the energy gradient criterion of rock burst tendency. The details are as follows: The plastic strain energy density of No. 25110 working face in Yuejin Mine is mainly distributed around the roadway, and the maximum plastic strain energy density is 3.30 × 10 7 J at 75 ms, all rock masses around the roadway are damaged, and the maximum damage variable value is 0.44.
The plastic strain energy density of No. 21032 return air uphill in Qianqiu mine is mainly distributed around the lower side of the roadway, and the maximum plastic strain energy density is 2.68 × 10 8 J at 80 ms, all rock masses around the roadway are damaged, and the maximum damage variable value is 1.76, indicating that the rock masses at the roof and two sides of the roadway have been damaged after the impact.
The energy gradient criterion of rock burst tendency was proposed according to the existing energy density criterion and the energy density analysis of surrounding rock.