Shock and Vibration

Volume 2014 (2014), Article ID 859648, 13 pages

http://dx.doi.org/10.1155/2014/859648

Research Article

Analysis of Seismic Damage of Underground Powerhouse Structure of Hydropower Plants Based on Dynamic Contact Force Method

Yang Yang,1,2 Juntao Chen,1,2 and Ming Xiao1,2

¹State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China

²Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering, Wuhan University, Ministry of Education, Wuhan 430072, China

Received 10 March 2014; Revised 4 August 2014; Accepted 16 August 2014; Published 2 September 2014

Academic Editor: Longjun Dong

Copyright © 2014 Yang Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Based on the characteristics of the dynamic interaction between an underground powerhouse concrete structure and its surrounding rock in a hydropower plant, an algorithm of dynamic contact force was proposed. This algorithm enables the simulation of three states of contact surface under dynamic loads, namely, cohesive contact, sliding contact, and separation. It is suitable for the numerical analysis of the dynamic response of the large and complex contact system consisting of underground powerhouse concrete structure and the surrounding rock. This algorithm and a 3D plastic-damage model were implemented in a dynamic computing platform, SUCED, to analyze the dynamic characteristics of the underground powerhouse structure of Yingxiuwan Hydropower Plant. By comparing the numerical results and postearthquake investigations, it was concluded that the amplitude and duration of seismic waves were the external factors causing seismic damage of the underground powerhouse structure, and the spatial variations in structural properties were the internal factors. The existence of rock mass surrounding the underground powerhouse was vital to the seismic stability of the structure. This work provides the theoretical basis for the anti-seismic design of underground powerhouse structures.

1. Introduction

The southwest region is the key area for hydropower development in China during the past few decades. A number of large-scale underground powerhouses of hydropower plants have been built in this region. It is also an earthquake-prone zone, with the seismic intensity usually above VII. Therefore, the underground powerhouses should have sufficient earthquake resistance capability which is vital to ensure normal operation of powerhouse and the safety of powerhouse personnel. Postearthquake investigations of the underground powerhouses in the epicentral region of Wenchuan earthquake, such as in Yingxiuwan, Yuzixi, and Gengda hydropower plants, reveal that underground powerhouses have, in general, stronger earthquake-resistant capability compared to their surface counterparts. Surrounding rock of the three underground powerhouses was generally stable, while the powerhouse concrete structure experienced obvious local cracking and failure. Concrete structure behaved as a weak link to affect the earthquake resistance capability of the underground powerhouses. Thus, it is very important to study the characteristics of seismic damage of underground powerhouse concrete structure.

It has been recognized that numerical method is an increasingly popular and effective strategy to solve the problem of structural dynamic response. A variety of works [1–5] have been devoted to the study of the dynamic response of underground powerhouse using the numerical method. Even though great achievement has been made in the field, the research on dynamic characteristic of powerhouse concrete structure is still lacking. From the microscopic point of view, the seismic damage and failure of underground powerhouse concrete structure are governed microscopic crack generation and propagation. In this way, the continuum method [6–8] (such as the finite element method, the finite difference method, and the boundary element method) and the discrete particle simulation method [9–22] (such as the discrete element method, the lattice-solid method, and the contact dynamics method) are the natural choice to simulate the damage and failure process. Among them, the finite element method is by far the most common and useful numerical method. The integration scheme it employs makes it appropriate for the widest variety of geologic and structural problems, and it can handle the most sophisticated constitutive relationships [23–25].

Underground powerhouse could be regarded as a system composed of powerhouse concrete structure and surrounding rock. Surrounding rock provides the external support for powerhouse concrete structure. As the concrete structure comes into direct contact with surrounding rock, seismic wave is propagated through rock mass to the powerhouse concrete structure. Therefore, the complex consisting of powerhouse concrete structure and surrounding rock undergoes a forced vibration. The simulation and analysis of the dynamic contact surface between surrounding rock and concrete structure is the key to the dynamic calculation of underground powerhouse concrete structure using the finite element method.

