#### Abstract

The power of the Distinct Element Method (DEM) in solving the stability of underground powerhouse caverns (UPC) is discovered, but when does the DEM need to be adopted? What problems can be solved? The flowchart is provided and applied to analyze the stability of UPC in this paper. With the guide of the flowchart, the damage index (*Di*) is used as a failure type (gravity-controlled or stress-induced) judgment indicator. Through the calculation of three typical engineering, the problems of random blocks stability, dynamic calculation, and support system evaluation are studied, respectively, with the help of the DEM code 3DEC. The method and results of this paper can give reference to engineering projects of its category.

#### 1. Introduction

Vast underground caverns are used for a variety of purpose in civil engineering. Stability of underground openings under different conditions is an essential issue in construction [1]. Rock mass is largely Discontinuous, Anisotropic, Inhomogeneous, and Not Elastic (DAINE) [2]. Numerical simulation has become an important method to solve the rock mechanics problem of this kind of engineering, while how to choose a suitable numerical method is a crucial stage.

In the continuous mechanics method for the analysis of the stability of underground caverns, three-dimensional stability analyses and displacement predictions of large UPC were carried out by using FLAC^{3D} [3]. For simulating the hydraulic-mechanical interaction in the process of cracking, a coupling method which based on the elastoplastic finite element method (FEM) is proposed [4]. A hybrid intelligent method is proposed for a large cavern excavated in alternate hard and soft rock strata, and the method is an integration of an evolutionary neural network and FEM analysis using a genetic algorithm [5]. Detailed performance-monitoring studies have been carried out for determining the deformations and stress distribution around underground powerhouse caverns in the nonhomogeneous rock mass, using the three-dimensional finite element method [6]. The stability of a large cavern group at great depth is discussed by large-scale three-dimensional geomechanical model tests and numerical simulations of FLAC^{3D} [7].

In underground excavated in jointed rock masses at relatively shallow depth, the most common types of failure are those involving wedges falling from the roof or sliding out of the sidewalls of the openings [8]. If the stability of rock mass has a close relationship with the discontinuity, the discontinuous method is indispensable, and the continuum method will not get a reasonable answer. Considering modeling in engineering practice, the interest has been placed on the adoption of discontinuum models which give a far more realistic and representative picture of rock mass behavior than equivalent continuum models [9]. The distinct element method (DEM) [10] and discontinuous deformation analysis (DDA) [11] are better suited than the finite element, boundary element, and finite difference methods to perform discontinuum analysis of underground excavations in jointed rock masses [12]. It is thought that DDA is a special type of discrete element method [13]. The relative advantages and shortcomings of the DEM and DDA are compared [14]. In 1988, 3DEC was developed for Falconbridge Mines in Sudbury, Ontario, to study rock bursting activity in deep mines. In 1993, the 3DEC models were constructed to analyse three parallel caverns in shallow jointed rock, Viikinmäki wastewater treatment plant, Helsinki, and Finland, to compare the cavern displacements predicted in the highly fractured zone with those observed [15].

Stability analysis with 3DEC was, thus, conducted to examine each of the potential issues (stress-inducing brittle failure, large deformation, and block instability) with various scenarios corresponding to each of the layout proposals [16]. Triaxial creep tests and back analysis of the time-dependent behavior of Siah Bisheh cavern by 3DEC [17]. The displacement-based back analysis using univariate optimization algorithm was applied, and numerical modeling results of 3DEC were in good agreement with measured displacements [18]. An expanded distinct element method (EDEM) was developed for simulating the crack generation and propagation due to the shear and tension failures in the matrix rock blocks [19].

There is no doubt that the application of the DEM has been widely recognized by the industry, especially in underground rock engineering. However, for different complex engineering problems, systematic and multifaceted research of how to use 3DEC should be proposed.

In this paper, we firstly present a practical flowchart of rock mechanics modeling for the steps in the UPC evaluation process. It is intended that these flowcharts should cover the principles and different aspects of the problems. Based on the flowchart, this paper is to conduct stability analysis of the caverns with DEM numerical modeling and study the mechanical behavior of the rock mass at three typical projects (Jurong, Dagangshan, and Qiongzhong; Figure 1), The three projects are all underground powerhouse cavern projects in large-scale hydropower projects. Jurong and Qiongzhong are pumped storage power plants, while Dagangshan is a conventional hydropower project.

