Abstract
The traditional Three-dimensional (3D) model for stability analysis of a large underground cavern group only includes the main powerhouse, the main transformer chamber, the tailrace surge tank, the main electrical wire hall, and part of the diversion tunnel and tailwater tunnel. Other caverns excavation are not taken into consideration, which is inconsistent with the actual project. In this study, a 3D Finite Element Method (FEM) model, including all the cavities, and a 3D simplified model containing only part of the cavities were established based on the actual underground cavern group project. We study the structural characteristics of the cavern group during excavation via 3D nonlinear FEM. At the same time, three numerical calculation schemes are designed for comparison. The results show that the simulated distributions of the stress field, deformation field, and plastic zone of the two numerical models are similar, but the magnitudes of the 3D FEM model, which considers all caverns are generally higher than the simplified model. The difference between the two simulated horizontal displacements of the main powerhouse reaches 53.72%. Therefore, from the perspective of engineering safety and actual construction, it is advisable to adopt the 3D FEM model, including all the cavities, when conducting a simulation analysis of the underground cavern group.
1. Introduction
Hydroelectric power which is an environmentally friendly, cheap, and safe way to obtain electricity can meet the growing electricity demand [1]. Water energy resources are mostly distributed in high mountains and valleys, not suitable for constructing ground buildings. In order to avoid adverse topographic, geological, and climatic factors on the surface, engineers often choose to construct underground powerhouse for the hydropower station. The underground powerhouse is generally composed of large or super-large caverns [2, 3]. A huge underground cavern group’s construction process is more complicated than that of a single cavern group due to the mutual influence of several caverns in the cavern group. However, the construction of an extensive underground cavern process is time consuming. Excavation will result in the redistribution of stress in the cavern. The stress field change around the cavern may lead to the deformation and instability of surrounding rock and rock collapse, which presents a significant threat to the structure, construction equipment and may cause the loss of life and property [4–6]. Therefore, during the excavation of underground caverns of hydropower stations, the stability of the surrounding rocks has always been the focus of the engineering and academic circles.
Many researchers have looked at the stability of the rock that surrounds underground caverns. One of the main reasons for the deformation of the surrounding rock is the ground stress. Before the excavation of the cavern, the rock and soil are generally in a state of natural stress balance, which is called the primary stress state or initial stress state (including self-weight stress and tectonic stress). After the cavern is excavated, the surrounding rock loses its original support and expands into the cavern space [7–9]. Zhang and Wang determined the damage characteristics of deeply buried caverns and the stress distribution characteristics of the cavern perimeter under different ground stress conditions by numerical model tests [10]. Ma and Zhang discussed the effects of high ground stress, deep unloading and construction processes under complex geological conditions. The comprehensive method combining experiment and simulation reveals law of deformation and failure of surrounding rock of large underground caverns [11]. Dai and Li et al. used a high-resolution micro-seismic monitoring system to assess the stability of underground powerhouses, and developed a comprehensive analysis method that can determine the damaged area and predict the macroscopic deformation of the surrounding rock of the underground powerhouse [12].
Under the action of high initial in-situ stress and excavation unloading of caverns, deep rock masses exhibit unique mechanical behaviors, such as ductility-brittle transition of rock mass, zonal cracking of surrounding rocks, accumulation and sudden release of deformation energy, brittle failure of unloading and so on [13–15]. Li and Wang et al. used the displacement of key points on the sidewall, the volume of the plastic zone and the volume of the crack to characterize the stability of the underground cavern. The optimal layout of underground caverns under different in-situ stress conditions and a modified formula for predicting sidewall displacement are proposed [16]. Zhang and Wang et al. proposed using index of damage degree to describe the Excavation damaged zone (EDZ) via the test results of 38 sets of drilling holes on the high side wall of the main powerhouse [17].
The evaluation of surrounding rock stability is an important research content, reflecting the stability degree of surrounding rock in the process of underground excavation and support and the comprehensive indexes of geological composition, environmental conditions, and engineering factors [18, 19]. Ren and Xu proposed a surrounding rock stability evaluation method that can analyze the deformation stability of surrounding rock and predict the collapse of a rock block quickly by combining DSPEM and RPPMCHE CWIS [20].
For the underground powerhouse of a hydropower station, the layout of the cavern group plays an important role in improving the stability of the surrounding rock of the cavern. Nowadays, great attention has been paid to the selection of the longitudinal axis of the cavern and the distance between the caverns in the design. Researches on the reasonable excavation process, support timing, support strength of the giant underground cavern group and related design theories have become major problems [21, 22]. Using the advantages of PSO and SVM algorithms, Jiang and Su et al. devised an intelligent optimization system for cavern excavation. This approach efficiently reduces the overall volume of the plastic zone and brittle failure of the surrounding rock, according to numerical simulation study [23].
Most of the underground caverns are located in deep canyons. Their geological conditions and geological structures are very complex, including various faults, joints, and weak interlayers [24, 25]. Duan and Feng et al. took the impact of weak interlayers into consideration and proposed an in-situ observation scheme based on geological monitoring data, revealing the development and failure mechanism of the rock mass failure process [26]. Hao and Azzam used numerical methods (UDEC) to establish induced plastic zones and evaluated the influence of faults on the failure and deformation behavior of rock masses around large underground caverns [27].
In the past, people believed that earthquakes did not damage underground facilities as severely as ground structures. Until the Kobe earthquake in the 1990s, little attention was paid to seismic research on underground facilities [28]. The Daikai subway station was the first urban underground building to be destroyed by seismic forces rather than ground instability [29]. Wang and Chen et al. investigated the seismic response of a huge underground powerhouse’s lining structure and suggested an explicit dynamic contact analysis method that took into account the surrounding rock’s bonding and damage characteristics as well as the lining interface [30]. Yang and Lu et al. used time-energy density analysis, amplitude spectrum analysis, and finite impulse response (FIR) digital filter methods, and pointed out that the low-frequency components in micro-seismic records were mainly caused by the instantaneous release of in-situ stress. The high frequency components are mainly caused by explosions [31].
Facts have proved that as the depth increases, the original stress will show a linear or nonlinear increase [32, 33]. The original ground stress caused by the overburden, site conditions, and excavation operations may lead to stress concentration, which can lead to cracks and damage [34, 35]. Rockburst hazards can cause significant damage to the foundation of structures and equipment, and may threaten the safety of workers [36, 37]. Fraadonbeh et al. used two self-organizing maps (SOM) and fuzzy c-means (FCM) robust clustering techniques to explore the relationship between rockburst related parameters [38].
Most of these studies on the stability of underground caverns use three-dimensional simplified models that include the main powerhouse, main transformation chamber, busbar tunnel, some diversion tunnels and tailrace tunnels, or only three main caverns, without considering the impact of the excavation of other caverns on stability.
This paper takes an underground cavern of a hydropower station as the research object, establishes a three-dimensional finite element model including all caverns. And studies the distribution of deformation, stress and plastic zone of the surrounding rock during the excavation of the cavern group. By comparing the calculation results with commonly used three-dimensional simplified models, the influence of different calculation models on the calculation results is analyzed.
2. Project Profile
The underground powerhouse system of an individual project adopts a tail layout, which contains the main and auxiliary powerhouse, the main transformer chamber, the main electrical wire hall, the main transformer transport tunnel, etc (see Figure 1). With an axial direction of N80°W, the main powerhouse is 169.00 m long, 25.40 m wide, and 54.425 m high. The auxiliary plant is on the right side of the main plant, and the installation space is between Unit 2 and Unit 3. The main transformer chamber is 139.40 m long, 20.20 m wide, and 21.20 m high. It runs parallel to the main powerhouse and is situated on the downstream side of it. The cutway drawing of Unit 4 is shown in Figure 2.