At present, the numerical calculation methods for dynamic contact problems are primarily the Lagrange multiplier method [26], penalty method [27, 28], and their improved versions [29, 30]. However, these methods tend to either increase the degree of freedom of the system or influence the time step of integration [31, 32]. These methods adversely affect the precision and speed of calculation when applied to the analysis of an underground powerhouse concrete structure, which involves a large number of contact elements and complex contact states. Liu et al. [33, 34] put forward the dynamic contact force method targeted at the dynamic response problem of the contact crack. The convergence and stability of this algorithm were easy to meet, making it suitable for large and complex contact systems. But it cannot reflect the bond-slip properties of the contact surface.

In this paper, a new dynamic contact force method, considering the bond-slip properties of the contact surface, is suggested using the fundamental integration formulation of the dynamic contact force method. The algorithm of the method considers the cohesive effect of the contact surface between concrete structure and the surrounding rock. It is capable of simulating the large slip phenomenon of the contact surface under dynamic loads. Based on the proposed algorithm, a finite element model considering dynamic constitutive properties of materials is built for the dynamic numerical analysis of the underground powerhouse structure of Yingxiuwan Hydropower Plant. Then the characteristics of seismic damage of the underground powerhouse structure under dynamic loads are studied by means of the numerical analysis and postearthquake investigations. The results will provide a theoretical basis for the antiseismic design of underground powerhouse.

2. Dynamic Contact Force Method considering the Bond-Slip Properties of Contact Surface

The contact model of underground powerhouse concrete structure and surrounding rock is shown in Figure 1. Before application of the dynamic loading, the nodes on the contact surface between concrete and surrounding rock belong to point-to-point contact. A certain amount of cohesive force exists between the nodes, which are in cohesive contact state. During dynamic loading process, the stress in some contact nodes would exceed the cohesive force and enter into the state of sliding contact or even separation. When a large relative sliding occurs between these contact nodes, they would come into contact with surfaces of the adjacent elements and then belong to point-to-surface contact.

Figure 1: Contact model of concrete structure and surrounding rock.

2.1. Fundamental Integration Formulation of the Dynamic Contact Force Method

According to the dynamic contact force method proposed by Liu et al. [33, 34], the problem of dynamic response of structure containing interface is discretized using the finite element method. Then, the differential equations of the system could be obtained as follows:

where , , and are the mass, damping, and stiffness matrix, respectively. is the displacement vector. is the velocity vector. is the acceleration vector. is the known vector of external forces. is the dynamic contact force vector.

The central difference method is used to solve the differential equation at a time . The time domain integral equation of the displacement and velocity of contact nodes containing the term of dynamic contact force could be obtained as follows:

where is the time step.

depends on the state of motion not only at time but also at time . Therefore, and cannot be obtained directly by using (2)~(4). In order to obtain the motion state of the contact node at , the contact force should be solved according to the contact conditions.

2.2. The Solution of Dynamic Contact Force under Point-to-Point Contact Condition

If relative sliding does not occur between contact nodes, or the relative sliding is small, the contact nodes are in or approximately in point-to-point contact condition. As shown in Figure 2, the node () on the contact surface of concrete and the node () on the contact surface of the surrounding rock are in point-to-point contact at time .

Figure 2: Relationship of point-to-point contact.

Suppose the contact node pair at time is in cohesive contact state; then in the normal and tangential direction, the contact nodes pair should meet the nonintrusive condition and the displacement compatibility condition with no relative sliding, respectively, as follows:

where is the unit normal vector of contact node.

Equation (2) is substituted into (5). According to the principle that a pair of dynamic contact force is equal in magnitude but opposite in direction, that is, , we have

where , are the normal and tangential components of , and , , , and satisfy the following equations, respectively:

In the above equations, and were obtained from the analysis of motion of the node pair. Therefore, they must satisfy the following inequalities.

(a)

If the value of does not satisfy (12), then the node pair would enter into the state of sliding contact. We have

(b)

If the value of does not satisfy (14), then the node pair entered into the separation state. We have

where and are the coefficients of static friction and kinetic friction, respectively. is the control area of the contact node to be calculated. is the cohesive force between the contact surface. Before time , if sliding or separation has not occurred between the node pair, , otherwise, .

2.3. The Solution of Dynamic Contact Force under Point-to-Surface Contact Condition

If larger relative sliding has occurred between contact nodes under dynamic loads, the contact nodes will be in the state of point-to-surface contact. Then there are no cohesive forces between the contact nodes and surfaces. As shown in Figure 3, at time , one node on the contact surface of concrete structure (or the surrounding rock) comes into contact with the surface of the surrounding rock (or concrete structure). Such a contact point on the contact surface is denoted as .