#### 2. General Flowchart of Modelling of UPC Using the DEM

The response of a rock mass is often dominated by discontinuities that cut through the rock because they are usually much weaker and more deformable than the intact rock. Cundall [10] described such a numerical model and applied it to the toppling failure of a rock slope. Subsequently, the method became known as the distinct element method or the discrete element method (DEM). The name “discrete element method” applies to a computer program only if it allows finite displacements and rotations of discrete bodies, including complete detachment; recognizes new contacts automatically as the calculation progresses. The term “distinct element method” was coined by Cundall and Strack [20] to refer to the discrete element scheme that uses deformable contacts and an explicit, time-domain solution of the original equations of motion (not the transformed, modal equations). It has been extended to deformable rock blocks and applied to such diverse systems as granular material, masonry structures, and hydraulic fracturing [21, 22].

Although rock discontinuities had previously been introduced (as specialized elements) into the finite element method [23], the DEM is different because it characterizes a joint as a nonlinear boundary condition, rather than an element, and it allows arbitrary displacement and rotation of rock blocks and unlimited freedom for any purpose to interact with any other object. Any viable DEM code requires an underlying process that continuously identifies pairs of neighboring blocks and the specific entities (corners, edges, and faces) that may interact between each pair. These tasks must execute in linear time for the code to be efficient for simulating systems with thousands of blocks and complex block shapes. One such scheme is described by Cundall [24], but there are other algorithms that also avoid the polynomial-time searches implied by a brute-force approach. The DEM is a powerful technique to perform stress analyses in blocky rock masses formed by discontinuities [12]. But, how is it used in plant design? When to use and what should be paid attention to when applying? There have been many earlier presentations on this subject: for example, the flowcharts developed by Hoek and Brown [25] and Hudson and Feng [26]. We provide an applicable flowchart in Figure 2; the work starts with a collection of geological data, and then, there is the division between types of instability—leading to the two rows in the flowchart considering instability due to adverse structural geology and excessively high rock stress. It, then, leads to decisions on the excavation processes, support evaluation, and operating maintenance of the site. It is noted that the flow chart is mainly aimed at how to choose suitable numerical simulation methods to solve different failure types of UPC projects.

Notes on the DEM evaluation strategy.(1)Requirement of structural designs and discontinuities data(2)Via precedent type analysis and consideration of the type of failures(3)Use of rock mass classifications, RMR, *Q*, GSI, and BQ to indicate required support(i)RMR: Bieniawski [27] published the details of a rock mass classification called the Geomechanics Classification or the Rock Mass Rating (RMR) system(ii)Q: Barton et al. [28] proposed a Tunneling Quality Index (*Q*) for the determination of rock mass characteristics and tunnel support requirements(iii)GSI: the Geological Strength Index (GSI), introduced by Hoek et al. [29](iv)BQ: basic quality index proposed by the Chinese standard for engineering rock mass classification-GB 50218-94, NMWR [30](4)Examination of earthquake disaster to the caverns(5)Use of numerical methods 3DEC to provide a detailed study

It provides the overall process as indicated in Figure 2. On this basis, we illustrate the application of the DEM in UPC from the following engineering examples.

#### 3. Structurally Controlled Failure Calculation

##### 3.1. Blocky Rock Mass and DFN

In the excavation of jointed rock masses at relatively shallow depth, the most common types of failure are those involving wedges falling from the roof or sliding out of the sidewalls of the openings. These wedges are formed by intersecting structural features, such as bedding planes and joints, which separate the rock mass into discrete but interlocked pieces. When the excavation of the opening creates a free face, the restraint from the surrounding rock is removed. One or more of these wedges can fall or slide from the surface if the bounding planes are continuous or rock bridges along the discontinuities are broken [8].

