# Field Detection and Numerical Simulation Study on Roof Asymmetric Fracture Mechanism of Large-Span Soft Rock Roadway with Discontinuity Surface

Authors: Liang Chen; Shouyang Gao; Bin Miao; Nan Yao; Shizhong Zhang; Wei Zhang
Journal: Geofluids
 (2022)

---

## Abstract

Existence of roof discontinuity surface is one of the extremely important factors causing the asymmetric fracture of roadway roof, especially for large-span soft rock roadway. In this paper, a crack density coefficient (D) is firstly defined by local thresholding-microcell segmentation method and used to analyze the evolution law of D with roof drilling depth (d) (the lower and upper of roof discontinuity are defined as regions I and II, respectively). Then, UDEC numerical simulation is used to investigate the asymmetry evolution law of roof total displacement, maximum principal stress, and crack density with stress release coefficient (α) considering the effect of discontinuity surface. Research results indicate that (1) the roof parameters (D) in regions I and II both show a negative logarithmic function decreasing trend with the increase of drilling depth. When d<4.5 m, the parameter (D) in region II is about 1.5 times that in region I; when d≥4.5 m, the parameter (D) in region I is almost zero, while the parameter (D) in region II maintains a slight decreasing trend. Roof failure presents asymmetric distribution characteristics along both sides of the discontinuous surface. (2) In the initial stage of open-off cutting excavation, the top left and right corners of roof as well as the bottom of discontinuous surface first occurred the tension-shear failure, and then as the α increases, the two side cracks gradually shift to the middle of roof in regions I and II with the appearance of stress concentration in two top corners. Meanwhile, the direction of maximum principal stress is transformed from “direction approximately parallel to the excavation surface” to “direction perpendicular to the discontinuous surface,” of which the location transfers to the deep anisotropically, forming an asymmetric stress loosening zones in the roof of regions I and II. The range of stress loosening zone in region II is significantly larger than that in region I. When the surrounding rock stress is completely released, the roof cracks in region I gradually transition from nonconnected to connected state and form a “quasi-right triangle” loosening zone. In addition, an “isosceles triangle-like” high crack density loosening zone with the roof middle as the axis is also formed in region II. The roadway roof presents a markedly asymmetric caving feature. (3) With the increase of discontinuity surface angles, the roof fracture range gradually decreases in region I and increases in region II as well as the structural features of roof pressure-bearing arch transform from “left-low right-high continuous asymmetric structure” (30°-60°) to “left-low right-high discontinuous asymmetric structure” (90°) to “unilateral partial pressure-bearing arch structure” (120°-150°). Research results can provide an extremely important reference for the optimization and design of support scheme with discontinuity-thick soft rock roof-large span roadways.


---

## Body

## 1. Introduction

Roof stability is an extremely important factor in ensuring safe and efficient mining. However, under the influence of geological tectonic stress, the discontinuity surface, such as faults, joints, and cracks, is widely distributed in rock stratum, which has a significant impact on the roadway-excavated roof stability, especially for the large-span soft rock roadway [1–4]. In the United States, the casualties caused by the roadway roof caving account for about 16% of coal mine accidents, while the proportion is as high as 33.2% in China [5]. Therefore, exploring the roof instability mechanism of large-span soft rock roadways with discontinuity surfaces can provide an important reference for the selection of support schemes and parameter optimization of such roadways.Previous scholars have carried out a large number of researches on roof instability mechanism of different strata structures by using the on-site monitoring, similarity simulation, and numerical simulation methods. The roof rock strata structures mainly included upper-soft and lower-hard strata, upper-hard and lower-soft strata, containing faults strata, and regenerated roof strata as well as laminated-jointed roof strata [6–11]. Above research results showed that the characteristics of roof strata structure have a significant impact on the roof instability mode of roadway, especially for large-span roadway. Gao et al. [12] investigated the shear failure phenomenon of a shale roadway roof by using an innovative numerical approach and confirmed the time sequence of marked microseismic activity, significant stress changes, and accelerated displacement during the process of a roof fall. Shen et al. [13] monitored the displacement, stress changes, and seismic activities during a roof fall by using an integrated roof monitoring system and found that the changes associated with roof falls follow a time sequence of initial seismicity, followed by stress changes and then displacement. Li [14, 15] researched the deformation and fracture characteristic of large-span open-off cut by applying theoretical analysis, numerical calculation, and site measurement. In addition, Yu et al. [16] studied the evolution laws of stress, displacement, and plastic zone of roadways with a compound roof and believed that the compound roof is vulnerable to abscission and instability with the separation of the layer caving. In addition, Wu et al. [17] analyzed the roof deformation and failure characteristics of large inclination roadway under the conditions of no weak structural plane and with the weak structural plane by using the FLAC3D software. The results showed that the slip of roof weak structural plane is an extremely important reason leading to the change of maximum deformation location and instability mode. Compared with faults, both bedding and joints are weak structural plane with smaller scale inside the stratum. Due to the existence of such structural planes, the slip-shear failure of roadway roof is easy to occur in the excavation process, which leads to asymmetric fracture phenomenon of roadway roof, and even induce roof collapse [18–20]. Fan and Jiang [21] focused on the role of weak structural bodies to investigate the roadway instability and found that the local weak structure planes inside the stratum can easily lead to asymmetrical instability of roadway, and the initial failure position is often near these weak structure planes. Niu et al. [22] believed that the large span roadway with the thin roof is prone to shear sliding failure along the weak structural plane of strata, which leads to shear failure of support system. Above scholars have carried out many studies on the roof instability phenomenon of roadway with different strata structures, but most studies only consider the effect of a single factor as the engineering background, such as fault, weak structure planes, and soft-hard strata feature. In practice engineering, the strength of surrounding rock on both sides of the weak structure plane is often different under the coupling effect of groundwater and discontinuity, which will lead to an upper-soft and lower-hard strata or upper-hard and lower-soft strata structure. However, comprehensive influence of such stratum structure and discontinuity surface on roof stability of coal mine roadway have rarely been considered in previous studies. This is critical for optimizing the roof support scheme and parameter.In this paper, the open-off cutting of B903 working face in No.4 coal mine of the Saier Group is taken as the engineering background. Firstly, a crack density coefficient (D) is defined by local thresholding-microcell segmentation method and used to analyze the evolution law of D with the roof drilling depth. Then, the UDEC numerical simulation method is used to investigate the asymmetry evolution law of roof total displacement, maximum principal stress, and crack density with the stress release coefficient (α) considering the effect of roof discontinuity surface. Finally, the roof asymmetric fracture mechanism of large-span soft rock roadways with discontinuity surface is revealed.


## 2. Field Detection on Asymmetric Fracture and Deformation Behavior of Roof

### 2.1. Engineering Geological Profile

Saier No.4 Mine is located in Tacheng area of Xinjiang Uygur Autonomous Region, as shown in Figure1(a). At present, B9# coal seam is mainly mined, of which the average thickness is 4.0 m and the average thickness of the direct roof is 24.9 m. Most of the coal seam is calcareous weakly cemented mudstone, which is gray and jointed, belonging to the thick mudstone roof coal seam. The floor is composed of mudstone, B10# coal seam, and clay rock from shallow to deep, in which the direct bottom mudstone is distributed in blocks with low hardness, locally mixed with plant debris and mostly plastic rock mass; the average thickness of B10# coal seam is 9.5 m, with relatively developed bedding and poor integrity; the average thickness of clay rock is 11.2 m, and the structure is dense, in which the expansion and disintegration characteristics are obvious, and the sandy mudstone bands are locally mixed. Figure 1(b) shows the comprehensive histogram of B9# coal seam. During the excavation process of open-off cutting in B903# working face, a weak structural plane running through the mudstone roof is exposed on the roof, with the structural plane moving towards 252° and inclination angle of about 60° from north to west. During the excavation of open-off cutting to stability, the right side of the structural plane is affected by the roof water, and the large deformation as well as roof failure problem occurs locally. The roof shows obvious asymmetric large deformation failure characteristics with the structural plane as the boundary point (Figure 1(c)). In order to facilitate the study, the left area of the roof discontinuity is defined as region I, and the right area is defined as region II.Figure 1
Location map (a), histogram of the strata (b), and fracture photos (c) of roadway [23].
(a)(b)(c)


