Numerical Modeling on Blasting Stress Wave in Interbedding Rheological Rockmass for the Stability of the Main Shaft of Mine

High rheological property and strong mining disturbance are prominent conditions in deep rock projects, and often present a coupled process to induce rockbursts, collapse and land subsidence. This paper aims to investigate the effect of surrounding rheological rockmass on the main shaft of mine with interbedding structure under coupled blasting stress wave condition. Based on elastic damage theory and the constitutive equation of rheological rockmass, considering the total strain rate in tensor form, a double-rock model was established. The model was then validated by comparing the numerical simulations with the test results. Matlab was used to develop the COMSOL software, allowing numerical simulation of the failure of shaft rock sequentially subjected to this complex loading path. The results from the rheology-impact numerical simulations show that the double-rock rheology was greatly affected by the dynamic disturbance. Under high creep stress and constant impact stress wave, the last failure of the double-rock with higher creep stress is more severe than that of double-rock with lower creep stress. The numerical simulation of the shaft in situ stress was used to predict the collapse of the rheological rockmass with interbedding structure. The preliminary results suggest that one contribution to this phenomenon is likely related to irreversible damage in the rock during the creep-disturbance process. Dynamic disturbance also promoted the failure of the rock near the interbedding structure, but also increased the strain and stress. This may reflect the initial compaction and elastic phase the rock near the interbedding structure related to the large compressive strength of the hard rock during the dynamic disturbance. The numerical results indicate that a circle-shaped spalling damage zone is prone to develop around the shaft with increased time. We discuss the reason for the decreased stress on point A2 near the interbedding structure of the shaft by visco-elastic theory. The results clearly showed insufficient stability of surrounding rockmass, thus the initial design of shaft is not reasonable. This study has important referential significance for main shaft design for similar mines.


