# Numerical Simulation and Analysis of Tunnel Failure Mode by Stochastic Fracture Model

Authors: Fukun Xiao; Lei Xu; Gang Liu; Zhiyuan Hou; Le Xing
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The dip angle, length, spacing, and fracture distance of rock fissure affect the morphology of roadway after collapse. The numerical simulation software CDEM is used to simulate the morphology of roadway collapse. The Monte Carlo model is used to simulate different types of crack models in two-dimensional plane and generate different crack models. The effects of crack angle, crack length, fracture distance, and spacing on the deformation of surrounding rock are analyzed. The influence of different rock burst on the failure strap-fall modes of fissure roadway and roadway in different sections is analyzed, and the stability law of roadway is studied. Under the condition of high stress, the roadway shape has little influence on the distribution of the principal stress difference of surrounding rock, but the equivalent excavation radius determines the distribution of the plastic zone of surrounding rock. The larger the ineffective reinforcement zone is, the larger the deformation around the roadway will be. The decrease of the angle between the structural plane and the vertical stress increases the failure range of the roadway under the gravity burst pressure. Under the horizontal tectonic stress type rock burst, when the structural plane inclination angle is 0°, the two-sided caving body fills the roadway and the roof caving range becomes smaller.


---

## Body

## 1. Introduction

Research on the development of joints and fractures is the basis of analyzing the impact of collapse roadway on emergency rescue. The complexity of geological conditions, unscientific roadway development, or improper supporting methods in coal mining will not only cause damage to the roadway to varying degrees, but also consume a large amount of manpower and material resources for secondary maintenance, which not only increases the cost of production, but also causes great difficulties for the next step of roadway excavation. It is more likely that the roof will fall in a large area when the surrounding rock of the roadway is affected by the rock burst in the incomplete state, which seriously threatens the safety of the roadway production.Wang [1] applied the principle of equivalent mechanics, position vector analysis, and virtual work stability to conduct a preliminary analysis on the mechanical properties and stability of blocks. Yan [2], based on the fabric characteristics of rockfill, used mathematical statistics as a tool to study the contact force between particles and the normal distribution of the contact force and obtained the probability distribution of the contact force between particles. Cong [3] analyzed the stability of the accumulation body in the reservoir area and obtained the stability of the accumulation body in the process of water level height change. Xu et al. [4] used digital image technology to study the mesomechanical properties of large accumulative bodies. Wang et al. [5] deduced the constitutive models of different dimensions of the caving ore granular body based on the constitutive relationship of the caving ore granular body, which clarified the relationship between the strength of the granular body and the flow generated during ore drawing. Barton [6] proposed classification and illustrated it by means of field examples and selected case records. Burnett and Holford [7] proposed an infinite element. Roscoe [8–10] proposed a new energy equation. Sun Guangzhong and Bai et al. [11, 12] determined a reasonable width for the corresponding narrow pillar. Bienlawski [13, 14] came up with a rock classification system. Kulatiab [15, 16] proposed a numerical method that combines analytical decomposition modeling with statistical simulation. Wang Li and Tang [17–19] have made an analysis of the stability of rock in underground engineering. Xing et al. [20] analyzed the geotechnical centrifugal model experiment. Yang et al. [21] have carried out a simulation test on the interaction of the construction of phase variable span tunnel on the ground subsidence. Tang et al. [22] introduced the numerical test of rock failure. Wang and Xing [23] introduced discrete elements in numerical simulation experiments. Wang [24] used discrete element method to conduct numerical analysis on the stability and mechanism of unsupported tunnels in jointed fractured rock mass at different buried depths. Goodman [25] gave a systematic overview of rock mechanics.The discrete element method is based on the continuum. This method can integrate the finite element algorithm with the discrete element algorithm, and the finite element method is used to perform operations in the model, while the discrete element method is used to process the boundary area of the model. By studying the model element, the internal force of the model element reaches the strain limit of the model and destroys it. It can not only deform or destroy the model of continuum shape under the action of external force, but also realize the destruction process of model material from continuum to discontinuum.


## 2. Numerical Simulation Model Establishment and Block Fracture Principle

In the semispring contact model, there is no need to determine what contact the block belongs to, so there is no need to calculate the contact area of the block, and the contact force can be directly calculated. In this method, the half spring is used as the main linked list to determine the corresponding target surface. When the half spring is inside the target surface or above the boundary, the contact can be made and the corresponding interpolation can be calculated. Because the semispring has its own characteristic area, there is no need to calculate the overlapping area between blocks in the process of measuring the contact force, thus reducing the calculation amount.