### 2.2. Characteristics of Asymmetric Strength Distribution

Figure2 shows the stress-strain curves of roof mudstone in regions I and II under the uniaxial compression. It can be seen from the figure that the uniaxial compressive strength (UCS) of the roof mudstone in regions I and II is 27.20 MPa and 15.86 MPa, respectively. The UCS of former is about 1.72 times that of the latter, indicating that the roof mudstone strength in region II is significantly lower than that in region I due to the influence of roof water and discontinuity. In addition, the postpeak curve of mudstone in regions II and I shows an obvious plastic failure and brittle failure characteristic, respectively, indicating that the mudstone in region II is more easily softened. With discontinuity as the boundary line, the mechanical properties of roof surrounding rock also show asymmetric distribution characteristics, which will exert a significant impact on the roof deformation characteristics and fracture mode of surrounding rock.Figure 2
Typical stress-strain curves of mudstone at different roof locations.


### 2.3. Characteristics of Asymmetric Fracture

#### 2.3.1. Crack Extraction Method of Borehole Peep Image

The borehole peep image can not only effectively observe the dynamic change of surrounding rock structure but also effectively characterize the fracture-loosening degree of surrounding rock. There are two main methods for characterizing borehole fractures. Firstly, the borehole peep image is expanded into a plane image, and then, the digital characteristics of fractures are extracted by using digital image processing technology; secondly, the relevant image processing software is used to directly process the high-definition borehole image taken by borehole peeping. Since the image quality is different under different equipment and photography environments, the large deviations are often presented in the extraction of fracture information [24]. Especially when far away from the light system, the difference between fracture images and background environments is small, and it is not easy to extract fracture characteristic. At present, for image segmentation technology, there are two widely adopted methods, namely, global and local threshold methods [25]. Among them, the adaptive histogram equalization method belongs to the range of local threshold method, which is a computer image processing algorithm to improve image contrast. This algorithm can play an important role in improving the local contrast of the image to obtain more image details.In order to effectively obtain the specific crack information in the local image of boreholes, the AHE algorithm is adopted in this paper to segment the image of the borehole peeping results (Figure3(c)). Compared with the global threshold algorithm, this algorithm can completely retain the image crack contour and improve the crack identification range, but there are local disturbance points. For weakening the influence of disturbance points on image processing results, the 8×8 corrosion structure element is used to reduce the corrosion noise of image binarization results, and then, the processed image is obtained as shown in Figure 3(d); then, the image is divided by equidistant concentric circles whose origin is in the middle of borehole image (Figure 3(e)); finally, the morphological region props function is used to count the area of each ring fracture, and the fracture size of each microelement is characterized. The specific digital feature extraction process of borehole peeping fractures is shown in Figure 3.Figure 3
Digital feature extraction process of borehole crack.Figure4 shows the crack identification results in the roof of regions I and II. It can be seen from Figure 4 that there are different degrees of fracture regions in the shallow surrounding rock of roadways, whether in the roof of region I or II; with the gradual increase of drilling depth, the degree of fracture development is gradually decreasing; the development degree of roof cracks in region I is lower, and the cross of cracks is less, mainly isolated cracks; the fracture density and opening degree of the roof in region II are larger than those in region I, and the fractures are mainly grid-shaped, with high degree of fragmentation and obvious instability trend.Figure 4
Fracture distribution morphology of roof in regions I and II under different borehole depths (note: the red represents identified cracks).
(a)
Region I(b)
Region II


#### 2.3.2. Definition of Crack Density

To quantitatively characterize the fracture characteristics of roof surrounding rock, the polar coordinate system is used to segment the fracture image by the MATLAB software. The coordinate origin is located in the middle of the borehole peep image. The schematic image segmented by polar coordinates is shown in Figure5. In the picture, ΔSABCD represents the microunit after image segmentation; red represents the cracks in cutting ring; Ri and Rj, respectively, represent the outer diameter and inner diameter of cutting ring. The schematic diagram of quantitative characterization principle of borehole fracture degree based on microelement method is shown in Figure 5.Figure 5
Schematic diagram of quantitative characterization principle of borehole fracture degree based on microelement method.According to the geometric relationship in Figure5, the ring area (Ri, Rj) where the microunit ΔSABCD is located can be expressed as follows:
(1)ΔSij=π×Ri2−Rj2,whereΔSij is the ring area (Ri, Rj) where the microunit ΔSABCD is located. The mean value of the roof fracture density D′ of the ring can be expressed as follows [23]:
(2)D′=∑1nΔSfπ×Ri2−Rj2,whereΔSf is the fracture area contained in the minimum statistical unit of the ring where the microunit ΔSABCD is located and n is the number of minimum elements with cracks in the ring where the microunit ΔSABCD is located. Thus, the average crack density coefficient of borehole peep image is as follows:
(3)D=∑1kD′Rmax−Rmin/J−1,whereRmax and Rmin are the maximum and minimum radius of the circle in the image recognition range, J is the ring spacing, and k is the number of rings.


## 2.1. Engineering Geological Profile

Saier No.4 Mine is located in Tacheng area of Xinjiang Uygur Autonomous Region, as shown in Figure1(a). At present, B9# coal seam is mainly mined, of which the average thickness is 4.0 m and the average thickness of the direct roof is 24.9 m. Most of the coal seam is calcareous weakly cemented mudstone, which is gray and jointed, belonging to the thick mudstone roof coal seam. The floor is composed of mudstone, B10# coal seam, and clay rock from shallow to deep, in which the direct bottom mudstone is distributed in blocks with low hardness, locally mixed with plant debris and mostly plastic rock mass; the average thickness of B10# coal seam is 9.5 m, with relatively developed bedding and poor integrity; the average thickness of clay rock is 11.2 m, and the structure is dense, in which the expansion and disintegration characteristics are obvious, and the sandy mudstone bands are locally mixed. Figure 1(b) shows the comprehensive histogram of B9# coal seam. During the excavation process of open-off cutting in B903# working face, a weak structural plane running through the mudstone roof is exposed on the roof, with the structural plane moving towards 252° and inclination angle of about 60° from north to west. During the excavation of open-off cutting to stability, the right side of the structural plane is affected by the roof water, and the large deformation as well as roof failure problem occurs locally. The roof shows obvious asymmetric large deformation failure characteristics with the structural plane as the boundary point (Figure 1(c)). In order to facilitate the study, the left area of the roof discontinuity is defined as region I, and the right area is defined as region II.Figure 1
Location map (a), histogram of the strata (b), and fracture photos (c) of roadway [23].
(a)(b)(c)


## 2.2. Characteristics of Asymmetric Strength Distribution

Figure2 shows the stress-strain curves of roof mudstone in regions I and II under the uniaxial compression. It can be seen from the figure that the uniaxial compressive strength (UCS) of the roof mudstone in regions I and II is 27.20 MPa and 15.86 MPa, respectively. The UCS of former is about 1.72 times that of the latter, indicating that the roof mudstone strength in region II is significantly lower than that in region I due to the influence of roof water and discontinuity. In addition, the postpeak curve of mudstone in regions II and I shows an obvious plastic failure and brittle failure characteristic, respectively, indicating that the mudstone in region II is more easily softened. With discontinuity as the boundary line, the mechanical properties of roof surrounding rock also show asymmetric distribution characteristics, which will exert a significant impact on the roof deformation characteristics and fracture mode of surrounding rock.Figure 2
Typical stress-strain curves of mudstone at different roof locations.


