Numerical Modeling on the Fracturing and Energy Evolution of Large Deep Underground Openings Subjected to Dynamic Disturbance

: The exploitation of deep resources is necessary for human development. At the same time, high-stress environments that are deep underground bring about great challenges vis- à -vis resource exploitation. A large deep opening is sensitive to high ground stress, and is easily inﬂuenced by external interference, which can lead to geologically hazardous occurrences. To investigate the evolution of fracturing and energy in large, deep stopes subjected to dynamic loads, we established a numerical model of a stope in the Gaofeng mine. Using ANSYS / LS-DYNA software, we implemented an implicit solution to initial static stress and an explicit solution for dynamic analysis. Based on our numerical results, we obtained the fracture behavior and energy evolution under coupled static and dynamic loads. To determine the response of ground pressure to mining activity, a 24-channel microseismic monitoring system was designed for the Gaofeng mine based on the numerical analysis.


Introduction
As human activities developed, rock engineering continued to occur underground via the exploitation of deep mines, the creation of hydropower station chambers, the production of underground oil and gas storage caverns, the utilization of underground carbon dioxide storage sites, and so on [1,2]. The high-stress environment of a deep rock mass threatens the safety of workers and machinery. Because deep resources are needed to develop human society, disaster control during the development of deep resources has become a challenge [3,4]. It is well known that the source of stress in deep mines is the superposition of original geo-stress and mining-induced stress [5]. Geo-stress continuously accumulates with the formation of overlying strata, and excavation disturbs the in situ stress in the rock mass. Therefore, the original stress can be redistributed around deep openings and high-energy regions can appear in critical areas due to the stress concentration effect [6].
Ore and rock masses store energy, and when energy is concentrated and disturbed by the appropriate dynamic amplitude, the stress state of those materials exceeds its strength, leading to collapse and rock burst disaster [7]. The main factors that lead to rock disaster are dynamic disturbances and high ground stress in deep rock engineering [8]. Although various new rock breaking technologies have been applied to rock breaking, drilling and blasting methods are still the main methods used in mining due to their high efficiency and easy operation. Because of its high frequency and high energy characteristics, the blasting load is the most significantly incurring factor that induces geotechnical A caving accident occurred in the No. 1 stope at the −179 m level, as shown in Figure 2. The structural parameters of the stope are as follows. It is perpendicular to the strike of the orebody. The width is 16 m and the length, which is the thickness of the orebody, is 37 m. In addition to the long pan, the surrounding rock mass is also broken, and both are disadvantageous for stope stability. The structural parameters of the caving region are as follows: the width is 4 m, the length is 10 m, and the average thickness is 2 m. The burst rocks were mainly in the shapes of slices and lenses. Furthermore, workers in the ramp heard clear cracking sounds and an ejection, which they said was similar to a gunshot. After a field investigation and analysis were performed, it was believed that the main cause of the accident was the long span of the roof without a pillar, which caused a stress redistribution and high geo-stress concentration in the roof; finally, caving was induced by the blasting stress wave from a nearby stope. For safety reasons and as a precaution, the workers cleared loose stones from the surface of the roof, and the cavity was filled quickly.

Boundary Conditions for the Static Solution
The orebody thickness in the study area varied between 20 m to 50 m, and the general axis of the stope was perpendicular to the strike of the orebody. Because of the air drifter's poor efficiency, the length of the caving step along the general axis of the stope was about 6 m. Under these circumstances, stope width and layer height played a significant role in stability and production capacity. Reasonable stope structural parameters were obtained via the Mathews stability graph approach, which provided a guide for the subsequent numerical simulation scheme. Then, A caving accident occurred in the No. 1 stope at the −179 m level, as shown in Figure 2. The structural parameters of the stope are as follows. It is perpendicular to the strike of the orebody. The width is 16 m and the length, which is the thickness of the orebody, is 37 m. In addition to the long pan, the surrounding rock mass is also broken, and both are disadvantageous for stope stability. The structural parameters of the caving region are as follows: the width is 4 m, the length is 10 m, and the average thickness is 2 m. The burst rocks were mainly in the shapes of slices and lenses. Furthermore, workers in the ramp heard clear cracking sounds and an ejection, which they said was similar to a gunshot. After a field investigation and analysis were performed, it was believed that the main cause of the accident was the long span of the roof without a pillar, which caused a stress redistribution and high geo-stress concentration in the roof; finally, caving was induced by the blasting stress wave from a nearby stope. For safety reasons and as a precaution, the workers cleared loose stones from the surface of the roof, and the cavity was filled quickly. A caving accident occurred in the No. 1 stope at the −179 m level, as shown in Figure 2. The structural parameters of the stope are as follows. It is perpendicular to the strike of the orebody. The width is 16 m and the length, which is the thickness of the orebody, is 37 m. In addition to the long pan, the surrounding rock mass is also broken, and both are disadvantageous for stope stability. The structural parameters of the caving region are as follows: the width is 4 m, the length is 10 m, and the average thickness is 2 m. The burst rocks were mainly in the shapes of slices and lenses. Furthermore, workers in the ramp heard clear cracking sounds and an ejection, which they said was similar to a gunshot. After a field investigation and analysis were performed, it was believed that the main cause of the accident was the long span of the roof without a pillar, which caused a stress redistribution and high geo-stress concentration in the roof; finally, caving was induced by the blasting stress wave from a nearby stope. For safety reasons and as a precaution, the workers cleared loose stones from the surface of the roof, and the cavity was filled quickly.

