Research on Zonal Disintegration Characteristics and Failure Mechanisms of Deep Tunnel in Jointed Rock Mass with Strength Reduction Method

Abstract

To understand the fracture features of zonal disintegration and reveal the failure mechanisms of circle tunnels excavated in deep jointed rock masses, a series of three-dimensional heterogeneous models considering varying joint dip angles are established. The strength reduction method is embedded in the RFPA method to achieve the gradual fracture process, macro failure mode and safety factor, and to reproduce the characteristic fracture phenomenon of deep rock masses, i.e., zonal disintegration. The mechanical mechanisms and acoustic emission energy of surrounding rocks during the different stages of the whole formation process of zonal disintegration affected by different-dip-angle joints and randomly distributed joints are further discussed. The results demonstrate that the zonal disintegration process is induced by the stress redistribution, which is significantly different from the formation mechanism of traditional surrounding rock loose zone; the dip angle of joint set has a great influence on the stress buildup, stress shadow and stress transfer as well as the failure mode of surrounding rock mass; the existence of parallel and random joints lead the newly formed cracks near the tunnel surface to developing along their strikes; the random joints make the zonal disintegration pattern much more complex and affected by the regional joint composition. These will greatly improve our understanding of the zonal disintegration in deep engineering.

1. Introduction

With the increasing development of global economy, the shallow resources and space are getting harder and harder to satisfy the growing demands of human beings for a better life. In recent decades, the deep-buried underground resources and transport infrastructure in mountain areas, especially those that are difficult to excavate or construct in the past, have become new growth points in energy exploitation, civil engineering and related fields. For example, the maximum depth of several tunnels on the Sichuan-Tibet Railway, China are over 1000 m, such as the Baxu Tunnel and Sangzhuling Tunnel [1,2], the maximum depth of the diversion tunnel at Jinping II Hydropower Station, China is around 2500 m [3] and the maximum depths of the Taimyrskii Mine in Russia [4] and Witwatersrand Mine in South Africa [5] are over 1000 m and 2000 m, respectively. In deep engineering, the complex geological conditions involving high ground temperature, high ground stress, high permeability and strong excavation disturbance leads to many characteristic geological problems, which are quite different from shallow buried underground engineering [6,7,8,9,10]. As one unconventional fracture phenomenon, the zonal disintegration phenomenon of deep surrounding rock has drawn a lot of experts’ attention [11,12,13]. In 1970, this phenomenon was firstly discovered in the roof of a gold mine in South Africa [5], which is characterized by the interval distribution of fractured zone and intact zone of roof rocks. In 2008, the zonal disintegration phenomenon in Dingji mine—910 m of Huainan mining area in China [12,14,15] was monitored by the borehole TV imager, as shown in Figure 1.

The formation mechanisms of this nonlinear deformation phenomenon are quite different from the traditional loosened zone of surrounding rock, which results in the difficulties of stability assessment and support design of deep tunnel. Therefore, it is necessary to reveal the mechanical mechanism and fracture characteristics of the zonal disintegration phenomenon. Through the theoretical analysis, numerical simulations and laboratory tests, some researchers [12,16,17,18,19] have discussed the formation conditions of zonal disintegration in deep engineering and suggested that when the tunnel axis is parallel to the direction of the maximum principal stress, the zonal disintegration is prone to occur. Based on this understanding, Zhang et al. [20] successfully reproduced the formation process of zonal disintegration through laboratory test using a true three-dimensional high-stress loading system, monitored the radial strain and displacement development and summarized the alternation distribution law of their peaks and troughs. Zhang et al. [21] conducted the zonal disintegration laboratory test using the soft-layered joint model with different spacing distances, and pointed out that the zonal disintegration is more obvious due to the weak interlayers considering the greater radial displacement and strain. Simultaneously, when the spacing of interlayers is denser, a greater number of fractured zones and larger damaged area will be observed. Pu and Xu [22] monitored the three-dimensional geological model with different uniaxial compressive strengths using the strain gauges and failure wires and found that during the loading process, both the radial tensile strain and the tangential compressive strain increased. According to their measurement results, the radius scale factor of the fracture area was also summarized. Factually, as a unique phenomenon in deep engineering [23], the zonal disintegration phenomenon generally appears in the complex environment of “high ground stress, high ground temperature, high permeability” [6,7,8], which directly leads to many difficulties in field and laboratorial tests, such as long-period, high-cost and hazardous operating environment and so on. Fortunately, the numerical simulation technology can effectively avoid these problems. By comparing the results of field monitoring, laboratory test and numerical simulation, the feasibility of the numerical simulation technology in studying zonal disintegration phenomenon has been verified [12]. For instance, Zhu et al. [24] used the finite element method to analyze the influence of rock physical and mechanical parameters, including internal friction angle, cohesion and Poisson’s ratio, on the fracture area and failure patterns of zonal disintegration phenomenon. Zhang et al. [25] studied the effect of confining pressure on the crack occurrence order and final failure range by computing the formation of zonal fractures using the 2D rock failure process analysis (RFPA) method. Wu et al. [26] applied the particle flow code to simulate the zonal fracture phenomenon of brittle rock mass, and carried out the sensitivity analysis of initial stress, lateral pressure coefficient and joint friction coefficient.

