Analysis of fractures of a hard rock specimen via unloading of central hole with different sectional shapes

Hard rock failure and rockburst hazards under high in situ stresses have been the subject of deep rock mechanics and engineering. Previous investigations employed cubic rock specimens with a central hole for simulation of rock fracturing around deep tunnels at a laboratory scale, while the failure characteristics and crack evolution behavior around different shapes of holes induced by excavation unloading response have been barely considered. A commercially combined finite‐discrete element method (combined FEM/DEM) was used to investigate the failure characteristics and crack propagation process of typical hard rock specimens (marble) via the unloading of central hole with different shapes. Rock heterogeneity was also considered in the model in combination with the engineering reality. The combined FEM/DEM approach was first validated by simulating uniaxial compression and Brazilian tests. Then, the parametrical analysis was conducted in detail on the basis of five different sectional shapes of central holes, including a circle, ellipse, U‐shape, trapezoid, and cube, and two lateral pressure coefficients. Analysis of crack propagation paths, released strain energy, displacement, and average velocity distribution of the monitoring points around the central hole suggests that the failure degree and destruction intensity are strongly related to the sectional shape and lateral pressure coefficients. Hard and brittle rock failure induced by the excavation unloading effect can be classified as stable failure (slabbing failure) and unstable failure (strain rockburst). The cubic, trapezoidal, and U‐shaped holes within the specimen are the most likely to induce unstable failure, while stable failure is the dominant failure pattern around circular and elliptical holes. The lateral pressure coefficient λ was also found to affect failure position and intensity (only for the axisymmetric section) around the central hole. The influence of rock heterogeneity on failure intensity and range around the central hole within the specimen was also discussed. This study emphasizes the importance and necessity of the excavation unloading effect when evaluating surrounding rock failure around deep tunnels.


