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

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


Introduction
The southwest region is the key area for hydropower development in China during the past few decades.A number of large-scale underground powerhouses of hydropower plants have been built in this region.It is also an 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 earthquake-resistant capability compared to their surface counterparts.Surrounding rock of the three underground powerhouses was generally stable, while the powerhouse concrete structure experienced obvious local cracking and failure.Concrete structure behaved as a weak link to affect the earthquake resistance capability of the underground powerhouses.Thus, it is very important to study the characteristics of seismic damage of underground powerhouse concrete structure.
It has been recognized that numerical method is an increasingly popular and effective strategy to solve the problem of structural dynamic response.A variety of works [1][2][3][4][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 2 Shock and Vibration powerhouse concrete structure are governed microscopic crack generation and propagation.In this way, the continuum method [6][7][8] (such as the finite element method, the finite difference method, and the boundary element method) and the discrete particle simulation method [9][10][11][12][13][14][15][16][17][18][19][20][21][22] (such as the discrete element method, the lattice-solid method, and the contact dynamics method) are the natural choice to simulate the damage and failure process.Among them, the finite element method is by far the most common and useful numerical method.The integration scheme it employs makes it appropriate for the widest variety of geologic and structural problems, and it can handle the most sophisticated constitutive relationships [23][24][25].
Underground powerhouse could be regarded as a system composed of powerhouse concrete structure and surrounding rock.Surrounding rock provides the external support for powerhouse concrete structure.As the concrete structure comes into direct contact with surrounding rock, seismic wave is propagated through rock mass to the powerhouse concrete structure.Therefore, the complex consisting of powerhouse concrete structure and surrounding rock undergoes a forced vibration.The simulation and analysis of the dynamic contact surface between surrounding rock and concrete structure is the key to the dynamic calculation of underground powerhouse concrete structure using the finite element method.At present, the numerical calculation methods for dynamic contact problems are primarily the Lagrange multiplier method [26], penalty method [27,28], and their improved versions [29,30].However, these methods tend to either increase the degree of freedom of the system or influence the time step of integration [31,32].These methods adversely affect the precision and speed of calculation when applied to the analysis of an underground powerhouse concrete structure, which involves a large number of contact elements and complex contact states.Liu et al. [33,34] put forward the dynamic contact force method targeted at the dynamic response problem of the contact crack.The convergence and stability of this algorithm were easy to meet, making it suitable for large and complex contact systems.But it cannot reflect the bond-slip properties of the contact surface.
In this paper, a new dynamic contact force method, considering the bond-slip properties of the contact surface, is suggested using the fundamental integration formulation of the dynamic contact force method.The algorithm of the method considers the cohesive effect of the contact surface between concrete structure and the surrounding rock.It is capable of simulating the large slip phenomenon of the contact surface under dynamic loads.Based on the proposed algorithm, a finite element model considering dynamic constitutive properties of materials is built for the dynamic numerical analysis of the underground powerhouse structure of Yingxiuwan Hydropower Plant.Then the characteristics of seismic damage of the underground powerhouse structure under dynamic loads are studied by means of the numerical analysis and postearthquake investigations.The results will provide a theoretical basis for the antiseismic design of underground powerhouse.

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

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 M, C, and K are the mass, damping, and stiffness matrix, respectively.U is the displacement vector.U is the velocity vector.Ü is the acceleration vector.F is the known vector of external forces.R 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.R  depends on the state of motion not only at time  but also at time  + Δ.Therefore, U +Δ and U+Δ 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.contact condition.As shown in Figure 2, the node () on the contact surface of concrete and the node (  ) on the contact surface of the surrounding rock are in point-to-point contact at time .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:

The Solution of Dynamic Contact
where n  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, R    = −R   , we have where N   , T   are the normal and tangential components of R   , and Δ 1 , Δ 2 , N   , and T   satisfy the following equations, respectively: In the above equations, T   and N   were obtained from the analysis of motion of the node pair.Therefore, they must satisfy the following inequalities: If the value of T   does not satisfy (12), then the node pair would enter into the state of sliding contact.We have If the value of T   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,  > 0, otherwise,  = 0.
Figure 3: Relationship of point-to-surface contact.

The Solution of Dynamic Contact Force under Pointto-Surface Contact Condition.
If larger relative sliding has occurred between contact nodes under dynamic loads, the contact nodes will be in the state of point-to-surface contact.
Then there are no cohesive forces between the contact nodes and surfaces.As shown in Figure 3, at time , one node  on the contact surface of concrete structure (or the surrounding rock) comes into contact with the surface of the surrounding rock (or concrete structure).Such a contact point on the contact surface is denoted as   .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 ‖Δ 3 ‖ < 0, node  is separated from the corresponding contact surface.N   , T   can be computed from ( 15).(b) If ‖Δ 3 ‖ ≥ 0, node  is in contact with the corresponding contact surface.It is assumed that (17) will result in ‖Δ 3 ‖ ≥ 0 for -number of nodes.Thus, ( 17) and ( 18) of the  nodes can be represented as follows: where [H] × is the coefficient matrix, relating Δ, , and .
Solving (20) and ( 21), the contact forces N  and T  can be obtained for the nodes in point-to-surface contact condition.If one node has ‖T  ‖ >   ‖N  ‖, then the node enters into the state of sliding contact.For this particular node, T  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 N  and T  of all the contact nodes from the above equations, the total dynamic contact force R  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.

Dynamic Constitutive Model of Concrete and Rock Mass
Under cyclic loading, the unloading stiffness of concrete and rock mass at the later yielding stage is lower than the stiffness at the initial linear stage.The plastic-damage model, proposed by Lubliner et al. [35], improved by Lee and Fenves [36,37], and generalized for 3D model by Omidi and Lotfi [38], could effectively simulate such a phenomenon.The model is suitable for the dynamic analysis of the quasibrittle materials, such as concrete and rock mass [3].According to the basic theory of plastic-damage model, the plastic-damage stress-strain relationship of rock mass or concrete can be expressed as follows: where  is the stress tensor. is the effective stress tensor.E 0 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. (0 ≤  ≤ 1) is the coefficient of restitution when the material shifts from tensile state to compressive state.

Shock and Vibration
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]. 1 and  2 are the first and second invariants of the effective stress tensor.σmax is the maximum effective principal stress.Concrete and the surrounding rock are used as friction material.The nonassociated flow rule can simulate the volume expansion properties under the compressive state.Therefore, the plastic-damage model employs nonassociated flow rule.The plastic potential function follows the Drucker-Prager hyperbolic function as follows: where  is the parameter which controls the potential function to approach the asymptote,   is the dilatancy parameter, and  0 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].

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

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  In order to guarantee the accuracy of dynamic calculation results, the maximum mesh size should be less than 7.5 m based on the cut-off frequency of seismic waves.As the maximum mesh size of the model is 7.0 m, the demand of FEM dynamic simulation can be satisfied.The profile of the model is shown in Figure 6.

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

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].

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

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.V   is the tensile or compressive damaged volume of each single element, and  is the number of damaged elements.

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  = 2.6 s.From  = 2.6 s to 20 s and from  = 30 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  = 60 s, the total damaged volume was about 460.3 m 3 .
Figures 9 and 10 show the plot of the damage coefficient when  = 0s and  = 60s, 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  = 0s and  = 60s 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.

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  = 0 s and  = 60 s was plotted as shown in Figures 13 and 14.When  = 0 s, a small amount of sliding or separation zones occurred on the lining.When  = 60 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.

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.

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.

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.

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 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.