3. Calculation Model and Excavation Scheme
3.1. Numerical Model
A schematic diagram of the full cavern model’s (R1) numerical calculation area is shown in Figure 3(a). The orientations of the three coordinate axes are as follows: the X-axis is the direction of the main plant, the direction pointing to the auxiliary plant is positive; the direction perpendicular to the axis of the main plant is the Y-axis, the direction pointing to the main transformation room is positive; the vertical direction is the Z-axis, and upward is positive. The origin is taken at the main center of Unit 1, with an elevation of 451.9 m. The numerical calculation area is 1320 m on the X-axis, 715 m on the Y-axis, and 718 m on the Z-axis. The calculation area includes the main plant, main transformer room, bus tunnel, main transformer transportation tunnel, connection corridor, blower room, exhaust fan room, access tunnel, main power plant delivery, exhaust tunnel, main transformer room delivery, and exhaust tunnel, Diversion tunnel, tailrace tunnel, high-voltage cable flat tunnel, high-voltage cable shaft, maintenance drainage gallery, exhaust shaft, construction branch tunnels 3, 4, 5, and 6, roof anchoring tunnels, and surrounding main powerhouses and main transformer rooms The outer three layers of drainage and anchoring holes while considering the geological structure of the plant area. A total of 649,612 gird nodes and 3,802,630 eight-node isoparametric elements are used in the 3D computational model (see Figure 3(b)).

