Microsimulation Study on Energy Release and Rock Block Ejection Force of Granite under Different Unloading Conditions

Rock burst is one typical dynamic disaster caused by excavation in deep underground engineering. High-stress unloading test is a common research method for rock bursts. Due to the limitation of laboratory test conditions, it is difficult to monitor the energy release and dissipation information during rock bursts in the unloading test. But the study of energy evolution law is more helpful to reveal the essential characteristics of rock burst. Therefore, the energy evolution process and ejection failure characteristics of granite after unloading were analyzed through the unloading simulation test in this paper, and the influence of unloading velocity, lateral stress, and axial stress were researched. The microstructure numerical model of the granite was established by using digital image processing technology and PFC2D software, aiming to match the real granite. The energy evolution process of unloading granite can be divided into three stages, namely the whole energy rapid release stage, sidewall energy slow-release stage, and rock block ejection stage. The area near the unloading sidewall is the main energy release and rock block ejection area. In the whole energy rapid release stage, the energy release velocity and dissipation velocity show similar law, i.e., a positive power function correlation with unloading velocity, a negative power function correlation with lateral stress, and a positive linear function correlation with axial stress. In the rock block ejection stage, with increasing the unloading velocity and axial stress, the rock block ejection force increases as a power function, while it decreases with increasing lateral stress. This research is an important supplement to the laboratory unloading test. It has theoretical guiding significance for rock burst hazard assessment during excavation in deep underground engineering.