Boundary Conditions for the Static Solution
The orebody thickness in the study area varied between 20 m to 50 m, and the general axis of the stope was perpendicular to the strike of the orebody. Because of the air drifter's poor efficiency, the length of the caving step along the general axis of the stope was about 6 m. Under these circumstances, stope width and layer height played a significant role in stability and production capacity. Reasonable stope structural parameters were obtained via the Mathews stability graph approach, which provided a guide for the subsequent numerical simulation scheme. Then,

Boundary Conditions for the Static Solution
The orebody thickness in the study area varied between 20 m to 50 m, and the general axis of the stope was perpendicular to the strike of the orebody. Because of the air drifter's poor efficiency, the length of the caving step along the general axis of the stope was about 6 m. Under these circumstances, stope width and layer height played a significant role in stability and production capacity. Reasonable stope Energies 2020, 13, 6102 4 of 18 structural parameters were obtained via the Mathews stability graph approach, which provided a guide for the subsequent numerical simulation scheme. Then, optimization of the structural parameters of the stope was implemented in the finite difference software FLAC3D. Finally, the optimal structural parameters for the Gaofeng mine were determined, and the width of the first mining step was 8 m [30]. Considering the capacity of the filling system, the layer height remained 4 m.
To reduce calculation costs, the numerical model was simplified to a quasi-plane strain model [31], as shown in Figure 3. The thickness of the numerical model equaled the side length of a single element. The sides of the square model were 100 m, and the model consisted of 111,512 hexahedral elements. The boundary conditions for the implicit solution are shown in Figure 3. A uniform stress load was applied to the left and right sides, as well as the top and bottom of the model. According to the measured geo-stress data at the −200 m level, σ x = 30.60 MPa and σ y = 20.22 MPa.
Energies 2020, 13, x FOR PEER REVIEW 4 of 17 optimization of the structural parameters of the stope was implemented in the finite difference software FLAC3D. Finally, the optimal structural parameters for the Gaofeng mine were determined, and the width of the first mining step was 8 m [30]. Considering the capacity of the filling system, the layer height remained 4 m.
To reduce calculation costs, the numerical model was simplified to a quasi-plane strain model [31], as shown in Figure 3. The thickness of the numerical model equaled the side length of a single element. The sides of the square model were 100 m, and the model consisted of 111,512 hexahedral elements. The boundary conditions for the implicit solution are shown in Figure 3.

Material Parameters
Material models appropriate for rock or concrete such as MAT_PLASTIC_KINEMATIC, MAT_BRITTLE_DAMAGE, MAT_JOHNSON_HOLMQUIST_CONCRETE, MAT_CSCM, and MAT_RHT were included in the LS-DYNA database [32]. Suitable material models were selected according to the requirements for different analyses, such as the strain rate effect, failure, and damage output. Although there were many material models available for rock in the material database, those that supported an algorithm for explicit and implicit conversion were limited. In this study, we needed to obtain the solution of the initial geo-stress. Then, the dynamic response of the rock mass under dynamic disturbance was obtained. The MAT_CSCM (continuous surface cap model) material model is a smooth or continuous surface cap model available for solid elements in LS-DYNA and it supported the explicit and implicit conversion algorithm. The material model included an isotropic constitutive equation, and hardening, damage, and rate effects were taken into account in this constitutive model. This model has been widely used in the dynamic analysis of reinforced concrete and rock.
The yield surface was defined in terms of three stress invariants: 1 J is the first invariant of the stress tensor, ' 2 J is the second invariant of the deviatoric stress tensor, and ' 3 J is the third invariant of the deviatoric stress tensor. All of these stress invariants can be written as:

