Excavation unloading‐induced fracturing of hard rock containing different shapes of central holes affected by unloading rates and in situ stresses

Slabbing failure and strain rock burst are the main failure patterns during the excavation and construction phases of deep tunnels in hard and brittle rock, which cause unexpected equipment damage and casualties. This study presents a numerical simulation with an unloading central hole within a hard rock specimen in the laboratory scale. A combined finite element/discrete element method (FEM/DEM) reproduces the crack initiation, propagation, and coalescence around the central hole during the entire failure process. The influence of sectional shapes, unloading rates of holes, and in situ stresses on the failure characteristics and mechanical response of typical hard and brittle rocks were investigated by analysis of the failure pattern, displacement distribution (average velocity) of monitoring elements, and released strain energy value. The numerical results indicate that the sectional shapes, unloading rates, and in situ stresses have a significant impact on the severity of destruction and failure range in hard rock under the excavation unloading conditions. Slabbing failure (stable failure) is always the dominant failure pattern around a circular hole, which shows higher bearing capacities and self‐stabilization. With the increase in unloading rates, more visible cracks are generated around the holes, and the displacement and average velocity of discrete blocks are further raised. This is particularly evident in rock specimens with holes having vertical walls, leading to intense unstable failure (strain rock burst) accompanied by a large amount of strain energy released. In situ stresses affect considerably the stability of the surrounding rock during the excavation unloading process. With constant in situ stress, the order of destruction severity around the central hole according to its sectional shape is cube > trapezoid > U‐shape > ellipse > circle. The destruction intensity is further aggravated with the increase in lateral pressure coefficient, and the failure regions are always observed in the roof and floor of the central hole. This study confirms that strain rock burst tends to be induced in hard and brittle rock tunnels with polygonal roadway section under high unloading rates and lateral pressure coefficients.