Figure 3: Relationship of point-to-surface contact.

It is assumed that node is in the state of cohesive contact with the corresponding contact surface at time . In the normal and tangential direction, the node should meet the nonintrusion condition and the displacement compatibility condition with no relative sliding, respectively, represented by (5). The displacements of the contact point at time and time are, respectively, given by

where is the shape function. is the node number of the contact surface.

Substituting (16) and (2) into (5), we have

where

In the above equations, two conditions should be discussed.

(a)If , node is separated from the corresponding contact surface. , can be computed from (15).(b)If , node is in contact with the corresponding contact surface. It is assumed that (17) will result in for -number of nodes. Thus, (17) and (18) of the nodes can be represented as follows:

 where is the coefficient matrix, relating , , and .

Solving (20) and (21), the contact forces and can be obtained for the nodes in point-to-surface contact condition. If one node has , then the node enters into the state of sliding contact. For this particular node, can be computed from (13), and (18) of the node should be removed from (21). Solving the new equation (21), the tangential contact force of other nodes can be known.

After solving and of all the contact nodes from the above equations, the total dynamic contact force can be obtained as follows:

Then, the displacement and velocity of each contact node for the next time step can be computed from (2) and (3). Figure 4 presents the flow chart for the proposed method.

Figure 4: The flow chart for the dynamic contact force method.

3. Dynamic Constitutive Model of Concrete and Rock Mass

Under cyclic loading, the unloading stiffness of concrete and rock mass at the later yielding stage is lower than the stiffness at the initial linear stage. The plastic-damage model, proposed by Lubliner et al. [35], improved by Lee and Fenves [36, 37], and generalized for 3D model by Omidi and Lotfi [38], could effectively simulate such a phenomenon. The model is suitable for the dynamic analysis of the quasibrittle materials, such as concrete and rock mass [3].

According to the basic theory of plastic-damage model, the plastic-damage stress-strain relationship of rock mass or concrete can be expressed as follows:

where is the stress tensor. is the effective stress tensor. is the initial stiffness of the material. is the strain tensor. is the plastic strain tensor. is the damage coefficient.

The structural damage is the result of the microcracks of the material. Under cyclic loading, the opening and closure of microcracks may happen, making the damage as a complex mechanism. When the state of stress especially changes from tensile to compressive, the stiffness weakened by the damage begins to recover. In order to simulate this phenomenon, the damage coefficient can be written as follows:

where and are tensile and compressive damage coefficients, respectively. and are the dimensionless constants as the functions of plastic strain. () is the coefficient of restitution when the material shifts from tensile state to compressive state.

The yield function of the model in form of effective stress is given as follow:

where and are the dimensionless constants and is a constant variable. For more details one can consult Omidi and Lotfi [38]. and are the first and second invariants of the effective stress tensor. is the maximum effective principal stress.

Concrete and the surrounding rock are used as friction material. The nonassociated flow rule can simulate the volume expansion properties under the compressive state. Therefore, the plastic-damage model employs nonassociated flow rule. The plastic potential function follows the Drucker-Prager hyperbolic function as follows:

where is the parameter which controls the potential function to approach the asymptote, is the dilatancy parameter, and is the maximum uniaxial tensile strength of the material.

An implicit cylindric anchor bar element method is adopted for the finite element implementation of the steel bars embedded in rock or concrete structures. The embedded steel bars are considered to improve the stiffness of the rock or concrete structures in this model. Therefore the stiffness of the steel bars can be superimposed onto the stiffness of rock or concrete element during numerical simulation. This method is detailed in [39].

4. Seismic Analysis of Underground Powerhouse Structure of Yingxiuwan Hydropower Plant

The dynamic analysis of underground powerhouse structure was performed with an in-house FEM numerical simulation platform, SUCED [40], which was designed for estimating seismic damage of large underground caverns group. The platform is based on the dynamic time-history method and uses explicit central difference method to solve the finite element equation.

4.1. Postearthquake Investigations of Underground Powerhouse Structure of Yingxiuwan Hydropower Plant