| INTRODUCTION
Mineral resource exploitation in the future has been prompted to access greater depths in the earth because of the exhausted reserves and energies at shallow depths. [1][2][3] Therefore, a large number of tunnels and chambers will be constructed and excavated in deep rockmass to meet the needs of underground engineering projects (Dammyr et al 4 ;Huang et al 5 ; Kang et al 6 ). Compared to the stress states of shallow buried depth, the mechanical characteristics of hard and brittle rocks at greater depth have been acknowledged to behave quite differently because of their complex stress states. 7,8 Because a large amount of energy is stored in the hard rockmass, the excavation unloading of deep tunnels and chambers will lead to stress redistribution and concentration, which create favorable conditions for rock fracturing and instability. 1,[7][8][9][10] In situ measurements have shown extremely unconventional failure phenomena in deep mining and tunneling engineering such as slabbing, zonal disintegration, and rockburst, which are of great threats to safe and efficient underground construction and activity. 9,[11][12][13][14][15][16][17] Therefore, it is necessary to gain better insight into the failure characteristics and mechanism response of hard rock around deep tunnels. 18,19 During recent years, some scholars have employed rock specimens with a central hole for simulation of failure behavior of deep tunnels at a laboratory scale and focused on the quantitative characterization of stress concentration and crack propagation behavior around the holes. [20][21][22][23][24][25][26][27] Li et al 28,29 investigated the failure properties of Iddefjord granite samples containing prefabricated holes under uniaxial compression tests at laboratory and numerical scales and rockburst phenomena such as rock block ejection, and slabbing failure similar to in situ measurements was observed around the prefabricated holes under high compressive stress. Zhou et al 30 analyzed the crack evolution law and failure behavior of marble specimens containing rectangular cavities under uniaxial loading. It was found that the number of holes and layout have a significant influence on the mechanical behavior of marble specimens. Yang et al 31 evaluated the strength, deformation, and crack evolution behavior of sandstone containing a single oval flaw under uniaxial compression using both experimental tests and numerical simulation. The results indicated that the oval flaw angle has a great effect on crack initiation and propagation and further explained the wing tensile crack propagation evolution law in detail. Zhu et al 32,33 investigated the effects of inclusions filled in prefabricated holes on the mechanical properties and fracture evolution of rock-like materials under uniaxial compression tests. They found both the inclusion type and shape are important factors affecting sandstone strength and deformation properties. Other researchers focused on slabbing failure and strain rockburst around deep tunnels under polyaxial stress states. Gong et al 34 investigated spalling rockburst characteristics using a true-triaxial test system, and the evolutionary process of rockbursts and V-shaped notches induced by spalling failure was fully monitored. Their experimental results indicate that initial stress states are of great significance in affecting rockburst tendency and failure degree. Zhong et al 35 performed biaxial compression tests to study the fracture properties of hard and brittle rock specimens containing a U-shaped opening. Crack initiation, propagation, and coalescence behavior were observed and analyzed in their tests.
Much effort has been dedicated to understanding the instability and failure of surrounding rock of deep roadway after unloading. [36][37][38][39] The mechanical response of rock around deep tunnels subjected to an unloading and loading mode should be quite different because of the disparate stress path between them. 7,8 When an underground opening is excavated, the accumulated energy in the hard rock will be released and deformation and fracture will occur near excavation boundaries. The excavation damaged zone (EDZ) is affected by multiple factors, such as unloading rate, rock lithology, tunnel shape, and in situ stresses. [40][41][42] It is widely confirmed that tunnel shape plays an important role in affecting the failure intensity and destruction area of tunnels at great depth. [23][24][25]43,44 This is primarily because the stress concentration, displacement, and deformation distribution are strongly related to the tunnel structure (straight line or curve) and cross-sectional area. Therefore, understanding unloading-induced failure around deep tunnels with different sectional shapes is necessary and urgent in reproducing the in situ underground engineering practice, which helps to provide a theoretical basis and technical guidance in supporting strategies and ground pressure control.
Because of the difficulties in reproducing the excavation unloading process during laboratory tests, numerical simulation has been considered an efficient and convenient approach for simulating the failure characteristics of hard rocks around deep tunnels. Manouchehrian and Cai 36,39 analyzed the role of weak planes around tunnels in rockburst occurrence and damage using the Abaqus two-dimensional (2D) explicit code. Li et al 38 proposed a mathematical physics model to characterize the unloading mechanisms of brittle rock under different stress paths in two dimensions using the universal discrete element code PFC2D. In their study, the influence of various in situ stresses and the unloading rate and path on the dynamic effects of a circular tunnel was fully investigated. Cao et al 37 developed a mathematical physics method for a two-dimensional circular opening to investigate the unloading vibrational mechanism. In their study, the unloading waveform subjected to great depth was evaluated for various lateral pressure coefficients and aspect ratios of rectangular tunnels. Wu et al 45  modeling techniques was proposed to capture the mechanical response of the entire TBM unloading process. Nevertheless, the entire failure process from continuum to dis-continuum in these studies was not captured. Therefore, it is of vital importance to employ a rational numerical approach that can better reflect crack initiation, propagation, and coalescence of rockmass 9,46 from the continuum state to dis-continuum domain. During recent years, the combined finite element method/ discrete element (FEM/DEM) approach, which integrates the finite element and discrete element technologies into a singular software, has been gradually confirmed and adopted by scholars in addressing issues related to rock mechanics and rock engineering. 47,48 For instance, Hamdi et al 49 adopted an integrated finite/discrete element method-discrete fracture network approach to analyze the Kiruna hanging wall surface subsidence with an emphasis on investigating the influence of discontinuity persistence and spacing. Trivino and Mohanty 50 assessed single hole-induced damage of a granitic outcrop through an experimental and FEM/DEM approach and the experimental and numerical results well agreed. Li et al 1 adopted an FEM/DEM approach to investigate the crack behavior of a splitting test of a circular ring and demonstrated that the failure mode and peak load are strongly related to the ratio of internal to external diameter.
Though an FEM/DEM approach has been used to model the rock failure process in rock mechanics, limited attention has been paid to the integrated response of the excavation unloading effect, sectional shape, and initial stress states (thus the lateral pressure coefficients) on hard rock failure around deep tunnels at a laboratory scale. In the present study, we adopted an FEM/DEM approach on the basis of a commercially available software to simulate typical hard specimen (marble) fracturing containing different shapes of central holes. The rock heterogeneity and excavation unloading process were realized by customizing material parameters and stress paths within the software. The entire failure process, released strain energy, maximum tangential stress, displacement, and average velocity stress of selected monitoring points around the excavation boundaries were analyzed and compared in detail. This study is intended to deepen our understanding of tunnel failure and potential dynamic hazards at great depth.