### 2.1. Establishment of Numerical Simulation Model

The geometric model is established as shown in Figure1. The model is 8 m in length and 8 m in height. The upper boundary of the model is set as in situ stress, the left boundary is horizontal tectonic stress, and the lower boundary is set as roller support. The roadway is built in the coal seam, and the roadway parameters are 2 m high and 2 m wide. The specific parameters for solving the model are mainly set in Table 1. Figure 1 shows the roadway stress model. The numerical model of 8 m ∗ 8 m is more consistent with the real tunnel shape.Table 1
Numerical model parameters.
ParameterThe mean modulus of elasticityE0 (MPa)Average compressive strengthE0 (MPa)Friction angle (°)Bulk density (kg/m3)Poisson’s ratioCoal seam18000253813400.24Rock72100903025300.25Figure 1
Force model of roadway.


### 2.2. Principle of Block Contact in the Model

Figure2 is a schematic diagram of contact pairs. The normal and tangential spring forces of the semispring are calculated as follows:(1)Fnj=−Knj×Δdjn,Fsj=−Ksj×Δdsj.Figure 2
Schematic diagram of contact pairs.Fnj and Fsj represent the normal and tangential forces of the J root; Knj and Ksj represent the normal and tangential stiffness of the J root; and Δdnj and Δdsj represent the normal and tangential displacement of the J root.In the process of failure calculation, the Mohr–Coulomb criterion is used to correct the spring force of the above equation, namely, the following equation, whereT represents the tensile strength, φ the internal friction angle, and C the cohesion force.(2)1If−Fnj≥T,Fnj=Fsj=0,C=0,T=0,2If−Fsj≥Fnj×tanϕ+C,Fsj=Fnj×tanϕ+C,C=0,T=0.


## 2.1. Establishment of Numerical Simulation Model

The geometric model is established as shown in Figure1. The model is 8 m in length and 8 m in height. The upper boundary of the model is set as in situ stress, the left boundary is horizontal tectonic stress, and the lower boundary is set as roller support. The roadway is built in the coal seam, and the roadway parameters are 2 m high and 2 m wide. The specific parameters for solving the model are mainly set in Table 1. Figure 1 shows the roadway stress model. The numerical model of 8 m ∗ 8 m is more consistent with the real tunnel shape.Table 1
Numerical model parameters.
ParameterThe mean modulus of elasticityE0 (MPa)Average compressive strengthE0 (MPa)Friction angle (°)Bulk density (kg/m3)Poisson’s ratioCoal seam18000253813400.24Rock72100903025300.25Figure 1
Force model of roadway.


## 2.2. Principle of Block Contact in the Model

Figure2 is a schematic diagram of contact pairs. The normal and tangential spring forces of the semispring are calculated as follows:(1)Fnj=−Knj×Δdjn,Fsj=−Ksj×Δdsj.Figure 2
Schematic diagram of contact pairs.Fnj and Fsj represent the normal and tangential forces of the J root; Knj and Ksj represent the normal and tangential stiffness of the J root; and Δdnj and Δdsj represent the normal and tangential displacement of the J root.In the process of failure calculation, the Mohr–Coulomb criterion is used to correct the spring force of the above equation, namely, the following equation, whereT represents the tensile strength, φ the internal friction angle, and C the cohesion force.(2)1If−Fnj≥T,Fnj=Fsj=0,C=0,T=0,2If−Fsj≥Fnj×tanϕ+C,Fsj=Fnj×tanϕ+C,C=0,T=0.


## 3. Calculation Model of Fracture Surface

### 3.1. Multiscale Computing Framework