INTRODUCTION
Rock burst is a common dynamic disaster in the deep hydraulic tunnel, traffic tunnels, and another underground rock engineering. With the expansion of the scale and depth of underground engineering construction, the frequency and intensity of rock burst disasters increase Feng et al., 2012;Zhao et al., 2018;Yin et al., 2019;Guo et al., 2019). Rock burst disaster has become a research focus in the field of deep rock engineering.
Rock burst mechanism and key disaster-causing factors identifications are the premise and foundation for prevention and control of rock burst. Many experts and scholars have carried out a lot of research work on rock burst mechanisms from the aspects of external stress load, mechanical properties of the rock, and dynamic load (Zhao et al., 2015;Xiao et al., 2016). He et al. (2012; divided the rock burst into two categories, namely, strain rock burst and impact-induced rock burst, and obtained their criteria. In order to understand the mechanism of rock burst and simulate the conditions of rock burst, Jiang et al. (2010) proposed a new energy index, namely local energy release rate, which can better evaluate the rock burst phenomenon in deep underground excavation. Li (2020) researched the principle of rock support to control rock bursts, and explained the energy source and conversion during rock burst in the conceptual model. Li et al. (2014) studied the influence of rock excavation unloading in the deep tunnel, and analyzed the influence of stress path and other factors on rock excavation unloading. Excavation unloading is an important cause of rock burst (Wu et al., 2009). True triaxial unloading test is often used to simulate rock burst process. For example, He et al. (2009) carried out one-side dynamic unloading test under true triaxial condition, analyzed the spectrum of full-wave acoustic emission data, obtained the frequency-amplitude relationship of acoustic emission signal in rock burst stage, and summarized the dynamic damage process and characteristics of rock under test conditions. Li et al. (2014) carried out true triaxial rock burst tests at different unloading velocities, and studied the damage evolution law of rock burst. Si and Gong (2020) carried out triaxial compression and biaxial compression tests under the same confining stress, and the effects of three-dimensional stress state and unloading velocity on rock burst mechanism was studied. Zhai et al. (2020) researched on rock burst process by using high-speed camera, acoustic emission system, and scanning electron microscope. Existing studies provide a lot of research basis for revealing rock burst mechanisms.
Rock burst is a dynamic failure phenomenon driven by energy. Studying the energy evolution law in the process of rock burst is more helpful to reveal the essential characteristics of rock burst (Wang et al., 2012). However, limited to the laboratory test conditions, the energy evolution process is less studied in the true triaxial unloading test, and more attention is paid to the study of disaster-causing stress conditions, acoustic emission characteristics, rock block fractal dimension, and ejection velocity, et al. Compared with laboratory tests, numerical simulation is an effective method to study rock burst energy (Sun et al., 2007). Therefore, this studytakes the granite which is commonly used in rock burst tests as the research object. Based on the micronumerical model of the granite, the energy evolution law and rock block ejection failure characteristics of granite during unloading were researched. The influences of unloading velocity and stress level on rock burst were studied.

Microstructure Characterization of Granite
Rock is a typical heterogeneous material, and its heterogeneity mainly includes two aspects. One is the heterogeneity of microparameters, and the other is the heterogeneity of minerals distribution. The weibull function is often used to characterize the heterogeneity of microparameters for rock materials with relatively homogeneous mineral structures, such as the sandstone. For the granite, its minerals distribution is very obvious, and the color difference between different minerals is also relatively large. Therefore, the microscopic structure is often extracted and characterized according to the color difference of rock minerals to establish the heterogeneous model for the granite (Chen et al., 2003;Chen et al., 2004;Zhu et al., 2008;Tan et al., 2016;Meng et al., 2018). The digital image processing technology is an effective method to obtain the microstructure of the rock. Its basic principle is to re-segment the digital image of rock into some groups according to the color of different minerals. The granite used in this simulation test (as shown in Figure 1A) is composed of black mica, transparent quartz (appearing gray), and white feldspar. The key to the characterization of the microstructure is to obtain the color segmentation thresholds among the three minerals. In this paper, the Otsu multithresholds segmentation method is used to obtain the segmentation thresholds of granite microstructures. The thresholds (t 1 , t 2 ) that satisfy the maximum gray variance S 2 among the minerals are taken as the optimal segmentation thresholds, and the formula is: where u is the whole gray value of the granite image, N i and u i are the number and average gray value of pixels of the three different minerals, respectively. The microstructure characterization image of the granite with a size of 100 × 100 mm obtained by the Otsu multi-thresholds segmentation method was shown in Figure 1B. The distribution and shapes of the three minerals in the characterization image are highly consistent with the real granite image.

Building of Microstructure Numerical Model
There are many numerical methods to simulate rocks' mechanical properties and behavior. The finite element method (FEM) and fast Lagrangian analysis of continue (FLAC) are continuum analysis methods. They can simulate Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 909371 the yield, plastic flow, softening, and large deformation of materials, especially in the field of elastic-plastic analysis, large deformation analysis, and simulation of the construction process. Peridynamics (PD) (Wang et al., 2016;Wang et al. 2017;Wang et al., 2018), field enriched finite element method (Zhou et al., 2021), and general particle dynamics (GPD) (Zhou et al., 2015;Bi et al., 2016) can simulate the initiation, propagation, and coalescence of microcracks as well as formation of macro-scale cracks. Particle flow code (PFC) (Potyondy and Cundall, 2014) is one discrete element method, and it can simulate rock macro mechanical behavior by analyzing the movement and interaction of particles. It is not limited by the amount of deformation and is convenient to deal with the problem of the discontinuous medium. These numerical methods above are very helpful for studying rock fractures. In this simulation, not only the energy evolution law in the process of rock fracture was studied, but also the rock block ejection characteristics after rock burst. Compared with other numerical methods, PFC can intuitively obtain the rock block motion behavior of rock failure (Yoon, 2007). Therefore, this paper used PFC 2D to build the microstructure model of granite and study the mechanical behavior during unloading. In the modeling process, firstly one granite specimen model with a size of 100 × 100 mm was generated by using PFC. In order to reduce the size effect , the radius of the particles was set to 0.4-0.5 mm, and the average ratio of specimen size to particle radius was 220. Then, the coordinates and gray information of each pixel in the characterization image were imported into PFC through the FISH function. According to the gray information of pixels, the color codes of particles with the closest coordinates were assigned respectively . The microstructure numerical model of the granite specimen as shown in Figure 1C. The black particles are mica, the gray particles are quartz, and the white ones are feldspar, as shown in Figure 1D. The three mineral particles were assigned with corresponding material parameters respectively, and the parameters of the cementation surface between minerals were set according to the parameters of low strength minerals. A parallel bond model was chosen to simulate the mechanical properties of rock materials. In the parallel bond model, the microparameters including the friction coefficient μ, contact stiffness k n and k s , parallel bonded stiffness k n and k s , parallel bonded strength σ n and σ s , and parallel bonded radius λ. Based on the work by Potyondy and Cundall. (2004), the microparameters can be determined by the following formulas: where E c and E c are Young's modulus of the particles and bond respectively; r is average radius of two contact particles. In general, k n /k s and k n / k s can be set to 2.5, λ is 1, μ is 0.5. E c and E c are proportional to the macro Young's modulus of the rock. σ n and σ s are proportional to the macro strength. Using Eqs 2, 3, the microparameters of the model can be estimated, but there are often some errors. Therefore, a trial-anderror method should be used to calibrate the suitable microscopic parameters (Peng et al., 2018). In the paper, the microstructure model parameters of the granite were selected from Zhang et al. (2019), as shown in Table 1.

Simulation Scheme
In thisarticle, numerical simulation test was carried out based on true triaxial loading and single-side unloading laboratory test of granite.
In the loading stage, the axial stress and lateral stress were applied to the granite specimen by the loading wall. When the stress was loading to the predetermined value, the displacement of the upper and the left wall was fixed, avoiding unreasonable loading displacement of the wall under stress constraint, which is different from the laboratory test. In the unloading stage, the right wall was applied with a certain velocity to unload the lateral stress. The simulation variables are unloading velocity, lateral stress, and axial stress, respectively. The specific schemes were shown in Table 2.

Energy Calculation and Monitoring Schemes
In PFC 2D , the work done by the external load on the granite specimen is the input energy, which can be calculated by the load and displacement of the loading wall, as follows: where E tot is the input energy, F ε is the load of the loading wall, and ΔU ε is the displacement increment of the loading wall. The input energy of granite under external load is divided into two parts, namely, strain energy and dissipation energy.
The strain energy includes two parts, one is the contact strain energy between particles, and the other is the parallel bond strain energy. The calculation formula of strain energy is as follows: where E pb str where E str is the strain energy, E c str and E pb str are the contact strain energy and the parallel bond strain energy, respectively. F n i and F s i are the normal force and tangential force at the contact point i of particles, respectively. F n i 、 F s i 、 M i is the normal force, tangential force, and bending moment respectively at the i parallel bond. A i is the parallel bond interface area, and I i is the moment of inertia. N c and N pb are the number of particle contact and the number of parallel bonds, respectively.
The calculation formula of dissipation energy E dis is: In the simulation process, in addition to monitoring the energy and crack of the model, the model was divided into five groups, and the energy data in different groups of the model were monitored too. Each group was 20 mm wide, and group No.1 was close to the unloading wall, as shown in Figure 1C.
The rock burst area and rock block ejection force are also important indicators for evaluating the rock burst intensity. A force measuring wall was arranged at 40 mm far away from the unloading wall to monitor the rock block ejection force.
In order to clearly and reasonably describe the process of energy release, the parameters need to be defined as follows: v rel energy release velocity-equals the ratio of released strain energy to time; vdi s energy dissipation velocity-equals the ratio of dissipation energy to time r dis energy dissipation ratio-equals the ratio of dissipation energy to released strain energy; F rock block ejection force-represents rock block ejection force on the force measuring wall; A rock burst area-represents the size of the rock burst pit.

ENERGY RELEASE LAW OF UNLOADING GRANITE Energy Release Process and Stage Division
Taking the test results of the granite with lateral stress of 30MPa, axial stress of 160 MPa, and unloading velocity of 0.1 mm/s as an example, the energy and stress evolution curves after unloading were shown in Figure 2. According to the energy release and dissipation law of the whole granite and each group, the energy release process can be divided into three stages: I-whole energy rapid release stage, II-side wall energy slow-release stage, and rock block ejection stage.

1)
Step I: whole energy rapid release stage. During the initial unloading stage, with decreasing lateral stress, the lateral deformation of the specimen is beginning to recover, and part of the strain energy accumulated in the specimen is released in the form of elastic energy, as shown in Figures  2A,B. The deformation recovery of the specimen gradually reduces from the unloading sidewall to the internal part ( Figure 3A), and the corresponding released energy is also gradually reduced, as shown in Figure 2C. According to the Mohr-Coulomb criterion, when the lateral stress decreases, the ultimate bearing strength of the specimen decreases. During unloading, the crack gradually generates, and part of strain energy is released in the form of dissipation energy.  With increasing crack number, the energy dissipation ratio increases gradually. Because cracks are mainly generated near the unloading sidewall of the granite, the energy release in this region is more obvious. 2) Step II: sidewall energy slow-release stage. After the lateral stress is reduced to zero and the lateral deformation of the specimen is recovered, the crack near the unloading side wall gradually connects, as shown in Figure 3B. Only the strain energy near the side wall still keeps releasing, and the energy release velocity is slower than that in Step I. The energy dissipation ratio gradually reduces due to the decrease of crack generating velocity.