(a)

(b)
The boundary conditions of the model during calculation: the normal displacement is constrained on the ±X surface of the mountain model, the normal displacement is constrained on the ±Y surface, +Z (surface) is the free boundary, and the −Z surface is constrained in X, Y, and Z displacements.
A schematic diagram of the Simplified model’s (R2) numerical calculation area is shown in Figure 4(a). The orientations of the three coordinate axes are as follows: the X-axis is the direction of the main plant, the direction pointing to the auxiliary plant is positive; the direction perpendicular to the axis of the main plant is the Y-axis, the direction pointing to the main transformation room is positive; the vertical direction is the Z-axis, and upward is positive. The origin is taken at the main center of Unit 1, with an elevation of 451.9 m. The numerical calculation area is 669 m on the X-axis, 550 m on the Y-axis, and 567 m on the Z-axis. The calculation area includes the main powerhouse, main transformer room, four bus tunnels, one main transformer transportation tunnel, one connection corridor, four water diversion tunnels, four tailrace tunnels, collection wells, and the machine room and shaft connected to the powerhouse And part of the access tunnels, wind tunnels, ventilation, and safety tunnels, and also consider the geological structure of the plant area. A total of 102,542 gird nodes and 594,312 eight-node isoparametric elements are used in the 3D computational model (see Figure 4(b)).

(a)

(b)
The boundary conditions of the model during calculation: the normal displacement is constrained on the ±X surface of the mountain model, the normal displacement is constrained on the ±Y surface, +Z (surface) is the free boundary, and the −Z surface is constrained in X, Y, and Z displacements.
3.2. The Modeling Details
3.2.1. Parameter Selection
The initial ground stress was obtained by using the Moore-Coulomb principal structure method. And according to the report of geological exploration, the rock mass is a medium bed∼thick layered structure with good integrity. Class III surrounding rocks and class IV surrounding rocks, mainly argillaceous siltstone, account for about 75% and 25%. Table 1 shows the suggested values of the physical and mechanical parameters of surrounding rocks.
3.2.2. Simulation Scheme of Construction Method
According to the construction experience of large underground cavern group at home and abroad, combined with the structural characteristics of underground cavern mentioned in this paper, three excavation calculation schemes are designed for two FEM models.
Excavation calculation Scheme 1 (T1): Model R1, Considering all caverns excavation, each excavation step is shown in Figure 5 below:

Excavation calculation Scheme 2 (T2): Model R2, The main and auxiliary powerhouse, main transformer tunnel, bus tunnel, main transformer transportation tunnel and some diversion tunnel and tailrace tunnel are taken into account for excavation. The excavation chamber of each excavation step is shown in Figure 6 below:

