# Numerical Study on Zonal Disintegration of Deep Rock Mass Using Three-Dimensional Bonded Block Model

Authors: Qingteng Tang; Wenbing Xie; Xingkai Wang; Zhili Su; Jinhai Xu
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

Zonal disintegration, a phenomenon of fractured zones and intact zones distributed alternately in deep rock mass, is different from the excavation-damaged zone of shallow rock mass. In this study, bonded block model of 3DEC was employed to study the fracture mode and origination condition of zonal disintegration. Initiation, propagation, and coalescence progress of fracture around the roadway boundary under different triaxial stress conditions are elaborated. Numerical simulation demonstrated that zonal disintegration may occur when the direction of maximum principal stress is parallel to the roadway axis. It is interesting to find that the fracture around the roadway boundary traced the line of a spiral line, while slip-line fractures distributed apart from the roadway boundary. The extent of the alternate fracture zone decreased as the confining pressure increased, and alternate fracture zone was no longer in existence when the confining pressure reaches a certain value. Effects of roadway shape on zonal disintegration were also studied, and the results indicated that the curvature of the fracture track line tends to be equal to the roadway boundary in shallow surrounding rock of the roadway, while the fractures in deep surrounding rock seems unaffected by the roadway shape. Those findings are of great significance to support design of deep underground openings.


---

## Body

## 1. Introduction

With the increasing depth of underground engineering, fracture zones distribution around the underground opening of deep rock mass is much different from the excavation damaged zone of shallow rock mass. Fissured zone, plastic zone, and elastic zone formed continuously in the shallow rock mass according to classic elastic-plastic theory [1]. While the zonal disintegration proposed by Shemyakin et al. [2] refers to the fracture zone and comparatively intact zone formed alternately in the deep rock mass (Figure 1(b)). It is significant to adjust the support design and excavation method for underground opening due to the emergence of the zonal disintegration phenomenon in deep rock mass.Figure 1
Fracture distribution around the tunnel: (a) based on elastic-plastic theory [1]; (b) field monitoring result [2].
(a)
(b)A radical debate on the existence of zonal disintegration was performed at Beijing [3] due to the great difference between zonal disintegration and classic fracture mode of excavation damage zone, and many studies have been carried out to validate the zonal disintegration by field monitoring. Zonal disintegration was initially observed in Witwatersrand gold mines which are located in South Africa [4]. Subsequently, zonal disintegration was monitored in several mines of the Norilsk minefield by using underground electric soundings and borehole periscope [2]. Furthermore, zonal disintegration was also observed in several deep coal mines of China [5, 6], which indicated that zonal disintegration was mainly caused by surrounding rock separating along with the weak structure. Similar phenomenon was captured at the Jinping underground laboratory [7].It is difficult to completely understand the generation and propagation process of fracture in surrounding rock mass around the underground opening owing to the fact that the monitoring result can only reflect the local failure of surrounding rock mass. So, scholars are dedicated to study the patterns and origination conditions of zonal disintegration by using a variety of research methods.Laboratory tests are seemingly the most direct method to validate the existence and fracture mode of zonal disintegration. Zonal disintegration was reproduced by loading along with the opening axis on a cylindrical similar-material-model with a circular opening which is confined at the sides and bottom [8]. However, the fractures are close to the excavation boundary (Figure 2(a)) because the loading process was executed after the excavation which is different from the in situ loading process that surrounding rock mass suffered. Subsequently, the laboratory test for zonal disintegration using triaxial loading equipment was carried out by Yuan et al. [9], a horse-hole cavern was excavated after in situ stress (σh:σH:σv = 1 : 3:1, the maximum horizontal stress is parallel to the cavern axis) applied. More obvious zonal disintegration phenomenon was obtained (Figure 2(b)). It is interesting to note that a spiral-line fracture could be captured around the cavern boundary. The laboratory test demonstrated that the high axial stress parallel to the opening axis is a key factor of zonal disintegration. However, there is still a lot of argument about the fracture mode of zonal disintegration.Figure 2
Laboratory test results of zonal disintegration: (a) load after opening excavation [8]; (b) excavation after in situ stress applied [9].
(a)
(b)Various numerical approaches were also adopted to carry out a comprehensive study of zonal disintegration due to the high cost and limits of the laboratory test. Qian et al. [10] adopted the 2D weak-element method which demonstrated that slip-line zonal fracture occurred when the in situ stress is greater than uniaxial compressive strength of surrounding rock mass, while annular zonal fractures were proved by Zuo et al. [11] using Realistic Failure Process Analysis Parallel system model (RFPA3D) of which axial stress parallel to the opening axis could be taken into account. Subsequently, zonal disintegration was deeply studied by using RFPA3D [12–15]. These studies show that zonal disintegration is a common phenomenon formed by shear failure zone under high stress, while the annular alternate fracture zone was generated under special stress condition that high horizontal stress is parallel to the tunnel axis. Meanwhile, the similarly result was obtained by using FLAC3D. Wang et al. [16] adopted the strain softening heterogeneous model with Mohr–Coulomb criterion with the tension cut-off failure criterion to study the zonal disintegration from the aspect of spatial strain localization and annular high and low shear strain zones alternately distributed around surrounding rock mass was obtained. Similar number and location of angular fracture zones were obtained by comparing the result of in situ monitoring by employing the maximum tensile stress criterion and strain energy density theory [17]. Subsequently, Zonal disintegration was also reproduced by using ABAQUS [18, 19]. Zonal disintegration and strain fluctuant law were discovered by using bilinear strain softening model and energy damage failure criterion. However, some scholars believed that rock mass around the excavation was subjected to dynamic loading during the excavation process. Annular zonal fractures were obtained by using the General Particle Dynamics (GPD) method [20, 21].Meanwhile, various theories have been proposed to explain the mechanism of zonal disintegration. The fluctuation law of radial and tangential stress in surrounding rock mass after underground opening excavated under hydrostatic pressure stress was deduced by the non-Euclidean continuum model [22]. Subsequently, the solution of non-Euclidean continuum model under nonhydrostatic pressure and the effect of cracks on non-Euclidean continuum model were studied [23]. Chen et al. [24] deduced that appearance of zonal disintegration is determined by the magnitude of in situ stress and mechanics properties of rock mass based on fracture mechanics. The density of fractures increases as in situ stress increases or cohesion and friction angle decreases.In a word, origination conditions and fracture mode of zonal disintegration have been comprehensively studied. However, radical debates on fracture mode of zonal disintegration have never stopped.Zonal disintegration may occur after the rock mass strength peaked. Thus, the pattern of zonal disintegration is determined by rock mass postpeak behaviors. However, most of the proposed numerical study was based on continuum theory, and this theory does not fully explain the spalling and slip failure of rock mass. Distinct element method was widely applied to study the discontinuous deformation of rock mass in recent years. Although a study of zonal disintegration using 2D distinct element method has also been reported [25], effect of the stress parallel to the excavation axis on zonal disintegration was not taken into account. As to theoretical approaches, existing theoretical approaches are almost based on the plane strain model, and stress parallel to the underground opening axis was not taken into account of which the significance has been proved by laboratory test. Furthermore, theories are mostly based on a variety of assumptions. Thus, there are numerous limitations to apply proposed theories to complex engineering site.This study employed the 3DEC bonded block model to study the fracture mode and origination conditions of zonal disintegration. Microcontact properties of the numerical rock sample were calibrated by the uniaxial compression test and Brazilian test so that the mechanical behavior of the numerical sample was consistent with the sandstone obtained from Zhujixi coal mine. Subsequently, a cuboid model with a circular roadway was constructed to study the zonal disintegration based on the mechanics properties of the numerical rock sample. Fracture mode of zonal disintegration around the roadway boundary under different triaxial stress conditions was investigated. Moreover, effect of roadway shape on zonal disintegration was also determined.