| INTRODUCTION
A large number of mines worldwide have gradually entered into the deep or ultra deep mining stage. 9,31,39 The dynamic hazards and safety problems related to deep rock mechanics, such as rock burst, zonal disintegration, and mine earthquake, have aroused widespread concern (Guzev et al). 125,11,17,27,28,4 3,48,49,55,58 A large amount of initial energy has been stored in the rockmass at great depth as it is usually subjected to higher in situ stresses. During the excavation unloading process, the stress states of the surrounding rock will be redistributed and concentrated, and the energy stored in the rockmass will be released and transferred. If the energy released by the unloading process under high-stress conditions is sufficient, dynamic hazards will arise and intense failure will take place thereafter. 47,51 In these incidents, brittle and hard rocks may exhibit distinct failure properties and mechanical responses compared with those under shallow buried conditions, 10,21,24,64 and the progressive generation of fractures parallel to the surface (spalling or slabbing) and sudden violent failure with a certain amount of blocks being ejected toward the free surface (strain or structural rock burst) are the dominant failure patterns. 9,13 Figure 1 presents the typical slabbing failure (stable failure) and strain rock burst (unstable failure) phenomena observed in several deep hard rock mines. 31,35,63 Typically, rock burst can be classified into three types, namely fault slip, pillar, and strain rock burst. 29,34 Of these categories, strain rock burst is the most common dynamic hazard encountered in underground engineering, especially near excavation boundaries of hard and intact rockmass. The severity of the destruction around the deep tunnel and chambers is strongly related to the strain energy because the failure is accompanied by energy dissipation and release. 52 Actually, slabbing failure is the precursor of strain rock burst in most cases. 11 Tan 44 found that the progressive failure process of strain rock burst can be summarized in three stages. At the very beginning, the surrounding rock is split into several slabs parallel to the excavation boundaries. Then, these thin slabs are cut into several blocks or rock plates with different lengths and sizes. Finally, these blocks and plates are ejected toward the free surface owing to the coupled buckling and pushing caused by the concentrated maximum tangential stress and specific radial stress inside the surrounding rock. In another case, only fractures parallel to the surface (slabbing failure) are observed, with gentle and moderate displacement response and energy evolution. Some valuable efforts have been devoted to understand the relationships between the slabbing failure and strain rock burst. Previous findings have confirmed that the failure intensity and occurrence of strain rock burst near the excavation boundaries are not only related to the stored strain energy in the rocks, but are also affected by the environments of strain energy accumulation and strain energy release upon failure. 6,18,25,62 Specifically, the initial stored strain energy value is mainly governed by in situ stresses, and the strain energy concentrated and released can be affected by the tunnel sectional shape and unloading rate, respectively. Here, the unloading rate refers to the advance rate of the excavation in mining and tunneling engineering, and the in situ stresses can be represented by lateral pressure coefficients. Previous researchers have examined the influence of these factors on the mechanical response of hard deep rocks from multiple perspectives. 7,16,19,26,37,63 Xu et al 53 found that an increase in the unloading strain rates and initial stress will lead to instantaneous violent failure (unstable failure) rather than progressive failure (stable failure) under true triaxial unloading tests, and the failure modes will change from a mixture of tensile to tensile failure coupled with more transgranular cracks. Zhao et al 62  It should be noted that although the laboratory tests, such as conventional triaxial and true triaxial tests, can well model the tunnel fracturing induced by the excavation unloading effect, the entire spatiotemporal evolution characteristics of the surrounding rocks were not successfully reflected. This is especially true when different sectional shapes are taken into account. Some scholars adopted rock specimens with one or multiple central holes to reproduce the failure behavior of deep tunnels in the laboratory scale and focused on the quantitative characterization of the stress concentration and crack propagation behavior around the holes. 4,50,56,60 However, most rocks in previous tests were subjected to a stress state of continuous loading, and this was manifested in that the specimens had compression-induced fractures instead of unloading-induced failure. These difficulties in reproducing the excavation unloading process in the laboratory scale are determined by the existing experimental techniques. Therefore, a numerical simulation might be an efficient and reasonable method for addressing this issue. Li et al 22 employed the unloading relaxation method to investigate the influence of the stress path on the excavation unloading response by using the PFC 2D software based on the discrete element method (DEM). Sun et al 42 adopted the rock failure process analysis (RFPA) and discontinuous deformation analysis (DDA) software to investigate the failure patterns of a circular tunnel under unloading conditions. Manouchehrian and Cai 29,30 simulated the dynamic rock burst induced by fault zones around deep circular tunnels by using the 2D Abaqus explicit code (a finite element analysis). However, owing to the selected numerical simulation techniques, the crack initiation, propagation, and coalescence paths could not be reproduced. 9 Therefore, a suitable numerical tool should be employed to describe the entire failure process from continuum to discontinuum states. Moreover, previous studies have mainly focused on circular tunnels, failing to provide a comprehensive comparison of the failure behavior between different sectional shapes. Moreover, the integrated response of unloading rates, in situ stresses, and sectional shapes was not clearly illustrated in the previous studies.
In the present study, we adopted the combined finite element method/discrete element method (FEM/DEM), a numerical code that has the advantages of both techniques and integrates the FEM and DEM into a singular tool to simulate the failure behavior of a typical hard rock specimen with central holes in the laboratory scale. The purpose is to model the real failure of the surrounding rock and the spatiotemporal evolution characteristics caused by tunnel excavation with various unloading rates, sectional shapes of central holes, and in situ stresses. Rock heterogeneity was also introduced in the model. The entire failure process, displacement distribution (average velocity) of the selected monitoring points, and released strain energy of the model were discussed in detail. This study is intended to deepen our understanding of tunnel failure at great depth and provide a theoretical guidance during the construction phase from the perspective of excavation method and tunnel section design, coupled with specific in situ stresses.

| Principles of the combined finite element method/discrete element method
The combined FEM/DEM was initially proposed by Munjiza 33 and has been widely used in the field of rock and geotechnical engineering in the past decade. This method can be used for simulation of laboratory tests, such as the Brazilian, uniaxial, and triaxial compression tests, hydraulic fracturing, 2,14,24,54 and in situ engineering problems such as those inherent to rock slopes, blasting, and block-caving induced surface subsidence. 1,15,32,45 The basic principle of FEM/DEM is to divide the solution domain into discrete elements along the existing discontinuous surface and further divide the finite element mesh within each discrete element. By introducing the damage mechanics and fracture mechanics theory into the FEM, the cracking of discrete elements is specified to take place either along the boundary or through the finite element. The motion and interaction of the discrete blocks are operated in the same way as in the DEM, in which the motion of a singular block is determined by its unbalanced force or moment in accordance with Newton's second law. The discrete blocks can move and rotate freely, and there is no need to satisfy the deformation coordination condition. For the fracturing simulations, for example, rock blasting, rock and concrete spalling, and structural collapse, the explicit algorithm was usually adopted owing to its ability in reducing memory consumption and optimizing the operating procedures when dealing with engineering problems. 54 In the current study, a commercial numerical tool named Elfen was adopted for modeling hard rock failure around the central hole. An obvious distinction from other numerical codes is its ability to allow the new cracks to cut across the original mesh (intraelement fracturing) rather than merely along or around the element boundary (interelement fracturing). Elfen has been verified for properly addressing rock engineering problems, and therefore, it will not be introduced in this study.