INTRODUCTION
The rheological properties of interbedding rock can cause deformation and failure of shaft wall rock. However, only the effects of instantaneous deformation convergence of the shaft are typically considered in shaft construction design. The shaft may experience a series of deformation failure phenomena such as significant convergence, deviation, cracking, squeezing spalling, and sag subsidence under long-term action of crustal stress (Bell and Genske, 2001). When subject to these changes such as Figure 1, the surrounding rock of shaft can suffer a large-scale crack and undergo transverse shear dislocation under the influences of blasting operation and mining disturbance, which can compromise mine safety and severely threaten the safety of operating personnel. In fact, the instability of the interbedding rock can trigger catastrophic landslide during the strong seismic force such as the 2008 Wenchuan earthquake (Cui et al., 2021). Sometimes there are liquefaction or subsidence on the multilayer geotechnical conditions (Li et al., 2022).
The current theoretical models and numerical simulation methods primarily establish a static model to analyze the longterm stability of the shaft wall rock but do not consider any longterm effects of the excavation stress wave disturbance on the rheological behavior of wall rock. In Reference (Chen et al., 2010), a three-dimensional finite element model was established to analyze the force and deformation characteristics of the wall rock of a long and large vertical shaft lined with supporting concrete based on the ideal elastoplastic material constitutive relationship and the Drucker-Prager yield principle. In Reference (Liu et al., 2020), the criteria for vertical shaft stability and the limit depth expression were derived to investigate the effects of differences in tension-compression strength of rock and intermediate principal stress on the limit depth of a vertical shaft. In Reference (Xie et al., 2019), The finite element method was used to perform numerical simulation of the effects of dynamic excavation of roadways in the middle section and large chambers in Jinchuan No. 3 mining area on the stability of the main and auxiliary shafts and analyze the effects of different excavation plans on the displacement of wall rock surrounding the main shaft. Sainoki et al. (2017) simulated the creep behavior of weak rock mass, with the validation of tertiary creep and the determination of model input parameters from field measurements. Barla et al. (Debernardi and Barla, 2009) used a viscoelastic-plastic model, an elastic-viscoplastic model, and an elastic-plastic-viscoplastic mode to analyze the squeezing problems of large size tunnels in rock masses of very poor quality. Sharifzadeh et al. (2013) simulated the timedependent behavior of the tunnel surrounding ground deformation in weak rocks considering the Burger-creep viscoplastic model. Fahimifar et al. (2010) obtained an analytical solution to predict time-dependent deformation of the excavation of circular tunnels in a visco-elastic Burger model under hydrostatic stress field. Zhang et al. (Zhang et al., 2016) studied the time-dependent behavior and stability assessment of double tunnels based on thermodynamics with internal state variable theory. Ma et al. (Ma and Liu, 2022) successfully used 3D discontinuous deformation analysis (DDA) to build a rockfalls model to investigate the movement characteristics of the collapsing rocks before the landslide. Zhou et al. (2021) first utilized the moment tensor inversion to present the fractures network with a combination of microseismic event locations. This is a great contribution to the stability of pillar in the underground mining.
It has long been accepted that rheological behavior of rock is one of the important causes of rock mass engineering deformation and instability. However, under deep mining conditions, the soft wall rock deformation arising from blasting vibration of hard rock mine may have different effects with time (Sun, 1997). Gao et al. (Gao et al., 2007), Fu Zhiliang et al. (Fu et al., 2010), and Zhang et al. (Zhang, 2010) have addressed the disturbance effects of dynamic disturbance on the rheological process of rock. Gao et al. (Gao et al., 2007) proposed the concept of rheological disturbance effects of rock, developed a testing machine to assess rheological disturbance effects of rock, used an experimental approach to assess the effects of dynamic disturbance on rock creep, and proposed a three-dimensional creep damage model that considered disturbance effects. Zhang et al. (Zhang, 2010) theoretically analyzed the time delay of deformation stability of plate or beam structures based on viscoelasticity theory and qualitatively defined the time delay of rockburst. Xu Zenghe et al. (Lyu et al., 2020) discussed the occurrence conditions and retardation time of jamb rockburst under rheological rock stratum, provided a theoretical explanation of the rockburst hysteresis, and derived the expression for rockburst hysteresis time.
This study focused on the stability of rheological interbedding rock of the main vertical shaft in the middle section of Xincheng Gold Mine Ⅺ# ore body-380 of Shandong Gold Mining Co., Ltd. Based on rheological damage theory, MATLAB programming was utilized to characterize and quantify the nonuniform damage of rock in the environment of COMSOL Multiphysics to obtain the numerical solution for theoretical models. The numerical solution and the test solution were compared to verify the rheological model. The on-site rock parameters and disturbance conditions were employed to analyze the variation rules of stress and the strain of the rheological interbedding rock surrounding the main vertical shaft under the joint action of the crustal stress and mining stress wave.

Numerical Model Setup
The rheology of rock can be characterized with a damage-based model, where the conservation equations for mass and momentum are derived on the macroscopic scale. All variables are averaged over the representative elementary volume (REV) of the rock.

Mechanical Equilibrium and Damage Evolution Equation
The rock is assumed to be elastic, with the constitutive relationship defined by a generalized Hooke's law. Here, a modified Navier equation, in terms of displacement under a combination of applied loading is expressed as where u i (i x, y, z) is the displacement (m), G is the shear modulus (Pa), ] is the drained Poisson's ratio, and F i represents the components of the net body force (N·m −3 ) in the i-coordinate. (i x, y, z) As illustrated in Figure 2, the damage of rock under tension or shear stress is initiated when the state of stress satisfies the maximum tensile stress criterion or the Mohr-Coulomb criterion, respectively, as expressed by Zhou Jingren et al. (Zhou et al., 2018).
where f t0 and f c0 are the uniaxial tensile and compressive strengths (Pa), respectively, θ is the internal frictional angle, and F 1 and F 2 are two damage threshold functions. According to elastic damage theory, the elastic modulus of an element degrades monotonically with the evolution of damage, and the elastic modulus of damaged material is expressed as follows: where D represents the damage variable, and E and E 0 are the elastic modules of the damaged and the undamaged material (Pa), respectively. In this kind of numerical simulation, each element and its damage is assumed isotropic, so the E, E 0 , and D variables are all scalars. Under stress conditions, the tensile stress criterion is applied preferentially. According to Figure 2, the damage variable can be calculated as where ε t0 and ε c0 are the maximum tensile principal strain and maximum compressive principal strain, respectively, when damage occurs, and n is a constitutive coefficient and is equal to 2.0. The damage variable calculated with Eq. 4 ranges from 0 to 1.0 for all kinds of damage. It should be noted that, in the numerical implementation of Eq. 4, tensile damage is always preferable to shear damage. Thus, the maximum tensile stress criterion is first used to evaluate if there is damage under tension, and then any elements that do not result in damage in tensile mode are then checked for shear damage potential with the Mohr-Coulomb criterion.