The choice of scale requires the following assumptions:(1)
Consistent with the principles of finite element method and related numerical method based on continuous model, the continuous model is dispersed in the form of grouping or grid, and the properties of the obtained material are determined by the properties of each internal unit, so it has the property of homogenization in a certain field. If the cracks in the unit cannot be ignored, the unit size cannot reach the standard, small size and high precision. When the size is refined and the stress field and displacement field in the field change little, it is the maximum allowable scale of the element.(2)
The damage in the nonuniform block is in a layered form, and the failure surface is in the middle of the failure layer. It has the property of equal strain, and the scale of failure surface is consistent with that of element. When the unit appears damaged, it represents the damage of the damage layer, which is its neutral surface.(3)
The thickness of the failure layer is confirmed according to the size of the tiny cracks, which may be the thickness of the shear band, and its length is higher than the “spacing” grade of the tiny cracks.(4)
In the failure surface, the material strength is not uniform and equal, and it has a certain distribution law. The failure surface is divided into many elements. According to the distribution law, the strength of different elements will be different. In the case of failure, because the failure surface has the characteristic of equal strain, the microelements with different strength will have different strains when the failure occurs, which leads to the failure surface showing a progressive form.(5)
The microelement strain strength in the failure surface is closely related to the material structure of small scale, which can be obtained by using the determined displacement boundary conditions. If the microunits can be represented by the lowest scale particles of a given material property, all multiscale calculations can be performed. If a smaller scale calculation is required, it is possible to redistribute the equal stress according to hypothesis (1) and eventually obtain smaller infinitesimals. Repeat the calculation until the unit scale of the given material properties is obtained.In exposition (1) and (2), the parameters of the material itself and some parameters of the corresponding failure surface are given in the same unit. Assumptions (3) and (4) refer to the comparison of low-level unit size of material properties and failure parameters.Figure3 explains the calculation problems of different size regions in the case of a certain unit. After the calculation process is completed, the damage process of the unit by external forces can be completed. It is the calculation process within the same size interval. Figure 4 shows the calculation process of the gradual destruction of each unit in an interval under the action of external forces. Repeating this calculation process will reduce the scope of calculation until the model we calculate is the same as the phenomenon encountered in the actual experiment or construction.Figure 3
Calculation of different scales within specified scope.Figure 4
Cell rupture and gradual scale computing block diagram update.


### 3.2. Strain Strength Distribution Criterion

Intensity of lower part of the first damage, which has the strong Coulomb friction law of high strain part to keep the material properties, in a unit area can be divided into two parts, namely, the fracture zone and continuous zone; with the increase of applied force, the strength of the lower half damage will happen and can produce destruction through failure surface gradually. The cell model has been completely destroyed, strain intensity distribution model is given as follows:(3)Fn=−αKnΔun,Δun>0,−KnΔun,Δun≤0,Fs=−αβFs′Fs′Fs′1−αβFntanφ,Δun≤0,αβFs′,Δun>0.There is a certain correlation between the tensile and shear integrity factorsα and β and the tensile strain ε_min and shear strain γ shown at the interface, as shown in equations (4) and (5). ε=Δun/L,γ=Δus/L,εmax, and γmax are the minimum strain values when the model material is damaged by tensile stress and shear stress, and εmin and γmin are the tensile strain and shear strain of the model material under the action of external force. Fε,Fγ are the probabilities of intact models directly related to strain, which are between 0 and 1. At the same time, Fε,Fγ can obey the more mainstream probabilistic models.(4)α=1,ε<εmax,Fε,εmin≤ε≤εmax,0,ε>εmax,(5)β=1,γ<γmin,Fγ,γmin≤γ≤γmax,0,γ>γmax.


## 3.1. Multiscale Computing Framework

The choice of scale requires the following assumptions:(1)
Consistent with the principles of finite element method and related numerical method based on continuous model, the continuous model is dispersed in the form of grouping or grid, and the properties of the obtained material are determined by the properties of each internal unit, so it has the property of homogenization in a certain field. If the cracks in the unit cannot be ignored, the unit size cannot reach the standard, small size and high precision. When the size is refined and the stress field and displacement field in the field change little, it is the maximum allowable scale of the element.(2)
The damage in the nonuniform block is in a layered form, and the failure surface is in the middle of the failure layer. It has the property of equal strain, and the scale of failure surface is consistent with that of element. When the unit appears damaged, it represents the damage of the damage layer, which is its neutral surface.(3)
The thickness of the failure layer is confirmed according to the size of the tiny cracks, which may be the thickness of the shear band, and its length is higher than the “spacing” grade of the tiny cracks.(4)
In the failure surface, the material strength is not uniform and equal, and it has a certain distribution law. The failure surface is divided into many elements. According to the distribution law, the strength of different elements will be different. In the case of failure, because the failure surface has the characteristic of equal strain, the microelements with different strength will have different strains when the failure occurs, which leads to the failure surface showing a progressive form.(5)
The microelement strain strength in the failure surface is closely related to the material structure of small scale, which can be obtained by using the determined displacement boundary conditions. If the microunits can be represented by the lowest scale particles of a given material property, all multiscale calculations can be performed. If a smaller scale calculation is required, it is possible to redistribute the equal stress according to hypothesis (1) and eventually obtain smaller infinitesimals. Repeat the calculation until the unit scale of the given material properties is obtained.In exposition (1) and (2), the parameters of the material itself and some parameters of the corresponding failure surface are given in the same unit. Assumptions (3) and (4) refer to the comparison of low-level unit size of material properties and failure parameters.Figure3 explains the calculation problems of different size regions in the case of a certain unit. After the calculation process is completed, the damage process of the unit by external forces can be completed. It is the calculation process within the same size interval. Figure 4 shows the calculation process of the gradual destruction of each unit in an interval under the action of external forces. Repeating this calculation process will reduce the scope of calculation until the model we calculate is the same as the phenomenon encountered in the actual experiment or construction.Figure 3
Calculation of different scales within specified scope.Figure 4
Cell rupture and gradual scale computing block diagram update.