Martin et al. [31, 32] pointed out that failure of underground openings in hard rocks is a function of the in situ stress magnitudes and the characteristics of the rock mass. At low in situ stress magnitudes, the failure process is controlled by the continuity and distribution of the natural fractures in the rock mass. However, as in situ stress magnitudes increase, the failure process is dominated by new stress-induced fractures growing parallel to the excavation boundary. This fracturing is generally referred to as brittle failure. It is indicated that the initiation of brittle failure occurs when the damage index (*D*_{i}) expressed as the ratio of the maximum tangential boundary stress to the laboratory unconfined compressive strength exceeds 0.4. Under the environment of brittle and hardness rocks (GSI > 40), the gravity-induced and structurally controlled failure are easy to happen, although stress-induced brittle spalling can be found in blocky and massive rock mass when *Di* > 0.15.

Fractures or faults are often modeled deterministically. Actual existing faults that have been mapped are represented explicitly in the model by specifying a location, dip, and dip direction. It is possible to consider hundreds of individual faults in this manner, but the procedure is quite labor intensive. The other alternative is a stochastic representation of faults. With this approach, which is called the discrete fracture network (DFN), specific faults are not modeled explicitly. Instead, a set of statistical parameters are defined, and a joint set is generated based on the statistical input. In this way, the joints do not represent specific mapped joints, and the joint set is nonunique. Different realizations will satisfy the statistical criteria [33].

##### 3.2. Overview of the Jurong Project

The Jurong pumped storage power station project is located in Jurong City, Jiangsu Province, 65 km apart from Nanjing City (Figure 1). The station is a daily regulating pumped storage power station with a total generating capacity of 1350 MW, aiming at modulating peak, filling the valley, frequency modulation, phase modulation, and emergency standby for the power system. The underground powerhouse complex consists mainly of the main and auxiliary powerhouse, the main transformer chamber, the bus tunnel, the traffic tunnel of the main transformer chamber, the cable tunnel of the traffic system, outlet shaft of 500 kV systems, the traffic tunnel into the factory, the ventilation and safety tunnels, drainage gallery, and drainage tunnel (Figure 3). The main and auxiliary powerhouse and the main transformer chamber are arranged in parallel with each other with a distance of 40 m. The underground powerhouse complex is located in the Northeast of Lun Mountain peak mountain slopes, which is magnificent. The main caverns are excavated in the Dengying Group of the Sinian Period, Mufushan Group , and Paotaishan Group of the Cambrian Period. The diorite porphyrite veins invasion also can be found. The design of the Jurong powerhouse complex was carried out by the HydroChina Huadong Engineering Corporation [34].

**(a)**

**(b)**

**(c)**

**(d)**

##### 3.3. Jurong Project Failure Mode

The maximum in situ stress of the Jurong plant is about 7.0 MPa, and the in situ stress regression analysis shows that the maximum principal stress is between 5.0 and 7.0 MPa. The uniaxial compressive strength of the intact rock is about 75 MPa. The engineering geological types of the surrounding rock mass are mainly Class III from the classification of BQ (GSI > 40). It can be considered that the damage index (*D*_{i}) of the plant is generally between 0.07 and 0.10. Accordingly, the numerical results show that stress-induced instability is not a common phenomenon, which only occurs at the position of small scale with the ratio of more than 0.15. From Figure 2, the conclusion can be deduced that the failure type is that the block will slide along the discontinuities. In other words, after excavation, the failure of the rock mass is controlled by discontinuity mostly. When the excavation of the opening creates a free face, the restraint from the surrounding rock is removed. One or more of these wedges will fall or slide from the surface. The DEM is an appropriate choice to simulate the behavior of the engineering rock mass.

##### 3.4. Jurong Project DFN

DFN modeling is increasingly being used to help provide solutions for many geotechnical and mining engineering problems. The DFN approach can be defined as the analysis and modeling process that explicitly incorporates the geometry and properties of discrete fractures as a central component controlling rock mass behavior. With the DFN module of 3DEC, the fracture population embedded into a rock mass is viewed as a set of discrete, planar, and finite-size fractures. By default, the discrete fractures are disk-shaped. The DFN modeling approach is stochastic: the geometrical features of a DFN model are only constrained through independent statistical distributions of its geometrical properties.