Yingxiuwan Hydropower Plant is one of the nine cascade hydropower plants of Minjiang River. This plant has an underground powerhouse, with three generating units. The seismic intensity of Yingxiuwan Hydropower Plant is up to XI, according to the distribution of seismic intensity of Wenchuan earthquake. Located just 8 km away from the epicenter, Yingxiuwan is one of the hydropower plants closest to the epicenter. A survey performed after the earthquake revealed some important findings. (a) Existing cracks on the arch of the underground powerhouse structure had propagated. A number of cracks appeared along the river and along the axis of the powerhouse after the earthquake (as shown in Figure 5(a)). (b) Part of sidewall lining cracked, but the crack width and the length were limited (as shown in Figure 5(b)). (c) The floor of turbine and generator layer showed closed cracks. The overlying floor was uplifted (as shown in Figure 5(c)). (d) The concrete of the inside wall of generator pedestal experienced a large number of cracks (as shown in Figure 5(d)).

(a)

(b)

(c)

(d)

Figure 5: Destruction of powerhouse concrete structure.

4.2. Numerical Calculation Model

Three-dimensional finite element model including surrounding rock and concrete structure of the powerhouse was established. The model consisted of 189,187 nodes and 171,432 hexahedron elements (eight nodes). The powerhouse is 52.8 m long, 17.0 m wide, and 37.2 m high. The model ranges along , , and -axes are 96.0, 88.2, and 147.7 m, respectively. In order to guarantee the accuracy of dynamic calculation results, the maximum mesh size should be less than 7.5 m based on the cut-off frequency of seismic waves. As the maximum mesh size of the model is 7.0 m, the demand of FEM dynamic simulation can be satisfied. The profile of the model is shown in Figure 6.

Figure 6: Profile of finite element model.

4.3. Mechanical Parameters of Material

The mechanical parameters of the rock mass and concrete material are shown in Table 1. The 3D plastic-damage model discussed in Section 3 was used for the numerical simulation of the surrounding rock and concrete structure.

Table 1: Mechanical parameter of rock mass and concrete.

Material
Deformation modulus (GPa)
Poisson ratio
Cohesion (MPa)
Friction angle (o)
Tensile strength (MPa)

Rock
10
0.25
2.18
41.6
1.97

Concrete (lining and inner structure)
25
0.17
2.0
46
1.30

4.4. Boundary Conditions

Free field boundary condition was applied to absorb the reflection waves along the four vertical boundaries. Viscoelastic artificial boundary condition was used to absorb the incident waves from the bottom of the model [41].

4.5. Dynamic Loads

Wolong was the nearest strong motion station to the epicenter to record strong ground motion from the Wenchuan earthquake. Its epicentral distance was 19 km. The seismic waves recorded in Wolong were concluded more representative than those recorded in other stations. The 20 s to 80 s sections of Wolong monitored acceleration data, which have high intensity and amplitude, were used as the seismic input load. The input acceleration time-histories used in the calculation model are shown in Figure 7.

(a) -direction

(b) -direction

(c) -direction

Figure 7: Input acceleration time-histories of calculation model.

4.6. Numerical Calculation Results

Static analysis of the surrounding rock excavation and concrete structure of the underground powerhouse was performed before performing a dynamic analysis. The results provided the initial conditions for a dynamic analysis.

A formula as introduced in (27) was used to calculate the damage volume of all damaged elements of underground powerhouse concrete structure:

where is the total damaged volume of the concrete structure. is the tensile or compressive damaged volume of each single element, and is the number of damaged elements.

4.6.1. Analysis of Evolving Process of Seismic Damage

The time history of the total damaged volume of powerhouse structure was plotted as shown in Figure 8. As can be seen, the total damaged volume was 63.2 m³ before seismic loading. Seismic damage on the powerhouse structure initiated when  s. From  s to 20 s and from  s to 40 s, the seismic waves entered into the two peaks. The total damaged volume increased sharply. At other times, the amplitudes of the seismic waves were small. So the total damaged volume increased slowly. When  s, the total damaged volume was about 460.3 m³.

Figure 8: Time history of the variation of the seismic damaged volume of the powerhouse structure.

Figures 9 and 10 show the plot of the damage coefficient when  s and  s, respectively. Before seismic loading, a small amount of damage zones distributed in the upper part of the lining. The damaged coefficient did not exceed 0.3. When the seismic loading was completed, a large amount of damage zones appeared in the lining. The damaged areas were mainly distributed in the arch and the upper sidewall. The damage coefficient at some locations was observed as high as 1.0. The concrete showed the risk of cracking failure.