3)
Step III: rock block ejection stage. After the crack near the side wall coalesced into the fracture zone, the rock block begins to eject from the specimen, and part of the strain energy is released in the form of kinetic energy, as shown in Figure 2D and Figure 3C. After the ejection of large size blocks, the bearing capacity of the specimen decreases sharply, and the fracture zone expands rapidly. A large number of small size blocks generate and eject in a short time. After rock blocks eject, there appears a triangular rock burst pit, which is similar to the engineering field example in Li (2020). The size of the rock burst pit equals the rock burst area.

Influence of Unloading Velocity on Energy Release
With the same lateral stress and axial stress, the energy evolution curves of granite under different unloading velocities were shown in Figure 4. With increasing unloading velocity, the duration time of energy release is significantly shortened, while the total energy release amount changes little. The farther away from the sidewall, the smaller the strain energy released. Both the energy dissipation ratio and total amount increase with increasing unloading velocity. In the whole energy rapid release stage, the strain energy release velocity and dissipation velocity increase with the unloading velocity in a power function relationship, as shown in Figure 5.

Influence of Lateral Stress on Energy Release
The energy evolution curves of granite under different lateral stress were shown in Figure 6. Under the same axial stress, with increasing lateral stress, the accumulated strain energy in granite before unloading decreases. This is because when the axial stress does positive work on the specimen, the lateral stress does negative work. The greater the lateral stress is, the greater the Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 909371 8 negative work is done, which causes smaller accumulated strain energy. The result is consistent with that of Huang and Li (2014).
When the lateral stress is 10 MPa, the energy is continuous released at a high velocity, and the specimen fractures whole in a short time after unloading. When the lateral stress exceeds 20 MPa, the three-stage characteristics of energy release of granite are relatively obvious. With the increase of lateral stress, the total amount and velocity of energy release and energy dissipation decrease gradually, which all show a power function relationship with lateral stress, as shown in Figures 7A,B. Also, the main area of energy release gradually decreases when the lateral stress increases. When the lateral stress is 30 MPa or 40 MPa, the energy release body is only within the area of 40 mm from the side wall, as shown in Figure 7C.

