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 aswhere 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 bewhere 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 aswhere , 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 aswhere 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 havewhere 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 havewhere is static compressive strength of concrete material and is amplification coefficient of dynamic compressive strength. The proposed equation by Euro-International Committee for Concrete iswhere .
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 aswhere 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 aswhere .
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.
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 aswhere 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 aswhere 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 havewhere 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 aswhere 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.
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.
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
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.
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].
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/s2 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)
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.
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.