Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 4106970, 14 pages
https://doi.org/10.1155/2017/4106970

Research Article

Seismic Response Analysis of Concrete Lining Structure in Large Underground Powerhouse

Xiaowei Wang,1,2 Juntao Chen,1,2 Ming Xiao,1,2 and Danqi Wu1,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, Ministry of Education, Wuhan University, Wuhan, China

Correspondence should be addressed to Juntao Chen; chenjt000@163.com

Received 24 March 2017; Accepted 7 June 2017; Published 17 July 2017

Academic Editor: Fabrizio Greco

Copyright © 2017 Xiaowei Wang 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 dynamic damage constitutive model of concrete material and seismic rock-lining structure interaction analysis method, the seismic response of lining structure in large underground powerhouse is studied in this paper. In order to describe strain rate dependence and fatigue damage of concrete material under cyclic loading, a dynamic constitutive model for concrete lining considering tension and shear anisotropic damage is presented, and the evolution equations of damage variables are derived. The proposed model is of simple form and can be programmed into finite element procedure easily. In order to describe seismic interaction characteristics of the surrounding rock and lining, an explicit dynamic contact analysis method considering bond and damage characteristics of contact face between the surrounding rock and lining is proposed, and this method can integrate directly without iteration. The proposed method is applied to seismic stability calculation of Yingxiuwan Underground Powerhouse, results reveal that the amplitude and duration of input seismic wave determine the damage degree of lining structure, the damage zone of lining structure is mainly distributed in its arch, and the contact face damage has great influence on the stability of the lining structure.

1. Introduction

More than a dozen large-scale underground powerhouses of hydropower stations are located in the high seismic area of Southwest China, and their stability under seismic action is essential for the safe operation of hydropower stations. Many scholars have done a lot of researches on the dynamic response of the underground powerhouse [1–4] and have made many great achievements. However, there are few studies on the dynamic response characteristics of the lining structure in underground powerhouse so far. The “Wenchuan” earthquake disaster investigation [5, 6] shows that lining structure is the weak part of the underground powerhouse, so it is significant to study its seismic response.

The seismic response analysis of the concrete lining structure in underground powerhouse mainly includes two aspects: one is the dynamic response of the concrete material and the other is seismic interaction of surrounding rock and lining structure. Dynamic constitutive model of concrete material is the basis for dynamic analysis of lining structure in underground powerhouse. Based on laboratory experiments or theoretical analysis, many scholars have established a variety of concrete dynamic constitutive models to study the dynamic response process of concrete material, mainly including dynamic elastoplastic constitutive model [7–9], nonlinear elastic dynamic constitutive model [10, 11], viscoplastic dynamic constitutive model [12, 13], and dynamic damage constitutive model [11, 14, 15], wherein the dynamic damage constitutive model is the most widely used. Because of the structural features of concrete material, minor defects are distributed inside concrete structure before loading and continue to develop under earthquake cyclic loading, so the failure process of concrete material is the process of damage accumulation. It is necessary to study the evolution process of internal defects of concrete material in order to fully understand its dynamic response characteristics [16], and this is why damage theory is widely used in the dynamic response analysis of concrete structures [17, 18].

Because the seismic interaction of surrounding rock and lining structure is of complex nonlinear characteristics, a large number of nonlinear calculations, as well as higher demands for both the seismic contact analysis model and solving algorithm, are needed. Therefore, it has been a difficult problem in the dynamic response analysis of lining structure in underground powerhouse. At present, the numerical method is mainly used to solve the dynamic contact problem, such as penalty function method [19], Lagrange multiplier method [20], numerical programming method [21], and contact element method [22]. However, these methods tend to have a large amount of computation, and the iteration does not usually converge. Therefore, it is difficult to be applied to the calculation of dynamic contact problems of large underground structures. Liu and Sharan [23] proposed a dynamic contact force method, which combines the contact constraints with the explicit center difference method. This method is widely used in large complex practical engineering calculation with better convergence and no iteration [2, 24] but it does not consider the bonding effects of the contact face itself. In fact, due to grouting effects, the lining structure is in strong bonding with the surrounding rock before the contact face cracking and the cohesion of contact face between lining structure and surrounding rock cannot be ignored.