Excavation calculation Scheme 3 (T3): Model R2, Considering the excavation of the main and auxiliary powerhouse and the main transformer chamber, the excavation chamber of each excavation step is shown in Figure 7 below:

4. Results
4.1. Characteristics of the Redistributed Stress after Excavation
Figures 8–10 shows the contour maps of the variation of maximum and minimum principal stress in surrounding rock of the characteristic section obtained after the completion of the three excavation schemes. After the cavern excavation, due to the existence of the initial stress field, the surrounding rock stress field is redistributed, which causes the radial stress of the cavern to be released and the tangential stress changes resulting in a tensile stress distribution around the cave.



The stress distribution of the three excavation schemes is similar. In the initial stage of excavation, the main tensile stress is distributed on the bottom plate because the high side wall is not formed. The top arch has good mechanical performance, and the bottom plate is equivalent to a flat beam, which is easy to be bent, so the maximum main tensile stress mostly occurs at the bottom plate. The high side wall is formed gradually with the excavation, the stress release degree near the side wall is increasing, and the maximum first principal stress is gradually transferred to the side wall. After the excavation of the underground powerhouse is completed, the range of the first principal stress of the cavern group are −5.24 to 0.77 MPa (T1), −5.08 to 0.56 MPa (T2) and −5.23 to 0.57 MPa (T3) respectively (“−” means compressive stress, “+” means tensile stress), and the compressive stress concentration areas are distributed in the downstream arch foot of the main plant, the intersection of the bottom plate and the upstream and downstream sidewalls, and the downstream arch foot of the main transformation chamber Place.
The compressive stress of the main powerhouse is concentrated at the intersection of the arch foot on the downstream side of the top arch and the upstream wall and the floor. The maximum values are −18.93 MPa, −15.13 MPa, and −15.80 MPa respectively. Similarly, the compressive stress of the main transformer room is concentrated at the arch foot on the downstream side of the top arch, and attention should be paid to protection. During the entire excavation process, the principal compressive stress did not exceed the surrounding rock’s compressive strength.
4.2. Characteristics of the Displacement Distribution after Excavation
Due to the release of in-situ stress, the surrounding rock produces spring back deformation, pointing to the cavern’s interior. It shows that under the action of the initial ground stress field, each side surface caused by cavern excavation has a “centripetal” movement trend, and the release displacement of surrounding rock excavation attenuates as the distance from the excavation surface increases.
After the cavern excavation, the displacement of the main powerhouse in the roof arch is within the range of −14 to −32 mm. The maximum downward displacement of the top arch lead is 32 mm, and the upward rebound of the bottom plate lead is 39 mm. Simultaneously, the displacements of the upstream and downstream sidewalls of the main powerhouse are within the range of −12∼50 mm and −36∼10 mm (see Figure 11). After the end of the excavation, the maximum horizontal displacement (U2) of the upstream wall perpendicular to the factory building’s axis is about 50 mm, which occurs in the upstream wall of the main factory building. The maximum horizontal displacement (U2) of the downstream wall perpendicular to the main factory building axis is about 36 mm, which occurs at the intersection of the bus hole and the main factory building. Moreover, the upstream side wall’s deformation is significantly greater than that of the downstream sidewall.

Figure 12 shows the cumulative value of vertical displacement (U3) of the roof and horizontal displacement (U2) of the sidewall in the middle section of unit 4 of the main workshop during each excavation step. It serves to show that the cumulative value of the vertical displacement (settlement of the arch) (see “P1, P2”) of the arch roof continues to increase as the excavation progresses. From the simulation results of the entire 7-story excavation, it can be seen that the excavation of the first and second floors has the most significant impact on the vertical displacement of the vault center of the cavern. When the first layer is excavated, the vertical displacement of the vault will increase by about 17 mm. When the second layer is excavated, the top arch’s vertical displacement is increased by about 3 mm. The subsequent excavation steps have little influence on the vertical displacement of the vault center. Each step from the third layer of excavation onwards has no more than 1 mm effect on the vault’s vertical displacement. The impact of the vertical displacement of the top arch is gradually decreasing. As the height of the sidewall continues to increase, the horizontal displacement (U2) (“P3∼P6”) of the sidewall also increases. The numerical simulation results also show that the fourth and fifth floors’ excavation will cause obvious deformation of the main building (“P5, P6” in Figure 12).