## 2. Microcontact Properties Calibration of Numerical Rock Sample

### 2.1. Bonded Block Model of Three-Dimensional Distinct Element Method

3DEC bonded block model is consisted of an assembly of tetrahedral blocks bonded at their contacts [26]. The bonded blocks can separate and slip along the contacts in microlevel so that the rock mass tensile and shear fracture can be simulated. Mechanical behavior of 3DEC Bonded block model is determined by the microcontact properties between tetrahedral blocks, and the contact behavior satisfies the Coulomb friction law.Mechanical interaction between two contacting blocks has a significant influence on the mechanical behavior of a DEM simulation in 3DEC. The equation between force and displacement of blocks was interpreted as follows [27].Force of contacts updates in a time step with the contribution of the displacement change. As Figure3 illustrates, in the normal direction of the contact plane, Fnmax is the maximum normal force and Tmax is the tensile strength of contact. If(1)Fn≤Fnmax,then(2)ΔFn=−Kn⋅ΔUn,where Kn is the normal stiffness of contact and ΔUn is the variation of normal displacement.Figure 3
Force equation ofith step and (i + 1)th step of contact.In the tangential direction of contact plane, the maximum tangential force of contact could be derived from the Coulomb friction law:(3)Fsmax=c+Fntanφ,where c is the cohesion of contact and φ is the friction angle of contact. Similarly, if(4)Fs≤Fsmax,then(5)ΔFs=−Ks⋅ΔUs,where Ks is the shear stiffness of contact and ΔUs is the variation of shear displacement.Else, if either normal force of contact is greater than the maximum normal force or tangential force of contact is greater than the maximum tangential force, then the tensile strength and cohesion of contact will decrease to zero.(6)Tmax=0,(7)Fsmax=Fntanφ.The effect of microcontact properties on the mechanical behavior of macronumerical model has been studied [28], which indicated that Young’s Modulus and Poisson’s ratio of 3DEC bonded block model are determined by normal stiffness and ratio of tangential stiffness to normal stiffness of microcontacts. Meanwhile, tensile strength, cohesion, and internal friction angle are calibrated by the contact tensile strength, cohesion, and internal friction angle, respectively.


### 2.2. Uniaxial Compression Test and Brazilian Test