Based on the above ideas, a dynamic constitutive model for concrete lining considering tension and shear anisotropic damage is presented in this paper, and the evolution equations of damage variables are derived. Combined with the proposed explicit dynamic contact analysis method considering bond and damage characteristics of contact face between lining and surrounding rock, a dynamic response analysis method of concrete lining structure in underground powerhouse is constructed. The method is applied to seismic stability calculation of Yingxiuwan Underground Powerhouse and used to study the seismic damage evolution process of concrete lining structure of underground powerhouse. The results provide references for seismic design of concrete lining in large underground powerhouse.

2. Dynamic Damage Constitutive Model for Concrete Lining

Based on a large number of experimental studies [25–28], it is found that the concrete material shows significant strain rate dependence and fatigue damage characteristics under dynamic cyclic loading. Its strain rate dependence is mainly reflected by the fact that its mechanical parameters increase with the increase of the strain rate, which has been proved by numerous experiments and theoretical studies [29–32]. Because concrete is a kind of heterogeneous material, its inherent micro cracks and voids continually develop and accumulate under seismic cyclic loading and finally cause damage on a macro level. In this process, the strength and stiffness of the concrete materials are continuously degraded with fatigue damage characteristics emerging. Based on the above facts, these two features of concrete material should be considered to build a dynamic constitutive model for concrete lining. Therefore, on the basis of damage theory, a concrete dynamic damage constitutive model considering the effects of strain rate is established in this section.

2.1. Dynamic Damage Model for Concrete Based on the Separation of Tension and Shear

The dynamic damage constitutive model for concrete needs to reflect the strain rate sensitivity as well as the generation and evolution of the damage; its general form can be written as

where means strain rate and means dynamic damage variable.

The rationality of the damage variable largely determines the accuracy of the damage model. The damage variable used in current engineering calculations is a scalar, which is simple and convenient for theoretical derivation and large-scale numerical calculation. However, as a kind of heterogeneous material, the damage evolution of concrete is directional under the action of load, especially the seismic cyclic load, so the measurement of damage should take the form of tensor [16]. Damage forms of concrete are mainly tensile damage and shear damage, and both are of different generation mechanisms, and so are their effects on strength and stiffness. Therefore, the dynamic damage variable can be divided into tensile damage and shear damage [33]:

where and are, respectively, tensile and shear damage variables, which are of the scalar form, and their specific evolution equations will be given in Section 2.2 and and are of fourth-order tensor form and are expressed as functions of the principal value of the stress tensor and the principal direction :

Combined with (2) and (3), the dynamic damage variable , which is of fourth-order tensor form, can be derived. According to the static constitutive relation of concrete, the dynamic form of concrete can be

where is dynamic Young’s modulus of concrete. In the absence of experimentation, it is advisable to take the static Young’s modulus; is fourth-order consistency tensor.

2.2. Evolution Equation of Tensile and Shear Damage Variables of Concrete

According to the derivation process above, in order to establish a complete dynamic damage model, it is essential to set up evolution equation of tensile and shear damage variables of concrete in dynamic action. There is difference between the tensile and shear damage evolution process under dynamic action and under static load, because damage evolution process under dynamic action relates not only to the strain but also to the strain rate. Therefore, the key point of the paper lies in setting up rate-dependent damage evolution equation.

Assuming the maximum tensile strain rule applying to the tensile damage of element in 3D (three-dimensional) stress state, article [34] substitutes the maximum tensile strain in three-dimensional stress state by equivalent strain, so the tensile damage evolution equation can be expressed as

where is the threshold of tensile strain damage and is the limit tensile strain, both considering strain rate, and is the residual strength coefficient, while is equivalent strain, which can be defined as

where , and are three principal strains, respectively. The function can be defined as

It is difficult to acquire the relation between the strain rate and the threshold of tensile strain damage and the limit tensile strain under dynamic action directly. According to a number of experiments [35], there is a high similarity between the static and the dynamic stress strain curves of concrete. According to the similarity,

where and are threshold of tensile strain damage and the limit tensile strain in static state respectively, while is dynamic strain amplification factor of concrete. Based on results of series of dynamic experiments on concrete, Euro-International Committee for Concrete proposed an equation:

where is static strain rate whose value is and is dynamic strain rate ranging from to .