Like the main powerhouse excavation, each free face of the main transformer room has a “centripetal” movement. The top arch and bottom plate of the main transformer room are mainly plumb displacements. With the gradual increase of the displacement during the excavation, the top arch’s vertical displacement after the excavation reaches 23 mm (see Figure 13). The main transformer room’s lateral wall is mainly horizontal displacement (U2) perpendicular to the axis of the cave. As the excavation progresses, the lateral wall’s height increases continuously, and the horizontal displacement of the lateral wall also increases continuously. After the main transformer chamber excavation, the horizontal displacement of upper and lower walls of the main transformer chamber is within the range of −23∼32 mm. The deformation characteristics of the surrounding rock of the main transformer chamber, in which the deformation of the upstream side is greater than that of the downstream side, are near related to the “group tunnel effect” caused by the excavation of the adjacent cavern (main powerhouse and the main electrical wire hall).

Figures 14 and 15 are schematic diagrams of the surrounding displacement of the main powerhouse and the main transformer chamber of the underground cavern after the excavation of T1 and T2 is completed, respectively. Comparing Figures 11 and 13, it can be seen that the displacement change laws are roughly similar. Table 2 shows the displacement values around the hole after the excavation of the three excavation schemes.

(a)

(b)

(a)

(b)
4.3. Characteristics of the Plastic Zone Distribution after Excavation
In the process of excavation, the plastic zone appears obvious, and the distribution of the plastic zone is local. The arch abutment of the main powerhouse appears plastic deformation first and spreads to the depth of surrounding rock with the main powerhouse’s excavation. Figures 16–18 are schematic diagrams of the distribution of plastic zones in the upstream and downstream walls of the main powerhouse and the main transformer room after the excavation of the three excavation methods. Table 3 shows the failure depth of plastic zone after excavation.



With the main powerhouse excavation, plastic zones of different degrees appear at the intersection of the upper and lower sidewalls and each cave. It is mainly distributed in the upper and lower sidewalls of the sixth and seventh floors of the main powerhouse, main electrical wire hall, main transformer transportation tunnel, traffic cavity, the horizontal hole underwater diversion, and the intersection with the main powerhouse, and other positions with abrupt shape changes. The main transformer chamber’s plastic zone is mainly distributed in the middle area of the sidewall, the roof arch of the main electrical wire hall, and the sidewall of the main transformer transportation tunnel. The plastic zone distribution on the upstream side of the main transformer room is larger than that on the downstream side due to the excavation of the main electrical wire hall and the main transformer transportation tunnel on the upstream side. There is no plastic zone connecting the main powerhouse, the main transformer chamber, and the workshop’s anchorage corridor.
4.4. Comparison of Three Excavation Calculation Schemes
Figure 19 shows the surrounding rock displacement cloud map of the workshop area completed by the cavern group’s excavation. It can be seen that the excavation of the auxiliary caves and chambers will affect the deformation of the rock around the main chamber and form a large chamber group effect on the main chamber. The results show that T1, compared with T2 and T3, increases the X-direction displacement of surrounding rocks by 7.59% and 13.91%. The displacement along Y-direction increased by 53.72% and 57.27%, respectively. The Z-direction displacement increased by 21.26% and 28.60%, respectively. It can be seen that the excavation of the surrounding cavity has the most significant influence on the Y-direction displacement, followed by the X-direction displacement and the Z-direction displacement. The excavation of all caverns will have a relatively obvious “group tunnel effect.”

The contour graphs of the maximum and minimum principal stress changes in the surrounding rock after excavation for the three excavation construction options are shown in Figure 20. In the ABAQUS software, it is ruled that tension stress is complimentary, while pressure stress is negative; therefore, in the figure, S-Max represents the maximum principle stress, while S-Min represents the minimum principal stress.