Figure 9: Damage coefficient distribution of the lining when  s.

Figure 10: Damage coefficient distribution of the lining when  s.

The damage coefficient distribution of the inner concrete structure of the powerhouse when  s and  s was plotted (Figures 11 and 12). As can be seen, when the seismic loading was completed, a wide range of damaged areas appeared at the floor of turbine and generator layer and the inside wall of generator pedestal. The damage coefficient was close to 1.0 at most of the damaged area. The concrete showed the risk of cracking failure. The damaged areas were also noticed at the crane beam corbel and the column of turbine layer. However, the damage coefficient was not large.

Figure 11: Damage coefficient distribution of the inner concrete structure when  s.

Figure 12: Damage coefficient distribution of the inner concrete structure when  s.

4.6.2. Analysis of the Sliding and Separation of Contact Surface

The distribution of the sliding and separation zone of the contact surface between concrete structure and surrounding rock when  s and  s was plotted as shown in Figures 13 and 14. When  s, a small amount of sliding or separation zones occurred on the lining. When  s, the sliding or separation zones of the contact surface were greatly expanded, which mainly occurred on the arch, the junction of the sidewall with the arch, and the junction of caverns. This is consistent with the distribution of the damaged areas of the lining structure.

Figure 13: Distribution of sliding and separation zone when  s.

Figure 14: Distribution of sliding and separation zone when  s.

5. Analysis of Seismic Damage Characteristics of Underground Powerhouse Structure

5.1. Comparison between Numerical Calculation and Postearthquake Investigation

The numerical calculation results showed that the damaged areas of underground powerhouse concrete structure mainly occurred on the arch, the upper sidewall, the floor of generator and turbine layers, the inside wall of generator pedestal, and the junction of caverns. This is basically consistent with the cracking failure areas noticed in the postearthquake investigations. Therefore, the dynamic calculation method for underground powerhouse concrete structure proposed in this paper is proved to be reasonable and effective. The results could truly reflect the damage characteristics of powerhouse structure and could provide theoretical basis for antiseismic design of such structures.

5.2. Influence of Seismic Waves on Seismic Damage of the Structure

The degree of damage of the underground powerhouse structure was closely related to the amplitude and duration time of seismic waves. The larger the amplitude and the longer the duration of seismic waves were, the more obvious the damage of the underground powerhouse structure was. In the numerical simulation of the underground powerhouse structure of Yingxiuwan Hydropower Plant, the distribution range of the damaged areas and the damage coefficient increased significantly at the time of the two obvious vibration processes of seismic waves. When the two obvious vibration processes were over, the amplitude of the subsequent seismic waves was smaller. However, due to the continued input of seismic waves, the extent of damaged areas and the damage coefficient increased to a certain degree. This shows that the amplitude and duration time of seismic waves determine the degree of the damage of an underground powerhouse structure. These are indeed the external factors causing the seismic damage of the underground powerhouse structure.

5.3. Influence of Structural Properties on Seismic Damage of the Structure

The structure of the underground powerhouse possesses significant spatial variation. The structure above the generator layer is mainly the lining, which is subjected only to the unidirectional constraint of the surrounding rock. The upper structure is relatively weak. Under seismic loads, the damaged areas were distributed extensively. The structure below the generator layer includes beams, plates, columns, and other massive concrete structures, which are subject to the constraints of surrounding rock on all sides. The proportion of the damaged areas was relatively smaller compared with that of the upper structure. Nevertheless, a large number of damaged areas occurred at the beams, plates, columns, and the inner wall of the generator pedestal, which have lower structural strength. It is apparent that the spatial variation of underground powerhouse structure led to the varying degrees of damage. This is an internal factor causing the seismic damage of the underground powerhouse structure.

5.4. Influence of Surrounding Rock on Seismic Damage of the Structure

The postearthquake investigations in many other plants showed that the characteristic seismic damage of surface powerhouses was predominantly the horizontal shear failure. These damages were serious and difficult to repair. On the other hand, the damage type of underground powerhouse structure was mainly the surface shedding and closed cracks. The degree of damage was generally lighter. The results of the numerical analysis showed that the damage of the powerhouse structure in the areas which could slide relative to or separate from the surrounding rock was more serious while that in the areas having good contact with the surrounding rock was lighter. Hence, the constraint of the surrounding rock remained influential to maintain the stability of the underground powerhouse structure under seismic loads. This influence of the surrounding rock is the main reason for the lesser damage of an underground powerhouse structure compared to its surface counterpart.