It is obvious that (5) to (9) compose a complete evolution equation of tensile damage variable. Similarly, evolution equation of shear damage variable can be defined as

where is the maximum compressive principal strain and is limit compressive strain of element. When Mohr-Coulomb criterion applies to the stress state of the element, then we have

where and are dynamic friction angle and Poisson’s ratio, which are usually assigned with static value, and is dynamic uniaxial compressive strength of concrete. According to the similarity of the static and the dynamic stress strain curves of concrete, we have

where is static compressive strength of concrete material and is amplification coefficient of dynamic compressive strength. The proposed equation by Euro-International Committee for Concrete is

where .

Similarly, (10) to (13) compose a complete evolution equation of dynamic shear damage variable. Obviously, tensile and shear damage variables in this paper are set on the basis of maximum strain criterion and Mohr-Coulomb criterion considering strain rate effect, which is brief and easy for calculation in finite element method.

2.3. Finite Element Simulation of Lining Element

Steel bars are distributed in the concrete lining of underground powerhouse, so their reinforcement effects on the lining structure should be considered in the process of finite element calculation of concrete lining. In this paper, it is assumed that the steel bars are uniformly distributed in all of the lining elements, and the stiffness of the lining element is the superposition of the stiffness of the concrete element and the stiffness of the steel element; its element stiffness matrix can be expressed as

where is stress transformation matrix and is equivalent elastic matrix for steel bar element; its expression can be seen in literature [36]. is equivalent elastic matrix for concrete element, and, based on the derivation of Sections 2.1 and 2.2, it can be expressed as

where .

3. Dynamic Contact Analysis Model of Lining Structure and Surrounding Rock

The dynamic contact analysis between large underground structure and surrounding rock is of numerous contact elements as well as complex contact states, and, due to the consolidation grouting and contact grouting, the contact face between the underground structure and the surrounding rock is in bonding state before the bonding failure with certain tensile strength and shear strength. Therefore, the suitable dynamic contact model of underground structure and surrounding rock should not only be able to support large-scale nonlinear calculation, but also reflect their dynamic contact characteristics. Based on the above facts as well as time-domain explicit integration method, a dynamic contact analysis model of surrounding rock and lining is established in this section, which combines the contact conditions and bonding properties.

3.1. The Explicit Time-Domain Integration Method for the Motion Equations of Contact System

The contact system model of surrounding rock and lining structure and the dynamic contact forces on their contact face is shown in Figure 1. According to the idea of Shu et al. [37], the two sides of the contact face are divided into surface of the surrounding rock and surface of the lining structure, the correspondent nodes (as and in Figure 1) on surfaces and are, namely, node pairs (pairs of nodes), the two nodes of each node pair occupy a common local coordinate system and the same global coordinates, and dynamic forces on the two nodes satisfy the interaction relationship; that is, . The main content of this subsection is to derive the time-domain progressive integration scheme for the nodes of the contact system according to the equation of motion.

Figure 1: Contact model and dynamic contact force on contact face.

After the finite element discretization, the contact system model can be regarded as a multi-DOF system consisting of multiple particles. The dynamic equilibrium equation of the nodes of contact system in seismic action can be defined as

where is the known external forces vector, and represent the normal and tangent contact force vectors, respectively, , , and are acceleration vectors, velocity vectors, and displacement vectors, respectively, and , , and are matrixes for mass, damping, and stiffness, where mass matrix is expressed by lumped mass matrix [38] and damping matrix is expressed by Rayleigh damping matrix [39]:

where is the matrix for function (its value is 1 within the region of node and 0 out of the region) and and are damping constants in proportion to mass and stiffness, respectively.

Solving the motion equation by central difference method, the explicit integral equation of contact nodes motion equation at time can be derived as

where is time step.

According to (18) to (21), the motion state of nodes at time can be calculated by the displacement, velocity, and contact force of nodes at time ,  and and depend on the motion state at time as well as at time . Therefore, to acquire the movement state of nodes at time , it is essential to calculate contact force and based on contact conditions.

3.2. Dynamic Contact Boundary Condition and the Solution of Dynamic Contact Force

The key to dynamic contact analysis is to solve the dynamic contact force on the contact face. The core idea of the solution is to predict first and then modify the dynamic contact force, that is, assuming the node pairs on the contact face is in bond state firstly, so the displacement inequality constraints of contact face can be transformed into equality constraints, and then correcting the calculated contact force according to the force constraints.