Discussion on the Antiseismic Measures of Underground Powerhouse Structure
According to the seismic damage characteristics, the constraint imposed by the surrounding rock remained as an important factor to the antiseismic stability of the underground structure.The powerhouse structure above the generator layer suffered from the problem of "underconstraint, " resulting in a large area of damage in the upper structure under seismic loads.For the structures below the generator layer, the problem of "overconstraint" made the action of compression caused by the surrounding rock deformation more obvious.
To ensure the antiseismic design, surrounding rock is recommended to be reinforced to improve its stability.In this way, the favorable constraint of the surrounding rock on the internal structure can be increased, while the unfavorable impact is minimized.Secondly, anchorage support or other support measures are recommended to strengthen the constraint of surrounding rock for the structure above the generator layer, so as to promote the antiseismic capacity of the upper structure.Finally, soft-seismic isolation layer can be added between the lower parts of the concrete structure and the surrounding rock.The soft-seismic isolation layer can not only weaken the compression caused by the surrounding rock deformation, but also reduce the seismic response of the structure to a certain extent.The enhancement in the antiseismic capability of the structure after incorporating the abovementioned measures was evident from the numerical simulation.For instance, Figures 15 and 16 show the distribution of the damaged areas in the powerhouse structure after adopting the abovementioned antiseismic measures.Significant improvement is seen in the structural performance when compared to the damage in Figures 10 and 12.The distribution range and the damage coefficient were greatly reduced.The stability of the powerhouse concrete structure was enhanced considerably.

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

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

Figure 4 :
Figure 4: The flow chart for the dynamic contact force method.
(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)).

Figure 6 :
Figure 6: Profile of finite element model.

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

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

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

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

Figure 16 :
Figure 16: Damage coefficient distribution of the inner concrete structure when  = 60 s.
Force under Point-to-Point Contact Condition.If relative sliding does not occur between contact nodes, or the relative sliding is small, the contact nodes are in or approximately in point-to-point

Table 1 :
Mechanical parameter of rock mass and concrete.