Mechanism of Rock Burst Revealed by Numerical Simulation and Energy Calculation

We propose a new computational method to calculate the storage elastic energy value of surrounding rocks based on numerical simulation and theoretical calculation. By calculating the difference value in energy under different force states and comparing them with the energy level when rock burst occurs, we get the mechanism of rock burst: when roadway and surrounding rocks are in the condition of large ratio bias force field, certain triggering stress causes mass release of the elastic energy of surrounding rocks around the roadway, and when the energy reaches a certain level, rock burst will occur. We also put forward the specific force field conditions and triggering stress values of rock burst, which is of great guiding significance for the mechanism disclosure, monitoring, and control of rock burst.


Introduction
At present, rock burst is a common dynamic disaster in underground coal mines, deep buried tunnels, and underground metal mines [1][2][3], mainly because its mechanism is not clear enough [4], which leads to unclear monitoring and early warning methods. In the past research, the mechanism of rock burst mainly includes strength theory [5], energy theory [6], stiffness theory [7], rock burst tendency theory [8], "three factors" theory [9], instability theory [10], "three criteria" theory [11], butterfly impact mechanism [12,13], dynamic and static combined load principle [14,15], and rock burst start-up theory [16]. ese theories reveal the occurrence mechanism of rock burst from various perspectives, but they either have their limitations or used too complex model and do not have better guidance for the monitoring and control of rock burst on-site. e famous "Ockham razor" principle [17] believes that the simplest explanation of the phenomenon is often more correct than the complex one. So, a feasible numerical model should be adopted for quantitative calculation to explain the occurrence mechanism of rock features [18]. And meanwhile, the real mechanism should be able to explain the occurrence principle of rock burst in various situations, rather than only explaining special cases. e failure of coal and rock mass is a change in state driven by energy [19]. erefore, studying the deformation and failure process of coal and rock mass from the perspective of energy will truly reflect its failure law [20]. According to the energy theory, when the energy given to the surrounding rock system by the outside exceeds the energy consumed before its failure, the rock burst will occur, and there is no doubt that this is correct, but there is no theoretical answer to how much energy will be released under what kind of mechanical state of system and triggering stress (ΔP) [21]. Here, the ΔP is often caused by roof breaking or leading abutment pressure [22][23][24].
e main monitoring methods of rock burst are microseismic (MS) monitoring [25][26][27], electromagnetic radiation [27], method of drilling bits [28], and stress monitoring [29]. Among these monitoring methods, MS monitoring is more representative and valuable, mainly due to its intuitiveness and application of vibration signal (MS event) monitoring. e main control measures are largediameter borehole pressure relief [30,31], blasting pressure relief [2,32], roadway support [33], and hydraulic fracturing [34], but the monitoring and control methods are heavily dependent on the mechanism; otherwise, the monitoring and control measures are not based at all, for example, what energy level of MS events can reach the level of giving an early warning of the rock burst, whether the large-diameter pressure relief measure really play a role, and how to play a role.
In view of the above facts, this paper takes numerical simulation and theoretical calculation as means and calculates the release energy of system under different mechanical force states and ΔP based on the simplest model. e calculation method and the revealed mechanism of rock burst are obviously of great significance for the prevention and control of rock burst.