Supposing that the surrounding rock and the lining structure are in bond contact state, then the contact nodes pair meets the constraint condition of no relative displacement in the normal and tangential direction:

where and are displacement vectors of contact nodes and in the global coordinate system, and (22) is also satisfied for other node pairs on contact face. The subscripts 1, 2, and 3 of   represent two tangential () components and one normal () component of nodal displacements along contact face (as shown in Figure 1). is the transformation matrix for the global coordinate system and local coordinate system of the contact node , and it is expressed as

Substituting (18) into (22), we can have

where are subvectors of tangential force along two directions , respectively, so we can have

Obviously, the normal and tangential contact forces in (24) and (25) are obtained on the assumption that the contact face between the surrounding rock and the lining structure is in bond state. However, the contact face is a weak part of the underground structure, and its damage under the earthquake cyclic load cannot be neglected. In fact, there are a variety of contact states between the surrounding rock and the lining structure under the action of seismic cyclic load, such as bonding contact, sliding contact, and separation. So it is necessary to check and correct the contact force of the contact face.

The main damage forms of contact face are the tensile cracking in the normal direction and shear slip in the tangential direction. This paper recalibrates the contact condition and corrects the contact force from the following aspects.

() If and , this indicates the contact face is in tension state and tensile cracks occurs on the contact face; the contact force should be corrected as

() If and , this means that contact face is in tension state but the tensile stress is smaller than tensile strength of the contact face; tangential contact force should be corrected as

() If , the contact face is in contact state, and calibrating tangential contact force is needed. Namely, if , which indicates there is shear slip on the contact face, the tangential contact force should be corrected as

where and are the cohesion and tensile strength of contact face, respectively. is the control area by node , while and are the static friction coefficient and kinetic friction coefficient, respectively. Note that when there is any tensile crack or shear slip on the contact face, correction of cohesion is needed. The contact forces and can be obtained by (24) to (28) and substituted into (18) to (21); the motion state of the contact nodes at time can be updated.

4. Project Example

4.1. Introduction

The Yingxiuwan Hydropower Station is located on the main stream of the Minjiang River, Yingxiu, Wenchuan, Sichuan Province. The underground powerhouse is buried at the depth of about 200 m, and its size is 52.8 m × 17.0 m × 37.2 m (length × width × height). The lining thickness of side wall and the arch in main powerhouse is 0.5 m and 0.9 m, respectively. The size of main transformed cavern is 59.4 m × 7.2 m × 27.9 m (length × width × height), and its lining is 0.5 m thick. The lining is made of C25 concrete. The powerhouse locates about 8 km from the epicenter of Wenchuan Earthquake, in the seismic zones of basic intensity VII. Because of the short distance to epicenter, the earthquake equals a prototype test for Yingxiuwan Underground Powerhouse. Therefore, the damage characteristics of its lining structure after the earthquake are representative. The proposed dynamic response analysis method of concrete lining is implemented into our research group’s dynamic finite element program SUCED [38] by a self-developed Fortran program, and, taking Yingxiuwan Underground Powerhouse as an engineering example, the seismic response characteristics of concrete lining structure are studied in this section.

4.2. Finite Element Model and Calculation Condition

The numerical model is composed of the main powerhouse and the main transformer cavern and discretized by hexahedron element with 8 nodes. There are 97652 elements in total, 7298 of which are lining elements. We define the longitudinal axis direction of main powerhouse and the vertical direction as the -axis and -axis in global coordinate system of calculation model, respectively, and the -axis is determined by the right-hand rule. The whole model is as shown in Figure 2; its range is , , and . Take plane as typical section, pick out 6 points on the plane as monitor points, and the arrangement is shown in Figure 3.

Figure 2: 3D finite element model.

Figure 3: Local finite element model of concrete structure and monitoring points.

Before dynamic analysis of lining structure, it is essential to simulate excavation and support of underground powerhouse; the initial geostress field is calculated by the weight of rock mass. Excavation calculation is implemented by 3D elastoplastic damage finite element method and the results are taken as the initial state of dynamic response analysis. The rock in the model area is mainly of class III, and mechanical parameters of materials used in dynamic calculation are shown in Table 1.