Material Parameters
Material models appropriate for rock or concrete such as MAT_PLASTIC_KINEMATIC, MAT_BRITTLE_DAMAGE, MAT_JOHNSON_HOLMQUIST_CONCRETE, MAT_CSCM, and MAT_RHT were included in the LS-DYNA database [32]. Suitable material models were selected according to the requirements for different analyses, such as the strain rate effect, failure, and damage output. Although there were many material models available for rock in the material database, those that supported an algorithm for explicit and implicit conversion were limited. In this study, we needed to obtain the solution of the initial geo-stress. Then, the dynamic response of the rock mass under dynamic disturbance was obtained. The MAT_CSCM (continuous surface cap model) material model is a smooth or continuous surface cap model available for solid elements in LS-DYNA and it supported the explicit and implicit conversion algorithm. The material model included an isotropic constitutive equation, and hardening, damage, and rate effects were taken into account in this constitutive model. This model has been widely used in the dynamic analysis of reinforced concrete and rock.
The yield surface was defined in terms of three stress invariants: J 1 is the first invariant of the stress tensor, J 2 is the second invariant of the deviatoric stress tensor, and J 3 is the third invariant of the deviatoric stress tensor. All of these stress invariants can be written as: (1) Energies 2020, 13, 6102 5 of 18 where S ij , S jk , and S ki denote the deviatoric stress tensors; P is pressure. These three invariants and the cap hardening parameter, κ are the basis for the yield function: where F f is the shear failure surface, denotes the Rubin three-invariant reduction factor, and F c represents the hardening cap. The Rubin three-invariant reduction factor determines the strength of the material for any state of stress relative to the strength under triaxial compression (TXC).
When the material is in a state of tension with a low confining pressure, the shear surface is used to determine the strength of the material: The values of α, β, λ, and θ are selected by fitting the model surface to strength measurements from TXC tests conducted on plain concrete cylinders. The strength of the material is modeled by a combination of the cap and shear surfaces in high confining pressure regimes. The isotropic hardening cap is a two-part function that is either unified or an ellipse: where L(κ) can be written as: where R is the cap ellipticity ratio (the ratio of its major to minor axes) and κ 0 is the value of J 1 at the initial intersection of the cap and shear surfaces before hardening occurs (before the cap moves). Strength under torsion (TOR) and triaxial extension (TXE) are modeled as Q 1 F f and Q 2 F f , respectively, where The parameters α 1 , λ 1 , β 1 , θ 1 , α 2 , λ 2 , β 2 , and θ 2 can be determined based on the experimental data from the TOR and TXE tests, respectively, when the pressure is compressive.
In the constitutive model, the damage parameter d is defined via the brittle damage and ductile damage, ranging from zero, for no damage, to one, for complete damage. Brittle damage accumulates under tensile stress that exceeds the brittle damage threshold. Brittle damage accumulation depends on the maximum principal strain, ε max , which can be written as: where τ b is an energy-type term and E is Young's modulus.
Ductile damage accumulates under compressive stress and depends on the total strain component, ε ij , as shown in the following Equations.
Energies 2020, 13, 6102 where τ d is an energy-type term and σ ij is the stress component. This model provides the user with some default parameters for normal strength concrete, and the parameters fit data reflecting unconfined compressive strengths between approximately 20 MPa and 58 MPa. For rock materials with a high strength, according to the experimental data, the user-defined parameters can be used [33].
To ensure that the simulation results fit well with the actual conditions of the stope in the Gaofeng mine, the numerical simulation parameters were checked against mechanical test results. Due to a large number of numerical simulation parameters, the laboratory test results were not sufficiently comprehensive, and some parameters were assigned the default value. The numerical simulation was accurate only when the numerical experimental data fit the laboratory test results well [34].
As the line with the red stars in Figure 4 shows, the experimental stress-strain curve of the ore has a uniaxial compressive strength of 117.84 MPa. After multiple verifications, the numerical results of the single-element test fit the experimental results well, and the final numerical parameters used for the ore are shown in Table 1. Due to the lack of fractures and pores in the numerical model, there was no compaction stage in the numerical curve.

Solution of Geo-Stress
The numerical simulation consisted of two stages: (1) implicit solution of geo-stress; (2) explicit solution of dynamic disturbance by introducing the initial stress state obtained from stage 1. The solution of geo-stress, shown in Figure 5, shows all regions around the stope in the compressive state under the stress boundary condition influence and stope geometry [35]. Figure 5 shows that the area of stress distribution in the vertical direction was larger than in the horizontal direction. Because the stope model in the numerical simulation was a rectangle, the stress concentration appeared in the corner of the stope, and the maximum stress rose to 52.98 MPa. The implicit solution only provided the initial geo-stress for the explicit solution, so the displacement data was deleted by modifying some commands [36].

Solution of Geo-Stress
The numerical simulation consisted of two stages: (1) implicit solution of geo-stress; (2) explicit solution of dynamic disturbance by introducing the initial stress state obtained from stage 1. The solution of geo-stress, shown in Figure 5, shows all regions around the stope in the compressive state under the stress boundary condition influence and stope geometry [35]. Figure 5 shows that the area of stress distribution in the vertical direction was larger than in the horizontal direction. Because the stope model in the numerical simulation was a rectangle, the stress concentration appeared in the corner of the stope, and the maximum stress rose to 52.98 MPa. The implicit solution only provided the initial geo-stress for the explicit solution, so the displacement data was deleted by modifying some commands [36].