However, the zonal disintegration is a 3D failure phenomenon involving the gradual process of stress buildup, stress shadow and stress transfer in the 3D space and the initiation, propagation and coalescence of 3D cracks. Currently, most of the simulated results are obtained by 2D numerical methods and cannot effectively reflect the influence of mechanical interaction of cracks in the 3D space on the formation of zonal disintegration by taking the rock heterogeneity and discontinuity into account. Therefore, in this study, to reveal the fracture characteristics and failure mechanisms of deep jointed rock masses, a series of 3D heterogeneous numerical models considering varying joint dip angles and random spatial distribution of joints are built up and the strength reduction method (SRM) is embedded to reproduce the zonal disintegration phenomenon of deep rock masses. The discontinuous deformation field, evolution rule of background stress field, alternate regional destructions and safety factor are therefore discussed, which provides theoretical basis of stability assessment, support design and safe construction in terms of zonal disintegration in deep engineering.

2. Materials and Methods

2.1. Basic Principles

The current numerical methods applied in rock mechanics can be categorized into three groups [27], i.e., continuum methods, discontinuum methods and hybrid continuum/discontinuum methods. The finite difference method (FDM) [28], finite element method (FEM) [29], boundary element method (BEM) [30], etc. are classical continuum methods. However, they generally cannot satisfactorily deal with the discontinuities widely existing in rock mass. The discrete element method (DEM) [31] and discontinuous deformation analysis (DDA) [32] are two mainstream discontinuum methods. However, the model block division for these methods often requires fine engineering geological survey and has a great influence on simulated results; the divided blocks generally cannot be broken during calculation and the much lower computing capacity and efficiency than continuum methods put limitations on the application of this kind of methods. The hybrid methods provide a feasible way of combining the strengths of both continuum methods and discontinuum methods. However, they still cannot handle the interaction of multiple joint sets in the 3D space. As a FEM-based method, RFPA is able to take the heterogeneity, nonlinearity and anisotropy of rock materials into account. Simultaneously, to simulate the evolution process of progressive failure of deep tunnels, the fundamental principles of strength reduction method are embedded into the RFPA3D code. Namely, the failure modes and safety factor can be gained without any assumption on the failure surface in a complex rock mass model.

To establish a 3D heterogeneous numerical model, a rock sample or structure will be discretized into hexahedral elements at first, which can be regarded as the basic mechanical units of the rock medium in the process of numerical simulation depending on the premise that the element size can reasonably reflect the characteristics of the rock medium.

In addition, with the aim of appropriately capturing the mechanical behaviors of rock mass, the non-uniformity of rock mass cannot be ignored. To fully reflect the heterogeneity of mesoscopic elements after discretization of rock mass model, the material properties of these elements, such as uniaxial-compressive strength and elastic modulus, are assumed to obey a statistical distribution function. As a common probability density function used in the field of rock mechanics, the Weibull distribution [33,34,35] is as shown in Equation (1).

$$\varphi \left(\alpha \right)=\frac{m}{{\alpha }_0}\cdot {\left(\frac{\alpha }{{\alpha }_0}\right)}^{m-1}\cdot e^{-{\left(\frac{\alpha }{{\alpha }_0}\right)}^m}$$

(1)