A cylindrical rock sample for the uniaxial compression test (Figure4(a)) and disc sample for the Brazilian test (Figure 4(b)) were created in 3DEC using bonded block model. Young’s modulus Poisson’s ratio and uniaxial compressive strength were calibrated by uniaxial compression test, the tensile strength was calibrated by Brazilian test according to the ISRM suggested method [29].Figure 4
3DEC bonded block geometry: (a) uniaxial compression test; (b) Brazilian test.
(a)
(b)Geometry of the 3DEC Uniaxial compression test model is shown in Figure4. The uniaxial compression test sample is composed of 36655 tetrahedral blocks which have a width of 0.6 m. The sample has a height of 16 m and diameter of 8 m with two rigidity equivalent platens. The 0.6 m width tetrahedral block was employed in consideration of computer performance limitations and macro qualitative analysis was conducted in this study. Effect of block size on the numerical study of zonal disintegration was studied in the subsequent section. In the same way, the Brazilian test sample includes a 16 m diameter 4 m thick disc which consists of 27078 tetrahedral blocks. The mechanics properties of the tetrahedral blocks and microcontact properties assigned to the 3DEC bonded block model are listed in Tables 1 and 2, respectively. The details of calibration procedure for microcontact properties could be found in Kazerani and Zhao [28].Table 1
Mechanics properties of the intact rock.
Material type
Density (kN/m3)
Young’s modulus (GPa)
Poisson’s ratio
Sandstone
2500
3.55
0.26Table 2
Calibrated microcontact properties assigned to 3DEC bonded block model.
Normal stiffness (GPa/m)
Shear stiffness (GPa/m)
Contact cohesion (MPa)
Contact friction angle (deg)
Contact tensile strength (MPa)
Residual friction angle (deg)
Residual cohesion (MPa)
Residual tensile strength (MPa)
17
6.8
23
18
9.2
13
1.3
0.52In the compression tests, loading was applied to the sample by applying a constant velocity of 0.1 m/s to the top platen while the bottom platen was fixed. Because one calculation step consumes about 10−7 s in 3DEC, the loading rate is slow enough so that the effect of loading rate on numerical test could be ignored [30]. Stress of the top platen and axial displacement of the rock sample were recorded during the loading process. Likewise, the same velocity in the vertical direction was applied at the top platen in the Brazilian test. The tensile strength of the rock sample was derived from the following formula:(8)σt=2PmaxπDt,where Pmax is the maximum load, D is the diameter of the rock sample, and t is the thickness of the rock sample. The simulation results of uniaxial compression test and Brazilian test are illustrated in Figures 5 and 6.Figure 5
Numerical test result of uniaxial compression test. (a) Frontal views at different view angles of fracture mode, red dash line for shear fracture and yellow dash line for tensile fracture. (b) Stress-strain curve of the numerical test and the laboratory test of sandstone.
(a)
(b)Figure 6
Numerical test result of Brazilian test: (a) fracture mode; (b) stress-strain curve of numerical test.
(a)
(b)Figure5(a) illustrates that tensile fractures and shear fractures form due to microcracks propagation and coalescence in the postpeak stage. Two main shear fractures represent macroshear fracture which can be distinctly observed generate at the right corner and middle part of the sample. Several vertical tensile fractures sporadically distribute at the surface of sample and connect to the shear fractures. As Figure 5(b) shows, stress-strain curve of the numerical rock sample under uniaxial compression is basically consistent with the laboratory test result of the sandstone obtained from Zhujixi coal mine. The uniaxial compressive strength of the numerical rock sample is 45 MPa. Stress radically decreases after uniaxial compression strength peaked and remains at a low level. The sample shows elastic brittle fracture pattern and postpeak behavior of the sample belong to the class I behavior proposed by Wawersik and Fairhurst [31].As Figure6 shows, a subvertical tensile fracture forms along the diameter of the disc sample. In addition, minor tensile fractures can also be seen from the vicinity of the platens. Tensile strength of the numerical rock sample is approximately 1.46 MPa.The rock mechanics test carried out using 3DEC bonded block model demonstrated that 3DEC can realistically simulate the complete process of rock deformation and spalling. Thus, zonal disintegration can be validly studied by employing this method.


## 2.1. Bonded Block Model of Three-Dimensional Distinct Element Method

3DEC bonded block model is consisted of an assembly of tetrahedral blocks bonded at their contacts [26]. The bonded blocks can separate and slip along the contacts in microlevel so that the rock mass tensile and shear fracture can be simulated. Mechanical behavior of 3DEC Bonded block model is determined by the microcontact properties between tetrahedral blocks, and the contact behavior satisfies the Coulomb friction law.Mechanical interaction between two contacting blocks has a significant influence on the mechanical behavior of a DEM simulation in 3DEC. The equation between force and displacement of blocks was interpreted as follows [27].Force of contacts updates in a time step with the contribution of the displacement change. As Figure3 illustrates, in the normal direction of the contact plane, Fnmax is the maximum normal force and Tmax is the tensile strength of contact. If(1)Fn≤Fnmax,then(2)ΔFn=−Kn⋅ΔUn,where Kn is the normal stiffness of contact and ΔUn is the variation of normal displacement.Figure 3
Force equation ofith step and (i + 1)th step of contact.In the tangential direction of contact plane, the maximum tangential force of contact could be derived from the Coulomb friction law:(3)Fsmax=c+Fntanφ,where c is the cohesion of contact and φ is the friction angle of contact. Similarly, if(4)Fs≤Fsmax,then(5)ΔFs=−Ks⋅ΔUs,where Ks is the shear stiffness of contact and ΔUs is the variation of shear displacement.Else, if either normal force of contact is greater than the maximum normal force or tangential force of contact is greater than the maximum tangential force, then the tensile strength and cohesion of contact will decrease to zero.(6)Tmax=0,(7)Fsmax=Fntanφ.The effect of microcontact properties on the mechanical behavior of macronumerical model has been studied [28], which indicated that Young’s Modulus and Poisson’s ratio of 3DEC bonded block model are determined by normal stiffness and ratio of tangential stiffness to normal stiffness of microcontacts. Meanwhile, tensile strength, cohesion, and internal friction angle are calibrated by the contact tensile strength, cohesion, and internal friction angle, respectively.


## 2.2. Uniaxial Compression Test and Brazilian Test