| Brief introduction of the combined finite element method/discrete element method approach
The FEM/DEM approach, which was proposed by Munjiza 51 , integrates the aspects of both finite elements and discrete elements to model multiple interacting deformable solids (Mietelman and Elmo 2014). The unique character of the FEM/DEM approach is that the transition from the continuum to dis-continuum state can be well reproduced through an entire failure process, which differs from the FEM (finite element method), DEM (discrete element method), and other approaches. The basic principal of the combined FEM/DEM approach is to divide the solution domain into discrete elements along the existing discontinuous surface and further divide the finite element meshes within each discrete element. By introducing damage mechanics and fracture mechanics theory into the FEM, the cracking of discrete elements is designated to occur on the boundary of the finite element. The motion and interaction of discrete blocks operate in the same manner 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 freely move and rotate, but there is no need to satisfy the deformational coordination condition. In this study, a commercial numerical software based on the combined FEM/DEM approach was adopted for modeling the failure process of a hard rock specimen containing a central hole with different sectional shapes induced by an excavation unloading effect. The Mohr-Coulomb constitutive model with rotating cracks was selected as it can simultaneously simulate rock fracturing resulting from both tension and compression. There are five plastic material parameters as follows: cohesion (c), friction angle (φ), dilatancy angle (ψ), tensile strength (σ t ), and fracture energy (G f ). The Mohr-Coulomb yield criterion is a generalization of the Coulomb friction failure law and is defined as follows: where τ is the magnitude of the shear stress, σ n is the normal stress, c is the cohesion, and φ is the friction angle. In a principal stress space, the yield surface is of a six-sided conical shape. The conical nature of the yield surface reflects the influence of pressure on the yield stress, and the criterion is applicable to rock, concrete, and soil problems. The Rankine tensile corner introduces additional yield criteria as follows: where σ i refers to each principal stress and σ t is the tensile strength. Although at present no explicit softening law is included for the tensile strength, indirect softening does result from the degradation of cohesion according to the following criterion: (1) = c − n tan , A conventional Mohr-Coulomb tensile yield surface is shown in Figure 1A. Because the Mohr-Coulomb criterion may overestimate the tensile strength of rocks, Paul 52 introduced the concept of a tension cutoff and modified the Mohr-Coulomb criterion, as shown in Figure 1B. To model the hardening or softening behavior, the material strength parameters are defined as functions of the effective plastic strain as follows: where ɛ −p is the effective plastic strain and In the present study, the linearized functions were adopted. In the Mohr-Coulomb with a rotating crack constitutive model, it is assumed that extensional failure can be considered as the dominant factor in rock brittle failure [18,22]. Under the circumstance of compression, the increment of inelastic strain Δ p 3 in the minimum principal stress direction can be associated with tensile strength degradation (damage) in the same direction, giving where σ t3 is the tensile strength in the direction of the minimum principal stress and n is the loading step. For the linear softening of the tensile strength, Equation 5 can be written as follows: where σ t is the initial tensile strength and H is the linear softening modulus. H is defined as c is the local element length scale (average dimension) and G f is the fracture energy. Once the yield surface represented by Mohr-Coulomb is satisfied at the elemental level, the tensile strength according to is updated. The fracture energy G f can be defined as follows: Actually, the specific fracture energy (usually referred to simply as the "fracture energy") G f is defined as the amount of energy needed to create a continuous crack in a unit area. The fracture energy, G f · A, can be obtained through a simple uniaxial tensile test. The total fracture energy is the area under the load-displacement curve. Note that A is the crosssectional area of the test specimen. Integrating over a localization bandwidth l c for a constant slope softening model provides the following: where in the finite element context the local control length l c = f(B e ), B e is the area of the element and E T is the tangential softening modulus. Note that the relationship between fracture toughness in mode I and fracture energy can be defined as follows: where K 2 Ic is the fracture toughness in mode I and E is the initial Young's modulus.
The Rankine rotating crack constitutive model is suitable for modeling the tensile failure of brittle materials, for example, rock, ceramic, and glass. The most important parameters in the Rankine rotating fracture model are fracture energy and tensile strength. The failure mechanism of brittle rocks is dominated by the formation of tensile cracks (mode I) when the tensile stress exceeds the tensile strength of the material. For mode I-dominated problems, the initial failure surface for both models is defined by a tension failure surface as follows (shown in Figure 2): where σ i is the principal stress invariant and f t is the tensile strength of the material. In the postinitial yield, the rotating crack formulation represents the anisotropic damage evolution by degrading the elastic modulus in the direction of the major principal stress invariant as follows: where ω is the damage parameter, nn is the local coordinate system associated with the principal stresses, E is the Young's modulus, and E d is the elastic damaged secant modulus. The damage parameter for the linear strain softening associated with microfracturing is given by the following: where ψ(ε) is a function of the total strain defined as follows: where E t is the tangential softening modulus. The advantages and applicability of using this constitutive model have been validated by previous literature. 49,53-55