Constitutive Law of Rock During Rheological Process
The total strain includes both elastic (ε e ) and creep strain (ε ce ): Where ε t is the total strain, ε e is the elastic strain, and ε ce is the creep strain. The elastic strain is defined by: where σ kk is the first stress tensor invariant, σ ij is the stress tensor, δ ij is the Kronecker function, ] is Poisson's ratio, and E is the elastic modulus.
The creep strain rate is defined by (Su et al., 2013) as follows: where A, n, m 0 , and n 0 is experimentally determined parameters. _ ε ce ij is the creep strain rate, S ij is the deviatoric part of σ ij , and σ e is the effective stress and is defined as Differentiating Eq. 5 with regards to time and considering a constant total strain, results in a total strain rate Combining the above equation with Eqs 6, 7 and Eq. 9, the total strain rate can be described by where _ ε t ij is the equivalent total strain rate tensor, _ σ ij is the equivalent relaxation stress rate tensor, and t is time. In Eq. 10, the stress is given in Pa and the time in hours.

Material Heterogeneity
Geo-material is assumed to be composed of many mesoscopic elements, and in this study the mechanical properties of these heterogeneous geo-materials are assumed to conform to a given Weibull distribution as defined in the following probability density function : where u is the mechanical parameter of the element (such as strength or elastic modulus); the scale parameter u 0 is related to the average of the element parameters; and the parameter m defines the shape of the distribution function. For the properties of the Weibull distribution, a larger value of m implies a more homogeneous material and a smaller value correlates with a more uniform material.

Numerical Model Description for the Interbedding Rock
As shown in Figure 3, granite diameter is 50mm, L50mm long, and sandstone diameter φ is 50mm, L50mm long, the overall double specimen size φ50mm×L100mm, simplified as plane stress model. There is a rotating cylinder boundary on the lower side of the double sample, and the granite, sandstone and the left and right of the sample are free boundaries. First, we apply a constant load P 0 on the upper side of the double sample for creep, and maintain the displacement for 12 h, and then the upper side of the double sample is dynamically disturbed. When loading the stress wave, the model loading is controled according to the time in Figure 3B, and the time increment of each step is 15 µs A total of 300 steps are calculated. The calculation model is divided into 1,256 triangular units and 22,583 nodes. The physical and mechanical parameters used in the numerical model are shown in Table 1.

Numerical and Experimental Results
As shown in Figure 4A, the strain time curve of multi-stage creep granite and sandstone composite sample after dynamic disturbance is obtained through test and numerical simulation. When the strain of the specimen reaches 3,117 με, the stress reached 15.18 MPa, the specimen began to creep, and after 12 h, the final strain increased to 3,212 με. The change is 2.95%. After 12 h, the sandstone and granite composite specimens were dynamically disturbed by drop hammer impact, resulting in the strain from 3212 to 3576 με (increased by 355 με). Therefore, during the first creep dynamic disturbance of the sandstone part of the combined sample, as shown in Figure 5A, uniformly distributed damage appears in the sample. According to the double specimen creep impact disturbance test process, the specimen was subjected to the second impact disturbance at 24 h.
After the second impact disturbance of the sample, the double sample entered the third creep stage, and the strain of the rock sample was 4,274 με at t = 36 h. Similarly, according to the double specimen creep impact disturbance test process, the specimen was subjected to the third impact disturbance at 36 h, and the strain caused by impact increased from 4,274 με Increased to 5,121 με. Then the sample entered accelerated creep. As shown in Figure 5A, Figure 6A, the sandstone sample showed an oblique crack at 30°to the horizontal direction at t = 36.6 h. Figure 4B showed the damage evolution process of the whole creep dynamic disturbance of the double specimen under the creep stress of 17.17 Pa. Firstly, the specimen is loaded at t = 8 min, and the specimen was less damaged. Then, in the creep stage, the damage points at t = 12 h increased. The sample was subjected to the first impact disturbance at 12 h. Under the action of dynamic load, the strain of the sample suddenly increased, resulting in the strain from 3,637 με Increased to 4,449 με (increased by 812 με). Then the sample entered the second creep stage at the creep stress of 17.17 Pa, and the sample strain increased to 4,942 με at t = 24 h. The sample was subjected to the second impact disturbance at t = 24 h, and the stress wave was repeatedly reflected inside the sample, so that the combined tensile stress after superposition was greater than the uniaxial tensile strength of the sample, at t = 24 h +400 μs caused the specimen to break in Figure 6B. The creep stress in Figure 4C was higher than that in Figures 4A,B, so one impact leaded to the failure of double specimens.

