Abstract
The stability of underground cavities is of increasing importance considering the predominant cavity locations built up in high mountain and canyon environments. Such cavity locations are characterized by a high initial in situ stress, which results in brittle fracture and deformation of the surrounding rock during cavity construction. This paper presents a numerical study of underground cavern stability considering four factors, namely, mechanical property of surrounding rock, cavern burial depth, lateral pressure coefficient in horizontal direction, and the angle included between plant longitudinal axis and horizontal principal stress. Analytical methods including the key point displacement in side wall, plastic zone volume, and splitting fracture volume are used to characterize the stability of underground cavern. A modified formula to predict side wall displacement is proposed based on prior work, which is applicable to 3D computation model by taking horizontal geostress in two directions into account. Eventually, the optimal layout of underground cavern is put forward under different conditions of geostress field.
1. Introduction
The large-scale hydropower projects regarding underground cavities are predominantly built in high mountain and canyon with tectonic complexity. In such cases, the excavation of surrounding rock mass during cavity construction typically gives rise to brittle fracture and rock burst under the condition of high in situ stress field. In a series of ongoing large-scale hydropower projects, oblique splitting and the development of large fissures in the high side walls have been extensively observed, which seriously undermines the stability of surrounding rock mass. As a consequence, to ensure the safety and security of underground projects, the optimization and refinement in the design of underground cavities are in great need for an improved understanding of splitting fracture.
In the design of underground caverns, it is critical to determine the layout of plant longitudinal axis, which plays a pivotal role in the stress distribution and stability of surrounding rock mass. The basic principle for determining the plant longitudinal axis is to choose the smaller angle included between the longitudinal axis and principal stress and the bigger angle included between longitudinal axis and fissure strike. However, in view of the increasing geostress for underground hydropower, it becomes important to elucidate whether the routine approach to determine the plant layout is most favorable or not. As a result, to characterize the stability of surrounding rock under multiple factors in particular the longitudinal axis of plant layout is of great practical significance and engineering application value.
This paper takes Shuang-Jiangkou, a hydropower station in southwest China, as example to characterize the stability of surrounding rock mass using numerical method FLAC3D. Four main factors are considered, that is, mechanical properties of surrounding rock mass, cavern burial depth, lateral pressure coefficient in horizontal direction, and the angle included between plant longitudinal axis and horizontal principal stress. The key point displacements of side walls and plastic zone volume are compared with criteria proposed to judge splitting fracture [1, 2].
A modified equation, based on key point displacement forecasting method, is proposed to enable the displacement prediction of side wall in 3D condition. The angle included between plant longitudinal axis and horizontal principal stress is yielded by the minimum principle of side wall displacement, plastic zone volume, and splitting zone volume. The most favorable layout of cavern longitudinal axis is acquired by the method mentioned above under different conditions of in situ stress.
2. Numerical Model
The underground cavities at Shuang-Jiangkou are chosen for analysis and five 3D finite difference models are built on the basis of angle included between the cavity axis and loading direction (Figure 1). In order to reduce calculation error, the models are divided into core area and peripheral area with the same amount of nodes and elements. The core area in each model has the same mesh while the peripheral area is meshed differently with the same amount of nodes and elements. By doing so, only the core area is influenced by the excavation unloading effect. The peripheral area merely plays a role in altering the principal stress direction. Hence, the effect of mesh difference in peripheral area is minimal. The boundary condition is applied to the peripheral area of the model. The bottom of the model is fixed and the top of the model has free boundary. The left and right sides of the model are applied by normal constraint. Horizontal stress is applied into the model to simulate the horizontal tectonic stress and vertical load is applied to simulate the influence of gravity field. The constitutive model adopted in this paper could be referred to in Appendix.
2.1. Model Scope
The underground cavities at Shuang-Jiangkou hydropower station consist of the main power house, the main transformer chamber, and tailrace surge chamber among others. As for the dimension of the underground cavern, the vertical length is 265 m and the horizontal length is 220 m. The excavated height of the main power house is 67.05 m and the excavated height of the main transformer chamber is 26.5 m. The layout of the excavations is tabulated in Table 1 and the excavation sequence of cavern groups is shown in Figure 2.
For simplicity, only main power house, main transformer chamber, and tailrace surge chamber are considered in the computational model. The overall model comprises 102048 nodes and 96444 elements, with 75858 elements from the core area (model shown in Figure 3).
3. Parameter Selection
Four factors are considered, that is, mechanical parameters of surrounding rock [3], cavern burial depth [4], horizontal pressure coefficient [5], and the angle included between the longitudinal axis of cavern and horizontal principal stress [6]. The stability of surrounding rock in underground caverns is thus investigated under multiple factors as mentioned above by means of numerical study. The selection of such four factors is deemed reasonable on the basis of scientific merit and wide application in engineering project.
The mechanical parameters of underground cavern in large-scale hydropower station in China are summarized in Table 2.
3.1. Mechanical Parameters of Surrounding Rock
See Table 3.
3.2. Horizontal Pressure Coefficient Selection
Initial geostress field could be achieved by the cavern burial depth with only gravitational stress field considered. However, it has been shown by extensive literature survey that the horizontal pressure coefficient [7] is usually larger than 1 due to the existence of tectonic stress. Hence, the lateral pressure coefficient is chosen as 1.0, 1.5, 2.0, and 2.5. The stability principle for surrounding rock is that the layout of plant longitudinal axis and horizontal principal stress is in small angle intersection. Thus, the chosen lateral pressure coefficients for two horizontal directions are tabulated in Table 4 according to the principle mentioned above.
3.3. Angle () Included between Longitudinal Axis of the Plant and Horizontal Principal Stress
Based on the principle with small intersection angle between longitudinal axis of the plant and horizontal principal stress, the angle is chosen as 0°, 10°, 20°, 30°, and 40°.
To reduce the computation error due to the geometry and meshing difference, the computational model is divided into two sections, the core section and peripheral section. The geometry and meshing are kept the same for the core area while the dimension and elements number are subject to the value of , the angle included between longitudinal axis of the plant and horizontal principal stress. In the initial stage of geostress, the normal stress is exerted on the peripheral section until equilibrium and then each surface of the model is confined to define the geostress. Figures 4 and 5 are the cross sections of computational models with the value of to be 20° and 0°, respectively.
4. Stability Approaches of Surrounding Rock
4.1. Key Point Displacement in Side Wall
This approach is to investigate the effect of multiple geostress factors on surrounding rock stability by measuring the relative displacement of key point in side wall [8]. The key point in main power house is shown in Figure 6. This point is chosen as representative since the downstream side wall of the main power house is most prone to suffer splitting fracture. Extensive study [9] has suggested that the displacement of the chosen point tends to be maximum in its vicinity.
4.2. Plastic Zone Volume
In this study, the element is defined as damaged once the deformation of the element becomes plastic [7]. The extent of damage for the surrounding rock is judged by the total volume of plastic zone due to excavation unloading. The main power house is taken as the main target since it is the most important construction in the underground cavern assembly. A stability index is proposed in the light of plastic zone volume of surrounding rock () and the exaction volume () for the main power house. It gives
5. Splitting Fracture Mechanism Revisiting
A splitting model of linear slippage crack groups has been given with a splitting criterion of cavern groups proposed [9]. A brief discussion of tensional crack deformation mechanism is given as shown below.
Rock removal by underground excavation gives rise to stress redistribution and concentration. The excavation alerts the compression state from triaxial condition to biaxial or even uniaxial condition. As a result, the elastic strain energy is dissipated by crack initiation, propagation, and coalescence, which finally leads to tensional failure as illustrated in Figure 7.
Due to the complexity of stress-strain state during the excavation, it is not suitable to evaluate the splitting fracture simply by the variation of stress or strain. Herein, energy principle considering crack propagation is adopted to elucidate the splitting deformation.
5.1. Splitting Model
A splitting model has been introduced by Martin and Chandler [10]. Figure 8 shows a linear sliding crack model, which was first utilized to study the effect of rock dilatancy. The crack body consists of an initial crack and wing crack. The wing crack is caused by the relative slippage of initial crack surface under compressive stress in far field.
Initial cracks are assumed to be distributed with even density and crack propagation is synchronous during loading. This simplified model is taken as the research object of fracture mechanism [11]. The dimension of initial cracks has been shown in Figure 8 and details about these cracks can be found elsewhere [12].
5.2. Splitting Criterion
The effective shear stress, , generated through and in a single linear sliding crack can be expressed as
Note that ;
Then, the stress intensity factor of linear sliding crack is
Note that c is half of the length of initial crack. For a single crack, the tensile force exerted on wing crack can be defined as the horizontal projection of effective shear force.
The resultant tensile force (F) across the whole crack surface is equivalent towhere A is the surface area of the tensional crack.
Substituting (2) and (4) into (5) gives
Crack coalescence is influenced by the stress intensity factor of the crack tip (). Thus, which results in the formation of a tensional crack can be expressed as
Since crack friction is not taken into account, the stress concentration caused by II-type crack can be ignored. So, with inequality , the inequality between and can be obtained:
Finally, stress inequality is produced after tensional failure occurs:
Inequality (9) can function as a fracture criterion for the surrounding rock, where is friction coefficient of crack while L is average length of ultimately formed splitting crack and h is the angle between initial crack and horizontal principal stress.
As for the splitting fracture analysis of surrounding rock, the stability of underground cavern excavation is computed by FLAC3D and then the splitting fracture criterion is used on the basis of simulation results. The distribution of splitting fracture zone and splitting depth are thus obtained by judging every element in the numerical model of underground cavern.
6. Numerical Results
Three indices for cavern stability, that is, key point displacement of side wall, plastic damage zone, and splitting failure zone, are evaluated to discriminate the effect of multiple geostress factors on surrounding rock stability in underground caverns. A modified formula to predict side wall displacement is proposed on the basis of work [13], which is applicable to 3D computation model and takes horizontal geostress in two directions into account. Eventually, the optimal layout of underground cavern is put forward under different conditions of geostress field [14].
6.1. Effect of Burial Depth on Surrounding Rock Stability
The burial depth is the main characteristic of geostress, which determines the value of vertical geostress. Normally, the cavern stability becomes adverse alongside the increase of burial depth.
Figure 9 exhibits an approximately proportional relationship between the key point displacement and burial depth of cavern under varying lateral pressure coefficient. It can be seen from Figures 10 and 11 that lateral pressure coefficient is the most significant factor in determining the plastic zone volume [13]. The plastic zone volume increases as the lateral pressure coefficient becomes larger. Moreover, the plastic zone is diminished with the increased burial depth when lateral pressure coefficient is larger than 2. The shallow burial depth results in smaller plastic zone when the lateral pressure coefficient is smaller than 1.5.
As indicated in Figure 12, the splitting volume is increased alongside the growth of burial depth. The smaller difference between two lateral pressure coefficients incurs smaller variation of splitting zone.
6.2. The Effect of Lateral Pressure Coefficient on Stability of Cavern Side Wall
Lateral pressure coefficient is an important indicator in elucidating the horizontal geostress. Four simulation cases are chosen as representative when the angle included between the longitudinal axis and the maximum principal stress in horizontal direction is zero. The effect of two-way lateral pressure coefficient is examined to investigate the surrounding rock stability of side wall in cavern. The parameters of four simulation cases are shown in Table 5.
6.2.1. Key Point Displacement in Side Wall
Figures 13 and 14 show the relationship between key point displacement in side wall and lateral pressure coefficient in four cases, respectively. Note that (1.0–2.5) refers to the lateral pressure coefficient 1.0 in x direction and the lateral pressure coefficients 1.0, 1.5, 2.0, and 2.5 in z direction, respectively. (1.0–2.5) z2.5 refers to the lateral pressure coefficient 1.0 in z direction and the lateral pressure coefficients 1.0, 1.5, 2.0, and 2.5 in direction, respectively.
It can be seen that the displacement in side wall increases with the growth of lateral pressure coefficient in x direction when the lateral pressure coefficient in z direction is kept constant and vice versa. The linear relationship is observed between the key point displacement and lateral pressure coefficient (x direction) with approximately quadratic function shown in lateral pressure coefficient (z direction). Nevertheless, the effect of lateral pressure coefficient in x direction is more pronounced since the slope of the curve is sharper under varying lateral coefficient pressure in x direction. This is attributable to distinctive excavation unloading in x direction since the horizontal stress is perpendicular to the side wall. However, the stress in z direction is parallel to the side wall with less leverage on the excavation unloading. Notably, the effect of lateral pressure ratio in z direction is marginal when it is less than 1.5, which indicates that only the lateral pressure coefficient in x direction could be considered [15].
As proposed in [13] by means of quasi-3D calculation, a displacement forecasting method was given to predict the key point displacement in side wall. It giveswhere is the height of the main power house (m); is the horizontal lateral pressure coefficient; is the volume weight ; is the burial depth of the cavern (m); is the deformation modulus (MPa); , , and are regression coefficients relevant to the geometric structure of cavern groups.
It should be noted that (10) is only applicable in quasi-3D state, which only takes the lateral pressure coefficient in x direction into account. Thus, the applicability of (10) is limited when 3D geostress is considered.
As suggested by above mentioned analysis, it may be construed that the key point displacement is not only related to the lateral pressure coefficient in both directions but also pertinent to the difference of lateral pressure coefficient in both directions. As a result, the modified lateral pressure coefficient λ is given below with considered:where is the lateral pressure coefficient in x direction. is the lateral pressure coefficient in z direction. is a coefficient pertinent to the difference between lateral pressure coefficients in both directions.
The approach to determine d is described as shown below. Firstly, choose different values of d. Then, (11) is used to determine the optimal value of d by least square method according to the relationship between and . Figures 15 and 16 are the plots of and for I and II surrounding rock, respectively.
As indicated in Figures 15 and 16, d is related to the rock property [16]. The value of d is assigned for the rock of same property to determine . Table 6 shows the value of four parameters obtained from regression. The fitted equation is expressed:
6.2.2. Plastic Zone
Figures 17–20 show the relationship between the lateral pressure coefficient and plastic volume. As indicated in Figures 17–20, the lateral stress in z direction plays a pivotal role in the size of plastic zone when the lateral stress in x direction is minimal [17, 18]. The influence of lateral stress in x direction becomes increasingly important alongside the growth of stress in x direction.
6.2.3. Splitting Zone
Figures 21 and 22 demonstrate that the splitting zone volume diminishes as the lateral pressure coefficient in z direction is increased. Nevertheless, the splitting zone volume is minimum with the lateral pressure coefficient 2.0 in x direction [17]. Before the excavation of the cavern, the increase of lateral pressure in x direction is likely to result in the decrease of splitting zone. However, due to the unloading from sidewall excavation, stress concentration is incurred since the surrounding rock is subject to tension by lateral pressure in x. The crack resulted from shear failure because of geostress. As the increase of horizontal pressure occurs, the failure pattern resulting in the crack formation is transformed from shear to tension. Within this process, the area of splitting zone is decreased until a threshold value is obtained. As suggested from Figures 21 and 22, the threshold value is believed to take place with the lateral pressure coefficient to be 2. When the pressure coefficient is increased to 2.5, the tension failure resulting in crack formation in the surrounding rock of sidewall begins to prevail. Therefore, the splitting zone volume increases correspondingly.
6.3. Effect of on Side Wall Stability
6.3.1. Key Point Displacement in Side Wall
It is revealed from Figure 23 that the displacement of side wall increases with the increase of when the lateral pressure coefficient is 1.0 in x direction. As the difference of lateral pressure coefficient between x and z directions is larger, the effect of angle included on the side wall displacement becomes greater. The minimum side wall displacement occurs when is 0°. As for the lateral pressure coefficient to be 1.5, the minimum side wall displacement happens under 10°. The optimum condition (minimum displacement for side wall) is observed under 20° with lateral pressure coefficient to be 2.0. The trend of case 2, of which the burial depth is increased by 500 m, is in line with that of case 1 (Figure 24) [19].
The mechanical properties adopted in case 3 (Figure 25) are higher than those in cases 1 and 2. Due to the fact that the deformation modulus in case 3 is double that in cases 1 and 2, the side wall displacement is diminished by 50%. The inflection point in case 3 is transferred to 10° with the lateral pressure coefficient being 2.0 as compared to the first two cases.
Overall, the optimum layout orientation () is primarily relevant to the lateral pressure coefficient in x direction under the criteria of minimum side wall displacement.
6.3.2. Splitting Zone
Figure 26 clearly shows that the splitting zone volume becomes smaller with increased lateral pressure coefficient. The optimum condition emerges at 0° or 40° under the same lateral pressure. The trend of splitting zone volume in Figure 27 is found to be the same as shown in Figure 26 although the burial depth in case 2 is increased by 500 m. The splitting zone volume in Figure 29 is relatively small presumably due to the enhancement of cohesive force and friction angle as compared to that in cases 1 and 2. However, the variation of splitting zone volume under lateral pressure coefficient and is roughly the same as the first two cases.
It may be concluded that the optimum condition emerges at 0° or 40°. The influence of on splitting zone volume is more distinctive under the condition of lower lateral pressure coefficient or less difference of lateral pressure between x and z directions. The splitting zone volume is reduced by 20% compared with the worst condition whereas the splitting zone volume is reduced by 8% under higher lateral pressure.
6.4. Optimum for Underground Cavern Stability
The variation of key point displacement, plastic zone volume, and splitting zone volume versus (Figure 28) is regressed by quadratic curve. The most favorable for the stability of underground cavern is summarized in Table 7 on the basis of respective criteria.
7. Conclusion
This paper presents a numerical analysis of underground cavern stability and considers four main parameters, that is, mechanical properties of surrounding rock mass, cavern depth, horizontal lateral pressure ratio, and angle included between longitudinal axis of the plant and horizontal principal stress. The following conclusions could be drawn from the current work.(1)The effect of burial depth: The displacement of sidewall increases as the burial depth increases. The plastic zone is associated with the lateral pressure coefficient. When the lateral pressure coefficient is larger than 2, the increase of burial depth could result in a decrease of plastic zone. When the lateral pressure coefficient is less than 1.5, the shallow burial depth is favorable for a small plastic zone. As for the splitting zone, the volume increases with the increase of burial depth. Besides, a smaller difference between the lateral pressure coefficients in two directions indicates a smaller splitting volume zone.(2)The effect of lateral pressure coefficient: Comparative analysis of simulation cases suggests that the key point displacement in the side wall is not only related to the lateral pressure coefficient but also related to the difference between the lateral pressure coefficients in two directions. As such, a modified formula to predict key point displacement is proposed under 3D stress condition. It is concluded that the geostress in z direction has a bigger influence in the plastic zone volume when the horizontal geostress in x direction is relatively small. Furthermore, the geostress in x direction has a bigger influence in plastic zone volume when the geostress in x direction is relatively large.(3)According to the minimum principle of key point displacement in side wall, plastic zone volume, and splitting zone volume, it is discerned that the most reasonable layout orientation for the surrounding rock stability is the angle included between the cavern axis and the maximum principal stress in horizontal direction.
Appendix
A. Constitutive Models
The synergic action of Zienkiewicz-Pande shear and tension failure yield criteria is adopted in this paper. Specifically, when the material is under elastic regime, Hooke’s law is applied. As the stress state meets the yield criteria, plastic deformation occurs. When Zienkiewicz-Pande shear criteria or tension failure criteria are met, a modified stress is proposed using a flow rule of relevant stressing condition.
A.1. Elastic Increment Equation
When the material is under elastic regime, Hooke’s law in terms of the increment of stress and strain gives
A.2. Synergic Yield Criteria and Flow Rule
The synergic yield criteria include the Zienkiewicz-Pande shear and tension failure principle. In FLAC3D, compressive stress is defined as negative. Hereby, the three principal stresses are expressed as
The Zienkiewicz-Pande shear failure criterion giveswhere ; indicates the relationship between Lode’s angle and the yield curve in plane π; .
The yield curve in meridional plane is expressed by hyperbola
The asymptotic line is expressed by Mohr-Coulomb envelope, which is situated in the plane of (Figure 29).
The tension failure criterion is given as shown below:
When plastic deformation occurs, the flow rules relevant to Zienkiewicz-Pande shear and tension failure principle are adopted. Hence, plastic potential functions of shear and tension failure are and , respectively.
The equation to identify shear region and tension region in the plane of iswhere
As shown in Figure 30, when , , the material has shear plastic deformation; when , , the material has tension plastic deformation.
Competing Interests
The authors confirm that there is no financial interest with any funding body except those mentioned below regarding the subject matter or materials presented in this manuscript. The funding mentioned below does not lead to any conflict of interests regarding the publication of this manuscript.
Acknowledgments
The study is jointly supported by grants from the National Natural Science Foundation of China (Grants nos. 41372294 and 51579142), Program for Changjiang Scholars and Innovative Research Team in University of China (no. IRT13075), and Opening Project Fund of State Key Laboratory of Mining Disaster Prevention and Control by Shandong Province and the Ministry of Science and Technology (no. MDPC201605) and China Scholarship Council.