## 2.3. Characteristics of Asymmetric Fracture

### 2.3.1. Crack Extraction Method of Borehole Peep Image

The borehole peep image can not only effectively observe the dynamic change of surrounding rock structure but also effectively characterize the fracture-loosening degree of surrounding rock. There are two main methods for characterizing borehole fractures. Firstly, the borehole peep image is expanded into a plane image, and then, the digital characteristics of fractures are extracted by using digital image processing technology; secondly, the relevant image processing software is used to directly process the high-definition borehole image taken by borehole peeping. Since the image quality is different under different equipment and photography environments, the large deviations are often presented in the extraction of fracture information [24]. Especially when far away from the light system, the difference between fracture images and background environments is small, and it is not easy to extract fracture characteristic. At present, for image segmentation technology, there are two widely adopted methods, namely, global and local threshold methods [25]. Among them, the adaptive histogram equalization method belongs to the range of local threshold method, which is a computer image processing algorithm to improve image contrast. This algorithm can play an important role in improving the local contrast of the image to obtain more image details.In order to effectively obtain the specific crack information in the local image of boreholes, the AHE algorithm is adopted in this paper to segment the image of the borehole peeping results (Figure3(c)). Compared with the global threshold algorithm, this algorithm can completely retain the image crack contour and improve the crack identification range, but there are local disturbance points. For weakening the influence of disturbance points on image processing results, the 8×8 corrosion structure element is used to reduce the corrosion noise of image binarization results, and then, the processed image is obtained as shown in Figure 3(d); then, the image is divided by equidistant concentric circles whose origin is in the middle of borehole image (Figure 3(e)); finally, the morphological region props function is used to count the area of each ring fracture, and the fracture size of each microelement is characterized. The specific digital feature extraction process of borehole peeping fractures is shown in Figure 3.Figure 3
Digital feature extraction process of borehole crack.Figure4 shows the crack identification results in the roof of regions I and II. It can be seen from Figure 4 that there are different degrees of fracture regions in the shallow surrounding rock of roadways, whether in the roof of region I or II; with the gradual increase of drilling depth, the degree of fracture development is gradually decreasing; the development degree of roof cracks in region I is lower, and the cross of cracks is less, mainly isolated cracks; the fracture density and opening degree of the roof in region II are larger than those in region I, and the fractures are mainly grid-shaped, with high degree of fragmentation and obvious instability trend.Figure 4
Fracture distribution morphology of roof in regions I and II under different borehole depths (note: the red represents identified cracks).
(a)
Region I(b)
Region II


### 2.3.2. Definition of Crack Density

To quantitatively characterize the fracture characteristics of roof surrounding rock, the polar coordinate system is used to segment the fracture image by the MATLAB software. The coordinate origin is located in the middle of the borehole peep image. The schematic image segmented by polar coordinates is shown in Figure5. In the picture, ΔSABCD represents the microunit after image segmentation; red represents the cracks in cutting ring; Ri and Rj, respectively, represent the outer diameter and inner diameter of cutting ring. The schematic diagram of quantitative characterization principle of borehole fracture degree based on microelement method is shown in Figure 5.Figure 5
Schematic diagram of quantitative characterization principle of borehole fracture degree based on microelement method.According to the geometric relationship in Figure5, the ring area (Ri, Rj) where the microunit ΔSABCD is located can be expressed as follows:
(1)ΔSij=π×Ri2−Rj2,whereΔSij is the ring area (Ri, Rj) where the microunit ΔSABCD is located. The mean value of the roof fracture density D′ of the ring can be expressed as follows [23]:
(2)D′=∑1nΔSfπ×Ri2−Rj2,whereΔSf is the fracture area contained in the minimum statistical unit of the ring where the microunit ΔSABCD is located and n is the number of minimum elements with cracks in the ring where the microunit ΔSABCD is located. Thus, the average crack density coefficient of borehole peep image is as follows:
(3)D=∑1kD′Rmax−Rmin/J−1,whereRmax and Rmin are the maximum and minimum radius of the circle in the image recognition range, J is the ring spacing, and k is the number of rings.


## 2.3.1. Crack Extraction Method of Borehole Peep Image

The borehole peep image can not only effectively observe the dynamic change of surrounding rock structure but also effectively characterize the fracture-loosening degree of surrounding rock. There are two main methods for characterizing borehole fractures. Firstly, the borehole peep image is expanded into a plane image, and then, the digital characteristics of fractures are extracted by using digital image processing technology; secondly, the relevant image processing software is used to directly process the high-definition borehole image taken by borehole peeping. Since the image quality is different under different equipment and photography environments, the large deviations are often presented in the extraction of fracture information [24]. Especially when far away from the light system, the difference between fracture images and background environments is small, and it is not easy to extract fracture characteristic. At present, for image segmentation technology, there are two widely adopted methods, namely, global and local threshold methods [25]. Among them, the adaptive histogram equalization method belongs to the range of local threshold method, which is a computer image processing algorithm to improve image contrast. This algorithm can play an important role in improving the local contrast of the image to obtain more image details.In order to effectively obtain the specific crack information in the local image of boreholes, the AHE algorithm is adopted in this paper to segment the image of the borehole peeping results (Figure3(c)). Compared with the global threshold algorithm, this algorithm can completely retain the image crack contour and improve the crack identification range, but there are local disturbance points. For weakening the influence of disturbance points on image processing results, the 8×8 corrosion structure element is used to reduce the corrosion noise of image binarization results, and then, the processed image is obtained as shown in Figure 3(d); then, the image is divided by equidistant concentric circles whose origin is in the middle of borehole image (Figure 3(e)); finally, the morphological region props function is used to count the area of each ring fracture, and the fracture size of each microelement is characterized. The specific digital feature extraction process of borehole peeping fractures is shown in Figure 3.Figure 3
Digital feature extraction process of borehole crack.Figure4 shows the crack identification results in the roof of regions I and II. It can be seen from Figure 4 that there are different degrees of fracture regions in the shallow surrounding rock of roadways, whether in the roof of region I or II; with the gradual increase of drilling depth, the degree of fracture development is gradually decreasing; the development degree of roof cracks in region I is lower, and the cross of cracks is less, mainly isolated cracks; the fracture density and opening degree of the roof in region II are larger than those in region I, and the fractures are mainly grid-shaped, with high degree of fragmentation and obvious instability trend.Figure 4
Fracture distribution morphology of roof in regions I and II under different borehole depths (note: the red represents identified cracks).
(a)
Region I(b)
Region II


## 2.3.2. Definition of Crack Density

To quantitatively characterize the fracture characteristics of roof surrounding rock, the polar coordinate system is used to segment the fracture image by the MATLAB software. The coordinate origin is located in the middle of the borehole peep image. The schematic image segmented by polar coordinates is shown in Figure5. In the picture, ΔSABCD represents the microunit after image segmentation; red represents the cracks in cutting ring; Ri and Rj, respectively, represent the outer diameter and inner diameter of cutting ring. The schematic diagram of quantitative characterization principle of borehole fracture degree based on microelement method is shown in Figure 5.Figure 5
Schematic diagram of quantitative characterization principle of borehole fracture degree based on microelement method.According to the geometric relationship in Figure5, the ring area (Ri, Rj) where the microunit ΔSABCD is located can be expressed as follows:
(1)ΔSij=π×Ri2−Rj2,whereΔSij is the ring area (Ri, Rj) where the microunit ΔSABCD is located. The mean value of the roof fracture density D′ of the ring can be expressed as follows [23]:
(2)D′=∑1nΔSfπ×Ri2−Rj2,whereΔSf is the fracture area contained in the minimum statistical unit of the ring where the microunit ΔSABCD is located and n is the number of minimum elements with cracks in the ring where the microunit ΔSABCD is located. Thus, the average crack density coefficient of borehole peep image is as follows:
(3)D=∑1kD′Rmax−Rmin/J−1,whereRmax and Rmin are the maximum and minimum radius of the circle in the image recognition range, J is the ring spacing, and k is the number of rings.