| Establishment of numerical model and boundary conditions
In the current study, the 2D geometry type (strain plane) analysis was adopted, in which the specimen thickness was assumed to be 1 mm. The height and width of the external boundary in the model were both 500 mm. The surrounding rock was represented by a single face, with a central hole indicating the tunnel to be excavated. The area of central hole was set to be 1963 mm 2 . The geometry size of each hole shape was derived by back-calculation of the constant area. In general, there is a common sense agreement that the external boundary of the model should be at least ten times the tunnel diameter or length to exclude the effect of the boundary on the stress redistribution around the tunnel and to avoid the influence on the numerical results of potential tensile stress waves induced by the excavation unloading, which are reflected back from the external boundaries toward the central hole. 9,29 Five commonly used sectional shapes were considered, which were the circle, ellipse, U-shape, trapezoid, and cube. By calculating the diameter or length of the central holes, it was found the ratio of external boundary height and width to the tunnel size (approximately 10 times for all the cases) that can guarantee the accuracy and rationality of the numerical results. Figure 2 illustrates the model geometry and boundary condition. Because visible cracks were always generated and gathered around the hole boundaries, a dense mesh was assigned in this area to obtain a precise numerical result. The mesh dimension of those elements situated far away from the hole boundaries should be sparser to improve the computation efficiency. The total number of unstructured triangular elements established by the model containing circular, elliptical, U-shaped, trapezoidal, and cubic holes was of 9700, 9640, 9678, 9750, and 9750, respectively.
Before excavation unloading, in situ stresses were applied to the external boundaries. A linear loading path over 0.01 seconds was employed in the present study to avoid unexpected fracturing due to transient loading. The in situ stresses were kept constant once reaching the predefined value throughout the simulation. Various in situ stresses were modeled with three lateral pressure coefficients k at 1.4, 1.7, and 2. Several monitoring elements around the hole boundaries were selected to obtain the distribution of displacement and average velocities. Elfen has a broad range of material modeling options that may be combined to create material models for different classes of problems. The Mohr-Coulomb constitutive model with a rotating crack was used in the present study because it can provide a realistic representation of both the tensile and shear failure. 14 This constitutive model requires three basic elastic parameters: Young's modulus (E, GPa), Poisson's ratio (μ), density (ρ, Ns 2 / mm 4 ), and five plastic parameters: cohesion (C, MPa), frictional angle (φ, °), tensile strength (σ t , MPa), dilation angle (ψ, °), and fracture energy (G f , N/mm). The Mohr-Coulomb yield criterion is a generalization of the Coulomb friction failure law and is defined by where τ is the magnitude of the shear stress, σ n is the normal stress, c is the cohesion, and φ is the friction angle. As the conventional Mohr-Coulomb criterion may overestimate the tensile strength of rocks, Paul 36 introduced the concept of a tension cutoff and modified Mohr-Coulomb criterion. To model the hardening or softening behavior, the material strength parameters are defined as functions of the effective plastic strain: where ɛ −p is the effective plastic strain, p 1 , p 2 , and p 3 are the extensional inelastic strains in the direction of maximum, intermediate, and minimum principle stress, respectively. A comprehensive introduction of the Mohr-Coulomb with rotating crack constitutive model is presented in Feng et al 11 and it will not be described in detail in this paper.
A typical hard and brittle rock based on marble extracted from the Hongling Zinc-Lead mine in Inner Mongolia, China, was selected for the numerical modeling, and the basic physical and mechanical parameters are listed in Table 1. Fracture energy was used for characterizing the difficulties in generating visible cracks in the FEM/DEM. The fracture energy G f can be obtained based on the following equation: where K IC is the fracture toughness in mode I, and E is the initial Young's modulus. The fracture toughness K IC can be estimated using an empirical relationship with the rock tensile strength 32 (Zhang) 59 : where σ t is the rock tensile strength determined from Brazilian tests. According to the above equations, the fracture energy was finally calculated to be 0.005 N/mm. The residual strength values (residual cohesion and frictional angle) were determined based on the classical cohesion weakening and frictional strengthening model (CWFS), which was proposed by Hajiabdolmajid et al 13 when modeling the hard failure of rock.