Table 1: Mechanical parameters of the rock and concrete.

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

Rock
10
0.25
2.18
42
1.97

Concrete
25
0.17
2.00
46
1.25

The free field artificial boundaries around the model are used to absorb the reflection waves, while the viscoelastic artificial boundary at the bottom is applied to absorb the incident waves. The seismic wave is the acceleration time history measured by Wolong Seismic Station near Wenchuan Earthquake, intercept 20 s to 40 s acceleration time history of it, which is of high amplitude and intensity, and then input seismic wave can be achieved after pretreatment [36]. The maximum acceleration of the input seismic wave is 4.22 m/ at , and its acceleration time history from three directions is shown in Figure 4.

(a)-direction

(b)-direction

(c)-direction

Figure 4: Acceleration time history curve of input seismic wave.

4.3. Dynamic Response Analysis of Underground Powerhouse Lining Structure

4.3.1. Time History Analysis of Stress and Displacement of Monitoring Points

Displacement time history curves of 6 monitoring points are shown in Figure 5. It can be seen that the forms and patterns of 6 curves resemble each other, with three obvious crests. And crests and troughs of 6 curves appear almost at the same time, which suggests all the monitoring points oscillating in synchrony. In 0 to 6 s, the displacements are small, and, after 6 s, the displacements fluctuate sharply, ranging from 1.5 to 5.5 cm. The three crests arise at around 10 s, 13.2 s, and 16.5 s, with the values of 4.8 cm, 5.4 cm, and 5.0 cm. After inputting seismic wave, the displacements increase from 0.5 to 1 cm more than before. The difference of the displacement time history curves lies in the amplitudes. Points B and C on side walls of main powerhouse have milder amplitudes, with a maximum relative displacement of 1.6 cm compared with Point A on the roof, indicating a large relative deformation of the lining in main powerhouse. The amplitudes of the three points in main transformer cavern are large, however with small differences between each other, indicating overall deformation of the lining in main transformer cavern.

Figure 5: Displacement time history curve of monitoring points.

The tensile strength of concrete is far less than the compressive strength, so tensile crack damage is common. Analyze the time history of the maximum principal stress of Points A, C, and D, whose tensile stress is positive, and the curves are shown in Figure 6. After excavation, the tensile stresses of the three monitoring points are 0.8 to 0.9 MPa. With seismic loading increasing, the amplitudes are small in 0~4 s, ranging from 0.8 to 1.0 MPa, and, in 4~15 s, the tensile stress has a sharp fluctuation, ranging from 0.45 to 1.64 MPa. The maximum tensile stresses of Points A and D are 1.64 MPa and 1.52 MPa, respectively, both over the tensile strength of concrete. In 15 to 20 s, the amplitudes start to decrease and finally reach a range of 0.8 to 1.0 MPa, still an increase compared with the initial stress. As shown in Figure 6, the time history curves of the maximum principal stress of points resemble the input seismic wave, and the fluctuation of tensile stress is mainly influenced by the input seismic wave. Points A and D have sharper amplitude compared with Point C, indicating the seismic wave has a greater influence on the arch of the lining, intensifying the tensile stress there and causing worse damage. It coincides with the fact that cracks were found on the arch of lining of Yingxiuwan Underground Powerhouse in earthquake damage investigation [31].

Figure 6: Maximum principal stress time history curve of monitoring points.

4.3.2. Damage Analysis of Lining Structure

In order to display the damage degree of lining more directly, this paper translates the damage tensor into a scalar damage coefficient , and is the average of principal values of damage tensor. Figure 7 shows the damage coefficient evolution process of lining structure under seismic action. Before the quake, namely, before (shown in Figure 7(a)), the damage areas scatter at the junctions of the arch and the side wall of lining of main powerhouse and at the junctions and the roof of lining of the main transformer cavern. The areas are relatively small. And the damage coefficient ranges from 0.1 to 0.2, with the maximum less than 0.3. According to Figure 5, the input seismic wave has a maximum amplitude of 4.22 m/s² at  s; the damage area and coefficient increase on a large scale at this time (shown in Figure 7(b)). The damage areas extend mainly to the arch and side wall of lining. Most areas at the arch of lining structure in main powerhouse are damaged in different degree, while the damage at the arch of lining structure in the main transformer cavern continues to develop, and damage even appears at some area of the side wall. There is serious damage at the junctions of the arch and the side wall, with the damage coefficient approximately equal to 1. The damage coefficient in the main transformer cavern ranges from 0 to 0.8, with maximum at the junctions.