The geometrical characteristics supported by the DFN module are the fracture size (diameter), orientation, and position distributions. The fracture size distribution is determined from a probability function that describes the distribution of disk diameters. In each case, the dip and dip direction are drawn from the specified distribution. From the DFN model, any number of different DFN realizations can be generated depending on a random seed. Numerous realizations of different DFN have been generated statistically. The statistical parameters associated with a DFN typically characterize the fracture size distribution, orientation distribution, and density of each fracture set [35]. The DFN used in this study consists of 3 fracture sets, each with a given statistical distribution (Table 1). Through the function of the DFN, we get the jointed rock model, as shown in Figures 3(a) and 3(b).

3DEC requires material properties for both the intact blocks and the discontinuities. Properties assigned to blocks are generally derived from laboratory testing programs (Table 2). Joint properties are conventionally derived from laboratory testing (e.g., triaxial and direct shear tests). The joint cohesion and friction angle correspond to the parameters in the Coulomb slip constitutive model (Table 3).

##### 3.5. Distributions of Different Displacement Blocks

To speed up the solution, a rigid blocks method is used for calculation. Based on the calculation result, we can divide the blocks into several types according to the displacement. It is a practical method to judge the stability through full displacement, although it is not perfect.

According to the absolute displacement, the surrounding rock mass of the Jurong underground powerhouse complex is generally divided into the following four categories: one is unstable blocks with a displacement greater than 5 cm, the second is fewer stable blocks with a displacement between 4 and 5 cm, the third is blocks with an obvious displacement within 3-4 cm, and the fourth is relatively stable or stable blocks with a displacement less than 3 cm. The results are shown in Figure 4, unstable blocks are mainly distributed in the main powerhouse, and the main powerhouse arch near the downstream side is the most dangerous position. It is necessary to fix a certain number of long anchor cables at the crown part of the main powerhouse cavern beside the systematic support. Also, there are a small number of unstable blocks in the upstream and downstream sidewalls of the main powerhouse, and the downstream blocks are more than the upstream. The volume of unstable blocks with a displacement greater than 5 cm is about 12 499.4 m^{3}, the volume of less stable blocks with displacement between 4 and 5 cm is about 2 844.4 m^{3}, and the volume of blocks with obvious displacement within 3-4 cm is about 3 941.5 m^{3}.

**(a)**

**(b)**

Simulation results indicate that the potential problem of stability will give priority to failure of structurally control. Also, the failure blocks have little relationship with the stress field after excavation. Figure 4 discovered, clearly, that instability blocks and latent instability blocks are not in accordance with the stress field distribution but have a close correspondence with the spatial location of the joints attitude, and the stress is not the primary control factor.

##### 3.6. Suggestion for the Design and Construction

Numerical analysis provides a useful tool to enhance understanding of the rock mass behavior after excavation. According to the absolute displacement, the surrounding rock mass of the Jurong complex and unstable blocks are mainly distributed around the main powerhouse, and the main powerhouse arch near the downstream side is the most dangerous position. It is necessary to fix a certain number of long anchor cables at the crown part of the main powerhouse cavern.

#### 4. Dynamic Calculation

##### 4.1. Principle of DEM Dynamic Calculation

The calculations performed in the DEM alternate between application of a force-displacement law at all contacts and Newton’s second law at all blocks. The force-displacement law is used to find contact forces from known (and fixed) displacements. Newton’s second law gives the motion of the blocks resulting from the known (and fixed) forces acting on them. If the blocks are deformable, the motion is calculated at the grid points of the triangular finite-strain elements within the blocks. Then, applying the block material constitutive relations gives new stresses within the elements.

The modeling of geomechanics problems involves media which, at the scale of the analysis, are better represented as unbounded. Numerical methods relying on the discretization of a finite region of space require that appropriate conditions be enforced at the artificial numerical boundaries. In dynamic problems, such boundary conditions cause the reflection of outward propagating waves back into the model and do not allow the necessary energy radiation. The alternative is to use quiet boundaries. Several formulations have been proposed. The viscous boundary developed by Lysmer and Kuhlemeyer [36] is used in 3DEC [33]. It is based on the use of independent dashpots in the normal and shear directions at the model boundaries. The free-field condition is applied to lateral boundary grid points. The lateral borders are coupled to the free-field grid by viscous dashpots to simulate a quiet boundary, and the unbalanced forces from the free-field zone are applied to the deformable-block boundary at the boundary grid points.