Method
e surrounding rocks of underground roadway or tunnel are in the state of three-dimensional stress, in order to simulate its surroundings more truly, and we take a certain range of bounded area and named it Ω (respect surrounding rock), set one hole (respect roadway or tunnel) in the center, and then apply an external force (P1, P2, P3) to the Ω area, where P 1 , P 2 , and P 3 are the maximum force, intermediate force, and minimum force applied to Ω, respectively. Under the external force, unit in the area will be subjected to the stress of (σ 1i , σ 2i , σ 3i ), where σ 1i , σ 2i , and σ 3i are the maximum principal stress, intermediate principal stress, and minimum principal stress of unit, respectively. So the elastic energy of unit (f(x, y, z)) is as follows [35]: where E i is the elastic modulus of unit and μ i is Poisson's ratio of unit. For the sake of expression, we define (P1, P2, P3) as the front state force (FSF). Under FSF, there is no plastic zone in the Ω region if it is an elastic model, the storage elastic energy of Ω (named U FSF ) can be expressed by Equation (2), while if it is an elastoplastic model, there will be some plastic units (Ω p ) exit in Ω, and others are elastic units (Ω s ); the stored elastic energy of Ω (named U FSF ′ ) is shown in Equation (3). e difference value (D-value) of elastic energy under FSF (D FSF ) should be expressed as Equation (4). Obviously, the D-value is the reduced energy of the system between elastic model and elastoplastic model under the FSF: Under the influence of ΔP, it is assumed that FSF changes to (P1 + ΔP, P2, P3); that is, the ΔP only increases to P 1 . For the sake of expression, we define (P1 + ΔP, P2, P3) as the latter state force (LSF). So, under LSF, there is no plastic zone in Ω if it is an elastic model, the storage elastic energy of Ω (named U LSF ) can be expressed by Equation (5), while if it is an elastoplastic model, there will be some new plastic units (ΔΩ p ) arising in Ω, and the total plastic units should be (Ω p + ΔΩ p ); elastic units should be(Ω e − ΔΩ p ), so the storage elastic energy of Ω (named U FSF ′ ) is shown in Equation (6). e D-value of elastic energy (D LSF ) under LSF should be expressed as Equation (7). e D-value is the reduced energy of the system between elastic model and elastoplastic model under the LSF: When rock burst occurs, the energy that causes huge damage to roadway and surrounding rock mass is mainly elastic wave energy [36]. In our method, all the released elastic energies will be converted into various forms of energy, such as dissipation energy, kinetic energy, and friction heat energy, only a part of which will convert into elastic wave energy, so we assume an energy conversion factor (β) of elastic wave. According to relevant research [37,38], the factor is generally between 1% and 10%. So, the release elastic wave energy (W) from FSF to LSF should have the difference value between the two D values and multiplied by β, as follows: where β is the energy conversion factor of elastic wave energy. e resulting calculation method is shown in Figure 1; the release energy caused by the change in mechanical state is finally obtained by four models. e three-dimensional stress of the unit mentioned in the above calculation method is very difficult to get in the theoretical calculation and field measurement, while FLAC 3D just can simulate and get the stress state of each unit in the system, so as to calculate the elastic energy of the unit. According to the above calculation method, we set up a calculation system, the length and width of the system are 200 m, the thickness is 1 m, and a hole with the diameter of 5.6 m is set in the center of the system, as shown in Figure 2. e system uses the same kind of coal as medium, its cohesion (C) is 3 MPa, and its internal friction angle (ψ) is 25°. e loading directions of P 1 , P 2 , and P 3 are up and down, front and back, and left and right, respectively. In order to simplify the calculation, we only change the size of P 1 to obtain the energy of the system when we simulate different force states. e initial force apply to the model is set to P 1 � P 2 � P 3 � 20 MPa, and then, only the value of P 1 increases 1 MPa each time until the system is fully plastic. When we calculate the energy under elastic model and elastoplastic model, we use the command of "model mech elastic" and "model Mohr" in the FLAC 3D program, respectively. In this way, we can get the storage elastic energy under different force states, so as to calculate the release energy of the system.