A cylindrical rock sample for the uniaxial compression test (Figure4(a)) and disc sample for the Brazilian test (Figure 4(b)) were created in 3DEC using bonded block model. Young’s modulus Poisson’s ratio and uniaxial compressive strength were calibrated by uniaxial compression test, the tensile strength was calibrated by Brazilian test according to the ISRM suggested method [29].Figure 4
3DEC bonded block geometry: (a) uniaxial compression test; (b) Brazilian test.
(a)
(b)Geometry of the 3DEC Uniaxial compression test model is shown in Figure4. The uniaxial compression test sample is composed of 36655 tetrahedral blocks which have a width of 0.6 m. The sample has a height of 16 m and diameter of 8 m with two rigidity equivalent platens. The 0.6 m width tetrahedral block was employed in consideration of computer performance limitations and macro qualitative analysis was conducted in this study. Effect of block size on the numerical study of zonal disintegration was studied in the subsequent section. In the same way, the Brazilian test sample includes a 16 m diameter 4 m thick disc which consists of 27078 tetrahedral blocks. The mechanics properties of the tetrahedral blocks and microcontact properties assigned to the 3DEC bonded block model are listed in Tables 1 and 2, respectively. The details of calibration procedure for microcontact properties could be found in Kazerani and Zhao [28].Table 1
Mechanics properties of the intact rock.
Material type
Density (kN/m3)
Young’s modulus (GPa)
Poisson’s ratio
Sandstone
2500
3.55
0.26Table 2
Calibrated microcontact properties assigned to 3DEC bonded block model.
Normal stiffness (GPa/m)
Shear stiffness (GPa/m)
Contact cohesion (MPa)
Contact friction angle (deg)
Contact tensile strength (MPa)
Residual friction angle (deg)
Residual cohesion (MPa)
Residual tensile strength (MPa)
17
6.8
23
18
9.2
13
1.3
0.52In the compression tests, loading was applied to the sample by applying a constant velocity of 0.1 m/s to the top platen while the bottom platen was fixed. Because one calculation step consumes about 10−7 s in 3DEC, the loading rate is slow enough so that the effect of loading rate on numerical test could be ignored [30]. Stress of the top platen and axial displacement of the rock sample were recorded during the loading process. Likewise, the same velocity in the vertical direction was applied at the top platen in the Brazilian test. The tensile strength of the rock sample was derived from the following formula:(8)σt=2PmaxπDt,where Pmax is the maximum load, D is the diameter of the rock sample, and t is the thickness of the rock sample. The simulation results of uniaxial compression test and Brazilian test are illustrated in Figures 5 and 6.Figure 5
Numerical test result of uniaxial compression test. (a) Frontal views at different view angles of fracture mode, red dash line for shear fracture and yellow dash line for tensile fracture. (b) Stress-strain curve of the numerical test and the laboratory test of sandstone.
(a)
(b)Figure 6
Numerical test result of Brazilian test: (a) fracture mode; (b) stress-strain curve of numerical test.
(a)
(b)Figure5(a) illustrates that tensile fractures and shear fractures form due to microcracks propagation and coalescence in the postpeak stage. Two main shear fractures represent macroshear fracture which can be distinctly observed generate at the right corner and middle part of the sample. Several vertical tensile fractures sporadically distribute at the surface of sample and connect to the shear fractures. As Figure 5(b) shows, stress-strain curve of the numerical rock sample under uniaxial compression is basically consistent with the laboratory test result of the sandstone obtained from Zhujixi coal mine. The uniaxial compressive strength of the numerical rock sample is 45 MPa. Stress radically decreases after uniaxial compression strength peaked and remains at a low level. The sample shows elastic brittle fracture pattern and postpeak behavior of the sample belong to the class I behavior proposed by Wawersik and Fairhurst [31].As Figure6 shows, a subvertical tensile fracture forms along the diameter of the disc sample. In addition, minor tensile fractures can also be seen from the vicinity of the platens. Tensile strength of the numerical rock sample is approximately 1.46 MPa.The rock mechanics test carried out using 3DEC bonded block model demonstrated that 3DEC can realistically simulate the complete process of rock deformation and spalling. Thus, zonal disintegration can be validly studied by employing this method.


## 3. Numerical Simulation of Zonal Disintegration

### 3.1. Numerical Test Scheme of Zonal Disintegration

A cuboid model (Figure7) with a circular roadway which includes an elasticity boundary and a synthetic rock mass was created based on the microcontact and intact rock properties mentioned in Section 2. The synthetic rock mass consisted of 240713 tetrahedral deformable blocks which have a width of 0.6 m which was identical to the width of the sample for rock properties test.Figure 7
Numerical model for zonal disintegration test: (a) cuboid model created in 3DEC; (b) sketch of the numerical model.
(a)
(b)As Figure7 illustrates, the cuboid model has a size of 50 m (width) × 50 m (height) × 4 m (thickness), and the synthetic rock mass has a size of 40 m (width) × 40 m (height) × 4 m (thickness) with a circular excavation which has a radius of 2 m. Theoretically, the larger the rock mass size is, the smaller the calculation error will be. Since stress disturbance range of tunnel excavation is about 5 times the diameter of the underground opening [32], considering the limitation of computer performance, the synthetic rock mass of 40 m (width) 40 m (height) was selected. Numerical simulations with increasing axial stress level were undertaken, which indicated that the influence of model size on the research purpose of this study could be ignored. If the in situ stress was directly applied to the boundary of synthetic rock mass, local stress concentration will occur in the model. Thus, a 5 m wide elastic boundary is constructed surrounding the synthetic rock mass for applying boundary conditions based on the fact that elastic boundary will undergo elastic deformation according to the stress disturbance range of roadway excavation. Frictionless roller boundary condition was assigned to the five faces except the top face of the cuboid model. Stress of 11.5 MPa which has a direction of the negative Z was applied to the top boundary as the weight of about 460-m overburden strata. A circular roadway was excavated after the model was equilibrated under in situ stress state (σx=12MPa,σy=60MPa,σz=12MPa).Previous studies indicated that zonal disintegration may occur when in situ stress is greater than uniaxial compression strength of rock mass [10]. Thus, a 60 MPa-stress of y direction was applied to the zonal disintegration test. And the magnitude of σx and σz equal to the weight of about 500 m overburden strata. The radial displacement, radial stress, and tangential stress were monitored along the monitoring line shown in Figure 7 which has an interval of 1 m with an increase of time step.