Under seismic loads, the surrounding rock deformed inward. The concrete structure was subjected to compression caused by the deformation of the surrounding rock. It was more pronounced in the lower part of underground powerhouse structure which would be subjected to constraints on all sides. For the underground powerhouse structure of Yingxiuwan Hydropower Plant, the floor of the generator layer was uplifted under the action of compression imposed by the surrounding rock. Therefore, the surrounding rock may have both advantages and disadvantages for underground powerhouse structure. However, in general, the advantages outweigh the disadvantages.

6. Discussion on the Antiseismic Measures of Underground Powerhouse Structure

According to the seismic damage characteristics, the constraint imposed by the surrounding rock remained as an important factor to the antiseismic stability of the underground structure. The powerhouse structure above the generator layer suffered from the problem of “underconstraint,” resulting in a large area of damage in the upper structure under seismic loads. For the structures below the generator layer, the problem of “overconstraint” made the action of compression caused by the surrounding rock deformation more obvious.

To ensure the antiseismic design, surrounding rock is recommended to be reinforced to improve its stability. In this way, the favorable constraint of the surrounding rock on the internal structure can be increased, while the unfavorable impact is minimized. Secondly, anchorage support or other support measures are recommended to strengthen the constraint of surrounding rock for the structure above the generator layer, so as to promote the antiseismic capacity of the upper structure. Finally, soft-seismic isolation layer can be added between the lower parts of the concrete structure and the surrounding rock. The soft-seismic isolation layer can not only weaken the compression caused by the surrounding rock deformation, but also reduce the seismic response of the structure to a certain extent. The enhancement in the antiseismic capability of the structure after incorporating the abovementioned measures was evident from the numerical simulation. For instance, Figures 15 and 16 show the distribution of the damaged areas in the powerhouse structure after adopting the abovementioned antiseismic measures. Significant improvement is seen in the structural performance when compared to the damage in Figures 10 and 12. The distribution range and the damage coefficient were greatly reduced. The stability of the powerhouse concrete structure was enhanced considerably.

Figure 15: Damage coefficient distribution of the lining when  s.

Figure 16: Damage coefficient distribution of the inner concrete structure when  s.

7. Conclusion

Based on the dynamic contact force method, which considers the bond-slip properties of the contact surface, a finite element model considering dynamic constitutive properties of materials is built for the dynamic numerical analysis of the underground powerhouse structure of Yingxiuwan Hydropower Plant. The numerical analysis results were compared against the findings from postearthquake investigations. The following conclusions can be made from this study.

(1)The dynamic contact force method fully considered the dynamic contact characteristics between the underground powerhouse concrete structure and the surrounding rock. It could simulate three contact states, namely, the cohesive contact, sliding contact, and separation of contact surface. It could be applied effectively in the dynamic calculation of large complex contact systems.(2)The numerical analysis and the postearthquake investigations revealed the evolution process and characteristics of the seismic damage of underground powerhouse structure. The amplitude and duration of seismic waves determined the degree of seismic damage. These are the external factors causing the damage of the underground powerhouse structure. The spatial variations of the structural properties led to the variation in the degree of damage. This is an internal factor causing the damage of the structure. The surrounding rock not only imposed favorable constraint but also caused unfavorable compression on the concrete structure. However, in general, the advantages outweigh the disadvantages.(3)In antiseismic design of the underground powerhouse structure, reinforced support system should be adopted to improve the stiffness of the surrounding rock and reduce the deformation. The upper part of the structure could be strengthened by anchorage support or other support measures to increase the constraining effect of the surrounding rock. The interstitial spaces between the lower parts of the concrete structure and the surrounding rock are recommended to be filled with soft-seismic isolation layer, so as to weaken the unfavorable compression due to surrounding rock and to reduce the seismic damage of the structure.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This study was supported by the National Key Basic Research Program of China (2010CB732005, 2015CB057900), the Major Program of the National Natural Science Foundation of China (91215301), the National Natural Science Foundation of China (51279136, 51209164), and the Research Fund for the Doctoral Program of Higher Education of China (20130141110015). These supports are greatly acknowledged and appreciated.