Results
After dozens of numerical simulations, we found that the farthest range of plastic zone reaches the model boundary when P 1 � 58.7 MPa under the elastoplastic model. We got the failure state of the elastic model and the elastoplastic model under each force state and discover that, under the elastic model, the system has not been damaged at all, but under the elastoplastic model, the stress concentration around the hole is obvious, the shape of the plastic area has experienced "circle," "ellipse," and "butterfly" successively, as shown in Figure 3. Finally, when the "butterfly leaf" extends to the boundary of the system, the system is completely plastic. e force state of the beginning of butterfly shape is P 1 � 40 MPa (the ratio of P 1 /P 3 is 2), the butterfly state is particularly obvious when P 1 � 50 MPa (the ratio of P 1 /P 3 is 2.5), the system is damaged in a large range when P 1 � 58.6 MPa (the ratio of P 1 /P 3 is 2.93), and the system is all plastic when P 1 � 58.7 MPa (the ratio of P 1 /P 3 is 2.935).
We calculate the elastic energy density of each unit, place them in the system according to the coordinates, and form the cloud chart to check the energy distribution; the cloud chart of a typical force state corresponding to Figure 3 is shown in Figure 4. With the increase in P 1 , the storage elastic energy of the system under elastic model and elastoplastic model is all increasing. When P 1 is small, the difference in energy distribution between the two models is small, but with the increase in P 1 , the energy change range of the elastoplastic model is larger than that of the elastic model, and the distribution of concentrated energy coincides with the position of butterfly leaf, indicating that there is a certain correlation between the distribution of energy and the distribution of plastic zone. Table 1 lists the calculation energy results of the two models under different force states, and we can find that, with the increase in P 1 , the storage elastic energy of the system increases from 10 9 J to 10 11 J. e storage elastic energy of the elastic model in each mechanical state is always greater than that of the elastoplastic model, and according to Method, the D-value is the reduced energy of the system under the corresponding force state. e D-value of elastic energy increased from 8.05E + 06 J to 9.32E + 08 J, indicating that with the increase in P 1 , the release energy increased gradually and correspondingly. Although the increment of P 1 is all 1 MPa, the release elastic energy (W) from FSF to LSF increases gradually from the beginning 3.71E + 05 J (P 1 increased from 20 MPa to 21 MPa) to the last 2.42E + 08 J (P 1 increased from 57 MPa to 58 MPa).
In order to see the trend of energy change under different force states more intuitively, we draw a change curve of the energy data, as shown in Figure 5, and we can clearly see that, with the increase in P 1 , the storage elastic energies of the two models are all increasing, but the increasing speed is different; the greater the stress is, the greater the difference between them is. e D-value of elastic energy shows different growth rates with the increase in P 1 ; before P 1 � 40 MPa, the rate is small, and the reduced energy is magnitude of 10 6 -10 7 J; between P 1 � 40 MPa and 50 MPa, the rate is increased, and the reduced energy is magnitude of 10 7 -10 8 J, but after P 1 � 50 MPa, the growth rate is obviously accelerated, and the reduced energy reach magnitude of 10 8 J.
e release elastic energy from FSF to LSF is shown in Figure 6. e magnitude of release energy is obviously different under different force states although the increment of P 1 is the same (1 MPa). e release energy curve should be divided into three stages, namely, "pregnant period," "growth period," and "upheaval period," respectively. In the pregnant stage, the release energy is mainly less than a magnitude of 10 6 J, and the corresponding P 1 is less than 49 MPa. In the growth stage, the release energy is mainly at a magnitude of 10 7 J, and the corresponding P 1 is between 50 MPa and 54 MPa, but in the upheaval period, when P 1 is big than 54 MPa, the release energy magnitude reaches 10 8 J.
Based on the above numerical simulation results, we can easily find that the main energy source of rock burst is the release of elastic energy of surrounding rock, which is the same as classical energy theory, but our results indicate that Shock and Vibration 3 the main source of energy released is the release energy of plastic zone, which are more specific. e more important conclusion is that ΔP has a greater probability of causing rock burst under the FSF condition of bias force field. e larger the P 1 /P 3 ratio of bias force field is, the greater the energy generated by the rock burst is, especially when the P 1 /P 3 ratio is above 2.5.

Discussion
We know that the magnitude of release energy during occurrence of rock burst is basically 10 7 -10 8 J [39,40]. e thickness of the system in our calculation method is 1 m, which means that our calculation result of the release energy is 1 m thick, while rock burst often causes tens to hundreds of meters of roadway damage. Taking 100 m failure length as an example, the elastic wave energy received by the MS system during the occurrence of rock burst is assumed to be 10 8 J, so the elastic wave energy per meter is 10 6 J. According to the calculation results in Results, the release energy is 10 6 J in the pregnant stage; considering β is 1%-10%, the elastic wave energy is 10 4 -10 5 J; similarly, the elastic wave energy of growth stage and upheaval stage is 10 5 -10 6 J and 10 6 -10 7 J, respectively. From the results, we can see that if ΔP is 1 MPa, the occurrence of rock burst is mainly in the growth stage and upheaval stage. e P 1 /P 3 ratios of growth stage and upheaval stage are all bigger than 2.5 (bias force field), and the shape of the plastic zone is "butterfly," which seems to be the precondition for the occurrence of rock burst. In order to better explain the issue, in this chapter, we mainly discuss the

