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

The dip angle, length, spacing, and fracture distance of rock ﬁssure aﬀect 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 diﬀerent types of crack models in two-dimensional plane and generate diﬀerent crack models. The eﬀects of crack angle, crack length, fracture distance, and spacing on the deformation of surrounding rock are analyzed. The inﬂuence of diﬀerent rock burst on the failure strap-fall modes of ﬁssure roadway and roadway in diﬀerent sections is analyzed, and the stability law of roadway is studied. Under the condition of high stress, the roadway shape has little inﬂuence on the distribution of the principal stress diﬀerence of surrounding rock, but the equivalent excavation radius determines the distribution of the plastic zone of surrounding rock. The larger the ineﬀective 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 ﬁlls the roadway and the roof caving range becomes smaller.


Introduction
Research on the development of joints and fractures is the basis of analyzing the impact of collapse roadway on emergency rescue. e 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][9][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][18][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. e discrete element method is based on the continuum. is 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.

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.

Establishment of Numerical Simulation
Model. e geometric model is established as shown in Figure 1 e 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. e roadway is built in the coal seam, and the roadway parameters are 2 m high and 2 m wide. e specific parameters for solving the model are mainly set in Table 1. Figure 1 shows the roadway stress model. e numerical model of 8 m * 8 m is more consistent with the real tunnel shape. Figure 2 is a schematic diagram of contact pairs. e normal and tangential spring forces of the semispring are calculated as follows: 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, where T represents the tensile strength, φ the internal friction angle, and C the cohesion force.

Multiscale Computing
Framework. e 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) e 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) e 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. e 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) e 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. Figure 3 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.

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. e cell model has been completely destroyed, strain intensity distribution model is given as follows: ere is a certain correlation between the tensile and shear integrity factors α and β and the tensile strain ε_min and shear strain c shown at the interface, as shown in equations (4) and (5). ε � (Δu n /L), c � (Δu s /L), ε max , and  c max are the minimum strain values when the model material is damaged by tensile stress and shear stress, and ε min and c min are the tensile strain and shear strain of the model material under the action of external force. F(ε), F(c) are the probabilities of intact models directly related to strain, which are between 0 and 1. At the same time, F(ε), F(c) can obey the more mainstream probabilistic models.

Using the Monte Carlo Method to
Produce Cracks e Monte Carlo model, which conforms to the distribution law of geometric parameters of rock mass structural plane, is generated by random sampling with computer. e number of fractures in rock mass usually follows Poisson distribution. e position of the crack center points obeys the uniform distribution in the research area. e joint occurrence (strike and dip angle) usually follows distribution, univariate or bivariate distribution, bivariate normal distribution or lognormal distribution, uniform distribution, etc. e 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. Figure 5 shows the establishment of the crack model.  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. e distribution characteristic values of fractures in the model are shown in Table 2. Figure 6 shows the fracture distribution with uniform distribution and different dip angles.

Numerical Simulation Analysis
To simulate the general situation of engineering geology, prestress loading of all complete models is P y /5, P x /5. When the model is stable, we excavate it in the form of command flow to produce mining roadway and adjust the stress to P y and P x to simulate the rock burst disaster of the roadway.
As can be seen from Figure 7, 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. e 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. e failure degree of the large dip angle model is 90°, and that of the small dip angle model is 0°. e 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. ere 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.
When the horizontal tectonic stress is greater than the vertical stress, the failure mainly occurs on both sides of the roadway. e 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.
e 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. is is shown in Figure 8. e 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. e 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 Figure 9.     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 Figure 10 that the fracture zone moves upward. e 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. e 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 Figure 11.

Conclusion
(1) e 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. e 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) e 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. e 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.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare no conflicts of interest.