References

Z. Z. Shen and H. C. Ren, “Dynamic response characteristics of underground powerhouse caverns for Sandaowan hydropower station,” Advanced Materials Research, vol. 382, pp. 80–83, 2012.

Y. Zhang, M. Xiao, and J. Chen, “Seismic damage analysis of underground caverns subjected to strong earthquake and assessment of post-earthquake reinforcement effect,” Disaster Advances, vol. 3, no. 4, pp. 127–132, 2010.

B.-Y. Zhao and Z.-Y. Ma, “Influence of cavern spacing on the stability of large cavern groups in a hydraulic power station,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 3, pp. 506–513, 2009.

W. Q. Sun, Z. Y. Ma, X. Yan, J. Qi, and X. Du, “Intelligent identification of underground powerhouse of pumped-storage power plant,” Acta Mechanica Sinica, vol. 21, no. 2, pp. 187–191, 2005.

J. Chen, Z. Zhang, and M. Xiao, “Seismic response analysis of surrounding rock of underground powerhouse caverns under obliquely incident seismic waves,” Disaster Advances, vol. 5, no. 4, pp. 1160–1166, 2012.

O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann, 2005.

G. González, M. Gerbault, J. Martinod et al., “Crack formation on top of propagating reverse faults of the Chuculay Fault System, Northern Chile: insights from field data and numerical modelling,” Journal of Structural Geology, vol. 30, no. 6, pp. 791–808, 2008.

G. E. Hilley, I. Mynatt, and D. D. Pollard, “Structural geometry of Raplee Ridge monocline and thrust fault imaged using inverse Boundary Element Modeling and ALSM data,” Journal of Structural Geology, vol. 32, no. 1, pp. 45–58, 2010.

D. O. Potyondy and P. A. Cundall, “A bonded-particle model for rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 8, pp. 1329–1364, 2004.

M. Xia and K.-P. Zhou, “Particle simulation of the failure process of brittle rock under triaxial compression,” International Journal of Minerals, Metallurgy and Materials, vol. 17, no. 5, pp. 507–513, 2010.

M. Xia and C. B. Zhao, “Simulation of rock deformation and mechanical characteristics using clump parallel-bond models,” Journal of Central South University, vol. 21, no. 7, pp. 2885–2893, 2014.

M. Xia, C. B. Zhao, and B. E. Hobbs, “Particle simulation of thermally-induced rock damage with consideration of temperature-dependent elastic modulus and strength,” Computers and Geotechnics, vol. 55, pp. 461–473, 2014.

J. S. Yoon, A. Zang, and O. Stephansson, “Simulating fracture and friction of Aue granite under confined asymmetric compressive test using clumped particle model,” International Journal of Rock Mechanics and Mining Sciences, vol. 49, pp. 68–83, 2012.

Z. Zhao, L. Jing, and I. Neretnieks, “Particle mechanics model for the effects of shear on solute retardation coefficient in rock fractures,” International Journal of Rock Mechanics and Mining Sciences, vol. 52, pp. 92–102, 2012.

P. A. Cundall and O. D. L. Strack, “A discrete numerical model for granular assemblies,” Geotechnique, vol. 29, no. 1, pp. 47–65, 1979.

Z. H. Zhao, “Gouge particle evolution in a rock fracture undergoing shear: a microscopic DEM study,” Rock Mechanics and Rock Engineering, vol. 46, no. 6, pp. 1461–1479, 2013.

P. Mora and D. Place, “A lattice solid model for the non-linear dynamics of earthquakes,” International Journal of Modern Physics, vol. 6, pp. 1059–1074, 1993.

F. Radjaï and F. Dubois, Discrete-Element Modeling of Granular Materials, Wiley-I STE, New York, NY, USA, 2011.

L.-J. Dong and X.-B. Li, “Three-dimensional analytical solution of acoustic emission or microseismic source location under cube monitoring network,” Transactions of Nonferrous Metals Society of China, vol. 22, no. 12, pp. 3087–3094, 2012.

L. J. Dong and X. B. Li, “A microseismic/acoustic emission source location method using arrival times of PS waves for unknown velocity system,” International Journal of Distributed Sensor Networks, vol. 2013, Article ID 307489, 8 pages, 2013.