SUBSIDENCE PROBLEMS DUE TO BLASTING WORKINGS IN THE MAIN SHAFT OF MINE 3.1 Numerical Model Setup
In the Xincheng Gold Mine, the XI # ore and the I # ore are connected. The XI # ore body reached the development stage, but was not explored. To analyze the effects of blasting disturbance on the long-term stability of the main vertical shaft during exploitation of the XI # ore body, a numerical model was established of exploitation and blasting disturbance of the XI # ore body versus the stability of rheological wall rock of the main shaft.
To establish the model, the main shaft with infinite length was simplified as a plane strain problem. The numerical model in this section has a size of 60 m × 60 m with a 6 m hole at the center, as shown in Figure 7. Considering the effects of soft ore and strong surrounding rock on the creep of rock mass, a hollow cylinder was embedded in the middle of the model. The size of the modelled specimen was maintained for all of the numerical simulations. The numerical model was discretized into 2,582 triangular elements and contains soft ore thickness of D = 0.5 m and overburden depth H = 380 m based on the site-specific geological conditions. This stress fluctuation on the boundary is considered to be the induced stress wave travels into the inhomogeneous rock. In this regard, the soft rocks are repeatedly compressed and the strong rocks are repeatedly stretched.
The static geo-stress condition governs the stress distribution and damage zone around the under-ground excavation. The direction of the maximum horizontal stress σ h, max is NWW -SEE in the Xincheng mine, and the maximum horizontal principal stress σ h, max and the minimum horizontal principal stress σ h, min increased with depth. Linear regression equations of σ h, max and σ h, min versus depth (H) were established by Cai et al. (2000) as σ h,max −0.44 + 0.0592H (12) σ h,min 0.44 + 0.0314H Where H represents the depth, σ h, max is the maximum horizontal principal stress, and σ h, min is the minimum horizontal principal stress. For rheological analysis, the boundary stresses in the vertical and horizontal directions

