Abstract
The analysis and simulation of progressive failure of surrounding rock is very important for analyzing the stability of surrounding rock in underground engineering. Size effect is also a key problem worth further study in engineering. Taking the underground powerhouse on the right bank of Baihetan as an example, the acoustic test results are collected and the relaxation and failure characteristics of the surrounding rock are summarized. Then, the numerical simulation of progressive failure of surrounding rock of underground powerhouse is carried out by using the finite discrete element method CDEM (continuum-based discrete element method). The results are compared with the acoustic test results of the surrounding rock relaxation layer, and the stress and displacement of surrounding rock characteristic points are analyzed. At the same time, the size effect of grid and mechanical parameter in the process of numerical simulation are discussed. The calculated fracture depth of surrounding rock is in good agreement with the acoustic test results, which shows the reliability of progressive failure simulation of surrounding rock of the underground powerhouse. When CDEM is used to simulate the excavation of tunnels with different tunnel diameters, the minimum grid size should be about 1% of the tunnel diameter. The mechanical parameters of rock mass have significant size effect, which needs to be analyzed in detail. The research results prove the superiority of the CDEM method in simulating the progressive failure of hard surrounding rock and its unique size effect characteristics, which can provide technical reference for the application of the CDEM method in other similar engineering problems.
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]. Some scholars have studied the failure mechanism of the surrounding rock in Jinping and Baihetan underground powerhouses [4, 5]. It is of great theoretical value and academic significance to study the progressive failure of the surrounding rock in underground engineering and explore the grid size effect and rock mass parameter size effect in numerical simulation analysis to ensure engineering safety and improve the theory of rock mass mechanics.
The numerical method to simulate the failure of the surrounding rock in underground engineering is mainly divided into three categories: continuous medium mechanics method, discontinuous medium mechanics method, and continuous discontinuous combination method. Among them, continuum mechanics methods include finite difference method (FDM), finite element method (FEM), and boundary element method (BEM). When simulating the mechanical behavior of the rock, these methods are difficult to reflect the strain softening phenomenon caused by local damage and failure of materials, and the characterization of macrocrack and fracture propagation is not intuitive. Therefore, new models or methods are introduced such as Martin’s high-order constitutive criterion [6], Cosserat’s micropolar model [7], nonlocal model [8], and meshless method [9]. Later, based on the addition of shape function with discontinuous property to the traditional finite element equation, the generalized finite element method (GFEM) [10] and extended finite element method (XFEM) [11–13] appeared. For example, Belyschko et al. [12] used XFEM to simulate a tunnel in fractured rock mass, and represented the crack with the discontinuity of internal displacement of rock mass.
Discontinuous mechanics method can directly simulate rock mass fracture and its interaction, but it is difficult to be used in the structural numerical analysis of large-scale complex projects due to the limitation of calculation ability and algorithm stability and the difficulty of defining mesoparameters. Discontinuous medium methods can be divided into two categories: discrete element method (DEM) [14–16], discontinuous deformation analysis method (DDA) [17], and numerical manifold method (NMM) [18]. For example, Duan et al. [15] used the DEM method to study the damage of hard rock caused by excavation unloading.
Based on the advantages and adaptability of the continuous method and discontinuous method, numerical methods combining continuous and discontinuous medium method (such as FEM/BEM, FEM/DEM, and FDM/DEM) are gradually developed. Among them, the coupling method of the finite element and discrete element is the most mature, and many achievements have been made [19–22]. For example, the author uses the CDEM method to simulate the progressive failure of the surrounding rock of Canadian URL test tunnel [19] and Baihetan exploration tunnel [18]. Besides, the coupling of finite element methods and peridynamics [23–25] also developed to solve similar problems.
Surrounding rock failure is a process from continuous to discontinuous, accompanied by multiple crack propagation and large displacement. The finite discrete element method combining the continuous method and discontinuous method is more suitable to study the gradual transformation process of the surrounding rock from continuous to discontinuous after excavation. However, in the finite discrete element method, the crack can only occur along the triangular boundary, so the mesh size and mesh direction have a significant impact on the simulation results. For the grid direction, unstructured grid generation is generally adopted [26, 27]. As for the grid size, there is still no good method to select the grid size around the rock test samples with different sizes and tunnels with different diameters [28, 29].
In addition, in the numerical analysis of underground engineering, the value of rock parameters is very important. Size effect is one of the inherent properties of rock materials. The size effect of rock mechanical parameters includes strength and deformation. The main research methods are tests and numerical analysis. Among them, the research on the size effect of strength is mainly on rock compressive strength [30, 31], less on shear strength [32], and whether there is size effect on rock deformation parameters is controversial. So far, a lot of scholars have carried out a large number of tests on different kinds of rock mechanical parameters, but no unified law of rock size effect has been summarized. And the size effect law of rock strength explained in the theory cannot correspond to the test results. Therefore, the research on the size effect of rock mechanical parameters is still an important research direction in geotechnical engineering.
The author has accumulated some achievements in using the finite discrete element method CDEM to simulate crack propagation, including the uniaxial and triaxial compression test of basalt, the URL test tunnel in Canada, and the Baihetan exploratory tunnel. In this paper, firstly, the same method is used to simulate the layer-by-layer excavation of the right bank powerhouse to show the reliability of the calculation results. Then, based on the previous simulation results of the indoor test and exploratory tunnel failure, the appropriate grid size and the variation law of mechanical parameters of basalt under different simulation scales are analyzed. The research results can provide technical reference for using CDEM to solve similar problems in other engineering.
2. Overview of Right Bank Powerhouse of Baihetan
2.1. Background of the Project
Baihetan Hydropower Station is one of the four step power stations in the lower reaches of Jinsha River, which will become the second largest giant hydropower station after the Three Gorges after completion.
The main caverns of the underground powerhouse on the right bank include the main and auxiliary powerhouse caverns, the main transformer caverns, the draft tube bulkhead gate chamber, and the tailrace surge chamber, which are arranged in parallel.
The powerhouse was excavated from top to bottom in ten layers. The total length of the underground powerhouse is 438 m, the height is 88.7 m, the top arch elevation is 624.6 m, the rock beam elevation ranges from 602.3 m to 604.4 m, and the widths below and above the rock beam are 31 m and 34 m, respectively.
2.2. Geological Conditions
The main and auxiliary powerhouses on the right bank are monoclinal rock strata, and the trend of rock strata intersects the powerhouse axis at a large angle of 60°∼70°. The exposed lithology is mainly aphanitic basalt, column jointed basalt, plagioclase basalt, almond basalt, breccia lava, and thin tuff from layer to layer. Among them, the exposed lithology of the right bank underground powerhouse is mainly massive basalt with breccia lava, which is hard, slightly new, and free of unloading. The rock mass structure is mainly submassive, with partial massive structure. The surrounding rock is mainly classified as Class III1, accounting for about 70%.
The surrounding rock of the right bank underground powerhouse caverns is developed with large soft interlayer staggered zone, small faults, random fractures, dense columnar joints, and other unfavorable structures, among which the interlayer staggered zone has the largest scale and the most prominent impact. The right bank underground powerhouse caverns are mainly affected by the interlayer staggered zones C3, C3-1, C4, and C5.
2.3. In Situ Stress
The initial stress of the right bank underground powerhouse is generally stable, except for the local area affected by large structures, such as interlayer belts, and is dominated by tectonic stress fields and denudation transformation. The horizontal stress is significantly greater than the vertical stress. The maximum principal stress direction of the right bank underground powerhouse area is N-S direction, horizontal. The intermediate principal stress is E-W direction, which inclines to the valley about 0°∼5°, and partially inclines to the mountain due to the influence of the interlayer zone. The minimum principal stress is nearly vertical. According to the measurement, the stress size meets the following relationship:
3. Relaxation Characteristics of the Right Bank Powerhouse
3.1. Acoustic Test Results
A total of 6∼9 acoustic test sections are arranged in the underground powerhouse on the right bank, and the typical section is shown in Figure 1. Acoustic detection holes are distributed at the top arch of 624.6 m elevation, arch shoulder of 622.6 m and 617 m elevation, and side wall of 608 m, 601 m, 591 m, 572.4 m, and 560 m elevation.
Multiple acoustic tests are carried out in different excavation stages, and the results are summarized in Tables 1–3. According to the test results:(1)Crown and spandrel: The relaxation depth of the surrounding rock is generally about 0.9∼2.5 m, with an average of about 1.54 m.(2)Side wall: The relaxation depth of surrounding rock at el. 608 m is about 2.6∼4.0 m, and the average depth is 3.50 m. The relaxation depth of the surrounding rock at el. 601 m is about 3.0∼5.6 m, with an average depth of 4.08 m. The relaxation depth of the upstream side wall at el. 591 m is about 1.0∼2.6 m, and the average depth is 2.20 m. The relaxation depth of the downstream side wall at el. 583 m is about 2.0∼2.6 m, with an average depth of 2.27 m. The relaxation depth of the upstream side wall at el. 572 m is about 1.4∼2.6 m, with an average depth of 1.99 m.
3.2. Failure Characteristics of Surrounding Rock
The rock mass of the right bank powerhouse is mainly brittle rock. And the average in situ stress is 26 MPa, belonging to the high ground stress level. Many typical stress-based failure phenomena occurred during the excavation of the powerhouse, such as spalling, fracture, relaxation collapse, and concrete spray layer cracking, as shown in Figure 2.
4. Numerical Simulation of the Right Bank Powerhouse
4.1. Geometric Model
Referring to the previous numerical simulation of URL test tunnel and Baihetan exploratory tunnel conducted by the author, the built-in software Gmsh is used for grid division, as shown in Figure 3.
The size of the whole model section is 600 m × 600 m. Area 1 is the powerhouse cavern to be excavated, and its grid size is 0.3 m. Area 2 is a square area close to the cavern with a range of 74 m × 133.9 m. Cracks first appear in this area during excavation, the grid size should be small enough to be 0.3 m. The range of area 3 is 300 m × 300 m, the grid size gradually changes from 0.3 m to 15 m. The range of area 4 is 600 m × 600 m, and the grid size changes from 15 m to 25 m. The total number of elements is 324244.
The right bank powerhouse is simulated according to the actual excavation process of eight layers, and the total excavation height is 73.9 m. The excavation heights from the first layer to the eighth layer are 13.6 m, 4.1 m, 11 m, 6.1 m, 11 m, 11 m, 4.9 m, and 12.2 m, respectively. The span of the top and bottom of the powerhouse is 34 m and 31 m, respectively.
4.2. Geostress and Boundary Conditions
The average burial depth of the powerhouse is about 500 m. The included angle between the first principal stress and the axis of the powerhouse is 10°, and the stress on the cross section of the powerhouse is mainly affected by the second principal stress. Ignoring the small dip angle of the second principal stress inclined to the valley, the in situ stress input by numerical simulation can be calculated, as shown in Table 4.
The boundary conditions adopted in the calculation are as follows: a normal displacement constraint is applied at the bottom. The normal stress is applied on the left and right sides, the shear stress is applied on the left and right, upper and lower sides, and both are applied under the condition of surface force. Then, the initial stress field of the model is obtained after calculation to the steady state.
4.3. Calculation Steps and Mechanical Parameters
As for the feedback analysis of the powerhouse tunnel group during the construction, a large number of numerical calculations have been carried out and rich results have been achieved. In their work, numerical analysis software FLAC3D or 3DEC is usually used for calculation, and the Hoek–Brown model and corresponding parameters are selected.
In this paper, the software CDEM is used for calculation. The constitutive model is based on the Mohr–Coulomb strength failure criterion. The acquisition methods of rock mass strength parameters and deformation parameters are as follows:(1)extract the rock mass density, uniaxial compressive strength UCS, geological strength index GSI, and the coefficient mi from the previous calculation report;(2)input the abovementioned parameters and tunnel buried depth into the RocData program to obtain the basic parameters of rock mass.(3)input the basic parameters of rock mass into the numerical calculation model as initial data, and the failure characteristics of surrounding rock with given parameters are calculated and compared with the failure characteristics observed.(4)In the process of feedback analysis, GSI and UCS shall be properly adjusted according to the difference of local rock mass quality until the calculation results are fully close to the observed. The final mechanical parameters adopted in the calculation are shown in Table 5.
The specific calculation steps are as follows: The initial in situ stress field is firstly calculated with a total of 1,000,000 steps. Then, the excavation simulation of the powerhouse start shall be carried out. Each layer of excavation shall be calculated in 3,000 steps until it is stable, which shall be carried out in 8 times. The total number of calculation steps is 1,024,000.
5. Results Analysis
5.1. Fracture Evolution of Surrounding Rock
The evolution process of internal cracks in layers I∼VIII of the powerhouse is obtained after the completion of excavation, as shown in Figure 4. At the same time, the damage of one-dimensional spring, i.e., contact element, is monitored during the calculation, and a total of three variables are statistically output, as shown in Figure 5. Among them, the total damage of spring is the proportion of elements that have been damaged. The total rupture of spring is the proportion of elements that have not only been damaged but also have the cohesion and tensile strength reduced to zero. The current damage of spring is the proportion of elements that are now in a damaged state.
As can be seen from Figure 4, after the excavation of layer I is completed, a circle of cracks appear in the surrounding rock of left and right bank spandrels and crown arches, and the cracks are evenly and symmetrically distributed. After the excavation of layer II, the crack propagation range increases, mainly along the exposed side wall of layer II, and the crack propagation depth increases slightly compared with that after the excavation of layer I. After the excavation of layer III, the number of cracks increases significantly. This is the design position of rock beam where the excavation surface is irregular. The cracks expand rapidly to the depth along the side walls on both sides, and the propagation depth is much greater than the crack depth of the top arch and left and right spandrels. After the excavation of layers IV∼VIII is completed, the law of crack propagation is relatively consistent. The cracks expand and connect along the exposed side wall, and the crack propagation depth is large. Under the action of cracks, the side walls on both sides tend to deform towards the free face, which can be reflected from the inclination of the side wall after the excavation of layer VIII. This phenomenon may be related to the high height of the side wall. As the excavation boundary of the pit in layer VIII is irregular, the crack distribution density is large at the turning point of the right footing, and the edge of the footing tends to rise.
There is a good correspondence between Figures 4 and 5. After the excavation of the powerhouse, the three fracture curves continue to grow with the advance of calculation steps. Taking the failure condition when reaching stability as an example, it can be seen from the total damage curve of spring that the proportion of cracks in stable state does not increase much during the excavation of layers I∼IV. From layer V to layer VIII, a large number of cracks are generated in the excavation process of each layer, and the proportion of cracks expands in an approximate S-shaped trend. At the end of each layer excavation, the crack propagation basically reaches a stable state. The total rupture curve of spring and the current damage curve of spring are the completely damaged elements and the elements now in the damaged state, respectively, and their proportion is far lower than that of the total elements ever damaged. The trend of growth is also relatively stable, and there is no significant S-shaped growth.
From the numerical simulation results, it can be seen that the failure simulation of the surrounding rock in the process of powerhouse excavation reflects the progressive failure characteristics of the surrounding rock. Firstly, after the excavation of powerhouse top arch, the stress of the surrounding rock begins to adjust, the tangential stress of the top arch increases continuously, and the vertical stress decreases continuously, resulting in circumferential splitting cracks in the top arch. Secondly, with the gradual increase of excavation depth, further stress adjustment of the surrounding rock begins. The local stress concentration begins to appear in the left and right side walls, and the cracks continue to sprout and expand, which provides favorable boundary conditions for the spalling and block falling of the side wall rock mass. And the upward bulging and cracking may occur in the bottom foot.
5.2. Comparison between Simulated and Observed Fracture
Taking the excavation of layer II, layer IV, and layer VI as an example, Figure 6 shows the feedback results of the relaxation depth of the crown arch and side wall after excavation of each layer. In order to facilitate analysis, the characteristic points at different elevations are numbered. Among them, 0 represents measuring point at crown arch of 624.6 m elevation; 1 and 1′ represent measuring points at upstream and downstream spandrels of 622.6 m elevation; 2 and 2′ represent measuring points at upstream and downstream abutment of 617 m elevation; 3 and 3′, 4 and 4′, 5 and 5′, and 6 and 6′ represent measuring points at upstream and downstream side walls of 608 m, 601 m, 591 m, and 583 m elevation, respectively. According to the feedback analysis results:(1)During the excavation of layer VI, the relaxation depth of the crown arch increases to a certain extent. The relaxation depth generally increases by about 0.0∼0.4 m, and the relaxation depth is 1.6∼3.2 m.(2)After the excavation of layer VI is completed, the relaxation depth of the upstream and downstream side walls at 608 m increases by about 0.2∼0.7 m, and the relaxation depth is 4.5 m and 4.6 m, respectively. The relaxation depth of upstream and downstream side walls at el. 601 m is 5.1 m and 6.0 m, respectively. The relaxation depth of upstream and downstream side walls at el. 591 m is 4.9 m and 6.0 m, respectively. The relaxation depth of upstream and downstream side walls below el. 583 m is 3.1 m and 3.9 m, respectively.
(a)
(b)
(c)
The simulation results are basically consistent with the deformation characteristics and relaxation depth measured in the field. The model adopted is simplified in the calculation in this paper, without considering any weak zone. The rock mass is considered to be homogeneous and isotropic, resulting in some errors in the calculation results. But on the whole, the numerical calculation can better reflect the failure of the surrounding rock of the underground powerhouse, and the results are reliable.
5.3. Stress Analysis of Surrounding Rock
In order to further analyze the failure of the roof arch and side wall of the surrounding rock, the corresponding measuring points are selected at different excavation stages around the tunnel. All measuring points are the points closest to the free face. The specific location distribution is shown in Figure 7. The evolution curves of the minimum principal stress of each measuring point on the left and right side walls with the calculation step are counted, as shown in Figure 8.
(a)
(b)
The evolution law of the minimum principal stress of the measuring points on the left side wall is basically the same as that of the measuring points on the right side wall. Take the measuring points on left side wall as an example. From layer I to layer VIII, the minimum principal stress of each measuring point first experiences a continuous growth, then drops rapidly at a certain time, and then oscillates and fluctuates near the lower value, and the maximum stress can reach about 70 MPa. Among them, the stress drop time of each measuring point from layer II to layer VIII is the excavation time of this layer, as shown in the figure, which is caused by strong unloading of the surrounding rock during excavation. The measuring points of layer I are slightly different. After the excavation of layer I, the stress does not decrease rapidly, but increases continuously. The possible reason can be closely related to the location of the measuring point, which is selected at the spandrels on both sides. After excavation, the stress concentration here is relatively significant. The stress does not begin to decrease until the excavation of layer II. In the excavation process of subsequent layers, the stress of the measuring point does not decrease to the minimum value, but fluctuates around 50 MPa.
5.4. Deformation Analysis of Surrounding Rock
The same measuring points are selected for the deformation analysis of the surrounding rock, and the evolution curve of the lateral displacement of the measuring points on left and right side walls with the calculation step is shown in Figure 9.
The lateral displacement of each measuring point on left and right side walls is symmetrically distributed. The left measuring point deforms to the right and the right measuring point deforms to the left, resulting in the deformation of the rock mass on both sides to the middle free face after the completion of excavation, which is easy to form spalling and block falling. From layer I to layer VIII, the displacement of measuring points of each layer begins to increase after the excavation of this layer is completed, and the growth trend is approximately linear. The displacement growth rate of the measuring points of layer I is the slowest. After the excavation of layer II is completed, the growth rate remains basically stable, approximately zero. The displacement growth rate of measuring points of layer II is the second. After excavation, the deformation increases slowly. When the excavation of layer III is completed, the displacement increases rapidly. Then, the deformation speed slows down gradually until the excavation is completed, and the final deformation is about 0.2 m. The displacement of measuring points of layers III∼VIII starts to increase linearly after the excavation of this layer, and the growth rate of measuring points on each layer is basically the same. At the end of the calculation, the lateral displacement of the measuring point of layer III is the largest, about 0.65 m. This layer is at the position of the rock beam, with the largest deformation.
It should be noted that, since the model for calculation in this paper is simplified, there may be large errors in the calculated data of displacement, and the evolution law of displacement should be paid more attention to.
6. Analysis of Size Effect
6.1. Size Effect of Grid
Because the crack of CDEM can only develop along the triangular element, in order to accurately simulate the stress field at the crack tip, the mesh is often required to be as small as possible, but it also brings the problem of increasing the amount and time of calculation. On the other hand, large mesh size will increase the yield displacement, resulting in too high calculation strength. Therefore, it is of great significance to study the appropriate grid size to improve the computational efficiency.
Table 6 summarizes the excavation diameter, grid size, and the ratio of grid size to excavation diameter of indoor test simulation, Canadian URL test tunnel simulation, Baihetan exploratory tunnel simulation, and Baihetan powerhouse simulation. It can be seen from the table that the ratio of the minimum grid size to the excavation diameter is about 1/100. At this time, the finite discrete element method CDEM can better simulate the progressive failure of the hard brittle rock and hard brittle surrounding rock.
6.2. Size Effect of Mechanical Parameters
Due to the existence of size effect, it is unreasonable to directly apply the mechanical parameters of a specific size of rock mass to engineering design and the establishment of the constitutive model. Similarly, the test results of indoor small-size rock samples cannot be directly extrapolated to the on-site rock mass. The indoor test, exploratory tunnel, and powerhouse are numerical simulations on three different scales. Table 7 summarizes the various parameters under different scales of numerical simulation. Among them, the indoor test simulation includes two types of basalt of almond basalt and cryptocrystalline basalt.
It can be seen from Table 7 that the mechanical parameters of almond basalt and cryptocrystalline basalt are slightly different, mainly in the value of deformation modulus, cohesion, and tensile strength. As far as the simulation results of the three scales are concerned, both the deformation parameters and the strength parameters show a certain degree of size effect. The deformation modulus first decreases and then increases from small scale to large scale. The cohesion and tensile strength decrease greatly. The internal friction angle increases gradually, but the increase range is small. The dilatancy angle and fracture energy remain unchanged.
6.3. Discussion
In the paper about URL test tunnel [19], the key factors of crack propagation simulation using the continuous discontinuous analysis method are discussed, including size effect, constitutive model, and excavation method. 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 previous paper focused on the grid size effect of the same engineering scale, while this paper focuses on engineering rock mass of different scales. The results show that, in order to obtain the simulation results as accurate as possible, the number of grids should not only be sufficient enough but also meet the condition that the ratio of grid size to the excavation diameter is about 1/100.
In addition, during the transition from rock to rock mass, due to the weakening of homogeneity and the addition of joints, the role of cohesion in the brittle failure of rock mass is gradually weakened, while the role of friction is gradually enhanced.
7. Conclusions
The progressive failure simulation of the surrounding rock in underground engineering is of great significance for analyzing and ensuring the stability of the surrounding rock. In this paper, the finite discrete element method CDEM is used to simulate the excavation of right bank powerhouse of Baihetan. The fracture evolution of the surrounding rock, the stress, and displacement of characteristic points are analyzed. The calculated results are compared with the relaxation layer depth measured by the acoustic test. At the same time, the size effect of grid and mechanical parameters are discussed when adopting CDEM for simulation. The main conclusions are as follows:(1)The fracture depth of the surrounding rock during the layer-by-layer excavation of the powerhouse is basically consistent with the relaxation layer depth measured by the acoustic test, which shows the reliability of the calculation results and reflects the rationality and practicality of the constitutive model and mechanical parameters adopted in this paper(2)When the finite discrete element method CDEM is used to simulate the fracture evolution of the surrounding rock in the process of tunnel excavation with different tunnel diameters, good simulation results can be obtained when the minimum size of the grid around the tunnel is controlled at 1% of the tunnel diameter(3)Under different scale simulating conditions, the mechanical parameters of rock and rock mass cannot be universal, and the variation laws of deformation parameters and strength parameters cannot be unified, so it needs to be analyzed in specific circumstances
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.
Acknowledgments
The work reported here was supported by China Postdoctoral Science Foundation with Grant no. 2021M690999.