L. J. Dong, X. B. Li, and G. Xie, “An analytical solution for acoustic emission source location for known P wave velocity system,” Mathematical Problems in Engineering, vol. 2014, Article ID 290686, 6 pages, 2014.

X. B. Li and L. J. Dong, “An efficient closed-form solution for acoustic emission source location in three-dimensional structures,” AIP Advances, vol. 4, no. 2, Article ID 027110, 9 pages, 2014.

D. M. Potts and L. Zdravkovic, Finite Element Analysis in Geotechnical Engineering, Thomas Telford, 1999.

M. Sharafisafa and M. Nazem, “Application of the distinct element method and the extended finite element method in modelling cracks and coalescence in brittle materials,” Computational Materials Science, vol. 91, pp. 102–121, 2014.

G. G. Gray, J. K. Morgan, and P. F. Sanz, “Overview of continuum and particle dynamics methods for mechanical modeling of contractional geologic structures,” Journal of Structural Geology, vol. 59, pp. 19–36, 2014.

N. Hu, “A solution method for dynamic contact problems,” Computers and Structures, vol. 63, no. 6, pp. 1053–1063, 1997.

D. Peric and D. R. J. Owen, “Computational model for 3-D contact problems with friction based on the penalty method,” International Journal for Numerical Methods in Engineering, vol. 35, pp. 1289–1309, 1992.

Y. Kanto and G. Yagawa, “Dynamic contact buckling analysis by the penalty finite element method,” International Journal for Numerical Methods in Engineering, vol. 29, no. 4, pp. 755–774, 1990.

J. C. Simo and T. A. Laursen, “An augmented Lagrangian treatment of contact problems involving friction,” Computers & Structures, vol. 42, no. 1, pp. 97–116, 1992.

T. A. Laursen and V. Chawla, “Design of energy conserving algorithms for frictionless dynamic contact problems,” International Journal for Numerical Methods in Engineering, vol. 40, no. 5, pp. 863–886, 1997.

G. D. Pollock and A. K. Noor, “Sensitivity analysis of the contact/impact response of composite structures,” Computers and Structures, vol. 61, no. 2, pp. 251–269, 1996.

T. Belytschko, W. K. Liu, and B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, Chichester, UK, 2000.

J. Liu and S. K. Sharan, “Analysis of dynamic contact of cracks in viscoelastic media,” Computer Methods in Applied Mechanics and Engineering, vol. 121, no. 1–4, pp. 187–200, 1995.

J. Liu, S. Liu, and X. Du, “A method for the analysis of dynamic response of structure containing non-smooth contactable interfaces,” Acta Mechanica Sinica, vol. 15, no. 1, pp. 63–72, 1999.

J. Lubliner, J. Oliver, S. Oller, and E. Oñate, “A plastic-damage model for concrete,” International Journal of Solids and Structures, vol. 25, no. 3, pp. 299–326, 1989.

J. Lee and G. L. Fenves, “Plastic-damage model for cyclic loading of concrete structures,” Journal of Engineering Mechanics, vol. 124, no. 8, pp. 892–900, 1998.

J. Lee and G. L. Fenves, “A plastic-damage concrete model for earthquake analysis of dams,” Earthquake Engineering & Structural Dynamics, vol. 27, pp. 937–956, 1998.

O. Omidi and V. Lotfi, “Finite element analysis of concrete structures using plastic-damage model in 3-d implementation,” International Journal of Civil Engineering, vol. 8, no. 3, pp. 187–203, 2010.

M. Xiao, “3-D elastoplastic FEM analysis of implicit cylindric anchor bar element for underground opening,” Chinese Journal of Geotechnical Engineering, vol. 14, no. 5, pp. 19–26, 1992.

Z. G. Zhang, M. X. Xiao, and J. T. C. Chen, “Simulation of earthquake disaster process of large-scale underground caverns using three-dimensional dynamic finite element method,” Chinese Journal of Rock Mechanics and Engineering, vol. 30, no. 3, pp. 509–523, 2011.

Z. Zhang, J. Chen, and M. Xiao, “Artificial boundary setting for dynamic time-history analysis of deep buried underground caverns in earthquake di saster,” Disaster Advances, vol. 5, no. 4, pp. 1136–1142, 2012.