Influence of Geostress Orientation on Fracture Response of Deep Underground Cavity Subjected to Dynamic Loading

Deep underground cavity is often damaged under the combined actions of high excavating-induced local stresses and dynamic loading. The fracturing zone and failure type are much related to the initial geostress state. To investigate the influence of geostress orientation on fracture behaviours of underground cavity due to dynamic loading, implicit to explicit sequential solution method was performed in the numerical code to realize the calculation of geostress initialization and dynamic loading on deep underground cavity.The results indicate that when the geostress orientation is heterotropic to the roadway’s floor face (e.g., 30 or 60), high stress and strain energy concentration are presented in the corner and the spandrel of the roadway, where V-shaped rock failure occurs with the release of massive energy in a very short time. When the geostress orientation is orthogonal to the roadway (e.g., 0 or 90), the tangential stress and strain energy distribute symmetrically around the cavity. In this regard, the stored strain energy is released slowly under the dynamic loading, resulting in mainly parallel fracture along the roadway’s profile.Therefore, to minimize the damage extent of the surrounding rock, it is of great concern to design the best excavation location and direction of new-opened roadway based on the measuring data of in situ geostresses.


Introduction
With the increase of excavated depth for underground rock engineering, numerous unconventional rock failures such as rock spalling [1][2][3], zonal disintegration phenomenon [4], and rockburst [5][6][7][8] have been often observed.One essential precursor for inducing these failures is the high geostress, which accumulates gradually in the surrounding rock during the excavation process.When underground roadway is situated in highly stressed rock mass, its mechanical behavior and failure type are influenced by various factors, including inner factors and external factors.The former include the rock mass's property, the roadway's type and size, and the spatial orientation of the roadway with geostresses, while the external factors that contribute to rock failure include the geostress state, the excavation method, and external dynamic disturbance [5,9].
Many studies have been devoted to investigating the rock failure mechanism of underground roadway, especially for the case when the roadway is situated in high geostress environment.Li et al. [10] investigated the influences of various in situ stresses, unloading rates, and paths on dynamic effects of circular opening, indicating that the increase of in situ stresses can greatly exacerbate the opening's dynamic effects.Lu et al. [11] also investigated the dynamic response of roadway due to transient release of initial geostress.While the excavation unloading process induces damage to underground cavity, dynamic loading is also a nonignorable cause to induce rock degradation or even failure [12,13].Zhu et al. [12] employed RFPA code to study the dynamic failure process of an underground opening under static geostresses and dynamic loading.Tao et al. [13] examined the multiple fracture zones around highly stressed cavities, and the results illustrated that the dynamic load and static stress gradient are two critical factors to induce multiple fracture.The previous researches indicated that the fracturing behavior and damage zone of underground roadways are generally different under various geostresses.
Meanwhile, underground rocks are always subject to complex stress fields, including gravity stress and tectonic stress [10].When excavating in complex stress fields, the new-opened cavity is often situated in an occasion when no former engineering experience can be referred from, due to the various vertical and horizontal stresses and their various spatial orientations.While there were some reports published focusing on the plastic zone and damage extent of underground opening under various lateral pressure coefficients [14], few were seen with regard to the influence of various geostress directions on the fracturing behavior of underground roadway.Pioneer work including this consideration can be found in [15] which investigated the effect of principal stress orientation on tunnel stability but did not involve the external factor of dynamic disturbance.
In deep underground excavation, due to the combined actions of geostresses and dynamic loadings, the surrounding rock around work face is usually induced to fail with the ejection of rock pieces.The rock failure process is widely accepted as a typical dynamic instability and energy release process.This paper explores the dynamic failure response of deep underground roadway subjected to dynamic disturbance under various orientations of geostress around the roadway.The lateral stress and strain energy distribution at the roadway's periphery are examined.Additionally, the strain energy releasing mechanism of the surrounding rock under dynamic loading is investigated, including its releasing magnitude and time.Finally, the influence of geostress orientation on the damage zone and fracturing pattern of the cavity was discussed.

Numerical Simulator Descriptions
The finite element method (FEM) program Ansys/Ls-dyna is widely used in nonlinear dynamical calculations to simulate sheet metal forming, bird strike and material failures, and so on.It is quite suitable for simulating rock failure due to large deformation and nonlinear dynamic loading.In this study, Ansys/Ls-dyna was explored to examine the dynamic behaviors of the underground opening subjected to dynamic loading.