## 3. Result Analysis

Figure6 shows the variation curve of the roof crack density coefficient with drilling depth in regions I and II. It can be seen from Figure 6 that the roof crack density coefficient shows a negative logarithmic decreasing trend as drilling depth increases gradually; when drilling depth≤4.5 m, the roof crack density in region II is about 1.5 times that in region I; when the drilling depth≥4.5 m, the roof crack density in region I is almost zero, while the roof crack density in region II basically maintains a slight downward trend. Therefore, the roof fracture degree and depth in region II are larger than those in region I; the roof overall presents asymmetric fracture characteristics.Figure 6
Average crack density coefficient varies with borehole depth in I and II regions.


## 4. Numerical Simulation on Asymmetric Progressive Instability Mechanism

### 4.1. Model Establishment and Validation

#### 4.1.1. Model Establishment

To effectively reveal the roof asymmetric fracture mechanism of large-span soft rock roadways with discontinuity, combined with the field geological conditions, the UDEC software is used to study the progressive evolution law of total displacement, maximum principal stress, and roof crack density of roadways. Figure7 shows the numerical simulation model of geometry and boundary conditions. In the figure, both the width and height of the model are 60 m. The roadway excavation area is set in the middle of the model. The excavation size of roadway is 8m×4m. The distance between the discontinuous structural plane and the upper left corner of the roof is 2.8 m, and the angle is 120°. The lateral pressure coefficient of the model is 1.5. Displacement constraints are applied to the lateral boundary and the bottom boundary of the model, and stress constraints are applied to the top boundary. The stress level is 15.6 MPa. In addition, in the simulation calculation process, the elastic constitutive model is adopted for the block, and the Mohr-Coulomb constitutive model is adopted for the discontinuous contact surface. In order to improve the calculation accuracy, the small hexagon block unit with an average length of 0.01 m is used in the surrounding rock near the roadway, and the hexagon unit with an average length of 0.05 m is used in other areas.Figure 7
Geometry and boundary conditions of numerical computation model.The process of roadway excavation is essentially the process of surrounding rock stress readjustment and gradual release to stability; the FISH function “ZONK” built in UDEC is used to simulate the gradual release of stress. In this simulation, the stress applied on the roadway surface is released in five stages, and the corresponding stress release coefficient (i) is gradually changed from 0 to 1, with an increment of 0.2 (the value indicates that 20% of the stress applied on the roadway surface is released). In each loosening stage, sufficient numerical calculation time steps are allowed to ensure the stable response of the model.As for accurately describing the difference in stress, deformation, and fracture distribution of deep surrounding rock in regions I and II on both sides of roadway roof, one measuring line and six monitoring points are arranged vertically in the middle of the two regions, respectively. The distance between the monitoring points of the same measuring line is 0.6 m, and the distance between the two measuring lines is 4.0 m. The first measuring point is 0.2 m away from the surface of cutting roof. Among them, the 1#-4# measuring points on the left side of the discontinuous structural plane are located in region I, and others are located in region II. Each measuring point on the right side of the discontinuous structural plane is located in region II. The specific arrangement of measuring points is shown in Figure7.


#### 4.1.2. Model Configuration

Taking the roof mudstone from regions I and II as the research object, uniaxial compression tests were carried out on cylindrical specimens with a size of1m×2m to check the accuracy of model parameters. In this simulation, the UDEC-Voronoi command is adopted to divide specimens into several discrete hexagonal elements, and the elements are set as the isotropic elastic model, and the joint plane is set as the Coulomb slip model. The physical and mechanics parameters are shown in Table 1. During the simulation, the specimen is fixed and loaded by loading plates on the upper and lower ends, and the mechanical parameters of loading plate are set as those of the steel plate. Among them, fixed constraints are applied on the lower loading plate, the x-direction speed of the upper loading plate is set to 0, the y-direction loading speed is set to 0.05 m/s, and the time step is 3×10−7 s/step; the calculated loading rate is equivalent to 1.5×10−8 s/step. Figure 8 shows the comparison of the axial stress-strain curve obtained from this simulation and experiment. It can be seen from the figure that the simulation results are basically consistent with the test results, indicating the rationality of the simulation parameter selection.Table 1
Physical and mechanics parameters for the simulated model.
LithologyMatrix propertiesContact propertiesDensity (kg/m3)E∗ (GPa)kn∗ (GPa/m)ks∗ (GPa/m)Friction angle (°)Cohesion (MPa)Tensile strength (MPa)Region I23503.416.86.8322.71.23Region II25007.038.415.4216.02.0Coal14002.88433.6301.40.5Mudstone25007.038.415.4216.02.0Claystone22006.0160649.5202.5Figure 8
Comparison of stress-strain curves obtained from simulation and test on roof rock in I and II regions.


### 4.2. Evolutionary Characteristics of Total Displacement

Figure9 shows the cloud diagram of total displacement under different stress release coefficients. The variation law of the total displacement at each measuring point of roof with the calculation time step is shown in Figure 10. The following can be seen from Figures 9 and 10:
(i)
With the increase of the stress release coefficient, the total roof displacement is asymmetrically distributed along both sides of the discontinuous surface. Wheni≤0.4, there is a small total displacement only in the shallow surrounding rock of region II, and the value is less than 20 mm, while there is no deformation in the surrounding rock of region I. When i increases from 0.4 to 0.8, the shallow surrounding rock in region I begins to deform slightly, but the maximum deformation is significantly lower than that in region II; at this time, the overall deformation of the roof gradually extends from the shallow to the deep, and the separation phenomenon begins to appear in the shallow surrounding rock. When the stress is completely released, the total displacement of surrounding rock surface in region II reaches 250 mm, while the total displacement of surrounding rock surface in region I is significantly lower than 75 mm, indicating that the rock structure is relatively complete(ii)
There are significant differences in the deformation amounts and rates of surrounding rock between regions I and II. When the distance from the measuring point to the roofsurface≤2.0 m, at the same interface, the deformation amounts and rates in region II are obviously greater than those in region I under any stress release degree, but the difference gradually decreases with the increase of depth. When the measuring point is more than 2.0 m away from the roof surface, the measuring points 5#, 6#, and 11#-12# are all in region II, and the deformation rates and amounts of the four measuring points are significantly lower, and the numerical differences are smaller, indicating that the surrounding rock at this part is relatively intact(iii)
The roof caving phenomenon has highly positive correlation with the stress release degree of surrounding rock. Wheni≤0.6, the overall deformation of the roof is smaller, and there is no obvious macro fracture zone on the roof surface; as i increases from 0.6 to 0.8, obvious shedding blocks begin to appear in the shallow part, but the deep surrounding rock is still relatively intact; when i=1.0, the fracture phenomenon of roof begins to extend to deep rock stratum, and the fracture height in region II is significantly greater than that in region I, and the risk of asymmetric roof caving is greatly increased.Figure 9
Cloud diagram of total displacement under different stress release coefficients.Figure 10
Variation law of total displacement at each measuring point with the calculation time step.


### 4.3. Evolutionary Characteristics of Maximum Principal Stress