Damping is used to describe the dissipation of system energy when the research object is under the action of dynamic loading. Local damping is chosen in this case. Local damping was originally designed to equilibrate static simulations. However, it has some characteristics that make it attractive for dynamic simulations. It operates by adding or subtracting mass from a grid point or structural node at certain times during a cycle of oscillation; there is overall conservation of mass because the amount added is equal to the amount subtracted. Thus, the use of local damping (local damping ratio = 0.03) is simpler than Rayleigh damping because we do not need to specify a frequency.

##### 4.2. Overview of the Dagangshan Project

The Dagangshan hydropower station dam site is in the upper part of the middle reaches of the Dadu River in Shimian County, Sichuan Province (Figure 1). The upstream relates to the tailwater of the Yingliangbao power station, and downstream relates to the reservoir of the Longtoushi power station. The dam site is about 40 km from Shimian County which located in the downstream of the Dadu River and about 72 km from Luding County which located in the upstream. The main task of the Dagangshan hydropower station is to generate electricity, and its control area is 62 727 km^{2}, accounting for 81% of the total drainage area of the Dadu River. The installed capacity of the power station is 2600 MW, the maximum water head is 178.0 m, the minimum water head is 156.8 m, and the rated water head is 160.0 m. The powerhouse cavern, transformer chamber, and tailrace surge tank in UPC are arranged in parallel (Figure 5).

The centerline distance between the main building and tailrace surge tank is limited in 141.20 m according to the requirement of guaranteed regulation, the centerline distance of the three caverns is 72.30 and 68.90 m, and the longitudinal axis of the main powerhouse is N°55°E with a size of 226.58 × 27.30 × 74.30 m. The elevation of the generator floor, turbine, and plant roof is 916.60 m, 944.00 m, and 991.80 m, respectively, and the span of the main cavern above the crane beam is 30.80 m, with the span of 27.30 m below the crane beam. The dimensions of the transformer chamber are 144.00 × 18.80 × 25.60 m. The rock mass of the UPC area is mainly composed of grayish-white and reddish biotite monzonitic granite, and there is diabase, granite aplite, diorite, and other vein rocks interspersed in granite, among which the diabase veins are widely distributed. The width of vein rocks is generally 0.5 m to 10 m, and the maximum width is up to 26 m, mainly with a steep dip angle. The rock classifications surrounding underground powerhouse caverns are mostly Class II and III from the classification of BQ (GSI > 40). The fresh rock around the cavern is intact. The rock structure is massive and submassive, and the rock blocks are firmly embedded. Besides, there is good fracture development of the top arch in local areas.

The underground cavern complex is constructed at a depth of approximately 390.00–520.00 m below the ground surface. The measured maximum principal stress value is 11.37–19.28 MPa. The laboratory intact rock uniaxial compressive strength is about 110 MPa. The damage index (*D*_{i}) of the plant is generally between 0.10 and 0.18. In addition, there are some auxiliary projects including the tailrace tunnel and busbar tunnels which are not considered in the calculation for shortening the time of calculation. According to the design scheme, the construction of the underground cavern complex is simplified as 9 steps, as shown in Figure 5. The typical cross joints are distributed in a rectangular region that contains the underground cavern complex. It is 400 m in the *X*-direction and 200 m in the *Z*-direction of the rectangular area. Each cavern consists of a circular arch and a vertical wall. Two group of typical joints (N12-25°E/NW 70–80°and N15°W/SW 65–80°) are included in the model.

##### 4.3. Dynamic Calculation Conditions