| Excavation unloading of central hole and material heterogeneity
The excavation unloading process was executed once the in situ stresses reached the predefined value. In the present study, the unloading relaxation method, which closely describes the actual excavation process, was adopted. A linear relaxation curve over 0.01 seconds was assigned to the central hole. Once the elements on the tunnel internal surface were deactivated after excavation, they were replaced by a series of forces equivalent to the internal forces of the removed elements. The forces were then released in accordance with a specific linear unloading path to simulate the excavation process. 22 This avoids unexpected fracturing induced by instantaneous unloading to 0, which would create a large unbalanced force. Various unloading rates representative of the unloading time at T v = 0.001, 0.01, and 0.02 seconds were considered. Taking T v = 0.02 seconds as an example, the initial unloading point starts from t = 0.01 seconds, and the terminal point of unloading should correspond to t = 0.03 seconds. The loading and unloading paths for all the cases are shown in Figure 3.
In the present study, material heterogeneity was also considered in the commercial software by specifying certain intervals for five independent material parameters: Young's modulus, Poisson's ratio, density, tensile strength, and fracture energy. A random number generator was adopted to create the values for each property of the element within the specified limits at either side of the mean value. 9 The upper limit and lower limit of each parameter are also listed in Table 1. A uniform distribution type was employed, indicating that these parameters were stochastically distributed within the marble specimen.
Normal and tangential penalties equivalent to the normal and shear stiffness were used to evaluate the contact forces of the newly generated cracks, which prevent them from overlapping. 2,15 The normal penalty is usually in the range of 0.5 E to 2.0 E, where E is Young's modulus. In the current study, the normal and tangential penalties were specified to be E and 0.1 E. 14 The smallest element size should be smaller than the original mesh element size to allow the meshes to cut across the newly generated cracks (intraelement fracturing), instead of merely along the mesh boundaries (interelement fracturing). By using the intraelement fracturing scheme, the single small fractures can be reasonably inserted, which leads to a more realist failure behavior and crack propagation path. The discrete contact parameters are listed in Table 2.