### 3.2. Analysis on Numerical Test Result of Zonal Disintegration

Initiation, propagation, and coalescence progress of cracks around the roadway boundary are illustrated in Figure8. It can be found from Figure 8 that shear fracture and tensile fracture form in surrounding rock mass under the high stress parallel to the axial of the roadway. Fissured zones and intact zones distributing alternately was captured. The trace line of fracture presents as a spiral line close to the roadway boundary which is similar to the laboratory result performed by Yuan et al. [9], while slip-line fractures densely distribute in the deep surrounding rock of the roadway.Figure 8
Fracture distribution around the roadway boundary at different time step. (a) 5000 steps. (b) 10000 steps. (c) 18000 steps. (d) 34000 steps.
(a)
(b)
(c)
(d)Crack initially generated at the vicinity of the roadway boundary which can be observed in Figure8(a). According to Figure 8(b), two tensile fractures formed as the crack propagated into deep surrounding rock of the roadway. Subsequently, shear fractures propagate into the deep surrounding rock of the roadway in all directions. At the same time, vertical and horizontal tensile fractures generate at different positions in surrounding rock mass and then connected with the shear fractures. Surrounding rock mass was divided into discontinuous zones by those shear fractures and tensile fractures. Tensile fractures are mainly located at shallow surrounding rock of the roadway, while shear fractures distribute in the deep surrounding rock of the roadway.The radial displacement, radial stress, and tangential stress around the roadway acquired by the monitoring line are illustrated in Figure9. Figure 9(a) indicates that stress redistributed in surrounding rock mass after the roadway was excavated. Stress around the roadway internal surface decreased radically due to the unloading effect of the excavation. The magnitude of radial stress increases gradually to in situ level as the distance to the roadway rib increases.Figure 9
Monitoring results of (a) the stress around the roadway and (b) the radial displacement around the roadway.
(a)
(b)In addition, tangential stress increases by 2-3 times of in situ stress and then decreases to in situ level as the distance to the roadway rib increases. The global variation law of radial stress and tangential stress of surrounding rock mass is consistent with the traditional elasticity-plasticity theory [33], while the variation of radial stress and tangential stress oscillates as the distance to the roadway rib increases which is a significant difference to the traditional theory. The change law of stress around the roadway is similar to the stress distribution of non-Euclidean geometry theory proposed by Guzev and Paroshin [22]. As Figure 9(b) shows, positive value of displacement represents that surrounding rock deformed toward to the internal of the roadway boundary. The displacement around the roadway boundary is about 473 mm. In general, the displacement decreased as the distance to the roadway rib increased. However, on account of fractures generated at different depths of surrounding rock, the rock mass on either side of the fracture was subjected to stress in opposite directions, and the deformation of deep side rock mass around the fracture was restrained. Thus, the displacement of the deeper rock mass could be greater than the deep side rock mass around the fracture. We can observe from Figure 9(b) that displacement fluctuates slightly at the deep surrounding rock of the roadway.


### 3.3. Effect of Block Size on Numerical Test of Zonal Disintegration

Block size is an inherent property of 3DEC bonded block model. The smaller the block size is, the more the model tends to be isotropic and homogeneous. Two numerical models, comprising tetrahedral blocks which have an edge length of 0.8 m and 1 m, were constructed in order to study the effect of block size on zonal disintegration. The same intact rock properties, microcontact properties, boundary condition, and loading procedure mentioned in Section3.1 were applied to the constructed models. Numerical test result which adopted different block size models is shown in Figure 10.Figure 10
Fracture distribution around the roadway boundary using different sizes of tetrahedral block: (a) block with 0.6 m edge length; (b) block with 0.8 m edge length; (c) block with 1 m edge length.
(a)
(b)
(c)As Figure10 shows, fracture modes of different width tetrahedral block models are similar. Tensile fracture, which is approximately parallel to the roadway boundary, forms in all three models. However, it is interesting to find that the extent of fractured zone decreases as the block size increases which may be due to that the confinement of block movement is stronger after contacts separation as the block size increases. The numerical test using different sizes of the block indicated that tetrahedral block which has a width of 0.6 m is small enough to study zonal disintegration.Numerical test of zonal disintegration indicated that the fracture zone and intact zone alternately distributed around the roadway boundary under high axial stress parallel to the roadway axis. It is worthy to note that zonal disintegration occurs with generation and connection of shear cracks and tensile cracks. The trace line of the fracture presented as a spiral line close to the roadway boundary which is different from the circular fracture zone in proposed studies.


## 3.1. Numerical Test Scheme of Zonal Disintegration