## 3.2. Strain Strength Distribution Criterion

Intensity of lower part of the first damage, which has the strong Coulomb friction law of high strain part to keep the material properties, in a unit area can be divided into two parts, namely, the fracture zone and continuous zone; with the increase of applied force, the strength of the lower half damage will happen and can produce destruction through failure surface gradually. The cell model has been completely destroyed, strain intensity distribution model is given as follows:(3)Fn=−αKnΔun,Δun>0,−KnΔun,Δun≤0,Fs=−αβFs′Fs′Fs′1−αβFntanφ,Δun≤0,αβFs′,Δun>0.There is a certain correlation between the tensile and shear integrity factorsα and β and the tensile strain ε_min and shear strain γ shown at the interface, as shown in equations (4) and (5). ε=Δun/L,γ=Δus/L,εmax, and γmax are the minimum strain values when the model material is damaged by tensile stress and shear stress, and εmin and γmin are the tensile strain and shear strain of the model material under the action of external force. Fε,Fγ are the probabilities of intact models directly related to strain, which are between 0 and 1. At the same time, Fε,Fγ can obey the more mainstream probabilistic models.(4)α=1,ε<εmax,Fε,εmin≤ε≤εmax,0,ε>εmax,(5)β=1,γ<γmin,Fγ,γmin≤γ≤γmax,0,γ>γmax.


## 4. Using the Monte Carlo Method to Produce Cracks

The Monte Carlo model, which conforms to the distribution law of geometric parameters of rock mass structural plane, is generated by random sampling with computer. The number of fractures in rock mass usually follows Poisson distribution. The position of the crack center points obeys the uniform distribution in the research area. The joint occurrence (strike and dip angle) usually follows distribution, univariate or bivariate distribution, bivariate normal distribution or lognormal distribution, uniform distribution, etc. The trace length and crack width of the structural plane obey the negative exponential distribution or lognormal distribution. A discontinuity network that is statistically completely equivalent to the real rock mass structure is generated. Figure5 shows the establishment of the crack model.Figure 5
Crack generation process.When the roadway is built in the mined geological body, the fissures will be distributed around the roadway due to the influence of mining. Based on the above uniform density function, plane fissures will be generated in the coal-rock mass through the secondary development of software combined with the strain strength distribution criterion. The distribution characteristic values of fractures in the model are shown in Table2. Figure 6 shows the fracture distribution with uniform distribution and different dip angles.Table 2
Geometric parameters of model.
Fissure angle (°)Crack length (m)Fracture spacing (m)Fracture slip (m)DistributionMeanVarianceDistributionMeanVarianceDistributionMeanVarianceDistributionMeanVarianceUniform—3Uniform0.80.3Uniform0.50.1Uniform0.070.1Figure 6
Fracture distributions for different dip angles in uniform distribution mode: (a) fissure angleα = 0°; (b) fissure angle α = 45°; (c) fissure angle α = 90°.
(a)(b)(c)


## 5. Numerical Simulation Analysis