| Numerical simulation validation
In this study, the failure behavior and mechanical response of the central hole within the hard rock specimen subjected to various unloading rates and in situ stresses were investigated. To verify the rationality of the selected model external dimension and internal hole size, we plot the maximum principal stress contour of the marble specimen containing a cubic hole with various unloading rates at t = 0.04 seconds when the lateral pressure coefficient k is 2 (the vertical stress is 25 MPa, and the horizontal stress is 50 MPa), as shown in Figure 4. The stress redistribution and concentration induced by the process of excavation unloading can be clearly observed around the central hole. When k = 2, the maximum tangential stress is gradually concentrated at the roof and floor of the hole, while it is reduced at both sidewalls. The higher stress concentration and accumulated strain energy finally lead to visible cracks and fracturing in these areas. Moreover, it can also be found that the stress redistribution area has not reached to the external boundaries of the model, and the stress states near the model boundaries is always kept at the predefined value (50 MPa). This shows that the stress redistribution and potentially reflected stress waves induced by the excavation unloading process would not affect the tunnel boundaries, therefore validating the rationality of the designed geometry model. Figure 5 also plots the tangential stress-time evolution curves of selected monitoring elements around the cubic hole with unloading rate of 0.02 seconds at k = 2. It can be clearly observed that during the 0-0.01 seconds period, the tangential stress is raised continuously in a linear manner. However, the tangential stress value at each monitoring point is somewhat diverse when t = 0.01 seconds (red circle). This indicates that the presence of material heterogeneity has already affected the stress distribution around the tunnel during the loading stage. As the excavation unloading continues, the tangential stress in the roof begins to increase, while there is a disparity in the stresses among these monitoring points. Figure 5A shows that the tangential stress value at t = 0.04 seconds is not located at the roof center but around the two corners, which is different to the case of a circular hole. 9 The tangential stress values around the corners have reached approximately 120 MPa, which is far beyond the material strength of marble (approximately 70 MPa). The maximum tangential stress presents a monotonically decreasing trend as the monitoring point approaches the roof center. This might be attributed to stress mutation around the corners, which weakens the degree of stress concentration at the roof center. After 0.03 seconds, the tangential stresses decrease rapidly, and the final values are always below 30 MPa. The reduction in tangential stress in these points demonstrates that the rockmass around the roof has entered into plastic behavior. Under such condition, the roof would fracture with visible cracks. The tangential stress-time evolution curves at the right side of the cubic hole are shown to be different from those at the roof points (see Figure 5B). The tangential stress presents a monotonically decreasing trend after the end of the initial loading point (t = 0.01 seconds), except for points 7 and 8. The tangential stress continues to decrease after the terminal unloading point (t = 0.03 seconds), and it finally decreases to approximately 0. Here, we call it the "stress reduction zone." As the tangential stresses were always lower than the rock strength, the rockmass around the opening is always in the stage of elasticity (considered to be stable). For the case of a circular hole, the tangential stress value at the center of the sidewall obtained by using Kirsch's solution (σ tan = P(1 + k) + 2P(1 − k) cos 2θ, where θ is the angle between the horizontal axis and the radius, k is the lateral pressure coefficient, and P is the in situ stress in vertical direction) should be 25 MPa. 3 The present study confirms and validates that the sectional shape may have a great impact on the mechanical behavior of a deep tunnel. Therefore, it is necessary to investigate the integrated response of sectional shapes, unloading rates, and in situ stresses.