Solution of Geo-Stress
The numerical simulation consisted of two stages: (1) implicit solution of geo-stress; (2) explicit solution of dynamic disturbance by introducing the initial stress state obtained from stage 1. The solution of geo-stress, shown in Figure 5, shows all regions around the stope in the compressive state under the stress boundary condition influence and stope geometry [35]. Figure 5 shows that the area of stress distribution in the vertical direction was larger than in the horizontal direction. Because the stope model in the numerical simulation was a rectangle, the stress concentration appeared in the corner of the stope, and the maximum stress rose to 52.98 MPa. The implicit solution only provided the initial geo-stress for the explicit solution, so the displacement data was deleted by modifying some commands [36].

Boundary Conditions for the Dynamic Solution
Influenced by the medium's propagation properties, the structure's geometry, and charge amount, among other factors, underground blasting vibration showed a feature of high uncertainty. Therefore, as shown in Figure 6, the real blasting load was usually expressed with a simplified equivalent blast loading curve [37]. A triangular loading function P d (t), which includes a rising stage and descending stage, has been widely used in the numerical simulation of blasting, which can be expressed as follows: where P max is the peak value of the dynamic disturbance and t e is the total duration of the triangular load. The rising stage duration varies from 0 to t r , and the peak value of the dynamic disturbance decreases from P max at t r to 0 at t e . For the triangular loading function, the rising stage indicates loading, and the descending stage indicates unloading.
where max P is the peak value of the dynamic disturbance and e t is the total duration of the triangular load. The rising stage duration varies from 0 to r t , and the peak value of the dynamic disturbance decreases from max P at r t to 0 at e t . For the triangular loading function, the rising stage indicates loading, and the descending stage indicates unloading. The distribution of openings in the underground mine is complicated; therefore, a dynamic disturbance could come from any direction for a certain stope. Because of the different charge amounts and distances to the explosive source, the amplitude of a dynamic load varied over a wide range. Considering the symmetry of the rectangular model, a dynamic load with different amplitudes was applied in the top and left boundaries of the model, as shown in Figures 7 and 8.  The distribution of openings in the underground mine is complicated; therefore, a dynamic disturbance could come from any direction for a certain stope. Because of the different charge amounts and distances to the explosive source, the amplitude of a dynamic load varied over a wide range. Considering the symmetry of the rectangular model, a dynamic load with different amplitudes was applied in the top and left boundaries of the model, as shown in Figures 7 and 8.

Energy Evolution around a Large Deep Stope
Notably, deep unexcavated rock masses can store energy. When the dynamic load propagates in a rock mass with a free surface, the rock mass experiences repeated loading and unloading [38]. If the stress exceeds the strength of the rock mass, the deep opening undergoes instantaneous collapse or stripping. Indicators of damage (e.g., stiffness, strength, strain energy density (SED), energy release rate, etc.) are used for risk assessment in deep engineering. To identify the transition and transformation of the energy in a rock mass and use the results and conclusions of laboratory tests to guide field engineering practice, the index of strain energy density is often used to study the rules of energy evolution in deep mines because it eliminates the influence of the size effect. As an energy analysis index without volume dependence, the strain energy density perfectly fits the quasi-3D model in this study [39] and can be written as: where W is the strain energy density, E is Young's modulus, and  is Poisson's ratio. 1  , 2  , and 3  denote the first, second, and third principal stresses, respectively.

SED Evolution around the Stope
Model I The storage of energy within a rock mass was expressed via the implicit solution for geo-stress, and a stable stress state formed after the implicit solution was applicable. In the process of finding an explicit dynamic solution, the dynamic disturbance on the top boundary with a load amplitude of 20

Energy Evolution around a Large Deep Stope
Notably, deep unexcavated rock masses can store energy. When the dynamic load propagates in a rock mass with a free surface, the rock mass experiences repeated loading and unloading [38]. If the stress exceeds the strength of the rock mass, the deep opening undergoes instantaneous collapse or stripping. Indicators of damage (e.g., stiffness, strength, strain energy density (SED), energy release rate, etc.) are used for risk assessment in deep engineering. To identify the transition and transformation of the energy in a rock mass and use the results and conclusions of laboratory tests to guide field engineering practice, the index of strain energy density is often used to study the rules of energy evolution in deep mines because it eliminates the influence of the size effect. As an energy analysis index without volume dependence, the strain energy density perfectly fits the quasi-3D model in this study [39] and can be written as: where W is the strain energy density, E is Young's modulus, and ν is Poisson's ratio. σ 1 , σ 2 , and σ 3 denote the first, second, and third principal stresses, respectively.