To simulate the general situation of engineering geology, prestress loading of all complete models isPy/5, Px/5. When the model is stable, we excavate it in the form of command flow to produce mining roadway and adjust the stress to Py and Px to simulate the rock burst disaster of the roadway.WhenPx=10MPa and Py=20MPa with the change of crack inclination angle: α=0°, α=45°, α=90°, the collapse mode of roadway is shown in Figure 7.As can be seen from Figure7, when subjected to the stress of rock burst, geological bodies with horizontal fractures have priority to produce horizontal fractures in the coal and rock mass, while when vertical fractures have geological bodies, vertical fractures have priority to develop, and so does inclined fractures, which lays a prerequisite for the different failure modes of the final instability of the roadway. The high stress area is mainly concentrated in the tunnel roof, roof, and the junction of the three conditions, and there are differences between the two types of failure. The failure degree of the large dip angle model is 90°, and that of the small dip angle model is 0°. The failure angle of the 0° tunnel roof is perpendicular to the stress, and its compressive strength is greater than that of the 90° tunnel roof dip angle model. There is a significant difference in the fall of the three tunnel models, indicating that when the inclination angle of the crack is consistent with the direction of the maximum load stress, the failure degree of the tunnel is the largest, while the failure degree is smaller when the crack is perpendicular to the direction of the stress. Figure 7 shows the influence of dip angle on the failure mode of roadway.Figure 7
Fracture mode of roadway in case ofPy = 20 MPa and Px = 10 MPa.
(a)(b)(c)When the horizontal tectonic stress is greater than the vertical stress, the failure mainly occurs on both sides of the roadway. The reason why the failure degree of the left and right sides is different is that there are more cracks on the left side than on the right side when the cracks are randomly distributed. The horizontal tectonic stress is in the same direction as the horizontal tectonic stress, so its failure mainly occurs on both sides of the roadway, cracks, and horizontal tectonic stress. When the roadway failure mainly occurs at an angle of 45° on both sides of the roadway, and it has a certain impact on the roof and floor. When the crack is under vertical and horizontal tectonic stress, it is 90°. Roadway damage mainly concentrates on both sides of roadway, but it has great influence on roof and floor. When the fracture distribution is consistent with the direction of the larger stress, the surrounding rock exhibits a weaker compressive strength, while when the fracture distribution is perpendicular to the direction of the larger principal stress, the compressive strength of the roadway is the largest, and the failure range is different to some extent. This is shown in Figure8.Figure 8
Failure mode of roadway in case ofPy = 10 MPa and Px = 20 MPa.
(a)(b)(c)The crack has a certain influence on the phenomenon of stress concentration, which is mainly reflected in the fracture scope of the roadway when there is no crack; the large carrier plate and two sides of the roadway are penetrated by the crack and collapse.The cross cracks in the rock mass are more in line with the actual working conditions. Under the influence of mining, cracks will be distributed around the roadway. Based on the uniform density function, the cross cracks in the coal and rock mass are generated. See Figure9.Figure 9
Fracture distributions for cross dip in uniform distribution mode.It can be observed from the displacement cloud image that the failure and collapse above the surrounding rock mass are dominant, and the model is close to the collapse without cracks. When the junction of cracks is distributed in the surrounding rock, the longitudinal influence range of the three zones of straddle falls is larger, and it can be found from Figure10 that the fracture zone moves upward.Figure 10
Failure mode of roadway in case ofPy = 20 MPa and Px = 10 MPa.The distribution of cracks is not uniform, so there is a certain difference between the left and right side of the damage; the right side is pulled off in the middle of the pressure, and the upper side is unstable and collapsed, while the left side is shown as a whole failure; and the whole surrounding rock is shown as a lateral failure range and is larger than the longitudinal failure range. The effect of fracture on stress concentration is limited at the initial stage, but it has a great effect on the range of tunnel failure at the later stage. See Figure11.Figure 11
Failure mode of roadway in case ofPy = 10 MPa and Px = 20 MPa.


## 6. Conclusion

(1)
The Monte Carlo model is used as a means to simulate different types of crack models on a two-dimensional plane, and the influence of crack angle on the deformation of roadway surrounding rock is discussed. The damage mode at the initial stage of roadway failure is basically the same as that of roadway without crack. In the tunnel model with cracks, part of the strain energy is used for the development and cracking of microcracks, and the damage range is smaller than that of the tunnel model without cracks. As the failure is most complete when the angle of the structural plane is consistent with the maximum load stress, it becomes more difficult to rebuild the emergency rescue passage in the collapsed roadway.(2)
The law of roadway caving deformation under different roadway sections is studied. Moreover, the failure mode and stress of roadway, as well as the influence of roadway section shape on roadway caving morphology, are calculated. Different underground rescue schemes should be made according to different tunnel shapes.(3)
When the fracture distribution is consistent with the direction of greater stress, the surrounding rock exhibits a weaker compressive strength, while when the fracture distribution is perpendicular to the direction of greater principal stress, the compressive strength of the roadway is the largest. The cross fissure in the surrounding rock increases the longitudinal influence range of the “three zones” of the roof, makes the bending zone move upward, and increases the influence range of the fracture zone and the caving zone. It should be considered as an important factor when making underground emergency rescue plan.



---

*Source: 3759301-2021-08-30.xml*