| Effects of unloading rates on hard rock failure with various hole shapes
The unloading rate is the key determinant of the unloading response when the in situ stress remains constant. 22 In For the case of a circular hole, when the unloading time is 0.02 seconds, very few cracks are observed around the central hole, indicating lower unloading disturbance. These visible cracks are mostly parallel to the excavation boundaries, which is a typical slabbing failure similar to that observed in laboratory tests and in situ monitoring in deep underground engineering. 20,41 With the increase in unloading rates (T v = 0.01 seconds), the failure extent and range are intensified to some extent, with more visible cracks generated. A further increase in unloading rates (T v = 0.001 seconds) will contribute to more cracks and an expanded failure zone near the excavation boundaries, as can be observed in Figure 6 (note that the larger failure zone around the floor should be due to the material heterogeneity, therefore leading to strength anisotropy in certain regions). This indicates that the severity of destruction and failure range will be aggravated with the increase in unloading rates. However, the failure intensity and range are still maintained to a relatively lower level. The dynamic failure phenomenon, such as slab splash or rock ejection, has never been observed around According to the description of Read 40 the multi-stage progressive failure process observed in the MBE test tunnel can be divided into four stages: initiation, dilation, slabbing and spalling, and then stabilization. Initially, a localized damage in the form of oriented flaws was generated at the boundary of the hole induced by concentration of tangential stress. Subsequently, extensive dilation occurs in these areas owing to shearing and crushing. With the development of the process zone, several thin slabs were formed by shearing, splitting, and buckling. Finally, the development of notch tips stops when the notch geometry provides sufficient confinement to stabilize the failure zone. These numerical results are in good agreement with the previous assumptions. 40 Figure 7 shows the failure process around an elliptical hole with different unloading rates. The development tendency of crack propagation affected by unloading rates is shown to be similar to what is observed for the case of a circular hole, but the severity of destruction and failure range become more violent and extensive, especially when the unloading time is reduced to 0.01 and 0.001 seconds. Under such conditions, slight strain rock burst is observed near the hole boundaries, with a few thin plates and blocks splashed toward the free surface. This might be attributed to the larger cross-sectional span at the roof and floor, which further weakens the stability of the surrounding rock. For the case of a U-shaped hole, when the unloading rate is lower (T v = 0.02 seconds), the cracks are mostly generated in the floor, whereas the roof presents good stability. This is expected as the arched structure has higher bearing capacity compared with a polygonal structure. With the increase in unloading rate, the severity of destruction and failure range tends to deteriorate. The typical strain rock burst phenomenon is observed when the unloading time is reduced to 0.01 and 0.001 seconds, with a large quantity of thin slabs and blocks ejected from the floor (see Figure 8). When looking at the fractures around the U-shaped hole, intense strain rock burst (unstable failure) occurs when the unloading time is equal or less than 0.01 seconds. Figure 9 shows that the floor is significantly damaged compared with the roof. In this regard, the severity of destruction and failure range of the model with an unloading time of 0.001 seconds is more obvious than that with 0.02 seconds. As for the cubic hole, when the unloading time is 0.02 seconds, only slabbing failure is observed around the hole boundaries, accompanied by roof caving and slight floor heave. With the increase in unloading rates, thin slabs are first generated in the roof and floor (see Figure 10). Subsequently, these slabs were cut into pieces and ejected toward the free surface due to the rapid and extensive release strain energy, which finally leads to intense strain rock burst. These results demonstrate that the slabbing failure is the precursor of strain rock burst. By observation of the entire failure process around the central hole within the marble specimen, it is indicated that the characteristic of the failure around deep tunnels is the integrated response of sectional shapes and unloading rates. When the unloading rate remains constant, the most favorable sectional shapes for crack propagation and failure development around the hole are the cube and trapezoid, followed by the U-shape, ellipse, and circle. To better illustrate the failure intensity of a rockmass induced by the unloading rates and sectional shape around the central hole, the displacement-time evolution curves of selected monitoring points at the roof and floor for the case of elliptical, U-shaped, and cubic holes are plotted. Given the fact that the potentially ejected elements may collide or contact mutually after a long unloading duration, we only plot the displacement-time evolution curves before t = 0.021, 0.03, and 0.04 seconds with an unloading time of 0.001, 0.01, and 0.02 seconds, respectively. As shown in Figure 11A (elliptical hole), when the unloading rate is lower (T v = 0.02 seconds), the displacement values at the roof that correspond to t = 0.04 seconds are 0.14, 0.13, and 0.14 mm, respectively. With the increase in unloading rates (T v = 0.01 seconds), these displacement values (when t = 0.03 seconds) reach 3.5, 0.21, and 0.25 mm, which are higher than those in the case of a lower unloading rate. When the unloading time is reduced to 0.001 seconds, the corresponding displacement values are raised to 10.9, 14.1, and 14.3 mm, indicating the severe displacement response induced by the rapid excavation unloading process under such condition. A similar variation trend could also be found in the floor. However, a careful comparison shows that the displacement values in the floor are always larger than those in the roof when the unloading rate remains constant. This should be due to the material heterogeneity that leads to asymmetry in the mechanical response. When looking at the U-shaped hole (see Figure 12), the displacement values at 0.021, 0.03, and 0.04 seconds present a monotonically increasing trend. Moreover, the displacement values at the floor are obviously larger than those at the roof, which is consistent with the description of the entire failure process in Figure 8. For the case of the cubic hole, the maximum displacement value at the roof and floor is approximately 30 mm (see Figure  13) when the unloading time is 0.001 seconds, which is evidently larger than that of the elliptical hole (15 mm), indicating a more severe failure intensity and stronger dynamic response. Figure 14 depicts the average velocities of three monitoring points for each hole shape with various unloading rates. The average velocities, a quantitative index for characterizing the failure intensity during the excavation unloading process, are calculated as the specific time quantum divided by the maximum displacement value. The numerical results indicate that for the case of a circular hole, the average velocities are always kept at a lower level, irrespective of the unloading rates. Under such condition, the severity of destruction and failure range are gentler, which is an indicative of stable failure. For the other four different shapes of holes, the failure characteristics of the surrounding rock exhibit a strong unloading rate dependency. In this regard, the severity of destruction and failure range are greatly intensified, with more thin slabs and blocks splashed violently toward the free surface, especially for the case of the cubic and trapezoidal holes (the maximum average velocity reaches approximately 2.5 m/s). Therefore, it is inferred that unstable failures, such as strain rock burst, are more prone to be induced when there is polygonal structure and higher unloading rates.
The rock failure is accompanied by transfer and release of energy. In this work, we also depict the released strain energy-time evolution of the model for each hole shape with various unloading rates. It can be observed in Figure  15 that during the initial loading stage (0-0.1 seconds), the released strain energy is always 0 for all the cases as the excavation unloading activity has not been executed; thus, the energy is gradually accumulated rather than being released. During the excavation unloading process, the strain energy will be released only when it approaches the terminal unloading point, which shows a similar trend with that of the process of crack initiation, propagation, and coalescence. The released strain energy remains constant once the entire unloading process is terminated. Figure 16 also plots the maximum released strain energy distribution for the model containing different hole shapes with various unloading rates. It can be observed in Figure 15 that the maximum released strain energy increases with the unloading rates regardless of the sectional shapes, showing that the failure process is accompanied by a greater release of energy if the selected excavation method has higher unloading rates, such as in the full-face drilling and blasting method. Under such condition, the rapid unloading will create large unbalanced forces and displacements near the excavation boundaries, and therefore, cause severe fracturing. 61,62 Contrarily, when tunnels are excavated by the jackleg, tunnel boring machine (TBM), and hydraulic rock-splitting methods, which have lower unloading rates, a gradual stress relaxation will be executed, which imposes a gentler and moderate release of strain energy. 22 As shown in Figure 16, when the unloading rate remains constant, the larger released strain energy always corresponds to the sectional shape with polygonal structure. For example, the released strain energy values of the trapezoidal and cubic holes with an unloading time of 0.001 seconds are 421.4 and 503.4 J, respectively, which are far higher than those in the cases of circular and elliptical holes, with maximum released strain energy of 157 and 216.5 J. Because the stress concentration is more prone to be induced around unexpected corners and vertical sidewalls, the higher energy intensity will be gathered near the excavation boundaries. Therefore, the failure depth will be further expanded owing to the sustained stress redistribution toward the interior of the surrounding rock. The failure zone will stop when the stress state reaches equilibrium.

