Abstract
The analysis and simulation of brittle failure of hard rock tunnels under high stress is essential for understanding and mastering the brittle failure characteristics of rock masses and for analyzing and regulating the stability of surrounding rocks in underground projects. The rupture and deformation of hard brittle basalt in the dam site area is one of the key problems faced by Baihetan Hydropower Station. In the paper, the spalling characteristics of the surrounding rock on the right bank of the exploratory tunnel are summarized. CDEM (continuum-based discrete element method) is used to carry out the numerical simulation and the fracture energy model considering cohesive force weakening and friction angle strengthening is adopted. The indoor test simulations are conducted first to verify the effectiveness of the model. Then the simulation of the right bank exploratory tunnel is conducted to study the brittle failure of the surrounding rock. The results are compared with the field exposure, and the stress and displacement of characteristic points of the surrounding rock are analyzed. The numerical results are in good agreement with the damage situation in the field and reflect the brittle failure characteristics of basalt under high-stress conditions, which helps to reasonably grasp the damage situation of the surrounding rock and take corresponding support measures, and also proves the superiority of CDEM method in solving hard rock fracture, providing a technical reference for similar hard rock brittle failure problems in engineering.
1. Introduction
Under high-stress conditions, the excavation of hard rock tunnels will cause rapid adjustment of secondary stress field of the surrounding rock, which may lead to damage destabilization phenomena such as local spalling, collapse, or even rock burst of surrounding rock, seriously threatening the project progress and personnel safety [1–3]. Only by fully understanding the damage mode of surrounding rock under different conditions and accurately grasping its damage mechanism, can corresponding control measures be formulated to ensure engineering safety. Therefore, it is necessary to gain better insight into the failure characteristics and mechanisms of hard rock tunnels [4, 5].
In recent years, many scholars have performed uniaxial compression [6], biaxial compression [7], and true-triaxial unloading tests [8, 9] on rock or rocklike specimens to simulate the spalling failure of the tunnel, but these studies focused on the material failure and did not consider the effects of the tunnel space structure. To analyze the influence of tunnel space structure on hard rock fracture, some scholars have employed rock specimens with a central hole to simulate the failure behavior of tunnels at a laboratory scale and focused on the crack propagation around the holes [10–12].
Because of the difficulties in reproducing the failure characteristics during laboratory tests and the rapid development of computational capacity, numerical methods have been increasingly used to simulate the failure characteristics of hard rock tunnels. The most representative work is the research related to the damage of the Canadian underground test tunnel. Martin firstly studied the cohesion loss and stress path effect of hard rock rupture in 1997 [13], based on which Hajiabdolmajid et al. [14, 15] proposed the CWFS (cohesion weakening and friction strengthening) model, Diederichs et al. [16, 17] proposed DISL (damage initiation and spalling limit) model to reproduce the formation of V-shaped damage zone. With reference to the CWFS model, Jiang et al. [18] proposed the hard rock intrinsic model, considering the effect of surrounding rock deterioration. All of the above are developed in the context of continuous media. However, because the displacement of adjacent interfaces must be coordinated in FEM, it is challenging to reach convergence in simulations. Some scholars adopted discrete-based methods such as PFC2D [19], DEM [20, 21], and DDA [22] to study related fracture issues in hard rock. However, the entire failure process from continuum to discontinuum in these studies was not captured. So, the combined finite discrete element method has been gradually confirmed and adopted by scholars to solve related problems [23–26]. For example, Zhong et al. [23] and Vazaios et al. [24] adopted the finite discrete element method to reproduce the V-shaped failure of the URL test tunnel.
Baihetan hydropower station in China is the largest power station under construction, and the underground powerhouse area is a high ground stress area with four types of basalt erupted at the dam site with internal crypto microfracture developed and elastic-brittle damage characteristics, determining that Baihetan is facing the critical rupture deformation problem of hard, brittle surrounding rocks [27]. Based on this, many scholars have carried out indoor uniaxial and triaxial tests to study the basic mechanical characteristics of basalt [28, 29]. Besides, during the excavation of the powerhouse, many typical hard rock damage phenomena such as spalling, rupture, and spray layer cracking have occurred [30, 31]. So, it is necessary to study the fracture deformation of hard, brittle basalt surrounding rocks.
In order to deeply study the brittle failure characteristics of basalt under high stress, several exploratory tunnels were arranged on each side of the left and right banks of the dam site area. In this paper, we first introduce the exploratory tunnel arrangement and the damage of the surrounding rock in exploratory tunnels. Then select the brittle intrinsic model and use CDEM [32] to compare indoor tests and numerical simulations for model validation. Lastly, simulate the exploratory tunnel excavation, analyze the brittle failure process of the surrounding rock, compare it with the field exposure, and analyze the stress and displacement of the surrounding rock in order to better recognize the brittle failure characteristics of basalt under high-stress conditions, and at the same time provide a reference for similar damage cases encountered in engineering.
2. Overview of Baihetan Underground Exploration Tunnels
2.1. Background of the Project
Baihetan Hydropower Station is one of the four-step power stations in the lower reaches of the Jinsha River, which will become the second-largest giant hydropower station after the Three Gorges after completion. In order to analyze the distribution of ground stress on the left and right banks, several exploratory tunnels are laid on each bank, and their location distribution is shown in Figure 1. The right bank tunnels consist of PD62 and its branch tunnels, which are mainly divided into four different orientations on the plane, among which PD62-2 is axially N80°E, orthogonal to the axis of the main tunnel, while PD62-4 and PD62-3 are axially N10°W, in line with the direction of the main tunnel. The left bank tunnels consist of PD61 and its branch tunnels; the orientations are not discussed here.
2.2. Spalling Characteristics of the Right Bank Exploration Tunnel
This paper takes the exploratory tunnels on the right bank as a research object. A detailed investigation of the spalling characteristics, including spatial distribution and intensity level, of the five exploratory tunnels, has been obtained, as shown in Figure 2.
The spalling is the direct manifestation of hard, brittle basalt under the action of large initial in-situ stress. The overall spalling on the right bank shows two basic patterns as follows:(1)The relationship between the spalling intensity of the exploration tunnel is(2)The spalling in PD62-2 is located in the top arch of the exploratory tunnel, while the spalling in the rest of the exploratory tunnels is generally inclined to the side arch on the upstream side of the river valley, as shown in Figure 3. The maximum spalling depth in PD62-2, PD62, PD62-1, and PD62-4 is 0.7 m, 0.7 m, 0.5 m, and 0.7 m, respectively. Regarding damage morphology, the spalling in both PD62-2 and PD62 shows typical V-shaped features, with the tip of the V shape extruding a distinct rupture shape and a large rupture depth. The spalling in PD62-1 also shows a typical V shape in general, but the tip morphology is poorer and slightly flatter compared to PD62-2 and PD62. Influenced by the gently inclined structural surface, the spalling of PD62-4 is significantly different from the other three tunnels. Although the most intense spalling part in PD62-4 also has a greater rupture depth, the morphology is much more soothing than the above tunnels. So the deeper spalling in PD62-4 can be considered as a special case, and will not be considered in the paper.
3. Mechanical Model of Basalt and Model Validation
3.1. Mechanical Model
The strain-softening model is commonly used to represent the postpeak weakening characteristics of rock mass in conventional continuum analysis. The strain-softening model assumes that the cohesion and frictional strength of the rock mass constitute its peak strength before plastic deformation occurs, and then both start to lose simultaneously and gradually decrease with increasing strain. However, Martin [13] found in their study of the URL test tunnel that only the cohesive force comes into play at the initial moment and gradually decreases with the development of damage, while the frictional strength comes into play and gradually increases with the development of damage, i.e., the CWFS model. Hajiabdolmajid et al. [14, 15] reproduced the V-shaped damage of the URL test tunnel based on CWFS, verifying the accuracy and feasibility of CWFS in simulating the deformation damage of the granite in the URL tunnel.
After summarizing the results of cyclic loading and unloading tests of amygdaloidal basalt, Chang [28] obtained the evolution of cohesion and internal friction angle with internal variables by fitting plastic internal variables based on the M-C yield criterion. The fitting results confirm that basalt also has similar properties to CWFS. This property must be considered in the calculation to describe the brittle damage behavior of basalt.
In CDEM, the intrinsic model is implemented in the following way: the finite element part uses a linear elastic model, and the elastic intrinsic incremental relations are
, refer to the average stress increment and average strain increment. is the volumetric strain increment. and are the bulk modulus and shear modulus of the material, respectively.
The fracture energy model is used for the discrete element part, which is essentially the maximum tensile stress model, and the M-C model, considering the linear softening effect of strength and the strengthening effect of friction angle. Firstly, the normal and tangential contact forces at the next time step on the interface are calculated using the incremental method as follows:
, are the normal and tangential force; , are the normal stiffness and tangential stiffness, which are obtained by inheriting from the element stiffness. , are the relative displacement increments in normal and tangential directions. Subsequently, equation (4) is used for the determination of tensile damage and the correction of normal contact force and tensile strength.
, , 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.
Equation (5) is used for the determination of shear damage and the correction of tangential contact force and cohesive force and internal friction angle.
, , 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. is the internal friction angle on the interface at the initial moment; n is the multiplication of the increase of internal friction angle, the value of which depends on the specific analysis.
3.2. Model Validation
In order to verify whether the above model can reflect the characteristics of the mechanical behavior of basalt under different confining pressures, the triaxial test curves of basalt with confining pressures of 0, 10, 30, and 70 MPa are compared with the numerical simulation curves, as shown in Figure 4. Among them, when the confining pressure is between 0 MPa and 10 MPa, the compaction section of the test curve is significant. Regardless of this compression deformation, the numerical simulation is compared with indoor test results, as shown on the right side of Figures 4(a) and 4(b). It can be seen from Figure 5 that numerical simulation curves of the fracture energy model can match well with test curves, indicating that the model can better reflect the brittle fall process of basalt under different confining pressures.
(a)
(b)
(c)
(d)
4. Numerical Simulation of Exploration Tunnel
4.1. Geometric Model
In discontinuous analysis, the good or bad calculation results are closely related to the element size. Referring to the numerical simulation of the URL test tunnel in Canada by Vazaios [24], the grid is divided as shown in Figure 5.
According to the actual excavated dimensions of the exploration tunnel, the simulated sidewall height of the exploration tunnel is 1.8 m, the top height is 2.2 m, and the bottom side width is 2.5 m. The size of the whole model section is 40 m × 40 m. Area 1 is the cave-shaped part to be excavated off, and its grid size transitions from 0.2 m to 0.03 m. Area 2 is a square area adjacent to the cavity, where the cracks appear first during excavation, with a section of 10 m × 10 m, and a small enough grid size of 0.03 m to accurately capture the crack expansion. The section of area 3 is 20 m × 20 m, with the grid size from 0.03 m to 1.5 m; area 4 is 40 m × 40 m, with the grid size from 1.5 m to 2.5 m. The total number of elements is 299662.
4.2. Geostress and Boundary Conditions
Taking PD62-2 exploratory tunnel as the simulation object, the geostress direction is approximated as follows [33]: the maximum principal stress orientation is in the N-S direction, the dip angle is nearly horizontal, and the size is 26 MPa; the middle principal stress is nearly horizontal, and the size is 18 MPa; the minimum principal stress is vertical, and the size is 10 MPa.
The boundary conditions adopted in the calculation are: a normal displacement constraint is applied at the bottom, a surface force of 26 MPa is applied at the left and right sides, and a surface force of 10 MPa is applied at the top.
4.3. Calculation Steps and Mechanical Parameters
Before the start of excavation calculation, the initial stress field of the model needs to be obtained first. The whole calculation process lasts for 800000 steps. Then the excavation of the exploratory tunnel is carried out. In order to eliminate the 3D effects caused by the actual excavation and to ensure the accuracy of the 2D calculation results, Vlachopoulos and Diederichs [34] proposed four methods to achieve the approximation of the 2D simulation to the 3D excavation. Curran et al. [35] introduced the face replacement method in detail.
In this paper, the face replacement method will also be used to approximate the 3D excavation effect. In brief, it is to continuously weaken the deformation modulus of the element in area 1, and after each modulus replacement, 5000 steps are calculated to reach a temporary equilibrium, and so on for five times. Eventually, the elements inside the exploratory tunnel are excavated one at a time, and then 10,000 steps are calculated to observe the fracture evolution of the surrounding rocks.
After several trial calculations for parameter adjustment, the final mechanical parameters used in the numerical simulation are shown in Table 1. The deformation modulus replacement is shown in Table 2.
5. Results Analysis
5.1. Fracture Evolution of Surrounding Rock
Figure 6 shows the evolution process of internal cracks in PD62-2, obtained by numerical simulation. The initial time step of formal excavation is 825100, and the results are output every 1000 steps, for a total of 10 times.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
As can be seen from the figure, the cracks first appeared at the bottom corners of the left and right sidewalls. Then, with the expansion of cracks at the bottom, small cracks began to appear at the top of the exploratory tunnel. With the advancement of the calculation step, cracks at the bottom and the top arch gradually expanded and connected. When calculating 10000 steps, that is, the total number of calculation steps is 835100, the crack propagation of the bottom and top arch basically reaches a stable state. At this time, the top arch forms a much more significant V-shaped failure, and the rock mass is broken and spalled. The cracks on both sides of the bottom also gather and connect to the middle depth, forming the shape of V-shaped cracks, but the rock mass is significantly less broken than the top arch.
During the calculation, the damage of the one-dimensional spring, i.e., contact element, is monitored. Three variables are statistically output, as shown in Figure 7. Among them, the total damage of spring is the proportion of elements that have been damaged. The total rupture of a spring is the proportion of elements that have not only been damaged but also have had their cohesion and tensile strength reduced to zero. The current damage of spring is the proportion of elements that are now in a damaged state. It can be seen from the figure that after the tunnel starts excavation, all three curves increase continuously with the advance of calculation steps. Before the 831500 calculation steps, the failure speed of the spring was faster; after that, the failure speed of the spring gradually decreased, and finally showed a convergence trend. Among them, the current damage ratio of the spring almost stops increasing, which has good correspondence with Figure 7.
5.2. Comparison between Simulated and Observed Fracture
We compare the failure of the exploratory tunnel finally obtained by numerical calculation with the field observation, as shown in Figure 8. As for the top arch, the numerical results are in good agreement with the field observations. The maximum depth of crack extension is about 0.7 m, and the failure depth of the top of PD62-2 is also about 0.7 m. The crack range of the crown arch calculated by numerical calculation is slightly larger than the field observation results, but under the combined action of surrounding rock stress and self-weight stress, the surrounding rock spalling and falling are concentrated in the core area of the crown arch, which is consistent with the field observation. In the numerical calculation, the bottom plate of the exploration tunnel also forms a penetrating crack, and the final shape is roughly V-shaped. This is because the stress boundary conditions applied in this paper are generally symmetrical, so the calculated failure area is also generally symmetrical. However, unlike the top arch, the bottom plate does not form a significant crushing failure.
In general, the numerical results can better reflect the spalling failure of the surrounding rock, and the calculation results are reliable.
5.3. Stress Analysis of Surrounding Rock
In order to further analyze the calculated damage to the top arch and bottom floor of the surrounding rock, 6 measurement points were taken at different depths of the top arch and bottom floor, respectively. Among them, the distance from the top arch 1#∼6# to the tunnel top is 0 m, 0.3 m, 0.4 m, 0.5 m, 0.7 m, and 1 m, respectively; the distance from the bottom floor 1#∼6# to the tunnel bottom is 0 m, 0.5 m, 0.6 m, 0.9 m, 1.0 m, and 1.1 m, respectively. The evolution curves of the minimum principal stress of each measurement point at the top and bottom of the tunnel with the calculation step after tunnel excavation are counted. The analysis results are shown in Figure 9.
Firstly, for the top arch, the trend of the minimum principal stress at measurement points 1#∼5# is to increase gradually with the excavation, and then to decrease at a certain point, and the deeper the measurement point is buried, the less significant the drop. At measurement point 6#, the minimum principal stress no longer has a drop process, but always grows slowly with the iteration step. The crack expansion pattern at the moment of stress drop from 1# to 5# measurement points is plotted on the top of the curve, and it can be seen from the figure that the moment when the stress starts to drop is the moment when the crack starts to expand at that location, and the position of each measurement point in the figure coincides with the crack boundary. We take the 1# measurement point as an example, which is a point close to the tunnel top. When the tunnel is first excavated, the minimum principal stress at this point gradually increases. It can be seen from Figure 6 that there is no crack at the top of the tunnel at the initial time, and the stress is accumulating. When the calculation reaches 826400 steps, the minimum principal stress begins to fall rapidly, and cracks are produced here due to the rapid release of energy. The cracks continue to expand and connect in the subsequent process, and the minimum principal stress has been at a low level and no longer rises. 2#∼5# measurement points can be analyzed similarly. Finally, when the calculation stops at 835100 steps, that is, the crack is in a stable state, the minimum principal stress at each measurement point also reaches a stable state, and the crack extends to the vicinity of measurement point 5#, the buried depth of this point is about 0.7 m, which is consistent with the measured failure depth, further verifying the accuracy of the calculation.
The minimum principal stress of each measurement points on the bottom floor can be roughly divided into three cases: (1) the principal stress of measurement points 1# and 2# always decreases slowly; (2) the principal stress of measurement points 3# and 4# increases first and then decreases; and (3) the principal stress of measurement points 5# and 6# always increases slowly. Among them, due to the unloading effect of tunnel excavation, cracks start to appear at the bottom of the tunnel at the beginning, and initially penetrate at 829100 steps, forming a V-shaped failure. 1# and 2# measurement points are located on the inner side of the V-shaped failure, so the principal stress decreases continuously with the unloading. While 3# and 4# measurement points are located on the outer side of the V-shaped failure, the crack gradually expands deeper as the iteration proceeds until 830500 and 831100 steps, when 3# and 4# measurement points begin to show stress drop one after another. Since the maximum depth of crack propagation is near measurement point 4#, no cracks are generated at measurement points 5# and 6#, and their principal stresses keep growing slowly.
5.4. Deformation Analysis of Surrounding Rock
Similar to the stress analysis of surrounding rock, 6 measurement points are taken at different depths of the tunnel top and bottom to analyze the surrounding rock deformation. Among them, the distribution of measurement points from 1# to 6# in the top arch is the same as that in the stress analysis, and the distance of each point from 1#∼6# in the tunnel bottom is 0 m, 0.2 m, 0.4 m, 0.5 m, 0.8 m, and 0.9 m, respectively. The evolution curves of the vertical displacement of each measurement point at the top and bottom of the tunnel with the calculation step after tunnel excavation are counted. The analysis results are shown in Figure 10.
Firstly, for the top arch, the vertical displacement of measurement points 1#∼6# shows an increasing trend in general, indicating that with the continuous sprouting, expansion, and connection of cracks, the surrounding rock is continuously spalling downward, resulting in the increasing vertical displacement. Among them, the vertical displacement of measurement point 1# is one order of magnitude higher than that of other measurement points, because the point is close to the free face of the tunnel top, and the point is most vulnerable to the unloading effect. The vertical displacement of measurement points 2#∼5# begins to increase rapidly at 827700 calculation steps, and the growth rate begins to slow down at 829400 calculation steps. It can be seen from Figure 6 that the interval between 827700 and 829400 steps is the stage of rapid crack propagation at the tunnel top. After 829400 steps, the V-shaped failure of the crack has been formed, the crack propagation speed gradually slows down, and the growth trend of vertical displacement is basically consistent with the propagation trend of the crack.
Similar to the minimum principal stress distribution law of the bottom floor, the vertical displacement of each measurement point of the bottom floor can be roughly divided into three cases: (1) the vertical displacement of measurement points 1#, 2#, and 3# always increases slowly; (2) the vertical displacement of measurement points 4# and 5# decreases first and then increases; and (3) the vertical displacement of measurement point 6# always decreases slowly. Among them, 1#, 2#, and 3# measurement points are in the inner side of the V-shaped failure of the bottom floor. With the tunnel excavation and unloading, the bottom floor crack penetration causes the bottom floor to bulge upward. Therefore, the vertical displacement of these three measurement points is increasing in the upward direction. Moreover, the displacement growth rate of 1# and 2# is larger than that of 3#, which conforms to the law that the closer to the free surface, the larger the displacement is. 4# and 5# measurement points are located in the area covered by the continuous expansion of subsequent cracks. When there is no crack, the overall displacement direction is downward under the action of self-weight stress. After the crack is generated, the displacement direction changes to upward and keeps growing. The 6# measurement point is located at the crack boundary during stable propagation, and its vertical displacement is less affected by crack propagation, and the displacement direction is always downward.
In general, the trend of stress and deformation at characteristic points of surrounding rock can better reflect crack propagation, which proves that the calculation results are reliable.
5.5. Essence of the Spalling
According to the stress and deformation analysis, the essence of spalling can be obtained. The underground tunnel produces stress unloading during excavation. Before the tunnel excavation, the rock mass only bears the action of the initial stress and is in a stable state of three-dimensional compression. After the excavation, the rock mass around the tunnel is unloaded, all the in-situ stress around the tunnel is released, and the rock mass is in a two-dimensional or one-dimensional stress state. The original equilibrium state of the surrounding rock stress field is broken, and the original approximately uniform initial stress field becomes an obvious uneven stress field. Stress concentration occurs in some surrounding rock areas, and the stress concentration areas continue to migrate to the interior of the surrounding rock, resulting in a certain high-stress area in the interior of the surrounding rock. When the adjusted stress state reaches the limit state of rock mass, the rock mass is damaged, and cracks begin to form in the surrounding rock of the tunnel boundary with the highest stress concentration and expand with the increase in pressure. It can be seen that the stress mutation caused by tunnel excavation and unloading, and local high stress, are the main causes of surrounding rock damage.
At the same time, the progressive failure of the surrounding rock is also reflected by the simulation of the roof spalling process. Firstly, the stress of the surrounding rock begins to adjust after the excavation of the exploratory tunnel, and the tangential stress of the top arch keeps increasing while the vertical stress keeps decreasing, resulting in circumferential splitting cracks in the hard, brittle basalt of the top arch. Furthermore, the deterioration of the stress condition of the top arch and the further deterioration of the surrounding rock make the splitting cracks propagate to the depths and the oblique intersecting shear cracks appear. Finally, the further propagation of shear cracks provides boundary conditions for the falling of flake rocks, resulting in the spalling of surrounding rock.
6. Conclusions
As the largest hydropower station under construction, the failure simulation of the hard, brittle surrounding rock in Baihetan underground exploration tunnel is of great significance to reveal the brittle failure characteristics of basalt. In this paper, the brittle constitutive model of basalt is studied in depth. The finite discrete element method CDEM is used to simulate the brittle failure of the surrounding rock, analyze the stress and displacement of the characteristic points, and compare them with the in-situ observed failure so as to test the rationality of the in-situ stress, constitutive model, and mechanical parameters. The main conclusions are as follows:(1)The finite discrete element method (CDEM) is used to simulate the indoor uniaxial and triaxial tests. The fracture energy constitutive model considering cohesion weakening and friction angle strengthening can better simulate the brittle failure behavior of basalt under different confining pressures.(2)In the numerical simulation of the exploratory tunnel, the maximum depth of crack propagation in the top arch of PD62-2 is about 0.7 m, which is consistent with the depth of surrounding rock spalling observed in the field. The stress and displacement of surrounding rock characteristic points further verify the accuracy of the calculation, which strongly illustrates the reasonableness and practicability of the constitutive model and mechanical parameters adopted in this paper.(3)Using CDEM, the brittle failure mechanism of hard surrounding rock under high stress can be obtained, which intuitively presents the evolution process of progressive failure of surrounding rock and provides scientific and technical support for the design and construction of tunnels.
Due to the problem of computational efficiency, numerical simulation is limited to two-dimensional calculations, which needs to be improved in the future.
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 reported here was supported by the China Postdoctoral Science Foundation with Grant no. 2021M690999.