A cuboid model (Figure7) with a circular roadway which includes an elasticity boundary and a synthetic rock mass was created based on the microcontact and intact rock properties mentioned in Section 2. The synthetic rock mass consisted of 240713 tetrahedral deformable blocks which have a width of 0.6 m which was identical to the width of the sample for rock properties test.Figure 7
Numerical model for zonal disintegration test: (a) cuboid model created in 3DEC; (b) sketch of the numerical model.
(a)
(b)As Figure7 illustrates, the cuboid model has a size of 50 m (width) × 50 m (height) × 4 m (thickness), and the synthetic rock mass has a size of 40 m (width) × 40 m (height) × 4 m (thickness) with a circular excavation which has a radius of 2 m. Theoretically, the larger the rock mass size is, the smaller the calculation error will be. Since stress disturbance range of tunnel excavation is about 5 times the diameter of the underground opening [32], considering the limitation of computer performance, the synthetic rock mass of 40 m (width) 40 m (height) was selected. Numerical simulations with increasing axial stress level were undertaken, which indicated that the influence of model size on the research purpose of this study could be ignored. If the in situ stress was directly applied to the boundary of synthetic rock mass, local stress concentration will occur in the model. Thus, a 5 m wide elastic boundary is constructed surrounding the synthetic rock mass for applying boundary conditions based on the fact that elastic boundary will undergo elastic deformation according to the stress disturbance range of roadway excavation. Frictionless roller boundary condition was assigned to the five faces except the top face of the cuboid model. Stress of 11.5 MPa which has a direction of the negative Z was applied to the top boundary as the weight of about 460-m overburden strata. A circular roadway was excavated after the model was equilibrated under in situ stress state (σx=12MPa,σy=60MPa,σz=12MPa).Previous studies indicated that zonal disintegration may occur when in situ stress is greater than uniaxial compression strength of rock mass [10]. Thus, a 60 MPa-stress of y direction was applied to the zonal disintegration test. And the magnitude of σx and σz equal to the weight of about 500 m overburden strata. The radial displacement, radial stress, and tangential stress were monitored along the monitoring line shown in Figure 7 which has an interval of 1 m with an increase of time step.


## 3.2. Analysis on Numerical Test Result of Zonal Disintegration

Initiation, propagation, and coalescence progress of cracks around the roadway boundary are illustrated in Figure8. It can be found from Figure 8 that shear fracture and tensile fracture form in surrounding rock mass under the high stress parallel to the axial of the roadway. Fissured zones and intact zones distributing alternately was captured. The trace line of fracture presents as a spiral line close to the roadway boundary which is similar to the laboratory result performed by Yuan et al. [9], while slip-line fractures densely distribute in the deep surrounding rock of the roadway.Figure 8
Fracture distribution around the roadway boundary at different time step. (a) 5000 steps. (b) 10000 steps. (c) 18000 steps. (d) 34000 steps.
(a)
(b)
(c)
(d)Crack initially generated at the vicinity of the roadway boundary which can be observed in Figure8(a). According to Figure 8(b), two tensile fractures formed as the crack propagated into deep surrounding rock of the roadway. Subsequently, shear fractures propagate into the deep surrounding rock of the roadway in all directions. At the same time, vertical and horizontal tensile fractures generate at different positions in surrounding rock mass and then connected with the shear fractures. Surrounding rock mass was divided into discontinuous zones by those shear fractures and tensile fractures. Tensile fractures are mainly located at shallow surrounding rock of the roadway, while shear fractures distribute in the deep surrounding rock of the roadway.The radial displacement, radial stress, and tangential stress around the roadway acquired by the monitoring line are illustrated in Figure9. Figure 9(a) indicates that stress redistributed in surrounding rock mass after the roadway was excavated. Stress around the roadway internal surface decreased radically due to the unloading effect of the excavation. The magnitude of radial stress increases gradually to in situ level as the distance to the roadway rib increases.Figure 9
Monitoring results of (a) the stress around the roadway and (b) the radial displacement around the roadway.
(a)
(b)In addition, tangential stress increases by 2-3 times of in situ stress and then decreases to in situ level as the distance to the roadway rib increases. The global variation law of radial stress and tangential stress of surrounding rock mass is consistent with the traditional elasticity-plasticity theory [33], while the variation of radial stress and tangential stress oscillates as the distance to the roadway rib increases which is a significant difference to the traditional theory. The change law of stress around the roadway is similar to the stress distribution of non-Euclidean geometry theory proposed by Guzev and Paroshin [22]. As Figure 9(b) shows, positive value of displacement represents that surrounding rock deformed toward to the internal of the roadway boundary. The displacement around the roadway boundary is about 473 mm. In general, the displacement decreased as the distance to the roadway rib increased. However, on account of fractures generated at different depths of surrounding rock, the rock mass on either side of the fracture was subjected to stress in opposite directions, and the deformation of deep side rock mass around the fracture was restrained. Thus, the displacement of the deeper rock mass could be greater than the deep side rock mass around the fracture. We can observe from Figure 9(b) that displacement fluctuates slightly at the deep surrounding rock of the roadway.


## 3.3. Effect of Block Size on Numerical Test of Zonal Disintegration

Block size is an inherent property of 3DEC bonded block model. The smaller the block size is, the more the model tends to be isotropic and homogeneous. Two numerical models, comprising tetrahedral blocks which have an edge length of 0.8 m and 1 m, were constructed in order to study the effect of block size on zonal disintegration. The same intact rock properties, microcontact properties, boundary condition, and loading procedure mentioned in Section3.1 were applied to the constructed models. Numerical test result which adopted different block size models is shown in Figure 10.Figure 10
Fracture distribution around the roadway boundary using different sizes of tetrahedral block: (a) block with 0.6 m edge length; (b) block with 0.8 m edge length; (c) block with 1 m edge length.
(a)
(b)
(c)As Figure10 shows, fracture modes of different width tetrahedral block models are similar. Tensile fracture, which is approximately parallel to the roadway boundary, forms in all three models. However, it is interesting to find that the extent of fractured zone decreases as the block size increases which may be due to that the confinement of block movement is stronger after contacts separation as the block size increases. The numerical test using different sizes of the block indicated that tetrahedral block which has a width of 0.6 m is small enough to study zonal disintegration.Numerical test of zonal disintegration indicated that the fracture zone and intact zone alternately distributed around the roadway boundary under high axial stress parallel to the roadway axis. It is worthy to note that zonal disintegration occurs with generation and connection of shear cracks and tensile cracks. The trace line of the fracture presented as a spiral line close to the roadway boundary which is different from the circular fracture zone in proposed studies.