After the earthquake of Wenchuan in 2008, an interim regulation has been issued by the Chinese authority (Hydropower and Water Resources Planning and Design General Institute) that a two-level seismic stability assessment shall be performed for the critical components of a large-scale hydropower plant. In this manner, the reference ground motions may be set into two levels, which, respectively, correspond to the seismic intensity of 63% exceedance probability during 50 years and 2% exceedance probability during 100 years . The two levels of ground motion are named as the term of operating basis earthquake (OBE) and safety evaluation earthquake (SEE), respectively. A basic earthquake intensity VIII is obtained, and the reference earthquake motion parameters are given, as shown in Table 4.

No seismic motion was recorded accurately in the history near the project site. Thus, a specified technique of artificial simulation of nonstationary seismic motion for a large-scale underground cavern complex was conducted. The simulated waveforms for OBE and SEE are shown in Figure 6. According to Chinese seismic codes, a discount factor of 2 is specified for the peak acceleration value in the seismic analysis of underground structures with an overburden depth greater than 50 m [37, 38]. Therefore, the peak values of OBE and SEE are set to be 168 and 279 gal in the following analyses, respectively. Also, the vertical component of the accelerograms is taken as two-thirds of the horizontal component. During the numerical simulation, the earthquake waveform was the input from the bottom of the model to simulate the propagation of the ground motion. As we mentioned earlier, the 3DEC can simulate the far-field rock mass through the viscous boundary and to eliminate the reflection effect of external waves. When 3DEC is used for seismic dynamic response analysis, the surrounding boundary of the model is set to viscous and free-field, and shear waves from the *X*- and *Z*-direction are applied to the bottom of the model.

**(a)**

**(b)**

##### 4.4. Dynamic Calculation Results

To study the stability of the UPC’s surrounding rock mass under seismic loading, the surrounding rock deformation is focused on under the OBE and SEE seismic condition. The corresponding numerical calculations are conducted, and the UPC’s surrounding rock deformation process diagram is obtained, which are shown in Figure 7. The figure shows that, at the beginning of the earthquake load, visible open cracks emerge in the right sidewall of the UPC, but the rock block does not fall off. Then, with the expansion of open cracks, obvious damages occur on the sidewall. Lately, the damaged area extends from the sidewall to the top arch gradually until a large area collapses eventually. The analysis of the deformation response of the UPC under dynamic earthquake load is a very complex problem, which is influenced by the model size, joint cutting shape, mechanical parameters of rock mass and joints, boundary conditions, and the effect of seismic loads.

##### 4.5. Suggestion for the Design and Construction

According to the simulations, the acceleration pulse is the primary cause of the destructive power associated with discontinuous effects. The failure mechanism of underground caverns is controlled by unfavorable discontinuities that are subjected to ground motions of different levels. The larger the PGA value, the larger the irrecoverable shear displacement in the geological discontinuity induced. This large shear displacement may induce cavern failure due to sliding or the formation of wedges with other intersecting discontinuities, thus threatening cavern stability. Under the action of two groups of discontinuity cutting, unstable blocks under earthquake conditions mainly appear in the main powerhouse cavern at this project. It should be noted that the support system is not considered in the analysis.

#### 5. Support System Calculation

##### 5.1. DEM Bolting Element

Rock bolts (anchors) have become a favorite technique for reinforcing rock masses all over the world [39, 40]. Rock bolts are installed to reinforce fractured rock mass by resisting dilation or shear movement along the fractures. Several rock bolts can create a reinforced zone of fractured rock mass to improve the self-supporting capacity of the rock. Rock bolts usually undergo tensile and shear loading in the field. When the bolted rock mass deforms, load transfer occurs between the bolt and the rock.

3DEC provides structural elements representing the rock support systems. Two types of the reinforcement model are provided in 3DEC: local and global. A local reinforcement model considers only the local effect of reinforcement where it passes through existing discontinuities. A global reinforcement model thinks the presence of the support along its entire length throughout the rock mass. From those, a global reinforcement cable element is selected in this study. Cable elements interact with the 3DEC grid through shear and normal coupling springs. Normal forces act perpendicular to the rock bolt element and shear forces parallel to the element. The coupling springs are nonlinear connectors that transfer forces and motion between the rock bolt elements and the grid at the rock bolt element nodes. The behavior of the shear coupling springs represents the shear behavior of grout while the normal coupling springs are primarily used to model the influence of the medium deformation. Cable elements in 3DEC allow the modeling of a shearing resistance along their length, as provided by the shear resistance between the grout and either the cable or the host medium. The cable is assumed to be divided into some segments of length with nodal points located at each segment end. The mass of each segment is lumped at the nodal points, as shown in Figure 8. Shearing resistance is represented by spring/slider connections between the structural nodes and the block zones in which the nodes are located.