eoretical Calculation of Plastic Zone.
Our scientific research team [12] has already obtained the theoretical calculation of plastic zone distribution under different force fields and applied to reveal the mechanism of earthquake [41][42][43].
e theoretical recessive formula of plastic zone boundary around a circular hole (model as Figure 7) is shown in the following equation: where a is the radius of the circular hole, r and θ are the polar coordinates of any point on the boundary of plastic zone, η is P 1 /P 3 , C is the cohesion of rock, and φ is the internal friction angle of rock. According to the theoretical formula, the calculated plastic zone results under different force states corresponding to Figure 3 are obtained as Figure 8. By comparing the two figures, we can find that the shape of plastic zone is basically the same in theoretical calculation and numerical simulation. e slight difference in the degree of smoothness is mainly due to the consideration of deformation compatibility equation in numerical simulation. Under the same basic conditions of numerical simulation, we get the change curve of the maximum radius (R max ) and area (S) of plastic zone with change in P 1 under the theoretical calculation, as shown in Figure 9. From the results, we can see that the change trend of R max and S is different, but the overall change trend is the same as the release energy curve obtained by numerical simulation, especially after the "butterfly" shape state (P 1 � 50 MPa and η � 2.5), and the theoretical calculation also shows that when η is 2.93, R max and S of the system will be very large. e occurrence of rock burst is often a moment of huge damage, and our results of energy release by numerical simulation and results of plastic zones calculated by theoretical calculation can explain the disaster obviously. at is to say, under the action of certain ΔP, the elastic energy of surrounding rocks is released instantaneously, which causes part of surrounding rocks to be thrown out. e occurrence of rock burst is more likely to occur under the condition of bias force field because small ΔP can produce huge energy release under this condition. Shock and Vibration

Triggering Stress.
Another key point of the problem is how large ΔP is, which is directly related to whether it can cause system to the energy standard of rock burst under certain FSFs. We have consulted some references and carried out a study by means of theoretical calculation and MS monitoring. Generally, the elastic wave is divided into P-wave and S-wave; P-wave is a longitudinal wave, and S-wave is a shear wave. e stress increment generated by P-wave (σ dP ) and S-wave (σ ds ) can be calculated by the following equations [34]:   180  160  140  120  100  80  60  40  20  180  160  140  120  100  80  60  40  20   180  160  140  120  100  80  60  40  20  180  160  140  120  100  80  60  40  20   180  160  140  120  100  80  60  40  20  180  160  140  120  100  80  60  40  Energy distribution of elastic model Energy distribution of elastoplastic model

Shock and Vibration
where ρ is the density of medium; C P and C S are the propagation velocities of P-wave and S-wave, respectively; and v pp and v ps are the peak particle velocity (PPV [44]) of P-wave and S-wave, respectively. Suppose μ � 0.29, ρ � 1400 kg/m 3 , C P � 3500 m/s, and C S � 2100 m/s (propagation speed of the S-wave is about 0.6 times of the P-wave [45]); the maximum PPV is 1 m/s according to a relevant study [40,46]. By calculation, the maximum of σ dP and σ dS can reach 3.74 MPa and 2 MPa, respectively. In fact, PPV measured by the MS system in a coal mine rarely reaches 1 m/s, the maximum value of PPV is 0.3 m/s according to the measured results in Poland [39], and the majority of PPV is about 0.0001 m/s to 0.1 m/s in the MS measurement results for half a year in a coal mine in Poland. If the maximum value of PPV is 0.3 m/s, the maximum of σ dP and σ dS can only reach 1.12 MPa and 0.88 MPa, respectively.
When σ dP reaches 3 MPa in the pregnant stage, it can trigger the release energy of 10 7 J, and the energy level may lead to the occurrence of rock burst, according to the analysis of Results. Similarly, when σ dP reaches 1 MPa in the growth period or σ dP reaches 0.1 MPa∼1 MPa in the upheaval period, it may lead to the occurrence of rock burst. erefore, we should monitor and forewarn the rock burst according to the FSF of roadway and the result of MS monitoring, instead of only relying on the MS system to reach a certain energy level. rough backcalculation, we get  Table 2.