SED Evolution around the Stope
Model I The storage of energy within a rock mass was expressed via the implicit solution for geo-stress, and a stable stress state formed after the implicit solution was applicable. In the process of finding an explicit dynamic solution, the dynamic disturbance on the top boundary with a load amplitude of 20 MPa was applied. When the dynamic disturbance propagated from the top boundary, the distribution of SED in the undisturbed region reflected the initial state, which was the same as the implicit solution results, as shown in Figure 9a. Meanwhile, due to the stress concentration effect, the energy concentration was also significant at the sharp corner of the stope, where the SED was at its maximum value. Low SED was distributed on the stope's sidewall. To track the evolution characteristics of the SED's free surface under dynamic disturbance, the positions of the four maximum values (P max1 , P max2 , P max3 , and P max4 ) and the four minimum values (P min1 , P min2 , P min3 , and P min4 ) were tracked. the middle of sidewalls, and the elements with a maximum SED returned to the four corners, a ure 9f shows. To make the contour lines of strain energy density at different moments comparabl e legend intervals were equal. Compared with the initial undisturbed state, there was an obviou istribution of SED under dynamic disturbance, and most regions showed energy release feature When the dynamic disturbance propagated to the vicinity of the stope, the region of the minimum strain energy density deviated from the sidewall and moved upward, which indicated that the rock mass elements in the area exhibited an unloading behavior and that the element energy was released during the dynamic load propagation. Furthermore, the high-SED area near the roof migrated to both ends of the floor, as Figure 9c shows. As the dynamic load continuously propagated downward, the four minimum SEDs fluctuated near the horizontal axis of the stope. The four corners of the rectangular cavity were still in a high-SED state under the influence of the stress concentration effect. The maximum SED was maintained on both sides of the floor when the dynamic load propagated from the roof to the floor. After the dynamic load passed through the floor of the cavity, the elements near the free surface underwent a self-regulation process. During that time, the four elements of minimum SED fluctuated around the horizontal axis of the stope, and the fluctuation decreased. When the dynamic disturbance ended, the elements with a minimum SED were observed at the middle of sidewalls, and the elements with a maximum SED returned to the four corners, as Figure 9f shows. To make the contour lines of strain energy density at different moments comparable, the legend intervals were equal. Compared with the initial undisturbed state, there was an obvious redistribution of SED under dynamic disturbance, and most regions showed energy release features.

Model II
Via the same research method mentioned above, the positions of the four maximum values (P max1 , P max2 , P max3 , and P max4 ) and four minimum values (P min1 , P min2 , P min3 , and P min4 ) for model II were also tracked, as described in this section. A dynamic load amplitude of 20 MPa was applied to the left boundary during the explicit solution. Before the dynamic stress wave reached the stope's vicinity, the SED around the stope maintained its initial state. When the dynamic disturbance began acting on the stope's left sidewall, the four elements of the minimum SED were distributed at the cavity's left sidewall and then migrated to the right with the movement of the dynamic load. As the dynamic stress wave propagated from the left sidewall to the right sidewall, the four elements of the maximum SED remained at the right two corners, as shown in Figure 10c-e. When the left dynamic disturbance ended, the four elements of the minimum SED were observed at the vertical axis of the model, and the four elements with the maximum SED returned to the four corners. ess wave propagated from the left sidewall to the right sidewall, the four elements of the maximum D remained at the right two corners, as shown in Figure 10c-e. When the left dynamic disturbanc ded, the four elements of the minimum SED were observed at the vertical axis of the model, an e four elements with the maximum SED returned to the four corners.
The legend of the contour lines in Figure 10 is the same as that in Figure 9, and the distributio ge of the SED contour line in Figure 10f is smaller than that in Figure 9f. The left dynamic loa d a stronger action of energy release than the top dynamic load. The stress environment aroun e stope was sensitive to the left-side dynamic disturbance.  Figure 11. To gain insight into the change in SED e values of the initial SED minus the final SED are listed in Figure 12. As Figure 11 shows, in an se, the dynamic response of the elements always started with SED decreasing, which indicates tha ring the dynamic loading of the energy storage region, the first response was to unload. For mod 1 P2, P3, and P4 showed a decreasing trend after the dynamic disturbance ended, and the maximum crease in SED was 52.50 KJ/m 3 , as Figure 12a shows. The elements of P5 on the sidewall absorbe ergy from the dynamic load, and their SEDs increased by 7.09 KJ/m 3 . Although the elements of P d P4 in the model had nearly equal initial SEDs, the changes in their SEDs were markedly differen er dynamic disturbances because the dynamic energy gradually weakened as the dynamic loa read through the rock mass. In model II, the SEDs of those five elements exhibited a simila ctuation law with that of the five elements in model I; the difference was that the elements in th fferent models exhibited different SED changes.
The different dynamic loading directions had different influences on the same position. The fir ample was that P1 and P3 in model I had the same positions as P3 and P5 in model II, respectivel The legend of the contour lines in Figure 10 is the same as that in Figure 9, and the distribution range of the SED contour line in Figure 10f is smaller than that in Figure 9f. The left dynamic load had a stronger action of energy release than the top dynamic load. The stress environment around the stope was sensitive to the left-side dynamic disturbance.