Influence of Axial Stress on Energy Release
The energy evolution curves of granite under different axial stress were shown in Figure 8. With the increase of axial stress, the total energy accumulated in the granite increases gradually before unloading and the duration of energy release after unloading increases too. The energy release velocity and dissipation velocity are linearly correlated with the axial stress, as shown in Figure 9A. After rock burst, the total amount of energy release and energy dissipation, and the energy dissipation ratio increase with the increase of axial stress, as shown in Figure 9B. Under low axial stress, such as 160 and 170 MPa, the main energy release area locates at the 40 mm range from the side wall. When the axial stress exceeds 180 MPa, the whole body is the main energy release area, as shown in Figure 9C.

Influence of Unloading Velocity
The unloading velocity has little effect on the rock burst area. However, due to the rapid release of energy at high unloading velocity, more energy is released in the form of kinetic energy, and the rock block ejection force is significantly increased, which shows one power function relationship with the unloading velocity, as shown in Figure 10.

Influence of Lateral Stress
Rock burst forms and rock block ejection force of the granites under different lateral stress were shown in Figure 11. When the lateral stress is 10 MPa, the rock burst area is about 3,000 mm 2 , and the rock block ejection force is about 3.21 kN. With increasing lateral stress, rock burst area, and rock block ejection force decrease, which shows one exponential function with lateral stress. When the lateral stress is 40 MPa, rock burst area reduces to 60%, and the rock block ejection force decreases to 44.6%, compared to that at 10 MPa.

Influence of Axial Stress
Rock burst forms and rock block ejection force of the granites under different axial stress were shown in Figure 12. When the axial stress is 160 MPa, rock burst area is about 1700 mm 2 , and the rock block ejection force is about 1.25 kN. When the axial stress is 200 MPa, rock burst area extends to 2.82 times, and rock block ejection force increases to 5.2 times, compared to that at 160 MPa. The rock burst area has one linear relation with axial stress, and the rock block ejection force has one power function relation with axial stress. CONCLUSION 1) Unloading process of the granite can be divided into three stages: I-whole energy rapid release stage, II-side wall energy slow-release stage, and III-rock block ejection stage. The release and dissipation of strain energy decrease with the increase of distance from side wall, and the main energy-releasing body and ejection body are near the side wall. 2) The unloading velocity mainly affects the duration of the three stages of energy release, the energy release velocity and the dissipation velocity. It has little effect on the total amount of strain energy released and rock burst area. The strain energy release velocity and the energy dissipation velocity are positively correlated with the unloading velocity by one power function, and the rock block ejection force is also positively correlated with the unloading velocity by one power function. 3) Under the same axial stress, with the increase of lateral stress, the released energy amount, dissipation energy amount, and the energy dissipation ratio of unloading granite gradually decrease. Moreover, the strain energy release velocity and the energy dissipation velocity are negatively correlated with lateral stress in a power function, and the rock block ejection force and rock burst area are also negatively correlated with lateral stress in a power function. 4) Under the same lateral stress, with the increase of axial stress, the released energy amount, dissipation energy amount, and the energy dissipation ratio of unloading granite gradually increase. Moreover, the strain energy release velocity and the energy dissipation velocity are positively correlated with the axial stress in a linear function. The rock burst area has one linear relation with axial stress, and the rock block ejection force has one power function relation with axial stress.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.