where $\alpha$ is a mechanical parameter of rock medium, such as elastic modulus, strength, Poisson’s ratio, weight, etc.; ${\alpha }_0$ is the mean value of the mechanical parameter; m is termed the heterogeneity coefficient reflecting the uniformity degree of rock material.

The strength reduction method is introduced into the solution process. At each reduction step after all loads are applied, the initial strength of one element will be reduced according to the following criterion [36]:

$$f_{\mathrm{s}}^{\mathrm{trial} }=\frac{f_0}{f_0^{\mathrm{trial} }}$$

(2)

where $f_0^{\mathrm{trial} }$ and $f_0$ are the test strength and the initial strength, respectively;$f_s^{\mathrm{trial} }$ is the test safety factor.

When the number of failure elements reaches the maximum at a step, it means that the model loses stability, and the corresponding test safety factor $f_s^{\mathrm{trial}}$ is the safety factor $F_S$ of the model. $F_S$ can be calculated using Equation (3):

$$F_S=\frac{1}{1-\left(k-1\right)\times \mathsf{\Delta }}$$

(3)

where k is the number of calculation steps until model failure; $\mathsf{\Delta }$ is the reduction coefficient of material strength, which should be set before simulation. When the instability of model occurs, $f_0^{\mathrm{trial} }$=$f_0-\left(k-1\right)\times \mathsf{\Delta }\times f_0$.

According to the deformation and damage state, mesoscopic elements in RFPA can be classified into three phases, i.e., the matrix phase, air phase and contact phase. The RFPA code is able to model the crack initiation, propagation and coalescence through the transformation of the three element phases. The Mohr-Coulomb criterion, with a tensile-off [33], is used as the strength criterion, as shown in Equation (4). Clearly, the tension failure mode is judged by the maximum tensile stress criterion and the shear failure mode is judged by the Mohr-Coulomb criterion.