Figure11 shows the cloud diagram of the maximum principal stress of roof under different stress release coefficients. The variation law of the maximum principal stress at each measuring point with the calculation time step is shown in Figure 12. The following can be seen from Figures 11 and 12:
(i)
The significant distribution difference of maximum principal stress(1)
The distribution range between the stress increasing zone and stress loosening zone of roof is obviously asymmetrical. Wheni≤0.4, the tensile-shear failure first occurs in the shallow surrounding rock of region II due to the low strength and large span, forming a shallow stress loosening zone; meanwhile, there is obvious stress concentration in the top left corner of the roof in region I, and the direction of maximum principal stress is transformed from “direction approximately parallel to the excavation surface” to “direction perpendicular to the discontinuous surface,” and the initial stress increasing zone is formed. When 0.4≤i≤0.8, the stress adjustment trend of roof surrounding rock in region II is obvious, and the stress increasing zone is constantly transferred to the deep, resulting in the formation of a wide range of stress loosening zone at the bottom, and the fracture phenomenon of surrounding rock is significant(2)
When the stress is completely released, the surrounding rock stress tends to be stable as the calculation time step adds; at the same time, the height of roof stress loosening zone in regions I and II is lower than 0.8 m and slightly larger than 2.0 m, respectively(ii)
The variation trend of the roof maximum principal stress is significantly different(1)
Wheni≤0.4, the surrounding rock stress is at the slow release stage, and the maximum principal stress at measuring points 1# and 7# is 14.71 MPa and 12.39 MPa, respectively, indicating that the stress release degree of surrounding rock in region II is obviously greater than that in region I after the roadway excavation. When 0.4≤i≤1, the surrounding rock stress is constantly adjusted, and the high stress is gradually transferred to the deep; when the α is 0.6 and 0.8, respectively, it can be seen that the stress decreasing degree in region II is faster than that in region I by matching the maximum principal stress of the monitoring point at the same depth; compared with those in region I, the stress loosening zone has a larger range and a wider failure depth in region II, which is consistent with the deformation monitoring results of numerical simulation(2)
When the stress is released completely, the surrounding rock deformation is at a stable stage as the calculated time step adds. At this time, the maximum principal stress at monitoring points 1# and 2# in region I is 0.77 MPa and 3.87 MPa, respectively, both in the postpeak residual strength stage, and the other monitoring points are in the prepeak stable bearing stage; the maximum principal stress at monitoring points 7#, 8#, 9#, and 10# in region II is 0 MPa, 1.16 MPa, 4.65 MPa, and 11.03 MPa, respectively, all in the postpeak residual strength stage, while the other monitoring points are in the prepeak stable bearing stage; the results show that the variation of surrounding rock stress is one of the important factors for the roof asymmetric instability of roadwaysFigure 11
Cloud diagram of maximum principal stress under different stress release coefficients.Figure 12
Variation law of maximum principal stress at each measuring point with the calculation time step: (a) region I; (b) region II.
(a)(b)


### 4.4. Evolutionary Characteristics of Fracture Density

To study the evolution characteristics and concentration distribution of cracks during the progressive failure of roof, the roof crack development degree is extracted by ARCGIS linear density analysis method, and the calculation principle of the extraction method is shown in Figure13. Linear density analysis is mainly used to calculate the density of linear elements in the field of each output grid cell. Its principle is to make a circle with the center of each grid cell as the center of the circle according to the search radius and calculate the sum of the length of each line falling into the circle and population and then divide by the circle area to obtain the grid cell density. The calculation equation is as follows:
(4)δ=∑k=1βLk.VkA,Figure 13
Schematic diagram of line density analysis method.whereδ is the linear density, β is the number of lines in a circle, Lk is the line length, Vk is the population, and A is the circle area.Figure14 shows the cloud diagram of roof crack density under different stress release coefficients. Figure 15 is the rose diagram of fracture orientation and frequency. The following can be seen from figures:
(i)
Wheni≤0.4, the locations with the high crack density are concentrated in the two top corners and discontinuous surfaces, which are consistent with stress concentration location. When 0.4≤i≤0.8, the crack density in the direction of connection between left corner (no. 3) and discontinuity bottom (no. 2) or middle (no. 4) is higher for the roof in region I; meanwhile, the isosceles triangle zone of roof in region II, which is composed of the connecting line between right corner (no. 1) and discontinuity bottom (no. 2), has a higher crack density. The crack development height of two regions gradually increase and extend to the deep as the stress release coefficient increases. When the stress is completely released (i=1.0), the roof cracks in regions I and II quickly expand to the deep with the increasing calculated time step; in addition, the fracture zone of “quasi-right triangle” is gradually formed in region I, in which the cracks are mainly concentrated in the range of 0-1 m, and an isosceles triangle-like fracture zone is formed in region II, where the fracture density is larger and the risk of roof collapse is more aggravated(ii)
The fracture development is mainly dominated by the separation cracks parallel to the excavation surface, supplemented by a small number of oblique cracks near the macroscopic fracture surface of the roof. The total number of surrounding rock fractures in regions I and II is 2351 and 11029, respectively, and the final height of fractures is about 1.2 m and 2.5 m, respectively, which is basically consistent with the field drilling peep resultsFigure 14
Cloud diagram of roof crack density under different stress release coefficients.Figure 15
Rose diagram of crack orientation and frequency: (a) region I; (b) region II.
(a)(b)


### 4.5. Roof Progressive Failure Mode of Long-Span Roadway with Discontinuities

Figure16 shows the gradual asymmetric instability process of the roof during the cutting excavation on B903 working face. Figure 16(a) shows that the surrounding rock is at the initial stress state before roadway excavation. From Figure 16(b), the tensile-shear failure occurs in the shallow surrounding rock at the corner of the roof and the bottom of the discontinuous surface in regions I and II at the initial stage of roadway excavation. Then, as the stress release coefficient increases, the crack gradually transfers to the middle of roof in regions I and II; in the transferring process, the stress concentration first occurs in top left and right corners of roof, and the direction of maximum principal stress is transformed from “direction approximately parallel to the excavation surface” to “direction perpendicular to the discontinuous surface.” The position of the maximum principal stress in the two regions extends nonisokinetically from the shallow to the deep and finally forms an asymmetric stress loosening zone of roof in regions I and II (Figure 16(c)).With the further increase of the stress release coefficient, the crack height of roof in two regions expands further; the roof cracks in region I gradually transition from nonconnected to connected state and form a “quasi-right triangle.” In addition, an “isosceles triangle-like” high crack density loosening zone with the roof middle as the axis is also formed in region II (Figure 16(d)); the loosening range of roof in region II is significantly larger than that in region I. The phenomenon of layer separation is prominent in two loosening zones; finally, the asymmetric caving instability phenomenon of roof is formed, which is characterized by “large scale caving occurs in region II, and partial layered block shedding occurs in region I.”Figure 16
Schematic diagram of progressive asymmetrical instability of roof of large-span roadway with discontinuous surface (Δh is the difference of roof fracture height between regions I and II).


## 4.1. Model Establishment and Validation

### 4.1.1. Model Establishment