| Effects of in situ stresses on hard rock failure with various hole shapes
In situ stresses have been confirmed to substantially affect the tunnel failure at great depth. 38,46,57 In this section, three lateral pressure coefficients (1.4, 1.7, and 2) representing various in situ stresses were modeled to investigate the rock fracturing around different shapes of holes within the marble specimen. Note that the vertical stress in this section is 30 MPa, and the horizontal stress is set to be of 42, 51, and 60 MPa, respectively. Figures 17-19 plot the entire failure process and crack patterns around circular, elliptical, and trapezoidal holes with various lateral pressure coefficients k with a constant unloading rate (T v = 0.01 seconds). The failure process and crack propagation path for the case of U-shaped and cubic holes is not discussed in detail owing to their similar variation trend. The numerical results show that the severity of destruction and crack distribution are strongly related to in situ stresses. It can be observed from these plots that the visible cracks are always distributed in the roof and floor. When k is greater than 1, the higher stress concentration and strain energy accumulation may take place around the roof and floor, which finally leads to visible fracturing and cracks in these areas. For the case of the circular hole, when k is 1.4, fewer cracks are generated around the excavation boundaries. Local slabbing failure can be found at the top and bottom of the circular hole. When k is raised to 1.7, the failure range is expanded, with more cracks generated in the roof and floor. The condition is further exacerbated when k is raised to 2. Under such condition, the fracturing degree at the floor is more severe than that at the roof, which should be caused by the material heterogeneity (similar phenomenon could also be found in the cases of elliptical and cubic holes). Actually, the asymmetrical crack distribution around a centrally asymmetrical tunnel is not only induced by the initial material heterogeneity, but is related to the progressive crack propagation, which affects the formation of subsequent ones. Cracking irrecoverably weakens the physical properties of the rock mass and further promotes the heterogeneities in nonirregular crack distributions. This has been confirmed in previous findings. 23 The analysis of failure behavior around circular holes suggests that the dominant failure pattern should be categorized as stable failure, although the failure degree and range are intensified to some extent. This indicates that the circular tunnel can maintain stability. As for the elliptical hole, when k is increased to 2, a slight strain rock burst phenomenon occurs where several thin slabs are ejected toward the free surface. It can be observed in Figures 17 and 18 that the severity of destruction and failure range around the elliptical hole is aggravated compared with the circular hole. When looking at the failure characteristics of the rockmass around the trapezoidal hole, an intense and violent strain rock burst can be observed in the surrounding rock when the lateral pressure coefficient k is 2. Under such condition, the thin slabs derived from the tensile cracks induced by the excavation unloading effect are cut into pieces of blocks with various sizes. The severity of failure and range are greatly intensified compared with the case of k = 1.4, as can be observed in Figure 19. It is also found that the favorable condition leading to rock failure should be in those holes with polygonal structure and higher in situ stresses.
To better illustrate the effect of in situ stresses on the failure characteristics of a deep tunnel, we plot the displacement-time evolution curves of selected monitoring points around the trapezoidal hole with various lateral pressure coefficients at a constant unloading rate (T v = 0.01 seconds), as can be seen in Figure 20. When k is 1.4, the displacement values corresponding to t = 0.04 seconds are in the range of 3.5-6.6 mm. With the increase in lateral pressure coefficient (k = 1.7), the maximum displacement value at t = 0.04 seconds is 22.5 mm, indicating that the failure intensity is significantly exacerbated. When k is raised to 2, the displacement values at t = 0.04 seconds are further increased (in the range of 23-42 mm). Figure 21 depicts the average velocities of each monitoring point around the trapezoidal hole affected by various k. It can be observed in the figure that the average velocities for six monitoring points remain at a low level when k is 1.4, indicating a moderate destruction intensity. With the increase in k, the average velocities are promoted, especially for the case of k = 2. This is an indication that the rockmass near the excavation boundaries is splashed in a violent and intense manner, which is the typical characteristic of strain rock burst (unstable failure). Figure 22 shows the released strain energy-time evolution curves for different sectional shapes with various lateral pressure coefficients, and Figure 23 plots the maximum released strain energy distribution of the model for different sectional shapes with various lateral pressure coefficients. On the one hand, for a constant k, the sequence of the maximum released strain energy is always On the other hand, for a certain shape of the central hole, the released strain energy presents a monotonically increasing trend with the lateral pressure coefficients. Taking the cubic hole as an example, the maximum released strain energy values for the cases of k = 1.4, 1.7, and 2 are 216.4, 396, and 525 J, respectively. This is primarily attributed to the fact that when the lateral pressure coefficient is high, the accumulated elastic strain energy stored in the specimen is elevated. Therefore, the strain energy will be released in large quantities once the excavation unloading activities are executed, which finally causes a violent and intense dynamic hazard phenomenon. The abovementioned analysis indicates that the unstable failure is more prone to be induced when there are larger in situ stresses (lateral pressure coefficients) and sectional shapes with polygonal structures.