$$\{\begin{array}{c}{\sigma }_1-\frac{\left(1+\sin \theta \right)}{\left(1-\sin \theta \right)}{\sigma }_3\ge {\sigma }_c,{\sigma }_1\ge {\sigma }_c\left(1-\frac{1+\sin \theta }{1-\sin \theta }\cdot \frac{1}{\lambda }\right)\\ {\sigma }_3\le -{\sigma }_t,{\sigma }_1\le {\sigma }_c\left(1-\frac{1+\sin \theta }{1-\sin \theta }\cdot \frac{1}{\lambda }\right)\end{array}$$

(4)

where ${\sigma }_c$ is the uniaxial-compressive strength, ${\sigma }_t$ is the tensile strength, and $\theta$ is the internal friction angle; $\lambda$ is the ratio of the compressive strength to tensile strength; ${\sigma }_1$ and ${\sigma }_3$ are the maximum and minimum principal stresses, respectively.

A matrix element represents the solid medium. Under uniaxial stress state, its mechanical behavior can be described by the constitutive relation of rock material shown in Figure 2. When the deformation of the element exceeds the ultimate tensile strain ${\epsilon }_{\max \_t}$ under uniaxial tensile state, the matrix element will be converted into an air element and its elastic modulus will decrease to a very low value so that it does not transfer stress any more. Simultaneously, when the deformation of the element exceeds the ultimate compressive strain ${\epsilon }_{\max \_c}$ under uniaxial compression state, it will enter the extrusion phase and be transformed into a contact element with a gradually increasing stiffness. Considering that most of the failure of rock materials is in the tensile mode at the mesoscopic scale, the priority is given to the maximum tensile stress criterion. Namely, it will be firstly judged if the stress state of an element satisfies the maximum tensile stress criterion. If not, the Mohr-Coulomb criterion will be checked next.

The strength reduction method is introduced into the RFPA method to realize the progressive failure progress simulation of the tunnel models. The uniaxial compressive strength and tensile strength of the mesoscopic elements will be reduced linearly in a certain proportion during calculation [36]. The stress, strain and failure states of mesoscopic elements will be analyzed after each strength reduction until the external and internal forces reach new equilibrium. Meanwhile, considering rock medium is a kind of quasi-brittle material and the strain energy will be released in the form of acoustic emission when new cracks occur, and the acoustic emissions caused by element failure will be monitored. Actually, when the stress state of an element reaches the failure criteria, the failure of the element will occur, and then the corresponding acoustic emission will be considered to happen. Namely, if one element gets damaged, an acoustic emission event will be generated. The location of the acoustic emission is the center of the element, and the energy magnitude of acoustic emission is calculated according to the stress–strain curve at the moment of failure. Therefore, it is reasonable to evaluate the stability of the model using the released count and energy of acoustic emissions. The reduction step with the largest number of acoustic emissions corresponds to the macro destruction of the numerical models [33,35]. Hence, it will be regarded as a criterion for judging model instability in this study, and the safety factor of the models can be therefore determined for comparing the stability of the tunnels with different joint sets.

2.2. Numerical Model Setting Up

In this paper, the feasibility and correctness of the combination of the strength reduction method and RFPA in studying the mechanical mechanisms of zonal disintegration are verified by comparing the simulated results with the experimental results of Gao et al. (2018) [37]. The prototype of the indoor test is a 910 m deep roadway tunnel in Dingji Coal Mine, Huainan, China. The similar materials are used to prefabricate a 0.6 m × 0.6 m × 0.6 m cube model and a circular tunnel is excavated at the center with a diameter of 0.1 m. The physical model is compressed along the tunnel axis with a load as high as 2.0 times of the material uniaxial tensile strength. Simultaneously, 1.5 times of the in situ stress load is applied in the horizontal direction, and 1.2 times of the material uniaxial compressive strength is applied in the vertical direction. A series of three-dimensional heterogeneous models are built up, whose material properties including elastic modulus, uniaxial compressive strength and tensile strength are assumed to be subject to the Weibull distribution. At the same time, according to Cai and He (2013) [38], the numerical models with the size of 8700 mm × 8700 mm × 4200 mm in the X, Y and Z directions, respectively, are established, and a circular tunnel with a radius of 700 mm is excavated at the center. The models whose axial direction is parallel to the Z axis are discretized into 2,500,000 hexahedron elements, as shown in Figure 3a. Besides, to reduce the end effect caused by the stiffness mismatch between the loading end and the model specimen [39], a backing plate with a thickness of 400 mm is set on every loading face of the tunnel models. Note that the elastic modulus of the backing plates is equal to the tunnel models but is considered to be homogeneous. Meanwhile, their strength is much greater than the tunnel models, as shown in the red parts of Figure 3b. The physical and mechanical parameters of surrounding rocks, mainly composed of hard rock granite, are listed in

Table 1. Physical and mechanical parameters of model material.

Parameter	Value
Elastic modulus (E)/GPa	80
Compressive strength $\left({\sigma }_c\right)$/MPa	150
Poisson ratio $\left(\mu \right)$	0.25
Friction angle $\left(\varphi \right)$/°	25
C/T coefficient	10
Residual strength coefficient	0.1
Heterogeneity coefficient	4
Reduction coefficient for strength	0.01

.

For the boundary conditions, the surfaces of x = 0, y = 0 and z = 0 are fixed along the normal direction and uniformly distributed loads are applied on the surfaces of x = 9.5 m, y = 9.5 m and z = 5.0 m, respectively. According to the previous research [13], the stress load along the axial direction of the tunnel is set as 2 times the confining pressure. Clearly, the confining pressure of 5 MPa is loaded on the surfaces of x = 9.5 m and y = 9.5 m, and the 10 MPa stress is loaded on the surface of z = 5.0 m in five steps before strength reduction, as shown in Figure 3c.

3. Results

3.1. Comparison of Numerical Simulation and Laboratory Test

In order to verify the correctness of the developed model in studying the zonal disintegration phenomenon, the comparison of the indoor physical test [37] and numerical simulation results of Model A is shown in Figure 4. For examining the numerical simulation results accurately, the different colors represent the different damage states of the elements in Figure 4a. Clearly, the elastic modulus of the elements corresponding to the red color is the smallest, which means the elastic modulus is reduced to the minimum value after those elements get fully failed; meanwhile, the elastic modulus of green elements is about 80 GPa, which has little change compared with the initial value. Namely, no failure and phase transformation occur to these elements.

The existence of scattered yellow and blue elements in the surrounding rock is caused by the different degrees of damage and the discreteness of the strength parameters of the elements in the heterogeneity model. Under the increasing external load, the stress and damage states of each element are closely related to its own properties, such as strength, elastic modulus, etc. Those elements with low phase transition threshold are easily damaged under a low load, such as yellow elements, but they may not grow and form cracks because of the unconnected distribution. However, the elements with high phase transition threshold, such as blue elements, may not get damaged even under high stress. From Figure 4a, it can be observed that the 3D RFPA method combined with the strength reduction method can successfully reproduce the zonal disintegration phenomenon characterized by the interval distribution of fractured zones and intact zones in surrounding rock. Furthermore, it can be found that the fractured ring farther away from the tunnel wall has thinner failure thickness, lower integrity and fuzzier boundary.

The elastic modulus evolution process of zonal disintegration obtained by RFPA numerical simulation is shown in Figure 5. It can be observed from Figure 5a that when the strength reduction is performed to the 94th step, the first fracture area appears near the tunnel wall, and the peak value of maximum principal stress ${\sigma }_1$ of surrounding rock is about 15 MPa, as shown in Figure 6. The failure area can be regarded as the traditional loose zone produced by the excavation of the tunnel. After that, when the model is calculated to the 98-2nd reduction step, the stress redistributes, the mean value of ${\sigma }_1$ rises up to 19.5 MPa, and the second complete circular fracture zone of surrounding rock is generated in an interval way, as shown in Figure 5b. Similarly, with the further reduction of the surrounding rock strength, the stress transfers to the deeper part of the surrounding rock again, and the third circular fracture area appears. The maximum value of ${\sigma }_1$ reaches 24.5 MPa now. In the same way, there are six fracture zones appearing in the surrounding rock of the model. Finally, when the sixth fracture zone located 2 m away from the tunnel wall is formed, the peak stress decreases to 27 MPa, and the stress caused by the excavation of the tunnel is basically released, which cannot cause more damage.

3.2. Numerical Simulation of Zonal Disintegration Phenomenon in Jointed Rock Mass with Different Inclinations

To discuss the zonal disintegration behavior of circular tunnel in jointed rock mass with different dip angles, the related models are established as shown in Figure 7, including three inclinations of 0°, 30° and 45°, respectively, i.e., Models B0, B30 and B45. The physical and mechanical parameters of these joints, composed of soft sandstone and other materials, are shown in

Table 2. Physical and mechanical parameters of joint material.

Parameter	Value
Elasticity modulus (E)/GPa	16
Compressive strength $\left({\sigma }_c\right)$/MPa	30
Poisson ratio $\left(\mu \right)$	0.25
Friction angle $\left(\varphi \right)$/°	25
C/T coefficient	10
Residual strength coefficient	0.1
Heterogeneity coefficient	200

. The size and boundary conditions of the models are the same as the jointless model A.

The elastic modulus diagrams shown in Figure 8, Figure 9 and Figure 10 depict the reduction steps corresponding to the formation of each fracture ring in Models B0, B30 and B45, through which we can intuitively understand the whole evolution process of the zonal disintegration phenomenon of jointed rock mass with different dip angles, from crack initiation, propagation to final coalescence.

Figure 8 shows the continuous damage evolution of the horizontal joint model B0 under the triaxial stress state, from which can be observed that with the rock mass strength weakening, the first fractured zone adjacent to the tunnel wall appears at the 96th reduction step. The elements in the corresponding region are damaged due to tensile failure and their elastic moduli gradually decrease with the developing damage. When the compressive/tensile strength continues reducing to the 97-7th step, the rock mass at a certain distance away from the first fracture ring forms an intermittent second fracture ring. Then, at the 97-15th step, a third fracture ring is observed at the interval of thicker intact rock mass. However, it is obvious that the new fracture ring owns a lower degree of integrity comparing with the two previous fracture rings. The rupture area of the fracture ring located above the tunnel develops along the horizontal direction due to the existence of the horizontal soft joint. After that, when the strength of the model drops to the 97-19th step, the fourth fracture ring appears. However, this fracture ring shows a strong intermittent feature and develops along the horizontal joint under the bottom of the tunnel. At the 97-21st and 97-23rd steps, the last two fracture rings with fuzzy boundary appears, and some of the damaged elements are even disorderly connected with previous one.

As shown in Figure 9 and Figure 10, compared with the horizontal joint circular tunnel Model B0, the zonal fracture modes of Models B0, B30 and B45 show similarities to a certain extent. For example, the number of fracture rings is the same, and the fracture rings are not as complete as the jointless Model A with a lower degree of integrity and fuzzier boundaries. Meanwhile, the differences between them are significant.

Clearly, the first fracture ring of Model B0 appears at the 96th reduction step, and the last fracture ring forms at the 97-23rd step. Under the same boundary condition, the first fracture ring of Model B30 advances to the 94-0th step, and the last fracture ring appearance postpones to the 98-31st step. As for Model B45, its first fracture ring appears earliest, being advanced to the 92nd step. However, its last failure ring is not complete until the 98-33rd step. Consequently, the simulated results demonstrate that the larger the angle between the joints and the horizontal direction is, the earlier the fracture ring appears. However, the last fracture ring displays the opposite formation rule, i.e., a model with larger joint dip angle needs more reduction steps to form the final failure mode, meaning several fracture rings distribute alternately, which also means that the corresponding surrounding rock is more stable.

Figure 11 displays the contrast diagrams of the accumulated acoustic emission energy and accumulated acoustic emission counts of each reduction step during the deformation and failure process of the different jointed Models B0, B30 and B45, respectively. The model with larger angle between the joints and the horizontal plane generally releases higher acoustic emission energy and more accumulated acoustic emission counts, which means a larger destruction area and more serious damage occurs before the final collapse. These results indicate that the model with larger inclination between the joints and the horizontal plane has higher ultimate load-bearing capacity and therefore higher stability, but the intensity of failure is higher before the final collapse. This conclusion also agrees with the above analysis of Figure 8, Figure 9 and Figure 10 from the perspective of energy.

In order to further analyze the zonal disintegration law of the jointed circular tunnels and reveal the influence of joints with different dip angles on the zonal disintegration phenomenon, the maximum principal stresses ${\sigma }_1$ of the elements on the x = 4750 mm and z = 850 mm line on both sides of the tunnels related to the occurred fracture rings are studied comprehensively. The detailed maximum principal stress ${\sigma }_1$ distribution along the broken line when the fracture rings occur in each model are shown in Figure 12. It can be observed from Figure 12a that the peak values of ${\sigma }_1$ in Model B0 are about 17~30 MPa. As illustrated by the black broken line with rectangular dots, when the first fracture zone is formed, the peak value is about 17.5 MPa. However, this value rises to 21 MPa at the 97-7th step when the second fracture ring is formed, and the spatial position of the peak value moves to the inner rock mass along the radial direction. The subsequent peak values related to each new fracture ring are 23.75 MPa, 25 MPa, 28 MPa and 28.75 MPa, respectively, and the peak positions of each curve are distributed outward, in turn. With the continuous process of stress buildup, stress shadow and stress transfer, the maximum principal stresses of the surrounding rocks gradually become stable at a relatively low level in the far-field rock masses, except for some slight fluctuations caused by the inhomogeneity of rock mass. The corresponding far-field rock masses are basically undamaged under the low excavation disturbance.

Obviously, the failure features of the zonal disintegration phenomenon show significant differences from the traditional excavation loose zone around shallow tunnels. The above analysis of the changes of surrounding rock stresses during the formation of each fracture ring at the key reduction steps suggests that the continuous process of stress buildup, stress shadow and stress transfer plays a critical role in promoting the surrounding rock masses of the circular tunnel to fail in the pattern of alternate regional destructions. Namely, the mechanical process of stress transferring and re-concentrating to the deeper part of the surrounding rock masses keeps the integrity of the rocks between two fracture rings and lead to zonal disintegration phenomenon in the view of mechanical mechanism. Moreover, when comparing the stress distributions of Models B0, B30 and B45 along the broken line, it can be concluded that the stress values show a downward tendency with the increasing of dip angle. In addition, the radius of the outermost fracture ring expands out gradually, i.e., 1.75 m, 1.9 m and 2.1 m away from the tunnel walls, respectively, which means that the destroyed range of surrounding rock masses extend wider if the joint dip angle increases.

3.3. Zonal Disintegration Evolution of Random Joint Model

The circular tunnel Model C with random joints distributed in surrounding rocks is shown in Figure 13. The physical and mechanical parameters of these joints are listed in Table 2 and the boundary conditions and strength reduction process are the same as Model A. The random joint model generates zonal disintegration phenomenon as shown in Figure 14, from which it can be observed that the basic fracture characteristics of Model C are similar to the previous models, i.e., there are also several fracture rings distributed at intervals near the tunnel wall. However, the fracture mode of the random joint model is significantly affected by weak joints and shows specific characteristics.

As shown in Figure 14, the first fracture ring appears at the 92nd strength reduction step, and then the second fracture ring is formed when the strength of the rock materials is reduced to the 96th step. Note that the rupture area D are not damaged along the radial direction as predicted, but forms a linear crack by connecting the upper end of the joint A and the lower end of the joint B. Similarly, when the rock strength is reduced to the 97-2nd step, the lower end of the joint B connects with the lower end of the joint C, forming a straight crack. Then, as shown in Figure 14d, the third fracture ring with low integrity degree is formed at the 97-8th step. After that, two new fracture rings with fuzzy boundaries and discontinuous shape occur at the 97-16th and 97-18th steps, respectively, as shown in Figure 14e,f. Furthermore, it can be observed from Figure 14 that during the failure process of the random joint Model C, the linear cracks formed by connecting close joints have priority over the damage path along the radial direction because of the stress concentration at joint ends and dense joint distribution. It means that the zonal disintegration phenomenon affected by random joints in surrounding rock masses has a low possibility to form the closed complete fracture rings because the existed joints may break the integrity of the potential fracture ring.

As shown in Figure 15, the maximum principal stresses of the randomly jointed rock around the circular tunnel are approximately 15~20 MPa, which are lower than the previous models. Additionally, the last fracture ring is located only 1.2 m away from the tunnel wall, which indicates the overall damage area of the surrounding rock is obviously smaller than the models discussed above.

As shown in Figure 16, as many as 20,488 acoustic emission events occur before macro instability, and there is a significant upward trend before fully researching the overall failure of the tunnel, demonstrating a salient failure precursory characteristic of the inhomogeneous surrounding rock. Besides, around 1.58 × 10⁷ J, acoustic emission energy is released by the damaged elements, which is also much lower than the above jointed rock models and jointless model because of its smaller overall damage area.

4. Discussion

Table 3. The safety factor of the circular tunnels with different surrounding rock masses.

Model	Strength Reduction Step at Failure	Critical Compressive/Tensile Strength (MPa)	Safety Factor
A	100	7.5/0.75	20.00
B0	96	13.5/1.35	11.11
B30	97	12/1.20	12.50
B45	98	10.5/1.05	14.29
C	96	13.5/1.35	11.11

lists the values of strength reduction times, critical compressive and tensile strengths and tunnel safety factors $F_S$ of the studied models obtained by the strength reduction method. It can be observed that the $F_S$ of Model A, i.e., the jointless model, is larger than the other models mainly because the jointless surrounding rock of the tunnel is the most stable, and the relatively closed complete fracture rings are produced in this model, showing the distinct zonal disintegration phenomenon. Moreover, after comparing the accumulated acoustic emission energy and the $F_S$ values of the jointed Models B0, B30 and B45, it can be concluded that with the dip angle between the joint set and the horizontal direction increasing, the surrounding rock of the tunnel will become more and more stable under the boundary conditions and material parameters in this study.

Furthermore, there is a positive correlation between the accumulated acoustic emission energy before final collapse and the safety factor of the tunnel, as shown in Figure 17. It is worth noting that the accumulated acoustic emission energy of each model compressed by initial stresses is very small, and the accumulated acoustic emission counts are basically negligible, as shown in Figure 11 and Figure 16. Therefore, it can be proven that the rock masses surrounding the tunnel are not damaged for the adopted parameters without reducing the strength. Besides, the safety factor of the random joint Model B0 is coincidentally equal to the horizontal joint Model C. However, the horizontal joints and randomly distributed joints play a completely different role in forming the zonal disintegration in terms of background stress evolution; i.e., the gradual process of stress buildup, stress shadow and stress transfer induced by the existed joints would result in different fracture patterns. Further comparison of the accumulated acoustic emission energy released in the two models demonstrates that at the same reduction steps, the latter release more energy than the former because of its highly fractured surrounding rock masses. Reasonably, Model B0 has a better stability than Model C under the influence of the distributed joints. Hence, the acoustic emission events provide more accurate evaluation information rather than the tunnel safety factor.

According to the above analysis, we know that in deep geological environments, the fracture characteristics of the zonal disintegration of surrounding rocks are very different from the traditional loose zone of surrounding rocks. Therefore, the support measures adopted for the two phenomena should also be different. Firstly, in terms of the failure mode, the traditional loose circle only produces fractures within a certain range away from the tunnel wall, while the zonal disintegration produces spaced failure areas within the interior of surrounding rocks, which leads to the different specific positions needing to be supported. Secondly, the primary process of stress redistribution inducing the loose zone of surrounding rocks can be finished immediately after tunnel excavation. However, the surrounding rocks will go through several times of stress redistribution when forming zonal disintegration. Hence, compared with the traditional support measures, the occurrence sequence of fracture rings should be considered for handling zonal disintegration, and the supporting schemes have to be determined. Finally, considering that the zonal fractures of surrounding rocks generally happen in high-stress environments in deep engineering, which may result in unimaginable accidents, it is necessary to further study more effective support methods.

5. Conclusions

With the aim of understanding the zonal disintegration phenomenon and revealing the intrinsic mechanical mechanisms, a series of 3D heterogeneous tunnel models considering varying joint dip angles are built up and tested by the rock failure process analysis method on the basis of the continuum mechanics, meso-damage mechanics and statistical strength theory. The strength reduction method is embedded to achieve the gradual fracture process, final failure mode and safety factor and to reproduce the characteristic fracture phenomenon of deep tunnels under high geo-stress level, i.e., zonal disintegration. The following conclusions can be reached:

The combined approach has been proven effective in capturing the zonal disintegration characteristics by the indoor physical test [37]. The detailed failure process of the jointed rock models demonstrates that the zonal disintegration is induced by the stress redistribution of surrounding rock masses. Namely, the continuous process of stress buildup, stress shadow and stress transfer play a critical role in promoting the surrounding rock masses to fail in the pattern of alternate regional destructions. It is the mechanical process of stress transferring and re-concentrating to the deeper part of the surrounding rock masses to keep the integrity of the rocks between two fracture rings and lead to zonal disintegration phenomenon in the view of mechanical mechanism, which is obviously different from the formation mechanism of the traditional surrounding rock loose zone.

The dip angle of the existed joint set shows great influence on the stress evolution and inner crack propagation of the surrounding rock mass. On the one hand, the fracture ring of the 45° joint model first appears at the 92-0th step, but the last ring forms at the 98-33rd step. Actually, a model with larger joint dip angle needs more reduction steps to generate the final failure mode. On the other hand, the acoustic emission energy released by Model B45 is the most, up to 4.77 × 10⁸ J, while Model B0 only releases 1.58 × 10⁷ J energy. This means that the model with larger inclination angle will be damaged more seriously before the final collapse. Besides, as the joint inclination increases from 0° to 45°, the safety factor of the model rises from 11.11 to 14.29. In addition, the model with larger joint inclination has higher ultimate load-bearing capacity and therefore higher stability. However, the intensity of failure is also higher before the final collapse.

Furthermore, the existence of random joints might significantly affect the integrity and regularity of the zonal disintegration pattern because the high stress concentrations at the ends of the existed joints could change the crack growth path and prevent the damage development along the radial direction. Namely, the zonal disintegration affected by randomly distributed joints has a low possibility to form the closed complete fracture rings. Meanwhile, the upward trend of acoustic emission energy released before nearly researching the overall failure of the tunnel could be regarded as a macro-failure precursor. Models C and B0 have the same safety factor of 11.11. The acoustic emission count of the former is 20,488, which is much higher than 12,801 of the latter. Hence, the acoustic emission events are able to provide more accurate evaluation basis rather than the tunnel safety factor.

These achievements will improve our understanding of the mechanical mechanism of zonal disintegration formation, provide insights into the characteristic fracture phenomenon of deep-buried tunnels, and establish the basis for design, construction and treatment of tunnels in deep engineering.