| Numerical validation of the combined FEM/DEM
In this study, uniaxial compression and Brazilian tests of a typical hard and brittle rock were simulated. The main purpose was to validate the accuracy and rationality of the combined FEM/DEM approach. For a 2D geometric type, the specimen thickness was assumed to be 1 by default. The specimen dimension of the uniaxial compression test and Brazilian test was 100 × 50 mm and φ 50 mm (diameter), respectively. Displacement controlled loading was imposed on the top-loading platen at a constant velocity. A ramp load style was applied with a total displacement of 3.0 and 1.0 mm at the top of the platen, respectively. Note that the specimen may have been completely fractured when the applied displacement was less than the critical value. The bottom platen was fixed in both the x and y directions in the two laboratory tests. Actually, the models with a lower fracture energy required less deformation before visible cracks were simulated; thus, hard rocks with brittle properties are suitable for such a case. The fracture energy can be obtained based on Equation 10, and the fracture toughness K IC can be estimated using an empirical relationship with the rock tensile strength as follows (Mitelman and Elmo 57 , Zhang 56 ): The fracture energy adopted in this study was calculated to be approximately 0.005 N/mm. Moreover, the hard and brittle properties of the selected specimen can be well reproduced following the cohesion weakening and frictional strengthening (CWFS) model proposed by Hajiabdolmajid et al 58 . To realistically represent the intrinsic properties of the underground engineering rockmass, 10,32,33 the material heterogeneity was fully considered and specified in the software. Rock heterogeneity was applied on each element via elemental basis, allowing each to have an independent value assigned for the five material parameters. A uniform distribution type was employed to characterize the material heterogeneity of the marble specimen. The distributions of Young's modulus, Poisson's ratio, density, tensile strength, and fracture energy within the for Brazilian test specimen are shown in Figure 3. The scale in the legend indicates the ratio of each parameter to the average value. The detailed physical and mechanical parameters of the marble specimen used in this study are listed in Table 1. Figure 4 shows plots of the typical failure modes of the marble specimen under the uniaxial compression test and Brazilian test. Under a uniaxial compression condition ( Figure 4A), a large number of cracks were generated in the specimen. These newly generated cracks were nearly parallel to the loading direction, which is a typical property of mesoextensional cracks. As the load increases, these cracks tend to propagate toward the diagonal direction within the specimen and finally form several macro shear bands. Therefore, the final failure mode under the uniaxial compression test was macroscopic shearing. Under the Brazilian test, the macro diametrical splitting can be clearly observed within the specimen. The macro diametrical cracks finally split the disk into two pieces, in which only vertical diametrical cracks initiate and propagate toward the upper and lower platens (see Figure 4B). The current numerical results were shown to be similar to what was observed in previous literature and findings, 28,29,[59][60][61][62] which validates the effectiveness and accuracy of the selected combined FEM/DEM approach. Figure 5A shows the typical stress-strain curve of the marble specimen under the uniaxial compression test. It can be seen that the nonlinearity before the peak strength point is barely observed. After the peak load, the curve suddenly collapsed with no obvious postpeak curve, which is an indication of the hard and brittle properties of the selected rock lithology. The loaddisplacement curve under the Brazilian test was also plotted as shown in Figure 5B and shows similar trends compared to the stress-strain curve of the uniaxial compressive tests. Note that the calculated uniaxial compression strength (60.4 MPa) and Brazilian tensile strength (2.55 MPa) in this study were slightly smaller than the anticipated theoretical strength. Actually, the testing results of the uniaxial compression strength and Brazilian tensile strength were 70 and 3.1 MPa, respectively. The distinction between the numerical and testing results can be attributed to the rock heterogeneity, which leads to the strength anisotropy of elements in the numerical model. Therefore, the stress concentration at certain locations in the marble specimen is further promoted. Another reason leading to weakening is because of the selected 2D computing mode used in the numerical simulation. According to the description of Hamdi et al 55 , 2D simulations provide lower strength values compared to those of three-dimensional (3D) simulations. These differences were interpreted to originate from an increased degree of complex interlocking between the elements in 3D simulations.