Time History of SED
Given the symmetry of the stope and dynamic loading boundary, five elements of each model were chosen to extract the data of the SED time history; the positions of the elements are shown in Figures 7 and 8, and the time history is shown in Figure 11. To gain insight into the change in SED, the values of the initial SED minus the final SED are listed in Figure 12. As Figure 11 shows, in any case, the dynamic response of the elements always started with SED decreasing, which indicates that, during the dynamic loading of the energy storage region, the first response was to unload. For model I, P 1 P 2 , P 3 , and P 4 showed a decreasing trend after the dynamic disturbance ended, and the maximum decrease in SED was 52.50 KJ/m 3 , as Figure 12a shows. The elements of P 5 on the sidewall absorbed energy from the dynamic load, and their SEDs increased by 7.09 KJ/m 3 . Although the elements of P 2 and P 4 in the model had nearly equal initial SEDs, the changes in their SEDs were markedly different after dynamic disturbances because the dynamic energy gradually weakened as the dynamic load spread through the rock mass. In model II, the SEDs of those five elements exhibited a similar fluctuation law with that of the five elements in model I; the difference was that the elements in the different models exhibited different SED changes.
Energies 2020, 13, x FOR PEER REVIEW 12 of 17 in the high initial energy regions, under the circumstance of a dynamic disturbance amplitude of 20 MPa, the initial stress state played a more significant role in the energy evolution of elements than the role of the dynamic disturbance.  Under all circumstances, the monitored SED results exhibited these stages during the dynamic disturbance, declining from sharply fluctuating to stable. In conclusion, elements on the central axis of the stope that were perpendicular to the dynamic loading boundary absorbed energy from the dynamic disturbance, and the remaining elements showed energy release characteristics.

Fracturing under Different Dynamic Stress Amplitudes
When external energy was input into a deep rock mass, energy was stored in the form of elastic energy. As the rock mass deformed under loading, the energy was consumed in the form of crushed energy and friction energy. If there was stored energy remaining in the rock, dynamic disaster occurred [40]. Because the uncertainty in the dynamic disturbance is high, to investigate the fracturing of large stopes, dynamic loads of 30 MPa, 40 MPa, 50 MPa, and 60 MPa were applied on the dynamic boundary. It is known that the compressive strength of rock is much greater than its tensile strength, and dynamic tensile failure is the primary cause for spalling around underground openings; therefore, the tensile failure zone was tracked based on the CSCM model [41]. All the hexahedral elements used in the numerical model were the same size, so the number of failure elements reflected the failure characteristics of the stope to some extent. The number of failure elements was automatically extracted from the numerical output file by a computer script. To illustrate the failure zone, as shown in Figures 13 and 14, the failure elements are shown in red, and the intact elements are shown in blue. Energies 2020, 13, x FOR PEER REVIEW 12 of 17 in the high initial energy regions, under the circumstance of a dynamic disturbance amplitude of 20 MPa, the initial stress state played a more significant role in the energy evolution of elements than the role of the dynamic disturbance.  Under all circumstances, the monitored SED results exhibited these stages during the dynamic disturbance, declining from sharply fluctuating to stable. In conclusion, elements on the central axis of the stope that were perpendicular to the dynamic loading boundary absorbed energy from the dynamic disturbance, and the remaining elements showed energy release characteristics.

