Research Article  Open Access
Analysis of Seismic Damage of Underground Powerhouse Structure of Hydropower Plants Based on Dynamic Contact Force Method
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 plasticdamage 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 antiseismic 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 largescale underground powerhouses of hydropower plants have been built in this region. It is also an earthquakeprone 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 earthquakeresistant 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 latticesolid 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 bondslip properties of the contact surface.
In this paper, a new dynamic contact force method, considering the bondslip 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 BondSlip 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 pointtopoint 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 pointtosurface contact.
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 PointtoPoint 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 pointtopoint 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 pointtopoint contact at time .
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 PointtoSurface Contact Condition
If larger relative sliding has occurred between contact nodes under dynamic loads, the contact nodes will be in the state of pointtosurface 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 .
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 pointtosurface 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.
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 plasticdamage 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 plasticdamage model, the plasticdamage stressstrain 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 plasticdamage model employs nonassociated flow rule. The plastic potential function follows the DruckerPrager 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 inhouse 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 timehistory 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)
4.2. Numerical Calculation Model
Threedimensional 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 cutoff 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.
4.3. Mechanical Parameters of Material
The mechanical parameters of the rock mass and concrete material are shown in Table 1. The 3D plasticdamage model discussed in Section 3 was used for the numerical simulation of the surrounding rock and concrete structure.

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 timehistories used in the calculation model are shown in Figure 7.
(a) direction
(b) direction
(c) direction
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^{3} 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^{3}.
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.
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.
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.
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, softseismic isolation layer can be added between the lower parts of the concrete structure and the surrounding rock. The softseismic 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.
7. Conclusion
Based on the dynamic contact force method, which considers the bondslip 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 softseismic 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. View at: Publisher Site  Google Scholar
 Y. Zhang, M. Xiao, and J. Chen, “Seismic damage analysis of underground caverns subjected to strong earthquake and assessment of postearthquake reinforcement effect,” Disaster Advances, vol. 3, no. 4, pp. 127–132, 2010. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 W. Q. Sun, Z. Y. Ma, X. Yan, J. Qi, and X. Du, “Intelligent identification of underground powerhouse of pumpedstorage power plant,” Acta Mechanica Sinica, vol. 21, no. 2, pp. 187–191, 2005. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
 O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, ButterworthHeinemann, 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 D. O. Potyondy and P. A. Cundall, “A bondedparticle model for rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 8, pp. 1329–1364, 2004. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 M. Xia and C. B. Zhao, “Simulation of rock deformation and mechanical characteristics using clump parallelbond models,” Journal of Central South University, vol. 21, no. 7, pp. 2885–2893, 2014. View at: Google Scholar
 M. Xia, C. B. Zhao, and B. E. Hobbs, “Particle simulation of thermallyinduced rock damage with consideration of temperaturedependent elastic modulus and strength,” Computers and Geotechnics, vol. 55, pp. 461–473, 2014. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 P. A. Cundall and O. D. L. Strack, “A discrete numerical model for granular assemblies,” Geotechnique, vol. 29, no. 1, pp. 47–65, 1979. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 P. Mora and D. Place, “A lattice solid model for the nonlinear dynamics of earthquakes,” International Journal of Modern Physics, vol. 6, pp. 1059–1074, 1993. View at: Google Scholar
 F. Radjaï and F. Dubois, DiscreteElement Modeling of Granular Materials, WileyI STE, New York, NY, USA, 2011.
 L.J. Dong and X.B. Li, “Threedimensional 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 X. B. Li and L. J. Dong, “An efficient closedform solution for acoustic emission source location in threedimensional structures,” AIP Advances, vol. 4, no. 2, Article ID 027110, 9 pages, 2014. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
 N. Hu, “A solution method for dynamic contact problems,” Computers and Structures, vol. 63, no. 6, pp. 1053–1063, 1997. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 D. Peric and D. R. J. Owen, “Computational model for 3D contact problems with friction based on the penalty method,” International Journal for Numerical Methods in Engineering, vol. 35, pp. 1289–1309, 1992. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 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. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 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. View at: Publisher Site  Google Scholar
 T. Belytschko, W. K. Liu, and B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, Chichester, UK, 2000. View at: MathSciNet
 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. View at: Publisher Site  Google Scholar
 J. Liu, S. Liu, and X. Du, “A method for the analysis of dynamic response of structure containing nonsmooth contactable interfaces,” Acta Mechanica Sinica, vol. 15, no. 1, pp. 63–72, 1999. View at: Publisher Site  Google Scholar
 J. Lubliner, J. Oliver, S. Oller, and E. Oñate, “A plasticdamage model for concrete,” International Journal of Solids and Structures, vol. 25, no. 3, pp. 299–326, 1989. View at: Publisher Site  Google Scholar
 J. Lee and G. L. Fenves, “Plasticdamage model for cyclic loading of concrete structures,” Journal of Engineering Mechanics, vol. 124, no. 8, pp. 892–900, 1998. View at: Publisher Site  Google Scholar
 J. Lee and G. L. Fenves, “A plasticdamage concrete model for earthquake analysis of dams,” Earthquake Engineering & Structural Dynamics, vol. 27, pp. 937–956, 1998. View at: Google Scholar
 O. Omidi and V. Lotfi, “Finite element analysis of concrete structures using plasticdamage model in 3d implementation,” International Journal of Civil Engineering, vol. 8, no. 3, pp. 187–203, 2010. View at: Google Scholar
 M. Xiao, “3D 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. View at: Google Scholar
 Z. G. Zhang, M. X. Xiao, and J. T. C. Chen, “Simulation of earthquake disaster process of largescale underground caverns using threedimensional dynamic finite element method,” Chinese Journal of Rock Mechanics and Engineering, vol. 30, no. 3, pp. 509–523, 2011. View at: Google Scholar
 Z. Zhang, J. Chen, and M. Xiao, “Artificial boundary setting for dynamic timehistory analysis of deep buried underground caverns in earthquake di saster,” Disaster Advances, vol. 5, no. 4, pp. 1136–1142, 2012. View at: Google Scholar
Copyright
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.