**(a)**

**(b)**

**(c)**

##### 5.2. Overview of the Qiongzhong Project

The Qiongzhong pumped storage power station is located in Qiongzhong City, Hainan Province, 164 km apart from Haikou City and 213 km apart from Sanya City (Figure 1). The station is a daily regulating pumped storage power station with a total generating capacity of 600 MW, aiming at modulating peak, filling the valley, frequency modulation, phase modulation, and emergency standby for the power system. The underground powerhouse complex consists of the main and auxiliary powerhouse, the main transformer chamber, the bus tunnel, the traffic tunnel of the main transformer chamber, the cable tunnel of the traffic system, the traffic tunnel into the factory, the ventilation and safety tunnels, drainage gallery, and drainage tunnel. The main powerhouse and the main transformer chamber are arranged in parallel with each other with a distance of 40 m. The main and auxiliary powerhouse has the size of 136.5 × 24.0 × 54.0 m and the transfer chamber of 98.87 × 19.00 × 20.6 m (Figure 9(a)). The engineering geological types of the surrounding rock mass are mainly Class II_{2} form the classification of BQ (GSI > 40). It can be considered that the damage index (*Di*) of the plant is generally between 0.11 and 0.12. The support system is demonstrated in Figure 9(b), and Table 5 gives the material property values used for the rock supports in the numerical model. The design of the Qiongzhong pumped storage power station project was carried out by HydroChina Zhongnan Engineering Corporation Limited [41].

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.3. Calculation of the Unsupported and Supported Scheme

The excavation of the cavern group is divided into 8 steps overall. The calculation results show that the monitoring points of the typical sections in each cavern increase gradually with the excavation process during the excavation process. After the cavern excavation ends, the deformation of the surrounding rock reaches its maximum and eventually stabilizes. The deformation of the crown is small. The deformation of the dome and abutment tends to be stable after 2-3 excavation steps. It has little change on the deformation of the crown in the subsequent excavation steps. With the excavation of the cavern, the deformation of the wall rock gradually increases, and the variation of the deformation of the monitoring points along the wall synchronize well with the excavation steps. The change of the deformation of the upstream and downstream sidewalls of the caverns is not different.