## 4. Sensitivity Analysis on Numerical Test of Zonal Disintegration

### 4.1. Effect of Stress Condition on Fracture Mode of Deep Rock Mass around the Roadway Boundary

It is well known that stress condition has a great effect on fracture modes of rock mass, such as tensile fracture tends to occur in uniaxial compression test while shear fracture forms in a triaxial compression test. It is necessary to study the effect of stress condition on the fracture mode of roadway surrounding rock mass. All conditions except stress condition of numerical simulation in this section are identical to the conditions informed in Section3.1 in order to investigate the fracture mode under different stress condition.Stress condition of the maximum principal stress is horizontal stress of which the direction is perpendicular to the roadway axis (σx=60MPa,σy=45MPa,σz=12MPa)was adopted in numerical study owing to the fact that horizontal stress is usually greater than vertical stress in the field measurement of geostress. The numerical simulation result is shown in Figure 11.Figure 11
Fracture distribution around the roadway boundary when the horizontal stress is maximum principal stress: (a) numerical result; (b) field observation result [34].
(a)
(b)As can be seen from Figure11, zonal disintegration was not captured under the stress condition of horizontal stress is the maximum principal stress. Fractures are merely formed around the roadway boundary. The trace line of V-shape fractures formed at the top and bottom of the roadway boundary. Sidewall of the roadway spalls due to generation of tensile fracture. The EDZ of the numerical result has a good agreement with the field monitoring result of Underground Research Laboratory [34].In addition, fracture mode varies under stress states of different confining pressure, which was comprehensively studied by Wu et al. with 2D DEM model [25]. The 3-dimension stress state was taken into account in this paper. Several numerical tests have been carried out under stress conditions of the stress parallel to the roadway axis is the maximum stress as confining pressure increases. The results of numerical test under different confining pressure (Table 3) are shown in Figure 12.Table 3
Stress conditions applied to the numerical model.
σ
x (MPa)
σ
y (MPa)
σ
z (MPa)
Stress condition A
12
60
12
Stress condition B
15
60
15
Stress condition C
30
60
30
Stress condition D
60
60
60Figure 12
Fracture distribution around the roadway boundary under different confining pressure: (a) stress condition A, (b) stress condition B, (c) stress condition C, and (d) stress condition D.
(a)
(b)
(c)
(d)As Figure12 shows, the extent of fractured zone decreases and fractures tend to be mainly distributed in shallow surrounding rock of the roadway as confining pressure increases. Only spiral-line-fracture could be observed around the roadway boundary when the confining pressure increases to 15 MPa. Furthermore, zonal disintegration cannot be observed when confining pressure increases to 30 MPa. The numerical simulation result indicated that zonal disintegration may occur when the stress condition satisfies a certain criterion and the direction of the maximum principal stress is parallel to the roadway axis. Zonal disintegration was also not captured under hydrostatic pressure. Rock mass deformed uniformly toward to the internal of the roadway from all directions. The displacement decreases as the distance to the roadway boundary increases.


### 4.2. Effect of Roadway Shape on Zonal Disintegration

Straight wall semicircular arch roadway and rectangular roadway model were constructed in order to study the effect of roadway shape on the fracture mode of zonal disintegration. All parameters except the roadway shape in this section are identical to the parameters used in Section3.1.As Figure13 shows, zonal disintegration phenomenon was obtained using different roadway shape model. The curvature of the fracture track line tends to be equal to the roadway boundary in shallow surrounding rock of the roadway. As we can see, the track line of fracture around the circular roadway boundary is approximately a spiral line, while fractures around the rectangular roadway boundary are orthogonal. However, fracture distribution in deep surrounding rock of the roadway was similar in different roadway shape model, arc fractures were all formed at a certain distance to the roadway boundary, which is basically consistence with the study carried out by Zuo et al. [13]. However, shear fractures which are parallel to each other formed in deep surrounding rock of the roadway, which is different from the previous study and it may be due to the different loading method.Figure 13
Fracture distribution around the roadway boundary using different roadway shape model: (a) circular roadway, (b) semicircular arch roadway, and (c) rectangular roadway.
(a)
(b)
(c)


## 4.1. Effect of Stress Condition on Fracture Mode of Deep Rock Mass around the Roadway Boundary