To effectively reveal the roof asymmetric fracture mechanism of large-span soft rock roadways with discontinuity, combined with the field geological conditions, the UDEC software is used to study the progressive evolution law of total displacement, maximum principal stress, and roof crack density of roadways. Figure7 shows the numerical simulation model of geometry and boundary conditions. In the figure, both the width and height of the model are 60 m. The roadway excavation area is set in the middle of the model. The excavation size of roadway is 8m×4m. The distance between the discontinuous structural plane and the upper left corner of the roof is 2.8 m, and the angle is 120°. The lateral pressure coefficient of the model is 1.5. Displacement constraints are applied to the lateral boundary and the bottom boundary of the model, and stress constraints are applied to the top boundary. The stress level is 15.6 MPa. In addition, in the simulation calculation process, the elastic constitutive model is adopted for the block, and the Mohr-Coulomb constitutive model is adopted for the discontinuous contact surface. In order to improve the calculation accuracy, the small hexagon block unit with an average length of 0.01 m is used in the surrounding rock near the roadway, and the hexagon unit with an average length of 0.05 m is used in other areas.Figure 7
Geometry and boundary conditions of numerical computation model.The process of roadway excavation is essentially the process of surrounding rock stress readjustment and gradual release to stability; the FISH function “ZONK” built in UDEC is used to simulate the gradual release of stress. In this simulation, the stress applied on the roadway surface is released in five stages, and the corresponding stress release coefficient (i) is gradually changed from 0 to 1, with an increment of 0.2 (the value indicates that 20% of the stress applied on the roadway surface is released). In each loosening stage, sufficient numerical calculation time steps are allowed to ensure the stable response of the model.As for accurately describing the difference in stress, deformation, and fracture distribution of deep surrounding rock in regions I and II on both sides of roadway roof, one measuring line and six monitoring points are arranged vertically in the middle of the two regions, respectively. The distance between the monitoring points of the same measuring line is 0.6 m, and the distance between the two measuring lines is 4.0 m. The first measuring point is 0.2 m away from the surface of cutting roof. Among them, the 1#-4# measuring points on the left side of the discontinuous structural plane are located in region I, and others are located in region II. Each measuring point on the right side of the discontinuous structural plane is located in region II. The specific arrangement of measuring points is shown in Figure7.


### 4.1.2. Model Configuration

Taking the roof mudstone from regions I and II as the research object, uniaxial compression tests were carried out on cylindrical specimens with a size of1m×2m to check the accuracy of model parameters. In this simulation, the UDEC-Voronoi command is adopted to divide specimens into several discrete hexagonal elements, and the elements are set as the isotropic elastic model, and the joint plane is set as the Coulomb slip model. The physical and mechanics parameters are shown in Table 1. During the simulation, the specimen is fixed and loaded by loading plates on the upper and lower ends, and the mechanical parameters of loading plate are set as those of the steel plate. Among them, fixed constraints are applied on the lower loading plate, the x-direction speed of the upper loading plate is set to 0, the y-direction loading speed is set to 0.05 m/s, and the time step is 3×10−7 s/step; the calculated loading rate is equivalent to 1.5×10−8 s/step. Figure 8 shows the comparison of the axial stress-strain curve obtained from this simulation and experiment. It can be seen from the figure that the simulation results are basically consistent with the test results, indicating the rationality of the simulation parameter selection.Table 1
Physical and mechanics parameters for the simulated model.
LithologyMatrix propertiesContact propertiesDensity (kg/m3)E∗ (GPa)kn∗ (GPa/m)ks∗ (GPa/m)Friction angle (°)Cohesion (MPa)Tensile strength (MPa)Region I23503.416.86.8322.71.23Region II25007.038.415.4216.02.0Coal14002.88433.6301.40.5Mudstone25007.038.415.4216.02.0Claystone22006.0160649.5202.5Figure 8
Comparison of stress-strain curves obtained from simulation and test on roof rock in I and II regions.


## 4.1.1. Model Establishment

To effectively reveal the roof asymmetric fracture mechanism of large-span soft rock roadways with discontinuity, combined with the field geological conditions, the UDEC software is used to study the progressive evolution law of total displacement, maximum principal stress, and roof crack density of roadways. Figure7 shows the numerical simulation model of geometry and boundary conditions. In the figure, both the width and height of the model are 60 m. The roadway excavation area is set in the middle of the model. The excavation size of roadway is 8m×4m. The distance between the discontinuous structural plane and the upper left corner of the roof is 2.8 m, and the angle is 120°. The lateral pressure coefficient of the model is 1.5. Displacement constraints are applied to the lateral boundary and the bottom boundary of the model, and stress constraints are applied to the top boundary. The stress level is 15.6 MPa. In addition, in the simulation calculation process, the elastic constitutive model is adopted for the block, and the Mohr-Coulomb constitutive model is adopted for the discontinuous contact surface. In order to improve the calculation accuracy, the small hexagon block unit with an average length of 0.01 m is used in the surrounding rock near the roadway, and the hexagon unit with an average length of 0.05 m is used in other areas.Figure 7
Geometry and boundary conditions of numerical computation model.The process of roadway excavation is essentially the process of surrounding rock stress readjustment and gradual release to stability; the FISH function “ZONK” built in UDEC is used to simulate the gradual release of stress. In this simulation, the stress applied on the roadway surface is released in five stages, and the corresponding stress release coefficient (i) is gradually changed from 0 to 1, with an increment of 0.2 (the value indicates that 20% of the stress applied on the roadway surface is released). In each loosening stage, sufficient numerical calculation time steps are allowed to ensure the stable response of the model.As for accurately describing the difference in stress, deformation, and fracture distribution of deep surrounding rock in regions I and II on both sides of roadway roof, one measuring line and six monitoring points are arranged vertically in the middle of the two regions, respectively. The distance between the monitoring points of the same measuring line is 0.6 m, and the distance between the two measuring lines is 4.0 m. The first measuring point is 0.2 m away from the surface of cutting roof. Among them, the 1#-4# measuring points on the left side of the discontinuous structural plane are located in region I, and others are located in region II. Each measuring point on the right side of the discontinuous structural plane is located in region II. The specific arrangement of measuring points is shown in Figure7.


## 4.1.2. Model Configuration

Taking the roof mudstone from regions I and II as the research object, uniaxial compression tests were carried out on cylindrical specimens with a size of1m×2m to check the accuracy of model parameters. In this simulation, the UDEC-Voronoi command is adopted to divide specimens into several discrete hexagonal elements, and the elements are set as the isotropic elastic model, and the joint plane is set as the Coulomb slip model. The physical and mechanics parameters are shown in Table 1. During the simulation, the specimen is fixed and loaded by loading plates on the upper and lower ends, and the mechanical parameters of loading plate are set as those of the steel plate. Among them, fixed constraints are applied on the lower loading plate, the x-direction speed of the upper loading plate is set to 0, the y-direction loading speed is set to 0.05 m/s, and the time step is 3×10−7 s/step; the calculated loading rate is equivalent to 1.5×10−8 s/step. Figure 8 shows the comparison of the axial stress-strain curve obtained from this simulation and experiment. It can be seen from the figure that the simulation results are basically consistent with the test results, indicating the rationality of the simulation parameter selection.Table 1
Physical and mechanics parameters for the simulated model.
LithologyMatrix propertiesContact propertiesDensity (kg/m3)E∗ (GPa)kn∗ (GPa/m)ks∗ (GPa/m)Friction angle (°)Cohesion (MPa)Tensile strength (MPa)Region I23503.416.86.8322.71.23Region II25007.038.415.4216.02.0Coal14002.88433.6301.40.5Mudstone25007.038.415.4216.02.0Claystone22006.0160649.5202.5Figure 8
Comparison of stress-strain curves obtained from simulation and test on roof rock in I and II regions.


## 4.2. Evolutionary Characteristics of Total Displacement