After all the caverns are excavated, the crown of the main powerhouse cavern deforms vertically downwards, and the deformation of the sidewall to the free surface and the rebound of the floor also occur. The maximum total displacement occurs near the intersection of the downstream sidewall of the main powerhouse cavern and the busbar hole in 1^{#} Set (Figure 10(a)), with a deformation of 28.49 mm. It is due to the existence of two discontinuities in the area near the right side of the sidewall which penetrates the busbar holes. The spacing of two the discontinuities is small, and the directions are approximately parallel and towards the main powerhouse cavern.

**(a)**

**(b)**

According to the design support scheme, the installation position of the prestressed penetration anchor cable is the upper of and between the downstream sidewall of the main powerhouse cavern and the upstream sidewall of the transformer chamber (Figures 9(b) and 9(c)). The total number of the anchors is 40. However, in the actual construction process, due to the geological conditions and the complexity of the construction, the installation of 8 anchor cables has been canceled (Figure 9(d)). When considering the support system, the final cumulative deformation values of all the critical monitoring points in the excavation process are lower than those in the unsupported condition. It shows that the current support system also has a good control effect on the deformation and stability of the surrounding rock mass. The maximum total displacement also occurs near the intersection of the downstream sidewall of the main powerhouse cavern and the busbar hole in 1^{#} Set (Figure 10(b)), with a deformation of 25.85 mm of the actual construction condition.

In order to explain, clearly, the deformation characteristics and reinforcement effects of various parts of the cavern under the conditions of system support, Figure 11 summarizes the displacement without or with the support of the key monitoring points of the main powerhouse and transformer chamber. The maximum total displacement occurs in the vicinity of the intersection of the downstream sidewall and the busbar hole in the main powerhouse 1^{#} Set (Section 1-1) with the discontinuity passing through.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.4. Cable Axial Force Distribution

This section analyses the stress of the bolts and anchors imposed by the system anchorage and demonstrates the reliability of the current construction support under the control of the discontinuities.

Figures 12(a) and 12(b) show the overall stress characteristics of the bolts system after excavation. In general, the stress of the anchor system is in good agreement with the deformation of the surrounding rock. The greater the deformation of the main powerhouse cavern and transformer chamber is, the higher the corresponding anchor bolt force is. The maximum position of the anchor bolt is located at the connection of the downstream wall of the main building and the busbar hole, which is consistent with the more substantial deformation at this position. The load of the standard mortar bolt of the crown of the main powerhouse of 1^{#} Set is 0–50 kN. The prestressed bolts in the arch seat are subjected to a more significant force, which reaches 100–150 kN. The standard bolt load of the upstream sidewall is 0–75 kN, and the value of the downstream sidewall is also 0–75 kN. Close values of bolt load are found in the chamber cavern. Throughout the model, the bolts of maximum load are in the upstream sidewall (crossed by the discontinuities) at 125–200 kN, part of the prestressed bolt stress close to the design strength (200 kN).

**(a)**

**(b)**

According to the design support scheme, the main place of the prestressed cable is the central of the downstream sidewall of the main powerhouse and the central of the upstream sidewall of the transformer chamber (Figure 8(b)). In the actual construction process, due to the complexity of geological conditions and constructions, the installation of 8 cables is canceled. The design load of the cable is 2000 kN, and the initial tensile lock average 1700 kN applies to the prestress cables. Figures 13(a) and 13(b) show the overall force characteristics of the final cavern cable system after excavation is completed under actual support conditions. The cable with maximum stress is located at middle rock mass cut by the discontinuity, between the main powerhouse and the transformer chamber, between 1^{#} Set and 2^{#} Set, reaching 2000 kN. The cable force of 1^{#} Set is between 1700 and 2000 kN. The cable force of 2^{#} Set is between 1700 and 1950 kN. The cables force of 3^{#} Set is between 1700 and 1810 kN.

**(a)**

**(b)**

##### 5.5. Suggestions for the Construction

Figures 14(a) and 14(b) show that the current support system has a good control effect on the deformation and stability of the surrounding rock mass although 8 noninstalled anchor bolts exist. By analyzing the stress of bolt and cable after the final excavation, it is considered that the anchoring system is safe under the construction support scheme. The overall stress state of the anchors is normal. The load of the most prestressed anchor bolts increases small compared with the prestressing force of the prestressed bolt, and a few of them appear close to the yield state, mainly distributed close to the discontinuity, at the downstream sidewall of the main powerhouse cavern, especially near the 1^{#} Set section.

#### 6. Conclusions

Each numerical method has its own advantages and disadvantages. The suitability and applicability of a numerical method must be ascertained for each individual case and on the objective of the study. Under the guidance of the flowchart, this paper studies the blocks stability and dynamic and anchoring problems of UPC through three typical projects. Based on the ratio of maximum initial stress to uniaxial compressive strength, the failure type can be judged. When the underground caverns are excavated in the hard and brittle rock mass, gravity-induced instability is more likely to happen. So, we all choose 3DEC to analyze the different problems of the three projects. With the DFN module of 3DEC, the fracture population can be embedded into a rock mass model. Because of its advanced DFN modeling and solving methods, 3DEC is a powerful DEM tool for research of UPC, which can provide useful suggestions for the design, construction, and operation. Furthermore, experience in the application of the DEM to practical situations of more cases is required so that an understanding may be developed of the flowchart and to answer where, when, and how the DEM is best applied.

#### Data Availability

Some data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request. Direct requests for these materials may be made to the provider as indicated in the Acknowledgements.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The work reported in this paper was financially supported by the National Natural Science Foundation of China (NSFC) under the Contract no. 51428902.