Practical Estimates of Rock Mass Parameters
The Generalised Hoek-Brown failure criterion for jointed rock masses as described by Hoek and Brown (Hoek and Brown, 1980) is where σ 1 and σ 3 are the maximum and minimum effective stresses at failure, respectively; σ ci is the uniaxial compressive strength of the intact rock pieces; s and α are constants that depend upon the characteristics of the rock mass; and m b is the value of the Hock-Brown constant m i for the rock mass. Lu et al.
The constant m i depends upon the mineralogy, composition, and grain size of the intact rock (Cai et al. 2004). The rock masses around the main shaft are coarse granite in the Xincheng gold mine, thus the value of m i equals 32.7 according to the rock mass geological characteristics listed in Table  2. Stereophotogrammetry is a method to extract the information of an area of interest by constructing a stereo-image from two or more photos (Huang et al., 2019). In this work, we applied stereophotogrammetry to obtain the joint density J v . of 4.9356 m −1 , which would indicate very good surface condition of the bedding planes and joints.
The strength of the rock mass σ cm is defined as where m b is the value of the Hock-Brown constant m i for the rock mass, and m b is calculated using Eq. 15. The Hock-Brown parameters s and α are defined by Eqs 16, 17, respectively. Laboratory tests should be performed to obtain σ ci , and the value of σ ci , is 100 MPa. When applying the Generalised Hoek-Brown criterion to jointed rock masses, the uniaxial tensile strength is lower than that obtained by the Mohr-Coulomb criterion (Yang, 2018). The uniaxial tensile strength (positive for tension) is then expressed as follows: The deformation modulus of a rock mass (E m ) is one of the significant parameters required to build numerical models for most rock engineering projects, such as open pit mining and tunnel excavations (Hoek and Brown, 1997). Based upon GSI in poor quality rock masses, the following modification to Serafim and Pereira's equation (Shen et al., 2012) is obtained for σ ci < 100MPa: The equivalent angles of friction and cohesive strengths are given by (Hussian et al., 2020) as follows: where, σ′ 3n is the ratio of σ 3max , the upper limit value of the minimum principal stress of the rock mass, to σ ci , the uniaxial compressive strength of intact rock: For a roadway project, the relationship between σ 3max , the upper limit value of the minimum principal stress of rock mass, and σ cm , the uniaxial compressive strength of rock mass (Hoek and Brown, 2019) can be described as: where, σ cm is the compressive strength of rock mass, MPa; γ is the weight of rock mass, kg/m 3 ; H t is the tunnel burial depth, m. When the horizontal stress exceeds the vertical stress, γH t is replaced by the horizontal stress. The generalized Hoek-Brown strength criteria were introduced to the geological strength indexes GSI for computing the strength parameters of rock mass, as shown in Table 4. The actual natural rock was nonuniform quasi-brittle material. Reference  introduced the Weibull statistical distribution function to describe a nonuniform distribution of the mesoscopic mechanical properties of rock. The values of material parameters were determined according to certain given Weibull distribution, and specific parameters are presented in Table 5. The soft ore material constants in the creep constitutive equation A 3 × 10 − 8s −1 are m 0 0.6 and n 0 1.2. Figure 8A presents the determined damage zone of the rheological interbedding structure surrounding the vertical shaft. The grey part represents tensile damage to the rock and the black part denotes intact rock. The damage values are continuously distributed from -1 (tensile damage) to 1 (shear damage). As can be seen in Figure 8A, there is nearly no damage to the sample interior for initial crustal stress σ x 10MPa or σ y 22MPa(t 1h). With the increase in time (t 40h), creep strain occurs and damage can sporadically occur in the concrete lined support of the vertical shaft. At constant crustal stress level but further increasing time (t 181h), the sporadically distributed damage points can aggregate, but no cracks occur at the periphery of the hole. Figure 8B present the curves for creep strain and damage law at Point A2 of the interbedding structure rock. Point A 2 of the wall rock is close to the ore area of the interior wall of the shaft, so the ore for the interior rock of the shaft restricts the deformation of the wall rock. When the strain of the test piece reaches 723 με, the test piece starts to experience creep. At 40 h, the strain of Point A 2 of the wall rock reached 805 με. At 181 h, the wall rock experienced creep due to rock damage and eventually the strain increased to 945 με.

Scenario II:damage Evolution Under Coupled Rheology-Dynamic Disturbance Conditions
As shown in Figure 9A, at t = 181 h and 210 μs, the stress along the radial direction for the heterogeneous rock model with a Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 930013 8 homogeneity index of 3.0 fluctuates near that calculated for the homogeneous elastic model, with no tensile stress induced around the interbedding structure. However, together with the propagation of stress wave, at t = 181 h and 1,050 μs, the tensile stress and the associated tensile damage ( Figure 9C) around the shaft. Figure 9G shows the strain and stress course curves of the A 2 points of the interbedding structure rock under dynamic disturbance. At 181 h, a blasting operation was performed near the shaft roadway. Both the stress and strain increased and the rock was in an elastic stage, indicating that no damage increased to the rock close to the soft ore for the shaft. As shown in Figure 9D, Figure 8E, the bursting stress wave was far away from the interbedding structure rock of the vertical shaft. After 840 μs, the stress changed from 45 to 104 MPa (increase of 59 MPa) and the strain increased by from 356 με to 1,309 με at t = 181 h and 1,050 μs under the influence of the blasting stress wave. But after t = 181 h and 1,050 μs, the observed increase in strain and decrease in stress may result from two aspects. On one hand, the rock sample is damaged, resulting in a decrease in the deformation modulus of the sample. With improved deformation performance, the stress decreases. On the other hand, dynamic disturbance leads to further compression of the rock and the strain will further increase. However, once the disturbance is removed, the rock maintains partial residual deformation without full restoration of elastically, so the rock stress decreases.