| Numerical model implementation
The geometry of the models used in this study is shown in Figure 6. Similarly, the model assumes a unit thickness to be 1 mm in 2D strain plane analysis. The rockmass is represented by a single face with a central hole representing the tunnel to be excavated. The area of the central hole was designated to be a constant of 1963 mm 2 , irrespective of the sectional shape. Five sectional shapes of deep tunnels of a circle, cube, U-shape, trapezoid, and ellipse were considered in the present study. Theoretically, the width and height of the outside boundary should be at least 10 times the tunnel diameter to exclude the boundary effect of the stress redistribution and concentration around the tunnel rockmass. 39 Additionally, for a 2D tunnel dynamic excavation problem, the unloading stress waves may be reflected back from the outside boundaries toward the central hole given the fact that the model dimension is finite. 9 To reduce the boundary effect induced by reflected stress waves and stress redistribution on the numerical results as much as possible, the outside model dimension was set to be 500 × 500 mm. Some typical elements around the roof, floor, and sidewalls of the central hole were selected to monitor the displacement and average velocity induced by the unloading effect. Because stress concentration resulting from the unloading process was expected around tunnel boundaries, very dense meshes were used in these nearby areas. The remainder of the regions that were far from the central hole was designated to be sparser to reduce computational duration. The total number of unstructured triangle elements established in the model containing circular, elliptical, U-shaped, trapezoidal, and cubic holes was 9700, 9640, 9678, 9750, and 9750, respectively.
The normal and tangential penalties equivalent to the normal and shear stiffness were used to evaluate the normal and tangential contact forces. 49 The normal penalty was usually within the range of 0.5E to 2.0E, where E is the Young's modulus. In this study, the normal and tangential penalties were set to be E and 0.1E. 53,55 The contact field defines the thickness of the contact layer, which corresponds to the maximum permissible penetration. The smallest element size should be smaller than the original mesh element size so that the generated cracks cut across the meshes instead of merely along the mesh boundaries. The discrete contact parameters are listed in Table 2.
In the current 2D plane analysis, in situ stresses were applied to the model boundaries. Two lateral pressure coefficients at k = 0.5 and 2 were selected. Vertical stress was applied on the top domain, while horizontal stresses were assigned on both sides of the model. The bottom domain was fixed in the y direction. The in situ stresses were applied by using a linear loading path with a total duration of 0.01 second. The purpose of employing such a loading path was to avoid unexpected fracturing when employing instantaneous loading. The in situ stresses were maintained constant once reaching the predefined value. The excavation unloading process was then executed after 0.01 second on the central hole until 0.02 second. During the excavation process, the internal surface force of the central hole was gradually released according to the linear relaxation curve. This avoided a sudden change in stress state, which would occur with instantaneous material removal. 9 The loading and unloading paths are shown in Figure 7.

| RESULTS AND DISCUSSION
In this study, the failure characteristics of the central hole within the marble specimen representing a deep underground tunnel subjected to two in situ stresses were investigated. For the condition of a lateral pressure coefficient k = 2, the vertical stress is 25 MPa and the horizontal stress is 50 MPa, while for the case of k = 0.5, the vertical and horizontal stresses are assumed to be 50 and 25 MPa, respectively.
To verify the rationality of the selected model dimension, we plotted the maximum principal stress contour of the marble specimen containing a circular hole at k = 2, as can be seen in Figure 8. The stress redistribution and concentration induced by the excavation unloading process can be clearly observed around the central hole. The maximum tangential stress appearing in the roof and floor of the hole began to increase while the maximum tangential stress at both sides of the circular hole gradually decreased. This figure shows that the stress redistribution and reflected stress wave had not reached the model boundaries and therefore did not influence the numerical results. Figure 9 shows the plots of the failure process and crack propagation around the central hole with different sectional shapes induced by the excavation unloading effect when k = 2. Crack initiation, propagation, and coalescence around the central holes are shown in detail. By observing these dynamic evolutionary processes, it is indicated that the crack initiation stage typically appears prior to the point at which the internal surface load is unloaded to 0 (corresponding to t = 0.02 second). As time passed, a visible crack was found in the roof and floor of the central hole near the excavation unloading boundaries. This was expected because of the predefined in situ stress states. According to the description by Feng et al, 9 when the lateral pressure coefficient is greater than 1, the maximum tangential stress is concentrated in the tunnel roof and floor and the radial stress is reduced to zero, which leads to higher strain energy accumulation in these areas. Figure 10 shows the plots of the maximum principal stress contour of the marble specimens with different central hole shapes when k = 2. It can be seen that the maximum principal stress in the roof and floor of the central hole reached approximately 110 MPa, irrespective of the sectional shapes. Under such conditions, the maximum tangential stress far exceeded the uniaxial compression strength of the marble, which led to macro failure in these areas. However, a more in-depth analysis of the crack propagation behavior shows that the failure intensity and damage range were strongly dependent on the sectional shapes of the central hole.

