Abstract
Accurate simulation of the failure process of hard brittle surrounding rockmass is very important for the analysis and control of the structural stability in deep underground engineering. In order to simulate the progressive failure process of the hard brittle surrounding rockmass, a continuous discontinuous deformation analysis method that couples the finite element and discrete element is adopted. Taking the URL test tunnel in Canada as an engineering case, the constitutive model of the contact considering the effects of cohesion weakening and friction strengthening is applied, and the 2D approximation to 3D excavation by applying elastic modulus reduction technology is adopted to simulate the range and depth of crack growth of the surrounding rockmass. Then, the comparison between simulated results and on-site monitoring results is performed, which shows good consistency. At the same time, the key factors in the numerical simulation of progressive failure in hard brittle rockmass are identified, including the number of elements, excavation effects, and constitutive models. The results show that the constitutive model determines the basic form of crack propagation, but in order to accurately simulate the progressive propagation of cracks, the number of elements must be sufficient enough and the effects of 3D excavation must be considered. The analysis accurately simulates the progressive failure characteristics of hard brittle surrounding rockmass under high stress, achieving the purpose of reasonably grasping the degree of damage to the surrounding rockmass, and provides technical reference and support on how to accurately simulate the failure of hard brittle surrounding rockmass using the finite discrete element method.
1. Introduction
Under the conditions of deep burial and high stress, the hard brittle rockmass exhibits completely different mechanical behaviors from those under shallow burial and low stress. Spalling and rock burst are the main failure modes of this type of rockmass. Studies have shown that, during the excavation of deep buried tunnels, the stress of surrounding rockmass redistributes and small cracks appear around the tunnel. These cracks continue to expand and penetrate each other under the action of external forces and eventually cause part of the rockmass to separate from the matrix rockmass, forming an excavation damage zone [1]. Tracking this process is helpful to understand the failure mechanism of hard brittle rockmass, so as to formulate corresponding control measures.
In recent years, numerical simulation technology has developed rapidly. When on-site monitoring or in-site observation is insufficient, numerical simulation methods can effectively track the failure process of surrounding rockmass. The methods are mainly divided into the continuous analysis method and discontinuous analysis method. For the failure mode and mechanism study of hard brittle rockmass under high stress, the most representative one is related work carried out by the URL test tunnel of Atomic Energy Canada Limited (AECL). The V-shaped spalling zone formed during the excavation process is a typical brittle failure of surrounding rockmass under high stress. Martin [2, 3] studied the deformation and failure mechanism of granite through conventional indoor triaxial loading and unloading tests. On the basis of Martin’s test, Hajiabdolmajid et al. [4, 5] put forward the CWFS (cohesion weakening and friction strengthening) constitutive model by deeply analyzing the evolution law of strength parameters of the rockmass with damage development. This constitutive model can well describe the brittle failure behavior of granite, and the V-shaped spalling zone is well reproduced through numerical simulation. Diederichs et al. [6, 7] proposed the DISL (damage initiation and spalling limit) model based on the Hoek–Brown constitutive model and also reproduced the formation of the V-shaped failure zone. Many scholars have also made many contributions in simulating hard brittle rock cracking with discontinuous analysis methods such as PFC3D [8–10] and UDEC [11–13], not elaborated one by one here.
In actual situations, we not only focus on the continuous and discontinuous mechanical properties of geological bodies but also often need to obtain the temporal and spatial transition process from the continuous state to discontinuous state. In order to solve this problem, the finite discrete element method that combines continuous and discontinuous analysis methods has been gradually developed to simulate the complex interactive failure process of deformed bodies [14–17]. At present, the method has been successfully used in the simulation of slopes [18, 19], tunnels [20–23], etc. Among them, Vazaios et al. [23] used the self-developed program to realize the progressive failure of the URL test tunnel and compared it with calculated results of the continuous analysis method. In addition, some scholars have combined extended finite element and discrete element [24], peridynamics, and finite element [25–27] for similar research, and the essence is also the coupling of the finite element method and other failure analysis method.
In this paper, the continuum-discontinuum numerical analysis method CDEM is used and the URL test tunnel is taken as an example to analyze the fracture response of the surrounding rockmass and to compare it with on-site monitoring, focusing on the progressive propagation of cracks. In the calculation process, the key effects of number of elements, excavation effects, and constitutive model on simulating crack formation are analyzed, which provide reference on how to correctly simulate the crack propagation of hard and brittle rockmass under high stress.
2. Overview of the URL Test Tunnel
2.1. Rock Characteristics and Mechanical Properties
The URL test tunnel of Canadian Atomic Energy Limited (ACEL) is located about 120 kilometers northeast of Winnipeg, Manitoba, Canada. The lithology is Lac du Bonnet granite. The circular Mine_by test tunnel is located at a depth of 420 m, which was built in 1989 to study the failure process induced by excavation disturbance and the progressive failure process of surrounding rockmass under high stress. The granite near this depth can be classified as hard, massive, and brittle rockmass, which can be assumed to be homogeneous and isotropic. According to indoor test results, the mechanical parameters of surrounding rockmass are shown in Table 1, which will be used to simulate the failure process of the circular test tunnel.
2.2. Failure Characteristics
The circular test tunnel has a length of 46 m and a diameter of 3.5 m. It is excavated by the technology of nonblasting hydraulic rock cracking. The excavation footage is 1 m. The axis of the tunnel is approximately parallel to the direction of the intermediate principal stress, making the ratio of the maximum principal stress to the minimum principal stress as the largest in the plane perpendicular to the tunnel axis, which is beneficial to promote the failure of the rockmass. The displacement, strain, and acoustic emission/microseismic (AE/MS) of the surrounding rockmass are monitored by advanced instruments and equipment.
During the excavation process, flake damage can be observed at the top and bottom of the test tunnel. With the advancement of the tunnel face, the damage develops radially. The spalling process mainly occurs in the range of about 2 times advancing diameters, and finally, a V-shaped notch is formed. The radial depth of the damage zone (measured from the center of the tunnel) is generally 1.3 to 1.5 times the tunnel radius, that is, the depth of the damage zone is between 525 mm and 875 mm, and the range angle is about 70 degrees, as shown in Figure 1.
According to the description of the failure area of the surrounding rockmass by Read [28], the damage is mainly concentrated in the notch area, that is to say, the rockmass outside the notch basically has no serious defects. This is because the progressive damage of the rockmass at the tunnel boundary will lead to spalling and then form a damage zone. With the gradual occurrence of spalling, due to the accumulation of constraints and the restriction of the end of the notch to the rockmass, no more plastic strain or damage will occur in the notch area, and the tunnel will gradually stabilize.
3. Finite Discrete Element Simulation
In order to obtain a better simulation consistent with on-site monitoring, multiple trial and error calculations have been repeated. The paper is introduced based on the final calculations, and then, the effects of some factors such as the number of elements, excavation effects, and constitutive models are discussed, which provides suggestions and directions for simulating the progressive failure of hard brittle surrounding rockmass using the finite discrete element method.
3.1. Method Introduction
The GDEM stress analysis system [29] is a high-performance finite element-discrete element calculation software based on CDEM (continuum discontinuum element method). This method couples the finite element and discrete element, performs finite element calculation inside the block, and performs discrete element calculation at the block and element interface. By introducing a breakable one-dimensional spring, as shown in Figure 2, and through the rupture of the block interior and the element interface, it cannot only simulate the deformation and movement characteristics of the material in continuous and discontinuous state but also realize the simulation of the progressive failure process of the material from continuous to discontinuous.
The explicit iterative calculations are adopted based on the time-history-based dynamic relaxation technology in CDEM. Therefore, the stress of elements and spring force of nodes at each time step can be obtained, and the strength judgment can be made according to the rupture criterion. The basic iterative process is shown in Figure 3.
3.2. Model Description
3.2.1. Geometric Model
Referring to the paper by Vazaios et al. [23], the software Gmsh is adopted to obtain a two-dimensional random mesh. In order to eliminate the boundary effects, the simulation range is about 17 times the excavation diameter, with the section size 60 m × 60 m. The model is divided into four regions, as shown in Figure 4. Among them, region 1 is the circular tunnel to be excavated, the element size of which transitions from 0.5 m to 0.03 m. Region 2 is the square area adjacent to the circular tunnel where cracks first appear during the excavation process. So, the element size should be small enough. The area range is 14 m × 14 m, and the element size is 0.03 m. Due to the large range of the entire model, in order to avoid the loss of computational efficiency caused by too many elements, two transition regions are specially set. Among them, the range of region 3 is 25 m × 25 m, and the element size varies from 0.03 m to 1.5 m. The range of region 4 is 60 m × 60 m, and the element size varies from 1.5 m to 2.5 m. The total number of elements is 574,828.
3.2.2. Boundary Conditions and Geostress
The boundary conditions of the model include displacement boundary conditions and stress boundary conditions. The measured in situ stress results show that the maximum principal stress in this area is 60 MPa, the intermediate principal stress is 45 MPa, the minimum principal stress is 11 MPa, the maximum stress ratio of which reaches 6 : 1. Considering the directions and converting the geostress from the principal stress space to the Cartesian space, the magnitudes of the normal stress and shear stress are shown in Table 2.
And, because the tunnel is deeply buried, the gravity-induced stress is not considered. Regarding the far-field boundary conditions, the displacement is fixed at the bottom of the model, and the corresponding normal stress and shear stress are applied to the upper surface and the left and right sides, respectively.
3.2.3. Constitutive Model
In the traditional continuous analysis method, in order to consider the reduction of the postpeak strength of the rockmass, the strain softening model and the CWFS model are commonly used. Among them, the strain softening model claims that both the cohesive force and frictional strength constitute its peak strength before plastic deformation occurs, and then, both begin to lose at the same time and gradually decrease as the strain increases. The CWFS model believes that only the cohesion of the rockmass plays a role at the initial moment, and the cohesion gradually decreases as the damage develops, and then, the friction strength begins to play a role and gradually increases as the damage develops. Hajiabdolmajid et al. [4, 5], Wu [30], and Liu [31] reproduced the V-shaped damage of the URL test tunnel based on the CWFS method, verifying the significance of friction strength strengthening in the postpeak stage.
Inspired by CWFS in continuum, this feature of rockmass is transplanted to the interface of discontinuous elements. The elastic constitutive model is adopted inside the element, and the input parameters include density, elastic modulus, and Poisson’s ratio. The fracture energy model is used on the element interface, and the input parameters include normal stiffness, tangential stiffness, cohesion, internal friction angle, tensile strength, tensile fracture energy, and shear fracture energy. The incremental relationship of the elastic constitutive model iswhere and are the average stress increment and strain increment, is the volumetric strain increment, and and are the bulk modulus and shear modulus of the material, respectively.
The fracture energy model is essentially the maximum tensile stress model and the M-C model, both of which consider a linear softening effect of the tensile strength and shear strength, respectively. The incremental method is used to calculate the normal and tangential force of the next step on the interface:where and are the normal and tangential force, and are the normal stiffness and tangential stiffness, which are obtained by inheriting form the element stiffness, and and are the relative displacement increments in normal and tangential directions. Then, following formula (3) is used to judge the tensile failure, and the normal force is corrected. That is, the tensile strength linearly weakens when the normal force exceeds the tensile strength:where , , and are the tensile strength of the interface at the initial time, this time, and the next time respectively, is the normal relative displacement on the interface at this time, is the tensile fracture energy, and is the area of the interface.
At the same time, formula (4) is used to judge the shear failure and tangential force is corrected. That is, the shear strength linearly weakens when the tangential force exceeds the shear strength:where , , and are the cohesion of the interface at the initial time, this time, and the next time, respectively, is the tangential relative displacement on the interface at this time, and is the shear fracture energy.
In the fracture energy model, the shear strength is composed of both cohesive force and frictional strength. In order to simulate the effects of CWFS, a small friction strength is given at the initial moment. At this time, the friction strength contributes little to the shear strength. Then, when the shear failure criterion is reached, the cohesion starts to lose gradually. At this time, the friction strength of the interface is set to increase by 10 times. Thus, the goal of cohesive weakening and friction strengthening is reached.
3.2.4. Calculation Steps
The first step is to calculate the initial geostress field of the model under the given boundary conditions, 832,800 steps were performed, and the initial result was saved for the following calculation. The next step is to simulate the excavation of the test tunnel, but how to simulate the three-dimensional effect with the two-dimensional model is a very important question. Because the tunnel is generally excavated in stages, the resistance of the surrounding rockmass and the excavation boundary are formed gradually not instantaneously. While in the two-dimensional model, excavation is simulated by emptying the tunnel elements, which means that the tunnel is excavated completely at the same time, resulting in a significant difference from the real situation. So, some measures should be taken to try to eliminate such excavation effects. Vlachopoulos and Diederichs [32] introduced four methods of using two-dimensional analysis to approximate three-dimensional tunnel excavation behavior, as shown in Figure 5. Curran et al. [33] specifically introduced the fourth method, that is, the face replacement method in the support designing in weak rocks. In short, the progressive excavation process is simulated by setting the elastic modulus of the to-be-excavated part to decrease gradually. The face replacement method, (d) in Figure 5, will be used in this paper.
(a)
(b)
(c)
(d)
In this paper, the specific performed process is as follows. After the initial elastic calculation is completed, the elastic modulus of the to-be-excavated area will be reduced, with other parameters keeping unchanged. And, the parameters of other areas also remain unchanged; then, 5000 steps are performed to achieve a temporary balance. Repeat the process 5 times until the 3D excavation behavior is almost approximated. The specific settings of elastic modulus are shown in Table 3. Then, the excavation process is performed by emptying the tunnel elements, and 3000 steps are performed to observe the crack growth.
3.2.5. Input Parameters
Due to the lack of experience in the simulation of crack propagation during tunnel excavation, it is difficult to directly solve the problem using the discontinuous analysis method, and the debugging process takes a lot of time. Therefore, the strategy we adopted is to first obtain good results in the category of continuum, then migrate to discontinuous analysis by considering contact elements, and obtain the final result through continuous test calculations.
In the continuous analysis, the strain softening model is adopted. The material parameters adopted are the same with Table 1. The goal is to make the failure state of the surrounding rockmass coincide with the on-site monitoring results. The software Gmsh is used for meshing, and the model region is the same with Figure 4, but the number of elements is 58336. By adjusting the value of tensile strength and dilatancy angle, the results of several schemes are shown in Figure 6. It is indicated that the increase in tensile strength can significantly reduce the number of elements in the tensile failure, and the increase in the dilatancy angle can make the V-shaped notch failure more obvious, and the range and depth of the damage zone are enlarged.
(a)
(b)
In the end, it is confirmed that the tensile strength of the surrounding rockmass is 10 MPa and the dilatancy angle is 30 degrees, which will be used in following calculations. The final input parameters in the model are shown in Table 4.
4. Simulated Brittle Response
4.1. Simulated Initial Stress Field
Firstly, the geostress field after calculating 832,800 steps is given in Figure 7, taking the first principal stress for example. As shown in the legend, the calculated stress of the whole model is between 53 MPa and 62 MPa, which is almost consistent with the nominal value 60 MPa. The direction of the first principal stress is indicated by the short red line, coinciding with the true direction of the first principal stress. The inclination angle is about 11 degrees. The accurate simulation of the geostress field guarantees the subsequent excavation calculation to be performed successfully.
4.2. Simulated Brittle Failure
After the 3D excavation approximation by repeating the elastic modulus reduction 5 times, the tunnel elements were removed completely, and the evolution process of the cracks of the rockmass is shown in Figure 8. Picture (a) shows the crack state at step 857,900 when there are only small cracks on the top and bottom of the tunnel. Picture (b) shows the crack state at step 858,000. At this time, the range of cracks at the top and bottom expands a little. And, at the same time, small tensile cracks appear in the vertical direction on the left and right sides of the tunnel. Picture (c) shows the crack state at step 858,900. The cracks at the top and bottom gradually extend to the deep part of the tunnel, and the tensile cracks on both sides also extend to the depth of the surrounding rockmass. The area of the tensile cracks has not expanded, being still a long and thin crack. Besides, at this time, the prototype of V-shaped failure has been basically formed. Pictures (d)∼(f) show the crack state at 859,400 steps, 859,900 steps, and 860900 steps separately. During these three stages, the cracks at the top and bottom expand to the depth further, the form of the V-shaped failure is more obvious, and the tensile cracks on both sides do not expand in a large area. After calculating 3000 steps, the crack growth has reached a relatively stable state.
(a)
(b)
(c)
(d)
(e)
(f)
The comparison between the failure state of the surrounding rockmass and the on-site monitoring results, when the crack no longer continues to grow, is further performed, as shown in Figure 9. The orientation of the top and bottom cracks is perpendicular to the direction of the first principal stress, and the direction of tensile cracks on both sides coincides with the first principal stress. After measurement, the range and depth of the V-shaped notch are in good agreement with the on-site monitoring results. The angle of the damage zone measured from the center of the tunnel is about 70 degrees. The depth of the relatively broken zone at the top is about 0.3 times the diameter, and the total depth of the cracks is about 0.5 times the diameter. The damage range and depth at the bottom of the tunnel in numerical calculation are similar to those obtained at the top. But in on-site monitoring, the actual damage zone measured at the bottom of the tunnel is smaller than the calculated value. There are probably two reasons, one is that the surrounding rockmass at the top falls off and may cover the ground; the other is that, under the combined function of self-weight stress and tectonic stress, the crack at the bottom is actually weaker than that at the top. But in the calculation, self-weight stress is ignored but symmetrical stress is adopted so that the same calculated results at the top and bottom are obtained.
In addition, compare the acoustic emission signals and the microseismic location around the tunnel with the crack propagation, as shown in Figure 9. It can be seen that, within the area where cracks appear, relatively dense acoustic emission signals and microseismic signals have been monitored, indicating that there is indeed crack appearing and energy releasing, which confirmed the correctness of the calculation results once again.
5. Key Factors Identification
This section mainly focuses on the key factors in accurately simulating the progressive propagation of cracks, which are all concluded from numbers of trials.
5.1. Effects of Number of Elements
Generally speaking, the element size has a great influence on numerical simulations [34]. In this calculation, the mesh and material parameters of the continuous model are used, and the contact is applied to the boundary of the elements to carry out continuous discontinuous simulation. The constitutive model of the contact element adopts cohesion weakening and friction strengthening, and the excavation approximation by elastic modulus weakening is considered. Therefore, compared with the final calculation scheme, there is no other variable but the number of elements.
The crack propagation process during 3000 calculation steps is shown in Figure 10, and three results at corresponding steps are chosen and analyzed. Compared with Figure 8, the location of the cracks around the tunnel is roughly right, but the range and depth of cracks extend more widely and deep, especially on the top and bottom. The cracks extend a long distance to both sides, and the V-shaped notch on the top is not completely formed. It can be inferred that more elements are needed to obtain ideal simulation effects.
(a)
(b)
(c)
5.2. Effects of Excavation Method
It can be seen from Section 5.1 that when there are not enough elements, the range and depth of the cracks expand widely although the crack location is roughly the same. The number of elements should be adjusted. At the same time, referring to the literature of Vazaios [23], the number of elements is finally added to 574828. The constitutive model of the contact considers cohesion weakening and friction strengthening, and the excavation effect is not considered but direct excavation is used for this calculation. In this way, compared with the final model, the only variable is the excavation method.
The crack propagation process during 3000 calculation steps is shown in Figure 11. Compared with Figure 8, cracks first appeared in the four corners of the tunnel instead of the top of the tunnel. As the calculation progressed, cracks began to appear on the top and bottom, but the cracks in the four corners still exist and get widely. When the crack propagation is almost in the stable state, the V-shaped notch has been roughly formed. At this time, the cracks on the top and bottom intersect with the cracks in the four corners. Compared with Figure 10, the increase of elements restricts the cracks from propagating farther and deeper, but new problems have appeared due to ignoring the excavation effects, which need to be further optimized and solved.
(a)
(b)
(c)
5.3. Effects of Constitutive Model
Finally, the influence of the constitutive model is compared. The material parameters used is obtained from the continuous analysis in Section 5.1. The number of elements is 573288, and the constitutive model of the contact element does not consider cohesion weakening and friction strengthening, but the excavation effect is considered. The only variable compared with the final model is the constitutive model.
The crack propagation process during 3000 calculation steps is shown in Figure 12. Compared with Figure 8, cracks initially appeared on the top and bottom of the tunnel. As the calculation progresses, the range of the cracks at the top and bottom expands rapidly and expands deeper and farther. After 3000 steps is finished, the crack expanding is still continuing, and the range is further expanded. Until this time, no V-shaped notch was formed on the roof of the tunnel and the cracks look messy. The reason could be that when the stress state of the element reaches the rupture criterion, the cohesive force begins to decrease, but the friction strength does not increase, resulting in the load-bearing capacity of the element not meeting the requirements; then, the rapid crack propagation happens. Therefore, choosing a reasonable constitutive model is very critical.
(a)
(b)
(c)
6. Conclusions
The brittle failure of hard brittle surrounding rockmass has always been a problem in engineering and academia. The coupling analysis method of the finite discrete element has therefore been developed rapidly, but how to accurately simulate the progressive propagation of cracks still needs to be studied in depth. In this paper, the URL test tunnel in Canada is taken as an example, the number of elements, excavation effects, and constitutive models are discussed through continuous adjustment of models and trial calculations, and finally, the result that is in good agreement with the on-site monitoring is obtained, and the following are summarized in conclusion:(1)The constitutive model of the contact element is the key factor that determines the basic form of crack propagation, which is generally obtained through laboratory tests. This paper draws on the experience summarized by the predecessors, applies the model considering cohesion weakening and friction strengthening to the contact, and achieves good results.(2)When dealing with crack propagation of hard brittle rockmass, the number of elements must be sufficient enough. Due to the element boundary being the contact surface, more elements can make the expansion of the stress path more accurate.(3)When using two-dimensional simulation to approximate three-dimensional excavation, the actual excavation method must be considered. Then, certain measures must be taken to simulate the real excavation situation as much as possible.
In addition, in this paper, the URL test tunnel is specified. When migrating to other similar projects with hard brittle rockmass, such experience can be copied, but it should be adjusted according to the specific situation. For example, the constitutive equation changes with different kinds of rocks or other methods need to be used to approximate the three-dimensional excavation effects.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The work was supported by China Postdoctoral Science Foundation with Grant no. 2021M690999.