(a)

(b)

(c)

Figure 7: Damage coefficient distribution of lining structure.

After inputting seismic wave, namely, after  s, the arch of lining is almost covered by the damaged area, while the damage area of the side wall also has a certain degree of expansion, and the damage coefficient is further improved (shown in Figure 7(c)). It can be seen from Figure 5 that the amplitude of the input seismic wave tends to decrease after  s, but the damage degree of the lining structure is further increased, which shows that the holding time of input seismic wave has a great influence on the seismic damage degree of the lining structure.

The tensile and shear damage distribution of lining structure after the earthquake is shown in Figure 8. The damage areas are mainly distributed at the arch and part of the side walls with a tendency to a further development. The tensile damage mainly distributes at the arch of the lining in main powerhouse, while shear damage mainly distributes at the arch and side walls of both the main powerhouse and the main transformer cavern. It can be seen from Figures 8 and 9 that the damage coefficient of tensile damage zone is relatively large, and damage situation is relatively serious, which means that measures should be taken to limit tensile damage in aseismic design.

Figure 8: Tension and shear damage distribution of lining structure.

Figure 9: Sliding and separation zone distribution of interface.

4.3.3. Slip and Crack Analysis of the Contact Face between Lining and the Surrounding Rock

The lining structure is adhered to the surrounding rock. If regarding the surrounding rock as a whole, it is obvious that the integrality of lining and the surrounding rock is not as good as that of rock itself. Therefore, in seismic cyclic action, it is difficult for the surrounding rock and lining structure to deform in synchrony, giving rise to slip and crack on their contact face, which is shown in Figure 9. According to the distribution of slip and cracks on the contact face, the damage areas are mainly at the arch and part of the side wall. Combined with Figures 8 and 9, it is obvious that there is an accordance of distribution between the slip and crack areas of contact face and the structure damage areas, suggesting that good contact state of the lining and surrounding rock is significant to the stability of the lining structure, and the constraint effects of the surrounding rock to the lining can efficiently reduce the damage of lining by seismic loading.

5. Conclusion

Based on dynamic damage characteristics of concrete lining structure in underground powerhouse, a dynamic constitutive model, which considers tensile and shear anisotropic damage, is established, and, combined with the proposed dynamic contact analysis model considering bond and damage of contact face between lining and the surrounding rock, a set of complete dynamic response analysis method of underground powerhouse lining structure is established. The method is used in the seismic damage analysis of lining structure in Yingxiuwan Underground Powerhouse; it concludes the following.

() Displacement time history curves and stress time history curves of monitoring points resemble each other. And all the points oscillate in synchrony. Compared with the side walls, the arch of lining in main powerhouse has relatively larger amplitudes of displacement and tensile stress. The tensile stress in the arch exceeds the tensile strength of concrete, with a maximum displacement difference of 1.6 cm compared to the side walls, which suggests a large relative deformation of lining structure in the main powerhouse. Displacements of points in the main transformer cavern are close to each other, suggesting a whole deformation of the main transformer cavern.

() Damage areas of lining structure are mainly distributed at the arch and parts of the side wall with a tendency to a further development, which coincides with the results of Yingxiuwan Underground Powerhouse in the earthquake damage investigation. Some areas at the arch of the lining in main powerhouse are damaged to the most severe degree. Tensile damage areas have a large damage coefficient with structure damaged severely. And measures should be taken to limit tensile damage of lining.

() During the vibration process with large seismic amplitude, the displacement and stress curves of the lining structure are obviously fluctuated, and the damage degree is significantly increased. After the peak acceleration of the input seismic wave, the amplitude of seismic wave tends to be gentle, but the damage degree of the lining structure continues to increase. These indicate that the damage degree of the lining structure is related not only to the amplitude of the input seismic wave, but also to the seismic wave duration.