Fracturing under Different Dynamic Stress Amplitudes
When external energy was input into a deep rock mass, energy was stored in the form of elastic energy. As the rock mass deformed under loading, the energy was consumed in the form of crushed energy and friction energy. If there was stored energy remaining in the rock, dynamic disaster occurred [40]. Because the uncertainty in the dynamic disturbance is high, to investigate the fracturing of large stopes, dynamic loads of 30 MPa, 40 MPa, 50 MPa, and 60 MPa were applied on the dynamic boundary. It is known that the compressive strength of rock is much greater than its tensile strength, and dynamic tensile failure is the primary cause for spalling around underground openings; therefore, the tensile failure zone was tracked based on the CSCM model [41]. All the hexahedral elements used in the numerical model were the same size, so the number of failure elements reflected the failure characteristics of the stope to some extent. The number of failure elements was automatically extracted from the numerical output file by a computer script. To illustrate the failure zone, as shown in Figures 13 and 14, the failure elements are shown in red, and the intact elements are shown in blue. The different dynamic loading directions had different influences on the same position. The first example was that P 1 and P 3 in model I had the same positions as P 3 and P 5 in model II, respectively, and all of these elements exhibited relatively low initial SED. The SED of P 1 in model I decreased by 2.25 KJ/m 3 , and the SED of P 3 in model II decreased by 9.94 KJ/m 3 . The SED of P 3 in model I decreased by 5.20 KJ/m 3 , and the SED of P 5 in model II increased by 6.27 KJ/m 3 , as shown in Figure 12. This result indicates that the dynamic load boundary played a vital role in the change in energy of elements with a low initial energy. As another example, both P 2 in model I and P 4 in model II were located in the top right of the stope, and they had a high initial SED. The final SED of P 2 in model I and P 4 in model II decreased by 49.89 KJ/m 3 and 62.14 KJ/m 3 , respectively. This result indicates that in the high initial energy regions, under the circumstance of a dynamic disturbance amplitude of 20 MPa, the initial stress state played a more significant role in the energy evolution of elements than the role of the dynamic disturbance.
Under all circumstances, the monitored SED results exhibited these stages during the dynamic disturbance, declining from sharply fluctuating to stable. In conclusion, elements on the central axis of the stope that were perpendicular to the dynamic loading boundary absorbed energy from the dynamic disturbance, and the remaining elements showed energy release characteristics.

Fracturing under Different Dynamic Stress Amplitudes
When external energy was input into a deep rock mass, energy was stored in the form of elastic energy. As the rock mass deformed under loading, the energy was consumed in the form of crushed energy and friction energy. If there was stored energy remaining in the rock, dynamic disaster occurred [40]. Because the uncertainty in the dynamic disturbance is high, to investigate the fracturing of large stopes, dynamic loads of 30 MPa, 40 MPa, 50 MPa, and 60 MPa were applied on the dynamic boundary. It is known that the compressive strength of rock is much greater than its tensile strength, and dynamic tensile failure is the primary cause for spalling around underground openings; therefore, the tensile failure zone was tracked based on the CSCM model [41]. All the hexahedral elements used in the numerical model were the same size, so the number of failure elements reflected the failure characteristics of the stope to some extent. The number of failure elements was automatically extracted from the numerical output file by a computer script. To illustrate the failure zone, as shown in Figures 13 and 14  When the stress wave propagated to the free surface, it was reflected as a tensile wave and caused fracturing. Due to the geometry of the stope, the roof was a larger free surface than the cavity's sidewall. Therefore, under a certain dynamic disturbance amplitude, the number of failure elements in the roof was always larger than the sidewall. The maximum number of failure elements induced by the upper and left dynamic disturbances were 1553 and 616, respectively.
We found that the spatial pattern of the failure elements under the low dynamic disturbance amplitude had little difference under the influence of the initial geo-stress, and the distributions of failure elements were geometrically similar. The failure zone under the top dynamic disturbance  When the stress wave propagated to the free surface, it was reflected as a tensile wave and caused fracturing. Due to the geometry of the stope, the roof was a larger free surface than the cavity's sidewall. Therefore, under a certain dynamic disturbance amplitude, the number of failure elements in the roof was always larger than the sidewall. The maximum number of failure elements induced by the upper and left dynamic disturbances were 1553 and 616, respectively.
We found that the spatial pattern of the failure elements under the low dynamic disturbance amplitude had little difference under the influence of the initial geo-stress, and the distributions of failure elements were geometrically similar. The failure zone under the top dynamic disturbance When the stress wave propagated to the free surface, it was reflected as a tensile wave and caused fracturing. Due to the geometry of the stope, the roof was a larger free surface than the cavity's sidewall. Therefore, under a certain dynamic disturbance amplitude, the number of failure elements in the roof was always larger than the sidewall. The maximum number of failure elements induced by the upper and left dynamic disturbances were 1553 and 616, respectively. We found that the spatial pattern of the failure elements under the low dynamic disturbance amplitude had little difference under the influence of the initial geo-stress, and the distributions of failure elements were geometrically similar. The failure zone under the top dynamic disturbance consisted of the main fracturing zone in the roof and horizontal cracks near the side walls, and the failure zone under the left dynamic disturbance consisted of the fracturing zone in the left side wall and vertical cracks near the roof and floor. The dynamic stress wave strength steadily reduced when spreading distance increased. Therefore, no failure occurred in the floor of model I and the right side wall of model II. With the increase in the top dynamic load, the failure zone in the roof became larger and interconnected with the horizontal cracks. Under the conditions of the dynamic loads from the left boundary, the failure zone was distributed on the left of the sidewall, and vertical cracks near the roof and floor. When the dynamic load exceeded 30 MPa in model II, the fracturing zone extended from the left sidewall to the roof and floor. With the continuous increase in the left dynamic load, vertical cracks near the roof and floor were gradually interconnected.
As the amplitude of the dynamic disturbance increased, the number of failure elements also increased, and the changes in the fracturing zone geometry were obvious. This finding indicates that the fracturing was governed by the dynamic load when the stope was subjected to a high dynamic load.