Modelling Layout.
A deep-buried roadway with straightwall-top-arch cross section was considered in this case.The specific geometries and loading conditions of the numerical model are shown in Figure 1.Provided that the mechanical characteristics in every cross section along the roadway axis are the same, the numerical model was assumed as a plane strain problem.Therefore, a single layer meshed model was established by employing three-dimensional eight-node solid element type.By using three-dimensional solid element to solve the plane strain problem, it is beneficial for improving the accuracy of the simulated results and reducing the computation time, as well as making the calculation process easier to converge.The calculating model was first prestressed by principal stresses, that is, the maximum principal stress  1 and the minimum principal stress  3 .In this case, the principal stress state with  1 of 25 MPa and  3 of 10 MPa was applied on the boundary faces.Then, a triangle stress wave (Figure 2) was applied at the element component which is 11.5 m away from  the roof.The total action time of the stress wave is 4.0 × 10 −4 s, of which the rising time is 1.0 × 10 −4 s and the peak stress is 40 MPa.All the faces are defined as nonreflecting boundaries to exclude the reflected stress waves that may be generated at the model boundaries.

Implicit to Explicit Sequential Solution Method.
The dynamic response of deep-buried opening involves two loadings of static geostress and dynamic disturbance, which is a coupled static and dynamic loading problem.The problem should be divided into two steps for numerical calculation: the static stress state of overstressed opening should be solved at first, and then the dynamic loading process will be computed on the basis of the static results.For this problem, the code provides a good solution method, that is, implicit to explicit sequential solution method.First, the implicit module is used to calculate the initial static geostress state of the opening.The strains, the displacements, and the stresses obtained from implicit calculation are imported into the explicit module, which is accomplished by creating a database file that updates the geometry and the stress history of the explicit element so that it matches the implicit static solution [16,17].The flow chart for implicit to explicit sequential solution process is presented in Figure 3.

Material Model for Rock
3.1.The Brittle Damage Model.Lots of material models have been developed to simulate the damage or fracture process for rock or rocklike materials, such as the Johnson-Holmquist model, the continuous surface cap model (CSCM), and the brittle damage model that are available in Ls-dyna.The models are designed for special purpose to take into account the erosion, strain rate effect and cracking, and so forth.The Johnson-Holmquist model is advantageous for considering compression damage but does not consider tensile damage as extensive [17].The CSCM includes an isotropic constitutive equation, yield, and hardening surfaces and is well used to simulate dynamic characteristics of rock failure [17,18], but it is very difficult to determine the input parameters.
In this paper, the brittle damage model was employed to simulate the dynamic fracturing characters of the roadway under high geostress and dynamic loading.The model is designed primarily for concrete though it can be applied to a wide variety of brittle materials [16], especially be suitable for simulating rock fracture under a tensile force [19].It also contains a minimal set of material constants which can be determined from the standard tests [20].More detailed demonstrations on this material model can be referred to from [16,21].

Determination of the Parameters for the Brittle Damage
Model.To validate the capability of the brittle damage model for simulating the fracturing behavior of rock, uniaxial compression tests (UCT) were conducted to compare the results from numerical simulation.The specimens for experiments were extracted from local surrounding rock buried in 500 m depth of Linglong gold mine where a violent bulking rockburst occurred in January, 2013 [22].The UCT were carried out on a servo-controlled testing system SANS-CHT4605.
Then, a meshed cylinder entity with the same size as the rock sample was built in the code.Uniaxial compression simulations were conducted when the boundary conditions and the loading scheme were the same with experimental conditions.The stress strain relation curves from simulation and experiment were compared once the numerical calculation was finished.Differences between the two curves were step-by-step narrowed by modifying the input parameters during the repeated simulation process.The experimental stress strain curve was compared with the simulated results, as shown in Figure 4.It is noticed that the simulated stress strain curve is basically close to the experimental curve, especially for the postpeak stage.Note that the peak strain of the simulations is less than experimental peak strain, which is mainly because the initial microcracks within rock are not taken into consideration in the code.Although there are differences in local strain value, the peak stress and overall strain are highly close, indicating that the material model is suitable for further simulation studies.The corresponding parameters as input were listed in Table 1.

Numerical Simulations and Results
The validated rock material was used in the numerical code to investigate the energy and failure characteristics of the cavity  in this section.First, the distributions of lateral stress and strain energy density were calculated when the underground roadway was excavated.Then, the strain energy evolution of surrounding rock due to dynamic loading was analyzed, and the fracturing behaviors of the opening were simulated under different geostress states.