| Failure process and crack propagation behavior
It can be seen from Figures 9 and 10 that the newly generated cracks around the cubic and trapezoidal holes most easily propagated, followed by the U-shaped, elliptical, and circular holes. For the case of a circular hole, only a few cracks were generated in the roof and floor. Typical slabbing failure could be found in these areas, which was quite similar to that observed in previous literature and in situ underground engineering. 1 As was stated by Ortlepp, 63 spalling or slabbing is generally defined as the formation of stress-induced slabs on the boundary of an underground excavation. It initiates in the region of maximum tangential stresses and results in a V-shaped notch that is local to the boundary of the opening. This study further illustrates that slabbing is the primary failure phenomenon for hard and brittle rock near underground excavation boundaries. For the elliptical hole, T A B L E 2 Discrete contact parameters adopted in the current study Normal penalty (P n /N/mm 2 ) 35810  the failure intensity and damage extent were found to be more severe compared to that of the circular hole, which might be attributed to the large cross-sectional span in the roof and floor that further weakened the self-stabilization of the surrounding rock. In this case, a slight strain rockburst phenomenon could be found in the floor, with several slabs splashing toward the free surface. It is also found from Figure 9B that the crack propagation was not asymmetrically distributed around the elliptical hole. This might have been because of the material heterogeneity that led to strength anisotropy and unbalanced strain energy accumulation around the elliptical hole. Figure 9C shows the plots of the entire failure process of the rock surrounding the U-shaped hole. It can be seen that when t = 0.0195 second (corresponding to the crack initiation stage), visible cracks could be generated in the floor, while the roof was not fractured. Compared to the linear section, the arch structure is considered to avoid obvious stress concentration and hence possess a higher bearing capacity. With further unloading, visible cracks were generated in the roof, but the failure extent and damage degree in the floor always prevailed. When t = 0.0304 second, one can clearly observe the strain rockburst phenomenon around the U-shaped hole, particularly near the bottom excavation boundaries. At this moment, a large number of slabs were ejected toward the center of the hole. According to the description by Tan, 64 the progressive failure process of strain rockburst can be summarized in three stages. During the first stage, the surrounding rock is split into several slabs parallel to the excavation boundaries affected by the stress concentration. Subsequently, these slabs are cut into blocks of different length and size. Finally, these blocks and slabs are ejected toward the free surface because of the coupled buckling (maximum tangential stress) and pushing (internal radial stress). The current numerical results agree well with previous speculation, which further verifies the strain rockburst mechanism. For the case of a trapezoidal hole, the failure intensity and extent were further aggravated compared to that of the U-shaped hole, as can be seen in Figures 9D and 10D. This might be a result of the polygonal structure facilitating a higher stress concentration around four undesirable corners and the large span in the floor further intensifying the instability of the surrounding rock; thus, a strain rockburst is more prone to be induced in these regions. Figure 9E shows the crack propagation process around a cubic hole. The crack simultaneously initiates in both the roof and floor. With further unloading, large numbers of cracks can be observed around the excavation boundaries. The strain rockburst is triggered once the accumulated strain energy exceeds the surface energy required for rock fracturing. 65,66 Because of the unwanted corners, the cubic hole will cause severe roof caving and floor heaving and even intense dynamic hazards. For the trapezoidal tunnel in deep underground engineering, the most serious areas should be around the floor (strain rockburst), while the roof always presents a stable failure state in the form of slabbing. For a cubic tunnel, the failure intensity and extent in the floor may be lower than that of the trapezoidal tunnel, while the roof is damaged to a great extent, which is unfavorable to roof maintenance and management, particularly in a deep coal mine. Figure 11 shows the plots of the entire failure process and crack propagation behavior around the central holes with different sectional shapes when the lateral pressure coefficient k = 0.5. It can be seen from Figure 9 that most visible because of the material heterogeneity, which leads to stress deviation and unbalanced energy distribution in certain regions. However, the anisotropy caused by material heterogeneity will not appreciably influence the numerical results because the initial material parameters were set to be uniformly distributed. For the cases of the U-shaped and trapezoidal holes, the failure intensity and destruction range tended to be uniform and analogous on both sides. This is different compared to the case when k = 2 (asymmetrical failure). For these centrally asymmetrical tunnels, the lateral pressure coefficient will not influence the failure intensity and damage destruction, but only the failure position. While for those tunnels with an axisymmetric cross section, such as the U-shaped and trapezoidal tunnels, not only does the failure position change, but the failure intensity and damage range are also affected to some extent. Zeng et al 44 investigated the failure characteristics of sandstone specimens containing different hole shapes under uniaxial compression by using laboratory tests and numerical simulation. Figure 12 shows the final failure pattern of specimens based on previous and current numerical results. It can be seen that the failure pattern and crack propagation paths were different between the two numerical results. In their tests, shear zones consisting of a series of tensile cracks were observed to run along the specimen from the hole to the corners for all the specimens. Moreover, tensile splitting cracks were initiated in the roof and floor of the central hole, except that containing an elliptical hole. However, the shear zones were not observed in the current numerical study (k = 0.5) and slabbing failure was the dominant failure pattern on both sides of hole walls. The different results between the previous and current research might be because of the following: 1. In the previous tests, the horizontal stress was not considered, and only vertical stress was applied to compress the specimen. Actually, deep rockmass is always subjected to a polyaxial stress state. Therefore, the tensile splitting cracks on the roof and floor and the shear zones along the diagonal direction are not fully restricted in their numerical results. 2. The excavation unloading process was not considered in previous research. This is important because deep tunnel failure is always induced by engineering disturbances such as mechanical excavation or blasting loads. Generally, the stress redistribution and concentration around deep tunnels and chambers is induced by the excavation unloading effect instead of monotonically compressive loads.
The aforementioned comparison and analysis suggest that the excavation unloading effect and polyaxial stress states should be considered when investigating the failure behavior and mechanical characteristics of deep tunnels.