() There is an accordance of distribution between the slip crack areas and the structure damage areas, suggesting that the lining structure has a higher possibility of being damaged in the areas which could slide relative to or separately from the surrounding rock, and the constraint of the surrounding rock contributes to the seismic stability of the lining structure. It is also confirmed that underground structure has better seismic stability than structure above the ground.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This study is supported by the National Key Basic Research Program of China (973 program) (Grant no. 2015CB057904) and the Major Program of the National Natural Science Foundation of China (Grant nos. 51579191 and 51279136). The supports are greatly acknowledged and appreciated.

References

H. Li, L. Zhu, T. Lu, J. Yang, X. Xia, and Y. Liu, “Seismic response analysis of an underground cavern groups in rock subjected to spatially non-uniform seismic ground motion,” Yanshilixue Yu Gongcheng Xuebao/Chinese Journal of Rock Mechanics and Engineering, vol. 27, no. 9, pp. 1757–1766, 2008.

Li. XJ and T. Lu, “Explicit finite element analysis of earthquake response for underground caverns of hydropower stations,” Journal of Hydroelectric Engineering, vol. 28, no. 5, pp. 1877–1882, 2009.

H. B. Li, X. D. Ma, and J. R. Li, “Study on influence factors of rock cavern displacement under earthquake,” Chinese Journal of Geotechnical Engineering, vol. 28, no. 3, pp. 358–362, 2006.

J. Y. Chen, Z. Q. Hu, and G. Lin, “Random seismic response analysis of large-scale underground caverns,” Journal of Hydraulic Engineering, vol. 2002, no. 1, pp. 71–75, 2002.

Y. Zhang, M. Xiao, and Y. Li, “Effect of Wenchuan earthquake on earthquake damage and dynamic response analysis of underground powerhouse of Yingxiuwan hydropower station,” Yanshilixue Yu Gongcheng Xuebao/Chinese Journal of Rock Mechanics and Engineering, vol. 29, no. S2, pp. 3663–3671, 2010.

Z. G. Zhang, Study on numerial simulation methods for seismic response of underground cavern complexs [Ph.D. thesis], Wuhan University, Wuhan , China, 2012.

J. Wu and J. Li, “Elastoplastic damage constitutive model for concrete considering strain rate effect under dynamic loading,” Tongji Daxue Xuebao/Journal of Tongji University, vol. 34, no. 11, pp. 1427–1440, 2006.

J. Y. Wu and J. Li, “Elastoplastic damage model for concrete considering damping effects,” Engineering Mechanics, vol. 23, no. 11, pp. 116–121, 2006.

A. Winnicki and C. Cichoń, “Plastic model for concrete in plane stress state. I: Theory,” Journal of Engineering Mechanics, vol. 124, no. 6, pp. 591–602, 1998.

J. W. Tedesco, J. C. Powell, C. A. Ross, and M. L. Hughes, “A strain-rate-dependent concrete material model for ADINA,” Computers & Structures, vol. 64, no. 5-6, pp. 1053–1067, 1997.

I. E. Shkolnik, “Influence of high strain rates on stress-strain relationship, strength and elastic modulus of concrete,” Cement and Concrete Composites, vol. 30, no. 10, pp. 1000–1012, 2008.

J. Li, X. D. Ren, and Q. P. Huang, “A viscoplastic damage constitutive model for concrete,” Chinese Journal of Theoretical and Applied Mechanics. Lixue Xuebao, vol. 43, no. 1, pp. 193–201, 2011.

S. Y. Xiao, H. N. Li, R. Q. Du, and G. Lin, “Four-parameter dynamic constitutive model of concrete,” Journal of Harbin Institute of Technology, vol. 38, no. 10, pp. 1754–1757, 2006.

Q. B. Li, C. H. Zhang, and G. L. Wang, “Dynamic damage constitutive model of concrete in uniaxial tension and compression,” Shui Li Xue Bao, vol. 12, pp. 55–60, 1994.

Q. B. Li, Z. C. Deng, and L. X. Zhang, “Dynamic damage constitutive model for concrete accounting for the difference of the initial elastic modulus,” Journal of Tsinghua University, vol. 43, no. 8, pp. 1088–1091, 2003.

G. Lin, J. Liu, and Z.-Q. Hu, “Current situation and progress of research on damage constitutive relation of concrete,” Dalian Ligong Daxue Xuebao/Journal of Dalian University of Technology, vol. 50, no. 6, pp. 1055–1064, 2010.