It is well known that stress condition has a great effect on fracture modes of rock mass, such as tensile fracture tends to occur in uniaxial compression test while shear fracture forms in a triaxial compression test. It is necessary to study the effect of stress condition on the fracture mode of roadway surrounding rock mass. All conditions except stress condition of numerical simulation in this section are identical to the conditions informed in Section3.1 in order to investigate the fracture mode under different stress condition.Stress condition of the maximum principal stress is horizontal stress of which the direction is perpendicular to the roadway axis (σx=60MPa,σy=45MPa,σz=12MPa)was adopted in numerical study owing to the fact that horizontal stress is usually greater than vertical stress in the field measurement of geostress. The numerical simulation result is shown in Figure 11.Figure 11
Fracture distribution around the roadway boundary when the horizontal stress is maximum principal stress: (a) numerical result; (b) field observation result [34].
(a)
(b)As can be seen from Figure11, zonal disintegration was not captured under the stress condition of horizontal stress is the maximum principal stress. Fractures are merely formed around the roadway boundary. The trace line of V-shape fractures formed at the top and bottom of the roadway boundary. Sidewall of the roadway spalls due to generation of tensile fracture. The EDZ of the numerical result has a good agreement with the field monitoring result of Underground Research Laboratory [34].In addition, fracture mode varies under stress states of different confining pressure, which was comprehensively studied by Wu et al. with 2D DEM model [25]. The 3-dimension stress state was taken into account in this paper. Several numerical tests have been carried out under stress conditions of the stress parallel to the roadway axis is the maximum stress as confining pressure increases. The results of numerical test under different confining pressure (Table 3) are shown in Figure 12.Table 3
Stress conditions applied to the numerical model.
σ
x (MPa)
σ
y (MPa)
σ
z (MPa)
Stress condition A
12
60
12
Stress condition B
15
60
15
Stress condition C
30
60
30
Stress condition D
60
60
60Figure 12
Fracture distribution around the roadway boundary under different confining pressure: (a) stress condition A, (b) stress condition B, (c) stress condition C, and (d) stress condition D.
(a)
(b)
(c)
(d)As Figure12 shows, the extent of fractured zone decreases and fractures tend to be mainly distributed in shallow surrounding rock of the roadway as confining pressure increases. Only spiral-line-fracture could be observed around the roadway boundary when the confining pressure increases to 15 MPa. Furthermore, zonal disintegration cannot be observed when confining pressure increases to 30 MPa. The numerical simulation result indicated that zonal disintegration may occur when the stress condition satisfies a certain criterion and the direction of the maximum principal stress is parallel to the roadway axis. Zonal disintegration was also not captured under hydrostatic pressure. Rock mass deformed uniformly toward to the internal of the roadway from all directions. The displacement decreases as the distance to the roadway boundary increases.


## 4.2. Effect of Roadway Shape on Zonal Disintegration

Straight wall semicircular arch roadway and rectangular roadway model were constructed in order to study the effect of roadway shape on the fracture mode of zonal disintegration. All parameters except the roadway shape in this section are identical to the parameters used in Section3.1.As Figure13 shows, zonal disintegration phenomenon was obtained using different roadway shape model. The curvature of the fracture track line tends to be equal to the roadway boundary in shallow surrounding rock of the roadway. As we can see, the track line of fracture around the circular roadway boundary is approximately a spiral line, while fractures around the rectangular roadway boundary are orthogonal. However, fracture distribution in deep surrounding rock of the roadway was similar in different roadway shape model, arc fractures were all formed at a certain distance to the roadway boundary, which is basically consistence with the study carried out by Zuo et al. [13]. However, shear fractures which are parallel to each other formed in deep surrounding rock of the roadway, which is different from the previous study and it may be due to the different loading method.Figure 13
Fracture distribution around the roadway boundary using different roadway shape model: (a) circular roadway, (b) semicircular arch roadway, and (c) rectangular roadway.
(a)
(b)
(c)


## 5. Conclusions

Zonal disintegration is extremely different from the excavation damaged zone around the underground opening based on elastic-plastic theory. The existence of zonal disintegration was initially validated by laboratory test which indicated that high stress parallel to the underground opening axis is the main trigger of zonal disintegration. Fracture mode of zonal disintegration has been intensely debated since it is difficult to obtain the fracture distribution of surrounding rock mass directly. Annular fractured zone and intact zone alternately distribute around the underground opening boundary as demonstrated in proposed numerical and theoretical studies. However, the previous study was almost based on continuum theory through which postpeak behavior of rock mass may not be well explained. The effect of the stress parallel to the underground opening axis was not studied in some papers which adopt 2D discontinuous medium theory.3DEC bonded block model was employed in this study. First, microcontact properties were calibrated so that the macromechanics properties of the numerical model were consistent with the sandstone sample. The uniaxial compression test and Brazilian test demonstrated that postpeak behavior of rock mass could be realistically simulated by using the 3DEC bonded block model. Subsequently, several numerical tests of zonal disintegration were conducted which indicated that zonal disintegration may occur when the stress parallel to the axial of the roadway is maximum principal stress. The trace line of fracture around the roadway boundary presented as a spiral line close to the roadway boundary which is different from previous numerical studies but has a good agreement with the laboratory test. Zonal disintegration is formed by generation and coalescence of shear fractures and tensile fractures. However, uncertainty of engineering geological conditions was not taken into account in this study. Reliability analysis method, like the first-order reliability method [35, 36], should be employed to extend this study to practical engineering. This still needs further study in the future.Finally, sensitive analysis was performed to study the effect of stress conditions and roadway shape on zonal disintegration. The following conclusions were drawn:(1)
Zonal disintegration may occur when the direction of maximum principal stress is parallel to the roadway axis and the magnitude of maximum principal stress is greater than the UCS of surrounding rock.(2)
A spiral-line fracture formed by coalescence of tensile cracks in shallow surrounding rock of the roadway, while slip-line fractures composed of shear cracks distributed in deep surrounding rock of the roadway.(3)
The extent of the fracture zone decreased as the confining pressure increased. Zonal disintegration will not occur when the confining pressure reaches a certain value.(4)
The curvature of the fracture track line tends to be equal to the roadway boundary in shallow surrounding rock of the roadway, while the fractures in deep surrounding rock of the roadway mass are similar as the roadway shape changes.



---

*Source: 3589417-2019-10-27.xml*