| CONCLUSIONS
The combined FEM/DEM was employed to investigate the failure characteristics around deep hard rock tunnels with five sectional shapes affected by various unloading rates and in situ stresses in the laboratory scale. The analysis of the failure process (crack initiation, propagation, and coalescence), displacement distribution (average velocity) of monitoring points, and released strain energy of the model indicates that the unloading rates and in situ stresses have a significant impact on the failure behavior and mechanical response of the rockmass around different shapes of holes. The severity of destruction and failure range will be intensified with the increase in unloading rates, although there is a clear distinction among the different shapes of holes. The intense and violent failure phenomenon was not observed for the circular hole, irrespective of the unloading rates, indicating that a stable failure in the form of slabbing failure is always the predominant failure pattern in this case. The possibility of strain rock burst (unstable failure) is increased for the case of larger unloading rate and for those tunnels with polygonal structure. The numerical results also demonstrate that the failure intensity presents a monotonically increasing trend with the increase in lateral pressure coefficients at a constant unloading rate, irrespective of the sectional shapes, which should be due to the higher stress level that precipitates the obvious stress concentration and strain energy accumulation induced by the excavation unloading process. The present study confirms that dynamic hazards are more prone to be induced in hard and brittle rock tunnels when they have a polygonal roadway section together with higher unloading rates and lateral pressure coefficients. In actual underground engineering, to reduce the risk of unexpected disasters with casualties as much as possible, it is recommended that tunnels with vertical wall configurations and rapid excavation processes (such as the blasting method) be avoided when they are concurrently subjected to high-stress states.