Stress Redistribution in Geostressed Environment.
When the underground roadway is created, the geostresses near the cavity workface will change: the tangential stress of the surrounding rock increases while the radial stress decreases gradually to zero at the opening's surface along the radial direction from outside rock mass.The tangential stress distribution at the periphery of the opening is different due to different initial stress states.Now that the closedform solution of tangential stress for a tunnel with straightwall-top-arch cross section is not straightforward derived, thus, the tangential stress around the roadway was calculated and plotted in the numerical program.Figure 5 presents the tangential stress distribution of the roadway under four different initial stress states, especially when the orientation angle of the maximum principal stress with the roadway's floor face () is 0 ∘ , 30 ∘ , 60 ∘ , and 90 ∘ , respectively.
It can be seen from Figure 5 that the compressive and tensile stress zone are alternatively distributed around the opening.Particularly, when the orientation angle  is 0 ∘ or 90 ∘ , the tangential stress distribution is symmetrical.High  compressive stress zones are shown in the places that are parallel to maximum principal stress, while lower compressive stress zones or even tensile stress zones are observed in the places that are vertical to maximum principal stress.In addition, obvious compressive stress concentration appears around the corner of the roadway for all the four cases.

SED Distribution in Static Geostressed
State.Based on the energy theory, the failure process of rock is actually a process of strain energy accumulation, dissipation, and release inside the rock.The excavating-induced stresses and the strain energy inside the surrounding rock keep changing continuously as the workface of the cavity moves forward.Some places of the rock mass accumulate energy, while others release energy.The rock will be induced to break once a certain failure criterion is satisfied.To express the energy characteristics of deep surrounding rock, strain energy density (SED) is used in this paper which is written as where  1 ,  2 , and  3 are the maximum, mediate, and minimum principal stress, respectively,  is the elastic modulus, and  is Poisson's ratio.The SED distributions of the roadway under the four cases of geostress orientation angles are shown in Figure 6.
From Figure 6, it can be seen that the SED distribution is generally similar to the distribution of compressive stress zone; that is, highly stressed zone promises high SED concentration region.The maximum values of SED for the four  cases are obviously various.The maximum magnitudes reach 87.5 kJ/m 3 and 83.6 kJ/m 3 when  is 0 ∘ and 90 ∘ , respectively.While  equals 30 ∘ or 60 ∘ , the SED at the corner is greatest, reaching 178.2 kJ/m 3 and 177.4 kJ/m 3 , respectively.This is because when the direction of maximum principal stress is parallel or vertical to the floor face (e.g.,  = 0 ∘ or 90 ∘ ), the induced stress state of the roadway is symmetrical, which is beneficial for moderate deformation of the surrounding rock.In this way, the strain energy distribution around the roadway is alleviated and balanced.
However, when the direction of maximum principal stress is heterotropic to the floor face (e.g.,  = 30 ∘ or 60 ∘ ), the stress state is extremely imbalanced, resulting in high stress concentration near the corner and spandrel of the roadway.In this regard, the surrounding rock masses accumulate a large amount of strain energy, which is disadvantageous for rock stability after excavation.