DISCUSSION
In the research of predicting the occurrence of temporal events, some scholars have made many outstanding contributions. Li et al. (2021a) used long-short-termmemory recurrent neural networks to predict the wind directions for a wind turbine generated largest energy. In Reference (Li et al., 2021b), the deep belief network method and the exponentially-weighted moving average (EWMA) control chart were utilized to predict the generator bearing failure on the condition of the dynamic and random stress. He and Kusiak (2018) applied the extreme learning machines to predict the performance of wind turbines at future. Numerical simulation is also a method to predict the occurrence of time events. We first compare the theoretical model with the numerical model on the rheological shaft with interbedding structure. The numerical model can predict the process of the rheological rockmass with interbedding structure for visco-elastic analysis. We assume that there is an infinite geologic body, under both in situ horizontal stress λp 0 and vertical stress p 0 . In this example, a circular tunnel with a radius r 1 is excavated with soft ore and a radius of r 1 + d on the excavation surface. In the model, the in situ stresses λp 0 and p 0 are applied at the external boundary in the X-and Y-directions, respectively, and the radial compressive stress σ r e can be obtained at the inner  boundary, which varies with the polar angle θ and can be obtained by application of the coordinate transformation formulas such as Eq. 25 (Yu et al., 1983). Immediately after the application of cavern wall rock displacement load, elastic effects occurred in the cavern wall rock. Under the action of crustal stress, the radial stress in the planar stress cavern wall rock can be expressed as σ e r p 2 ( VERY BLOCKY-interlocked, partially disturbed rock mass with multifaceted angular blocks formed by four or more discontinuity sets 3<J v ≤10 ----BLOCKY/DISTURBED-folded and/or faulted with angular blocks formed by many intersecting discontinuity sets 10<J v ≤30 ----DISINTEGRATED-poorly interlocked, heavily broken rock mass with a mixture or angular and rounded rock pieces  GH + G c (3κ + 1) (r 1 + d) 2 − r 2 1 3 (28)  H (r 1 + d) 6 (κ c + 3) + 3(r 1 + d) 4 r 2 1 (3κ c + 1) + 3(r 1 + d) 2 r 4 1 (κ c + 3) + r 6 1 (3κ c + 1) Where r is the distance from the excavation center to some point of rock mass, d is the thickness of the soft ore, E is the elastic modulus of the hard rock mass, γ is Poisson's ratio of the hard rock mass, E c is the elastic modulus of the soft ore, and γ c is Poisson's ratio of the soft ore. This rheological model represents the strain and displacement caused by long-term in situ stresses, in which the stresses of each point are time-dependent, as shown in Figure 10A.
In general, the rock creep process can be influenced by a series of factors including stress level and loading time. The creep strain can be described as (Kraus, 1980): ε ce Aσ m0 e t n0 (33) used the maximum tension strain and the Mohr-Coulomb criteria as the damage assessment criteria to establish a rheological numerical model of the double-rock. The numerical result suggests that the interbedding structure rock exhibits significant creep occurs under multi-stage creep test with dynamic disturbance. The damage of the double-rock creep was an oblique crack at 30°to the horizontal direction under the stress = 15.18 MPa and impact high = 30 cm. 2) The rheological numerical simulation of mining disturbance of the interbedding structure surrounding the main shaft at -380 m reveal that the strain and stress increase simultaneously and exhibit consistency after the ore rock close to the shaft is subject to stress wave disturbance. After the interbedding structure near the shaft is subjected to stress wave disturbance, the strain increases and the stress decreases, causing unloading due to the residual deformation produced by the rock before and after dynamic disturbance.