Figure9 shows the cloud diagram of total displacement under different stress release coefficients. The variation law of the total displacement at each measuring point of roof with the calculation time step is shown in Figure 10. The following can be seen from Figures 9 and 10:
(i)
With the increase of the stress release coefficient, the total roof displacement is asymmetrically distributed along both sides of the discontinuous surface. Wheni≤0.4, there is a small total displacement only in the shallow surrounding rock of region II, and the value is less than 20 mm, while there is no deformation in the surrounding rock of region I. When i increases from 0.4 to 0.8, the shallow surrounding rock in region I begins to deform slightly, but the maximum deformation is significantly lower than that in region II; at this time, the overall deformation of the roof gradually extends from the shallow to the deep, and the separation phenomenon begins to appear in the shallow surrounding rock. When the stress is completely released, the total displacement of surrounding rock surface in region II reaches 250 mm, while the total displacement of surrounding rock surface in region I is significantly lower than 75 mm, indicating that the rock structure is relatively complete(ii)
There are significant differences in the deformation amounts and rates of surrounding rock between regions I and II. When the distance from the measuring point to the roofsurface≤2.0 m, at the same interface, the deformation amounts and rates in region II are obviously greater than those in region I under any stress release degree, but the difference gradually decreases with the increase of depth. When the measuring point is more than 2.0 m away from the roof surface, the measuring points 5#, 6#, and 11#-12# are all in region II, and the deformation rates and amounts of the four measuring points are significantly lower, and the numerical differences are smaller, indicating that the surrounding rock at this part is relatively intact(iii)
The roof caving phenomenon has highly positive correlation with the stress release degree of surrounding rock. Wheni≤0.6, the overall deformation of the roof is smaller, and there is no obvious macro fracture zone on the roof surface; as i increases from 0.6 to 0.8, obvious shedding blocks begin to appear in the shallow part, but the deep surrounding rock is still relatively intact; when i=1.0, the fracture phenomenon of roof begins to extend to deep rock stratum, and the fracture height in region II is significantly greater than that in region I, and the risk of asymmetric roof caving is greatly increased.Figure 9
Cloud diagram of total displacement under different stress release coefficients.Figure 10
Variation law of total displacement at each measuring point with the calculation time step.


## 4.3. Evolutionary Characteristics of Maximum Principal Stress

Figure11 shows the cloud diagram of the maximum principal stress of roof under different stress release coefficients. The variation law of the maximum principal stress at each measuring point with the calculation time step is shown in Figure 12. The following can be seen from Figures 11 and 12:
(i)
The significant distribution difference of maximum principal stress(1)
The distribution range between the stress increasing zone and stress loosening zone of roof is obviously asymmetrical. Wheni≤0.4, the tensile-shear failure first occurs in the shallow surrounding rock of region II due to the low strength and large span, forming a shallow stress loosening zone; meanwhile, there is obvious stress concentration in the top left corner of the roof in region I, and the direction of maximum principal stress is transformed from “direction approximately parallel to the excavation surface” to “direction perpendicular to the discontinuous surface,” and the initial stress increasing zone is formed. When 0.4≤i≤0.8, the stress adjustment trend of roof surrounding rock in region II is obvious, and the stress increasing zone is constantly transferred to the deep, resulting in the formation of a wide range of stress loosening zone at the bottom, and the fracture phenomenon of surrounding rock is significant(2)
When the stress is completely released, the surrounding rock stress tends to be stable as the calculation time step adds; at the same time, the height of roof stress loosening zone in regions I and II is lower than 0.8 m and slightly larger than 2.0 m, respectively(ii)
The variation trend of the roof maximum principal stress is significantly different(1)
Wheni≤0.4, the surrounding rock stress is at the slow release stage, and the maximum principal stress at measuring points 1# and 7# is 14.71 MPa and 12.39 MPa, respectively, indicating that the stress release degree of surrounding rock in region II is obviously greater than that in region I after the roadway excavation. When 0.4≤i≤1, the surrounding rock stress is constantly adjusted, and the high stress is gradually transferred to the deep; when the α is 0.6 and 0.8, respectively, it can be seen that the stress decreasing degree in region II is faster than that in region I by matching the maximum principal stress of the monitoring point at the same depth; compared with those in region I, the stress loosening zone has a larger range and a wider failure depth in region II, which is consistent with the deformation monitoring results of numerical simulation(2)
When the stress is released completely, the surrounding rock deformation is at a stable stage as the calculated time step adds. At this time, the maximum principal stress at monitoring points 1# and 2# in region I is 0.77 MPa and 3.87 MPa, respectively, both in the postpeak residual strength stage, and the other monitoring points are in the prepeak stable bearing stage; the maximum principal stress at monitoring points 7#, 8#, 9#, and 10# in region II is 0 MPa, 1.16 MPa, 4.65 MPa, and 11.03 MPa, respectively, all in the postpeak residual strength stage, while the other monitoring points are in the prepeak stable bearing stage; the results show that the variation of surrounding rock stress is one of the important factors for the roof asymmetric instability of roadwaysFigure 11
Cloud diagram of maximum principal stress under different stress release coefficients.Figure 12
Variation law of maximum principal stress at each measuring point with the calculation time step: (a) region I; (b) region II.
(a)(b)


## 4.4. Evolutionary Characteristics of Fracture Density

To study the evolution characteristics and concentration distribution of cracks during the progressive failure of roof, the roof crack development degree is extracted by ARCGIS linear density analysis method, and the calculation principle of the extraction method is shown in Figure13. Linear density analysis is mainly used to calculate the density of linear elements in the field of each output grid cell. Its principle is to make a circle with the center of each grid cell as the center of the circle according to the search radius and calculate the sum of the length of each line falling into the circle and population and then divide by the circle area to obtain the grid cell density. The calculation equation is as follows:
(4)δ=∑k=1βLk.VkA,Figure 13
Schematic diagram of line density analysis method.whereδ is the linear density, β is the number of lines in a circle, Lk is the line length, Vk is the population, and A is the circle area.Figure14 shows the cloud diagram of roof crack density under different stress release coefficients. Figure 15 is the rose diagram of fracture orientation and frequency. The following can be seen from figures:
(i)
Wheni≤0.4, the locations with the high crack density are concentrated in the two top corners and discontinuous surfaces, which are consistent with stress concentration location. When 0.4≤i≤0.8, the crack density in the direction of connection between left corner (no. 3) and discontinuity bottom (no. 2) or middle (no. 4) is higher for the roof in region I; meanwhile, the isosceles triangle zone of roof in region II, which is composed of the connecting line between right corner (no. 1) and discontinuity bottom (no. 2), has a higher crack density. The crack development height of two regions gradually increase and extend to the deep as the stress release coefficient increases. When the stress is completely released (i=1.0), the roof cracks in regions I and II quickly expand to the deep with the increasing calculated time step; in addition, the fracture zone of “quasi-right triangle” is gradually formed in region I, in which the cracks are mainly concentrated in the range of 0-1 m, and an isosceles triangle-like fracture zone is formed in region II, where the fracture density is larger and the risk of roof collapse is more aggravated(ii)
The fracture development is mainly dominated by the separation cracks parallel to the excavation surface, supplemented by a small number of oblique cracks near the macroscopic fracture surface of the roof. The total number of surrounding rock fractures in regions I and II is 2351 and 11029, respectively, and the final height of fractures is about 1.2 m and 2.5 m, respectively, which is basically consistent with the field drilling peep resultsFigure 14
Cloud diagram of roof crack density under different stress release coefficients.Figure 15
Rose diagram of crack orientation and frequency: (a) region I; (b) region II.
(a)(b)


## 4.5. Roof Progressive Failure Mode of Long-Span Roadway with Discontinuities