All the three excavation methods can redistribute the stress of surrounding rock, and the stress is concentrated in the vault and bottom of the three caves. The maximum pressure stress induced by T2 and T3 is −22.20 MPa and −22.28 MPa, respectively, while the maximum pressure stress induced by T1 is −40.21 MPa. They cause tensile stress zones in the surrounding rock, which are concentrated in the upper and lower sidewalls and end walls of the main workshop and main transformer room. T1 tensile stress was the highest at 0.90 MPa. Then T2, 0.88 MPa; The minimum stress provided by T3 is 0.87 MPa. The results show that the excavation of all caverns will affect the stress distribution and stress extreme value of the surrounding rock of the underground powerhouse.
Figure 21 for cavities completed the excavation after the workshop section typical cloud cave surrounding rock plastic zone distribution. Figure 14 shows the T1, T2, and T3 model of plastic zone distribution, plastic zone, mainly in the main factory building abutment, 6 and 7 of the main workshop layer of upstream and downstream of the upstream and downstream sidewall and main transformer chamber second upstream and downstream sidewall, plastic zone than the downstream and upstream wall plastic zone distribution range is significant, the main plant and there is no breakthrough of the plastic zone area between main transformer room.

There is little difference between T2 and T3 in the plastic zone, in which the maximum depth of the plastic zone in the upstream wall of the T2 model is 14 m, and the depth of the plastic zone in the upstream wall of the T3 model is 13 m. The T1 plastic zone was significantly higher than T2 and T3, and the maximum plastic zone depth was 20 m. It indicates that R1 model considers more caves, and its excavation disturbance area is more extensive, which has adverse effects on the surrounding rocks of the main powerhouse and the main transformer chamber. We should pay attention to and strengthen the support strength at the intersection of the arch, the side wall and the opening. During the construction stage, timely support measures should be taken to strengthen random support to prevent large damage and deformation.
5. Conclusions
In this paper, a 3D finite element model including all caverns and a 3D simplified model containing only part of the caverns such as the main powerhouse and the main transformer chamber are established respectively, and the numerical analysis of the underground cavern group is carried out by using the nonlinear finite element method. In addition, a three-dimensional simplified model is established and three excavation calculation schemes are designed for comparison, and the following conclusions are obtained.(1)The numerical results using the model of all caverns show that the distribution law of the rock displacement field, stress field and plastic zone around the cave is similar to the results of the numerical model used in the stability analysis of the traditional underground cavern group.(2)The stress, displacement, and plastic zone of the model containing all the caverns are larger than the model that only contains the main plant, main transformer room, bus tunnel, etc. The maximum horizontal displacement difference of the plant area is 53.72%, and the maximum deviation of the third stress is 39.57%. The maximum difference in the depth of the plastic zone in the workshop area is 30.77%. And the position of the extreme value of stress and displacement also changed.(3)After the excavation of the auxiliary cavern is completed, there will be places where stress is concentrated. The stress concentration parts should be strengthened to prevent cracking damage. Thereby reducing the probability of rockburst and collapse events and ensuring construction safety.(4)The excavation of the small caverns around the powerhouse will weaken the integrity of the surrounding rock and have a certain impact on the deformation of the surrounding rock. Three-dimensional finite element models of all caverns.
Based on the above conclusions, it is reasonable to use the numerical simulation analysis of the 3D finite element model including all the caverns. Future work should consider supporting situations on this basis.
Data Availability
The codes used in this paper are available from the author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This study was funded by the National Natural Science Foundation of China (grant no. 51579089). Data curation, C. S. and H. Z.; Formal analysis, S. H. and Y. X.; Methodology, H. Z. and S. H.; Writing-original draft, S. H. and R. Y.; Writing-review & editing, Y. X. and S. H.; Supervision, R. Y. and E. C. All authors have read and agreed to the published version of the manuscript.