| Displacement and average velocity of monitoring points around the central hole
In this study, we plotted the displacement-time evolution curves of selected monitoring points around the central holes when k = 2, as can be seen from Figures 13-17. The purpose was to investigate the failure intensity and extent of deep tunnels affected by sectional shapes induced by the excavation unloading process. Because the displacement value before t = 0.02 second was nearly 0, the abscissa starts at 0.018 second (except for the case of the circular hole). Given that the ejected elements may collide or come into mutual contact after a long unloading duration, we only plotted the displacement-time evolution curves before t = 0.03 second. The monitoring points were selected near the roof and floor of the central hole because of the possible fracturing and visible cracks in these areas. Figure 13 shows the displacementtime evolution curves of the monitoring points around the circular hole. It can be seen that the displacement values in the roof present a linear increasing trend during the period 0-0.01 second because of the applied in situ stress that leads to initial squeeze and settlement. The displacement values continue to increase during the initial unloading stage. The displacement presents a rapid increasing trend when the excavation unloading process is completely terminated. After t = 0.02 second, the displacement values for all the points remained constant though there is a disparity among them. The displacement-time evolution curves in the floor show a distinct variation trend compared to those in the roof. The displacement values corresponding to point 7, 278, and 275 present the trend of initial increase, then decrease, and finally stabilization. This is because the opposite motion trend in the y direction, induced by the excavation unloading effect, offsets a portion of the displacement. Moreover, the maximum displacement values of these monitoring points are always maintained at a low level (within the range between 0.25 and It can be seen that the displacement values at t = 0.03 second greatly increased in the case of the trapezoidal and cubic holes, indicating that the ejected distance of these elements is too far from the original position. Under such a condition, the failure is extremely intense and could be treated as an unstable failure (eg, strain rockburst).
To further illustrate the failure intensity and damage extent affected by the integrated response of the sectional shapes and excavation unloading process, the average velocity of selected monitoring points was calculated on the basis of a corresponding maximum displacement value, as can be seen in Figure 18. For the circular hole, the monitoring points around the floor were not considered in this study because of their complex displacement-time evolution. It can be seen from Figure 18A that the average velocities for all the monitoring points were maintained at a very low level, though the roof was fractured with visible cracks. In this regard, slabbing failure was the dominant failure pattern around the circular hole. For the case of the elliptical hole, the average velocity of point 283, 277, 280, and 8 exceeded 0.7 m/s, indicating that the failure intensity is somewhat aggravated. Under such a condition, a slight rockburst phenomenon is observed near the excavation boundaries. Figure 18C shows that the average velocities around the U-shaped hole are further promoted, particularly for those points in the floor. Compared to the case of the circular elliptical hole, the rockmass in the roof  36,39 investigated the average velocity distribution around a deep circular tunnel subjected to excavation unloading using Abaqus. They found that when the maximum average velocity was 0.14 m/s, rock failure can be considered stable failure in the form of spalling, spitting, or shallow slabbing. When the average ejected velocity increased to 3.4 m/s, the failure transformed to unstable failure. The current numerical results show similarity when compared to their findings. Note that the classification of rockburst grades (slight rockburst and intense rockburst) and failure intensity (slabbing failure and strain rockburst) is merely through a relative comparison of the monitored maximum displacement and average velocity values among each case, which is not based on the exact border in this study.

| Released strain energy of specimens with different sectional shapes of central holes
Rockmass instability and fracturing can be induced by underground engineering activities, such as deep mining and tunneling projects, because of the release and transfer of energy in the surrounding rock. 67 Therefore, it is of great significance to describe the failure and damage around deep tunnels or chambers from an energy perspective (Sanchidrian et al 67  termination of the unloading point (t = 0.02 second). The released strain energy remained constant upon reaching the peak value. The maximum released strain energy values of the marble specimens with different shapes for different lateral pressure coefficients are also shown in Figure 21.
It can be observed that the released strain energy induced by the excavation unloading process largely depends on the sectional shape and lateral pressure coefficient k. When k = 2, the released strain energy for the case of the circular, elliptical, U-shaped, trapezoidal, and cubic holes was 90.1, 116, 224.6, 274, and 345.6 J, respectively. The more strain energy that is released, the more severe the failure intensity and damage extent of the surrounding rock. In the case of k = 0.5, the sequence of the released strain energy is shown to be different from the case of k = 2. The sequence of the released strain energy for the case k = 0.5 changed to E cube > E U-shape > E trapezoid > E circle > E ellipse . Under such a condition, the minimum released strain energy does not correspond to the circular hole but the elliptical hole. The corresponding energy value of the specimen containing the elliptical hole is 75 J, which is less than the energy value of the circular hole of 103.6 J. Because fracturing is mainly concentrated on both sides of the hole walls, the limited cross-sectional span of the elliptical hole on both walls will prevent visible cracks from propagating and coalescing. Therefore, the failure intensity and damage range are restricted to some extent. Conversely, the U-shaped hole presents more severe fracturing on both holewalls compared to that of the case of the trapezoidal hole when k = 0.5, as can be seen in Figure 21. This is attributed to the tilt walls always having a higher bearing capacity compared to that of the vertical walls. The released strain energy of the specimen containing a U-shaped hole when k = 0.5 is greater than that of the case of k = 2. Similar results could also be found in the specimen containing a trapezoidal hole. Table  3 lists the normalized released strain energy increase for different central hole shapes. The normalized released strain energy increase for the case of circular and cubic holes is less, at only 13% and 7.76% (the lower disparities were mainly influenced by the material heterogeneity in the marble specimen), respectively. The normalized increase for the remaining three shapes is somewhat higher, particularly for the case of the elliptical hole (54.7%). This was attributed to the transformed failure position being affected by different lateral pressure coefficients and the asymmetry of the structure in the horizontal axis that leads to the distinct failure characteristics and mechanical behavior of the surrounding rock around the central hole within the specimen. Current numerical results suggest that the failure intensity and damage range around deep underground engineering tunnels and chambers should be evaluated on the basis of the integrated response of the sectional shape, excavation unloading process, lateral pressure coefficient, and rock heterogeneity when evaluating deep underground excavation engineering.

| CONCLUSIONS
In this study, the entire failure process (including crack initiation, propagation, and coalescence behavior) of a typical hard rock specimen (marble) with unloading and different central hole shapes was reproduced by employing a commercially combined finite/discrete element method. The main purpose was to explore the mechanical response properties and crack evolution behavior of a deep tunnel affected by tunnel shapes and lateral pressure coefficients during the excavation unloading process at a laboratory scale. The linear excavation unloading process and rock heterogeneity were fully considered to better reflect the practical underground engineering background. Five commonly used tunnel shapes and two lateral pressure coefficients were selected. Numerical results show that when the lateral pressure coefficient k = 2, failure was always observed in the roof and floor of the central hole, while the failure position was transformed to the sidewalls when k decreased to 0.5. The visible cracks are always parallel to the excavation boundaries, which validates the essence of the extensional failure in hard and brittle rock. Analysis of the crack propagation behavior and displacement (average velocity) distribution indicates that dynamic hazards are more likely induced around the U-shaped, trapezoidal, and cubic holes, which might be because of the higher stress concentration around undesirable corners. The slabbing failure (stable failure) was always the dominant failure pattern around the circular and elliptical holes. The failure intensity and damage range show a distinct variation trend in terms of different lateral pressure coefficients. The sequence of released strain energy when k = 0.5 is E cube > E U-shape > E trapezoid > E circle > E ellipse but E cube > E trapezoid > E U-shape > E ellipse > E circle when k = 2. It is the noncentrosymmetric structure of the U-shaped, trapezoidal, and elliptical holes that leads to the variation in failure intensity and damage extent of the marble specimen. Therefore, the initial stress states should be considered when determining appropriate tunnel shapes in deep underground engineering to avoid potential failure as much as possible. This study confirms the differences in crack evolution behavior between compression-and unloading-induced fracturing around tunnels with various sectional shapes and also