SED Dissipation due to Dynamic Loading.
To further examine the strain energy releasing response of the opening subjected to dynamic disturbance, four special elements which are located at roof, left sidewall, right sidewall, and floor, as marked in Figure 2, were chosen to trace the SED change law. Figure 7 presents the SED evolutional curves for the four locations under different geostress orientations.
When  is 0 ∘ , the strain energy accumulates at the roof of the roadway in a very high level.The stored strain energy declines sharply to zero in just 1.16 × 10 −3 s, indicating that the roof surrounding rock fails fiercely due to the dynamic loading, whereas little influences on the change of strain energy are observed at the opening's sidewalls and floor under the action of dynamic stress wave.
When  equals 30 ∘ , it is noticed that the SED at the roof is much greater than other locations of the opening and fully released within 1.48 × 10 −3 s, indicating a violent dynamic failure occurs in the rock mass.But there is little change on SED at the sidewalls and floor for this case.
For the case of  = 60 ∘ , the strain energy distribution of the roadway is obviously different.Under the dynamic disturbance, the SED of the right sidewall experiences a long time of punctuation and declines to a lower level eventually.The SED of other places, such as the left sidewall and the roof, decreases slowly till zero with a long releasing period of 4.17 × 10 −3 s.When  reaches 90 ∘ , the excavating-induced stresses are symmetrically distributed around the roadway.The change law of the strain energy at the left sidewall and the right sidewall is generally the same, reducing from 14.13 kJ/m 3 to 0 after 4.81 × 10 −3 s.On the contrary, there is little strain energy accumulated in the roof and floor, where spalling failure is induced subjected to the dynamic disturbance.
In terms of the energy releasing time, it can be observed that the strain energy releasing time for the roof when  = 0 ∘ and 30 ∘ is much less than that when  = 90 ∘ .When the geostress orientation angle is 0 ∘ or 30 ∘ , the roadway's roof performs obviously high stress concentration, gathering considerable strain energy there.Under the direct impact of dynamic loading from the top rock mass, the strain energy within the roof is released sharply, indicating a fierce rockburst happens in this occasion.Differently, it takes a longer period of time for the SED to be released at sidewalls when  = 90 ∘ .In this regard, the high strain energy concentrates mainly at the sidewalls, where less impact will act upon the surrounding rock when subjected to the dynamic loading, so the energy is released relatively slower.

Fracturing Zone and Failure Type.
To understand how the geostress orientation influences the fracture zone around the opening under dynamic loading, further simulations were extended in the numerical code.The plastic region and fracture pattern are presented in Figure 8, where the fringe levels with various colors indicate the magnitudes of plastic damage.The fracturing zone is equivalently represented by deleting the failure elements automatically.
From Figure 8, the damage zone and extent of the surrounding rock mass have distinct differences under different orientation angles of the principal stress with the floor face.In general, the fracturing zone mainly performs around the surrounding rock in which it is parallel to the maximum principal stress.This is because the failure place was previously subject to high compressive stress, and it is quite prone to collapse under the superposition action of dynamic disturbance, with releasing a large amount of strain energy as well.
When the angle is 30 ∘ or 60 ∘ , the damage zone and extent are much severer, presenting V-shaped breakout in high accumulated energy places.In this occasion,  the excavating-induced stresses distributed around the surroundings are nonuniform, resulting in high stress concentration in the corner and spandrel of the roadway.The highly stressed surrounding rock will be largely fractured when experiencing dynamic loading.When the principal stress direction is vertical or parallel to the floor face, that is, when  is 0 ∘ or 90 ∘ , the failure of the surrounding rock is characterized as parallel fracture along the roadway's profile.

Conclusions
An implicit to explicit sequential solution method was employed in Ansys/Ls-dyna to realize the calculation of dynamic loading on highly stressed underground roadway.A validated material model was used to explore the fracture response of deep-buried roadway due to static geostress and dynamic loading, especially to investigate how the geostress orientation influences the strain energy release evolution and failure pattern of the underground roadway subjected to dynamic loading.The numerical results indicate that when the geostress orientation angle  is 30 ∘ or 60 ∘ , tangential stress and strain energy concentrate mainly in the corner and the spandrel, where V-shaped breakout with releasing massive energy is observed when experiencing dynamic loading.When the geostress orientation is 0 ∘ or 90 ∘ , the stress and strain energy distribute symmetrically around the roadway, leading to parallel fracture along the roadway's profile with slower release of strain energy.This study can contribute to optimal selection of excavation location and direction of new-opened cavity according to the in situ geostress state, as well as to mitigating the fracture and damage extent of the surrounding rock when proper support strategy is carried out.

Figure 1 :
Figure 1: Numerical model with specific geometries and loading conditions.

Figure 2 :
Figure 2: Triangle stress wave used as dynamic disturbance.

Figure 3 :
Figure 3: Flow chart for implicit to explicit sequential solution process.

Figure 4 :
Figure 4: Comparison of the experimental stress strain relation with numerical result.

Figure 6 :
Figure 6: SED distributions under four different orientation angles of the maximum principal stress with the floor face (unit: kJ/m 3 ).

3 )
sid d d d d d d d d d d d d d d d d dewa e e e e e e e e e e e ll Roof Floor SED (KJ•m − = 90 ∘

Figure 7 :
Figure 7: SED evolutional curves for the four locations under various geostress orientations.

Figure 8 :
Figure 8: The plastic region and fracture pattern of the roadway under various geostress orientations.

Table 1 :
Input parameters for the rock model.