Figure16 shows the gradual asymmetric instability process of the roof during the cutting excavation on B903 working face. Figure 16(a) shows that the surrounding rock is at the initial stress state before roadway excavation. From Figure 16(b), the tensile-shear failure occurs in the shallow surrounding rock at the corner of the roof and the bottom of the discontinuous surface in regions I and II at the initial stage of roadway excavation. Then, as the stress release coefficient increases, the crack gradually transfers to the middle of roof in regions I and II; in the transferring process, the stress concentration first occurs in top left and right corners of roof, and the direction of maximum principal stress is transformed from “direction approximately parallel to the excavation surface” to “direction perpendicular to the discontinuous surface.” The position of the maximum principal stress in the two regions extends nonisokinetically from the shallow to the deep and finally forms an asymmetric stress loosening zone of roof in regions I and II (Figure 16(c)).With the further increase of the stress release coefficient, the crack height of roof in two regions expands further; the roof cracks in region I gradually transition from nonconnected to connected state and form a “quasi-right triangle.” In addition, an “isosceles triangle-like” high crack density loosening zone with the roof middle as the axis is also formed in region II (Figure 16(d)); the loosening range of roof in region II is significantly larger than that in region I. The phenomenon of layer separation is prominent in two loosening zones; finally, the asymmetric caving instability phenomenon of roof is formed, which is characterized by “large scale caving occurs in region II, and partial layered block shedding occurs in region I.”Figure 16
Schematic diagram of progressive asymmetrical instability of roof of large-span roadway with discontinuous surface (Δh is the difference of roof fracture height between regions I and II).


## 5. Discussion

Gao et al. [12] studied the failure mode of roadway roof without the discontinuous surface, and they believed that shear failure of the roadway roof initiates at the roadway corners and then progressively propagates deeper into the roof, finally forming a large scale roof failure. But, the roof failure model is symmetrical as shown in Figure 17(a). In fact, due to the influence of diagenesis and geological structures [26, 27], the discontinuous surface is widely distributed in rock strata. Under the erosion of local groundwater, the strength of rocks will also be reduced, and eventually, surrounding rocks with different strength will be formed on both sides of the discontinuous surface. Under the coupling influence of these two aspects, it is easy to cause roof asymmetric failure of roadway (Figure 17(b)).Figure 17
Roof fracture and instability characteristics of roadway with or without discontinuity: (a) fracture characteristic without discontinuity; (b) fracture characteristic with discontinuity.
(a)(b)To further reveal the influence of discontinuity angles on the asymmetric failure mechanism of roadway roof, numerical simulation models of roadway as discontinuity angle of 30°, 60°, 90°, 120°, and 150°are established, as shown in Figure18. Figure 19 presents the vector cloud diagram of the maximum principal stress of roadway roof under different discontinuity angles. Figure 20 shows the monitoring curves of the maximum principal stress at each monitoring point under different discontinuity angles. It can be seen from Figures 19 and 20 that:
(i)
When discontinuityangle<90°, with the gradual increase of the dip angle, the stress increasing zone of roof in region I gradually shifts from the top left corner to the middle, and the stress concentration decreases gradually, but the failure depth has no significant increase trend. With the increase of dip angles, the effective bearing capacity provided by the shallow surrounding rock in region I to the surrounding rock in region II is gradually reduced, so that the stress increasing zone of roof in region II gradually shifts to the deep, the stress loosening range is gradually enlarged, and the overall stability of roof is getting worse and worse(ii)
When discontinuityangle=90°, as the span of region I is smaller than that of region II, the value of the stress increasing zone in region I is larger than that in region II, but the height is smaller than that in region II, and the roadway failure is distributed in the shallow part of region I and below the middle part of region II(iii)
When90°≤discontinuityangle≤150°, under the influence of reaction force, the roof bearing capacity of region II to region I gradually increases, making the roof stability of region I gradually be improved, and decreasing the roof stability of region II, resulting in that only a small range of separation failure occurs of roof in region I, and a large range of collapse and caving occurs in region II. Above all, with the increase of discontinuity surface angles, the roof fracture range gradually decreases in region I and increases in region II as well as the structural features of roof pressure-bearing arch is transformed from “left-low right-high continuous asymmetric structure” (30°-60°) to “left-low right-high discontinuous asymmetric structure” (90°) to “unilateral partial pressure-bearing arch structure” (120°-150°). The structure feature of pressure-bearing arch is closely related to the sliding movement of discontinuity. When discontinuity angle<90°, the stratum structure shows upper-soft and lower-hard characteristic, which is not conducive to the sliding movement of discontinuous surface, so it is easy to form continuous bearing arch; when discontinuity angle=90°, the vertical movement of discontinuity causes the uneven stress distribution of two side roof and then leads to the occurrence of discontinuous stress concentration zone; when 90°≤discontinuityangle≤150°, the stratum structure shows upper-hard and lower-soft characteristic, which is not conducive to the stability of the lower strata, so that the area of stress increasing zone in region II is transferred to the deep, resulting in the occurrence of unilateral partial load-bearing arch structureFigure 18
Schematic diagram of discontinuity angle arrangement and monitoring point distribution.Figure 19
Vector illustration of maximum principal stress under different discontinuity surface angles.Figure 20
Variation law of maximum principal stress at different measuring points with the calculation time step.


## 6. Conclusions

(i)
A crack density coefficient (D), which can effectively characterize the development degree of roof cracks, is defined by local thresholding-microcell segmentation method and used to investigate the evolution law of roof fracture degree with drilling depth. The roof parameters (D) in regions I and II both show a negative logarithmic function decreasing trend with the increasing drilling depth. When d≤4.5m, the parameter (D) in region II is about 1.5 times that in region I; when d≥4.5 m, the parameter (D) in region I is almost zero, while the parameter (D) in region II maintains a slight decreasing trend. The roof failure presents asymmetric distribution characteristics along both sides of the discontinuous surface(ii)
In the process of the open-off cutting excavation to stability, the total displacement and maximum principal stress of roof are asymmetrically distributed along both sides of the discontinuous surface. Wheni≤0.4, the tensile-shear failure first occurs in top left and right corners of roof as well as the bottom of discontinuous surface with the appearance of stress concentration in two top corners. When 0.4≤i≤0.8, as the stress release coefficient increases, the direction of maximum principal stress is transformed from “direction approximately parallel to the excavation surface” to “direction perpendicular to the discontinuous surface” with the location of maximum principal stress gradually transfers to the deep. Meanwhile, the stress loosening zones are formed in the shallow surrounding rock, and the total displacement of roof in regions I and II also increases significantly. When the stress is completely released, the height of roof stress loosening zone in regions I and II is less than 0.8 m and slightly larger than 2.0 m, respectively. The roof asymmetric caving phenomenon occurs(iii)
In the process of the open-off cutting excavation to stability, the crack distribution of roof is asymmetrically distributed along both sides of the discontinuous surface. Wheni≤0.4, the locations with high crack density are concentrated in the two top corners and discontinuous surfaces, which are consistent with stress concentration location. When i increases from 0.4 to 0.8, the crack density in the direction of connection between left corner and discontinuity bottom or middle is higher for the roof in region I; meanwhile, the isosceles triangle zone of the roof in region II, which is composed of the connecting line between right corner and discontinuity bottom, has a higher crack density. The crack development height of two regions gradually increases and extends to the deep with the increase of stress release coefficient. When the stress is completely released, the roof cracks in region I gradually transition from nonconnected to connected state and form a “quasi-right triangle.” In addition, an “isosceles triangle-like” high crack density loosening zone with the roof middle as the axis is also formed in region II. The loosening range of roof in region II is significantly larger than that in region I. The phenomenon of layer separation is prominent in two loosening zones, and the risk of the asymmetric roof collapse along discontinuities increases(iv)
With the increase of discontinuity surface angles, the roof fracture range gradually decreases in region I and increases in region II as well as the structural features of roof pressure-bearing arch are transformed from “left-low right-high continuous asymmetric structure” (30°-60°) to “left-low right-high discontinuous asymmetric structure” (90°) to “unilateral partial pressure-bearing arch structure” (120°-150°)



---

*Source: 8945239-2022-06-21.xml*