Discussion
Mining operations under high stress conditions at depth are often faced with rock bursts and other geological disasters, both of which challenge the safety and efficiency of mining. From the above research results, a large area of the fracturing zone in the Gaofeng mine was formed under strong dynamic disturbance, which will eventually cause huge losses of life and property. Furthermore, when the mining exploration went to the −250 m level, the ground pressure significantly changed due to mining activity. In engineering practice, pretreatment measures for eliminating potential hazards is required when workers and machinery enter a large opening. However, due to the complexity of the mining environment, geological disasters in deep mines exhibit high uncertainty. It is difficult to manually realize on-site workplace hazard screening. In addition, such a screening is limited to disaster-inducing factors in mining operation design. From the numerical results presented in this study, we observed that the large size introduced high risks to stabilize the deep opening under dynamic disturbance. In addition, reasonable stoping layouts also had significant effects on deep engineering stability [42,43]. The stress redistribution induced by an unreasonable mining scheme could partly contribute to stress concentration in the sharp corner of openings, and rock burst may be subjected to a very low dynamic disturbance [44,45]. In such cases, in order to keep the safety of the working environment, the maximum detonation charge must be reduced, which is not acceptable due to the low mine production capacity. Although some geotechnics have been proven to be effective to prevent rock burst, lower personnel exposure is the first line of defense for a mining operation [46].
It is known that the redistribution of geo-stress induced by mining operations may lead to the initiation, development and propagation of a microcrack [47]. All signals of elastic waves generated by rock deformation and failure can be captured by sensors and recorded with a data acquisition instrument [48]. As an effective early warning approach, microseismic monitoring techniques have become a common means to monitor the stability of engineering in deep mines [49,50]. To determine the response of the geo-stress in mining activities, a 24-channel microseismic monitoring system ( Figure 15) was designed for the Gaofeng mine and is now under construction.
initiation, development and propagation of a microcrack [47]. All signals of elastic waves generated by rock deformation and failure can be captured by sensors and recorded with a data acquisition instrument [48]. As an effective early warning approach, microseismic monitoring techniques have become a common means to monitor the stability of engineering in deep mines [49,50]. To determine the response of the geo-stress in mining activities, a 24-channel microseismic monitoring system ( Figure 15) was designed for the Gaofeng mine and is now under construction. In total, 16 microseismic sensors were installed at the −200 m and −175 m levels, including 12 uniaxial geophone sensors and four triaxial geophone sensors. The data acquisition instruments were arranged in a special chamber at the −200 m level. An optical cable was used for data transfer, and In total, 16 microseismic sensors were installed at the −200 m and −175 m levels, including 12 uniaxial geophone sensors and four triaxial geophone sensors. The data acquisition instruments were arranged in a special chamber at the −200 m level. An optical cable was used for data transfer, and real-time data viewing was realized in the control center at the 690 m level. The microseismic data was also sent to the remote console at Central South University via wireless communication technology. If the fiber optic system was out of order, the technician still collected data at the −200 m level. A real-time microseismic monitoring system is expected to help in disaster early warning.

Conclusions
In this study, a numerical model for a large deep underground opening (8 m × 4 m) was built based on the optimization of a stope structural parameter. An implicit solution of the initial static stress and an explicit solution for dynamic analysis were implemented into ANSYS/LS-DYNA software. Combining the SED, fracturing zone extent, and number of failure elements, the numerical investigation of fracturing and energy evolution of a large, deep underground opening subjected to dynamic stress was presented. The main conclusions are as follows: (1) Considering energy distribution, the main dynamic response of the stored energy in a rock mass under dynamic disturbance is energy release. When the dynamic disturbance finished, the distribution range of the SED decreased obviously. When the direction of dynamic loading corresponded to the maximum principal stress, the effect of energy release was even more significant.
(2) In every case, the dynamic response of the monitored elements always began by decreasing SED, which indicated that when the dynamic loading affected the energy storage region, the first response was to unload. The elements far from the loading boundary absorbed dynamic energy, which led to their SEDs increased.
(3) For large-sized deep rock engineering, the main failure zone is governed by the size of free surface. In a certain initial stress environment, the geo-stress plays a significant role in the fracturing behavior of rock mass. When the amplitude of the dynamic disturbance is relatively low, the initial geo-stress dominates the fracturing of the rock mass. When the peak value of the dynamic load is high, the number of rock mass failure elements increases, and the failure mode is governed by the dynamic disturbance.
(4) In engineering practice, pretreatment measures for eliminating potential hazards are required when workers and machinery enter a large underground opening. To realize long-term effective ground pressure control in the Gaofeng mine, an early warning of surrounding rock mass deformation and failure should be provided in a timely manner.