Comparison of Nonbias Force Field Release Energy.
In order to better explain that the bias force field is the main factor of rock burst, we carried out numerical simulation under the condition of uniform force field and formed Figure 10. Under the condition of uniform force field, the storage elastic energy of the elastic model and the elastoplastic model is all increasing, which is the same as the bias force field, but the D-value of elastic energy is obviously different from the bias force field, and it does not show any increase in the rate of change, which shows that the energy released from FSF to LSF does not increase under the uniform force field. e release energy level of 10 7 J can be achieved only when the    forces in three directions are all increased by 5 MPa, but this kind of ΔP basically does not exist. Figure 11 shows the statistics of release energy under different FSF and ΔP conditions. If FSF is P 1 � P 2 � P 3 � 40 MPa and ΔP � 1 MPa, the release energy is 1.86E + 06 J; if FSF is small ratio bias force (P 1 /P 3 � 1.5) and ΔP � 1 MPa, the release energy is 1.94E + 06 J, and the energy magnitude is not met in with the occurrence of rock burst. While if FSF is big ratio bias force (P 1 /P 3 � 2.75) and ΔP � 1 MPa, the release energy is 1.43E + 08 J, which will probably lead to rock burst; if FSF is huge ratio bias force (P 1 /P 3 � 2.92), only ΔP � 0.2 MPa is needed, which will result in the energy release of 4.78E + 07 J and probably lead to rock burst. is fully proves a preliminary fact that high force field is not terrible, and the real danger is bias force field.
From the above analysis, it is fully confirmed that FSF is the main reason for the occurrence of rock burst, while ΔP is 6      the secondary factor. When the combination of the two factors causes the release energy to reach a magnitude of 10 7 J, rock burst will be caused. is important discovery is of great significance to the monitoring of rock burst, indicating that the monitoring should not only focus on ΔP, and the more important is FSF. It is also of great value to the governance of rock burst, which should aim to reduce the D-value energy of surrounding rock so that the energy can be released in advance; in other words, if the governance measures do not trigger the release of energy, it indicates that the governance measures are invalid. e only drawback is that our theory cannot be verified in the field at present (but a study [47] has confirmed the existence of butterfly damage in the laboratory) because the FSF of the roadway is difficult to measure, but we believe that, with the promotion of the theory, there will be corresponding stress measuring equipment for monitoring.

Conclusion
We put forward a method to calculate the release energy of surrounding rock system from one mechanical state to another. rough calculation, we mainly get the following conclusions: (1) e mechanism of rock burst is that, under a certain FSF, ΔP causes the instantaneous release of elastic energy of surrounding rock system and brings about huge damage; the release energy mainly comes from units of elastic to plastic. erefore, FSF and ΔP are the two factors that lead to the rock burst, and FSF should be the main factor. (2) Under different FSF conditions, the same ΔP will lead to different energy releases, which is closely related to the P 1 /P 3 ratio of FSF. Different FSFs are divided into three stages according to the ratio: pregnant period, growth period, and upheaval period. e P 1 /P 3 ratio dividing points of the three stages are 2.5 and 2.75.
(3) FSF and ΔP should be considered as the main factors in the monitoring and early warning of rock burst, and FSF should be the main factor. e ΔP of 3 MPa, 1 MPa, and 0.1 MPa of the example case is the minimum value to cause rock burst of the three stages, respectively. It is speculated that most of the rock burst is produced in the growth period and upheaval period, unless the PPV is superlarge.

Data Availability
is paper is a basic theoretical study, using numerical simulation and theoretical calculation methods. If necessary, we can provide numerical simulation program.
Shock and Vibration 13