C. B. Du and Q. Z. Su, “Dynamic damage of concrete gravity dams under earthquake excitation,” Engineering Mechanics, vol. 20, no. 10, pp. 170–175, 2003.

J.-Y. Chen, J. Li, and G. Lin, “Seismic response analysis of concrete gravity dam based on strain rate-dependant damage model,” Harbin Gongye Daxue Xuebao/Journal of Harbin Institute of Technology, vol. 37, no. 6, pp. 786–789, 2005.

D. Perić 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, no. 6, pp. 1289–1309, 1992.

K.-J. Bathe and A. Chaudhary, “A solution method for planar and axisymmetric contact problems,” International Journal for Numerical Methods in Engineering, vol. 21, no. 1, pp. 65–88, 1985.

S. Tangaramvong and F. Tin-Loi, “An FE-MPEC approach for limit load evaluation in the presence of contact and displacement constraints,” International Journal of Solids and Structures, vol. 49, no. 13, pp. 1753–1763, 2012.

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.

B.-Y. Zhang and H.-Q. Chen, “A new LDDA iterative algorithm for dynamic contact forces,” Gongcheng Lixue/Engineering Mechanics, vol. 24, no. 6, pp. 1–6, 2007.

J. M. Raphael, “Tensile Strength of Concrete,” ACI Journal, vol. 81, no. 2, pp. 158–165, 1984.

G. Lin, J. Y. Chen, and S. Y. Xiao, “Dyanmic behavior of concrete and nonlinear seismic response of arch dam,” Shui Li Xue Bao, vol. 6, pp. 30–36, 2003.

S. Y. Xiao and G. Lin, “A Drucker-Prager consistent rate-dependent model,” Engineering Mechanics, vol. 20, no. 4, pp. 147–151, 2003.

X. R. Ge, Y. Jiang, and Y. D. Lu, “Testing study on fatigue deformation law of rock under cyclic loading,” Chinese Journal of Rock Mechanics and Engineering, vol. 22, no. 10, pp. 1581–1585, 2003.

J.-H. Yon, N. M. Hawkins, and A. S. Kobayashi, “Strain-rate sensitivity of concrete mechanical properties,” ACI Materials Journal, vol. 89, no. 2, pp. 146–153, 1992.

J. Eibl and B. Schmidt-Hurtienne, “Strain-rate-sensitive constitutive law for concrete,” Journal of Engineering Mechanics, vol. 125, no. 12, pp. 1411–1420, 1999.

D. Lu, G. Wang, X. Du, and Y. Wang, “A nonlinear dynamic uniaxial strength criterion that considers the ultimate dynamic strength of concrete,” International Journal of Impact Engineering, vol. 103, pp. 124–137, 2017.

X.-L. Du, Y. Wang, and D.-C. Lu, “Non-linear uniaxial dynamic strength criterion for concrete,” Shuili Xuebao/Journal of Hydraulic Engineering, vol. 41, no. 3, pp. 300–309, 2010.

J. Li and X. D. Ren, “A review on the constitutive model for static and dynamic damage of concrete,” Advances Mechanics, vol. 40, no. 3, pp. 284–297, 2010.

S. Y. Wang, L. Sun, C. Yang, S. Q. Yang, and C. A. Tang, “Numerical study on static and dynamic fracture evolution around rock cavities,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 5, no. 4, pp. 262–276, 2013.

X. Q. Zhou and H. Hao, “Modelling of compressive behaviour of concrete-like materials at high strain rate,” International Journal of Solids and Structures, vol. 45, no. 17, pp. 4648–4661, 2008.

J. J. Jiang, X. Z. Lu, and L. P. Ye, Finite Element Analysis of Concrete Structures, Tsinghua University Press, 2005.

L. Shu, L. Jingbo, S. Jiaguang, and C. Yujian, “A method for analyzing three-dimensional dynamic contact problems in visco-elastic media with kinetic and static friction,” Computers and Structures, vol. 81, no. 24-25, pp. 2383–2394, 2003.

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.

Y. Yang, Seismic response analysis of surrounding rock and structure for underground powerhouse of hydropower station [Ph.D. thesis], Wuhan University, Wuhan, China, 2012.