Abstract
Layered rock mass generally has adverse effects on the seismic stability of deep-buried underground powerhouses. Aiming at the complexity of seismic wave field in deep-buried underground powerhouses, the obliquely input methods of P and SV waves considering the incident direction and multi-incident surfaces are constructed in this study, respectively. It can convert the seismic wave input problem into the problem of solving equivalent nodal force acting on the artificial boundaries. Based on the characteristics of dynamic interaction between interlayers in layered rock mass, an explicit dynamic contact force method considering the seismic deterioration effect and bond-slip characteristics of interface is also presented to simulate various contact states such as bond, separation, and sliding. The combined application of the above methods to analyze the seismic response of underground powerhouse at Azad Pattan hydropower station shows that they could accurately simulate the seismic damage evolution process of deep-buried underground powerhouses in layered rock mass. The numerical results indicate that the obliquely incident seismic waves contribute to a larger seismic reaction of underground powerhouses, which largely lies in the amplitudes of the displacement and stress fluctuations. After considering the dynamic contact, the obvious seismic deterioration effect and interlaminar dislocation displacement occur between soft and hard rock, and the sidewalls suffer more severe deterioration degree than the top arch. The seismic damage area of lining structure is mainly distributed in the place where the soft rock strata pass through and the upper structure with larger free surface. Additionally, two major damage modes of underground powerhouses in layered rock mass, namely, damage due to structural deformation and interlaminar dislocation, are derived from the numerical results and can be reasonably explained by the corresponding damage mechanisms.
1. Introduction
Due to the abundant water power resources in Northeast Pakistan, more than a dozen hydropower stations have been planned and constructed to alleviate the severe power shortage of Pakistan's State Grid, such as Karot, NJ, SK, and Azad Pattan hydropower stations, thus forming a large number of underground cavern groups. These hydropower stations are located in the southern foothills of the western Himalayas, adjacent to the alpine and gorge region in Southwest China. Strong earthquakes occur frequently and a variety of stratigraphic lithologies are widely distributed in this area, especially for the soft-hard alternant strata, resulting in the extremely complicated engineering geological conditions [1]. The survey performed after the great “5.12” Wenchuan earthquake revealed that the underground structure located in poor geological conditions was vulnerable to seismic damage including lining spalling, collapse, rupture, dislocation, invert uplift, and the reinforcing rebar distortion, especially at interface between soft and hard interbed layered rock mass [2]. In addition, the underground powerhouses located in the alpine gorge area of high earthquake intensity are generally buried deep. Consequently, it is of great significance to study the seismic response characteristics and failure mechanisms of deep-buried underground powerhouses in layered rock mass for the safe operation of the cavern.
The seismic damage of underground powerhouses in layered rock mass induced by earthquakes is a dynamic failure process, affected by the shaking of the seismic waves and the distortion of the interlayer staggered zone [3]. Therefore, the seismic response analysis of a deep-buried underground powerhouse in layered rock mass mainly includes two aspects: one is the near-field seismic motion input for deep-buried powerhouses, and the other is the dynamic contact analysis between interlayers. In the existing dynamic analysis of underground powerhouses, it is traditionally acknowledged that the seismic waves are vertical incidence from the bottom of the model, and the top of the model is built directly to the ground surface [4]. In order to reflect the influence of the ground reflection wave on the seismic response of underground structures, the above method is feasible for the dynamic calculation of shallow-buried underground powerhouses. But for deep-buried caverns under near-field earthquakes, when the epicenter is close to the engineering site, the number of finite elements has multiplied, inevitably leading to excessively large model scale and low computational efficiency. In addition, the seismic motion may propagate from any direction to the underground cavern, leading to the nonuniform spatial deformation of the underground structure [5]. Du et al. [6] studied the seismic response characteristics of tunnel structure under the obliquely incident seismic waves, and their achievement showed that the dynamic responses of tunnel structure differ from those under vertical incidence. Zhao et al. [7] analyzed the seismic response of large underground rock caverns under different incident directions of P waves, and the results indicated that the obliquely incident seismic waves dramatically influence the stability of underground rock caverns. Zhao et al. [8] and Huang et al. [9] investigated the effect of the obliquely incident seismic waves on the dynamic response of subway station structure and rock-fault-tunnel system. Wong et al. [10] developed a mathematical model to study the responses of pipes under plane seismic waves with different incident angles. Stamos and Beskos [11] proposed a special direct boundary element method to analyze the seismic responses of long tunnels under obliquely incident plane harmonic waves. Naggar et al. [12] studied the impacts of the incident angles of seismic waves on the internal forces in composite and jointed tunnel linings. Kouretzis et al. [13] presented a 3D thin shell strain analysis to study the earthquake-induced strains in long cylindrical underground structures subjected to seismic shear wave. In conclusion, some achievements have been made in the studies about the seismic response of shallow-buried underground structure under the obliquely incident seismic waves. However, restricted by the alpine canyon terrain, most of underground powerhouses in Northeast Pakistan are deeply buried. There are few studies on the influence of the near-field obliquely incident seismic waves on the seismic response characteristics of deep-buried underground powerhouses so far.
The dynamic interaction at the interface between soft and hard layered rock under seismic action possesses sophisticated contact nonlinear characteristics and is accompanied by the vibration deterioration effect. Consequently, it has always been a complex problem in the seismic response analysis of underground powerhouses in layered rock mass. Liu and Sharan [14] proposed a dynamic contact force method based on the explicit central difference method by introducing the contact conditions. This method has been widely applied to the dynamic calculation of the complex contact system because of its high computational efficiency and good convergence, compared with the Lagrange multiplier method [15], penalty method [16], and linear complementarity method [17]. However, the bond-slip properties of contact surface and the impact of the seismic cyclic loading on the vibration deterioration of interface are not considered in the method. In fact, the deterioration phenomenon of the mechanical properties of rock mass structural plane under dynamic loading has been confirmed by a large number of rock dynamic tests [18]. Lee et al. [19], Jafari et al. [20], and Li et al. [21] revealed the fact that the peak shear of rock joints is adversely influenced by the cyclic loading and the loading rate through a series of cyclic shear tests.
Aimed at the seismic wave field characteristics of deep-buried underground powerhouses, the input methods of the obliquely incident seismic waves are established based on 3D viscoelastic artificial boundary conditions. A dynamic contact force algorithm considering the seismic deterioration effect and the bond-slip characteristics of interface is also built based on the dynamic interaction characteristics of the contact surface. The methods constructed in this paper are applied to the seismic stability calculation of an underground powerhouse of the Azad Pattan hydropower station in Pakistan, in order to analyze the effects of layered rock mass on the seismic response of the underground powerhouse under the obliquely incident seismic waves. These results are expected to provide valuable references for the antiseismic design of underground powerhouses deeply buried in layered rock mass.
2. Input Method of Obliquely Incident Seismic Waves in Deep-Buried Underground Powerhouse
2.1. Basic Ideas
Compared with the ground and underground buildings of the city, the particularity of underground powerhouses lies in that they are completely buried in the infinite domain mountain medium. In view of the problem that the finite element model is inconvenient to be built to the ground surface in the dynamic time-history analysis of deep-buried underground powerhouses, an input method of 3D seismic waves that are obliquely incident to the model is set up based on the viscoelastic artificial boundary conditions. Using this method, the top direction of the model can be artificially truncated without being built to the surface, and the 3D viscoelastic artificial boundary is set at the top. Then, the seismic wave field reflected from the ground surface to the top artificial boundary is calculated by the analytical method, plotted in Figure 1, which is incident into the model through the viscoelastic artificial boundary conditions. According to the wave field decomposition theory, the total wave field around the model can be decomposed into internal and external wave fields. The external wave field is the wave field transmitted from the inside of the model to the infinite domain foundation through the artificial boundary, which can be solved by the explicit finite element step-by-step integral method. Consequently, the key to realizing the accurate input of obliquely incident seismic waves for deep-buried underground powerhouses is to obtain the obliquely incident wave field in the infinite domain foundation and the corresponding surface reflected wave field. Based on the 3D viscoelastic artificial boundary [22], the infinite domain seismic wave field is converted into equivalent nodal forces acting on the nodes of artificial boundaries, realizing the interaction between the internal and external wave field of the model.
2.2. The Solutions for Equivalent Nodal Force of Obliquely Incident Seismic Waves
The differences between obliquely and vertically incident seismic motion for deep-buried underground powerhouses include three aspects: the first is that the obliquely incident seismic waves will generate a complex wave-type conversion at the free surface of the half-space [23]. For example, the incident P waves could be reflected by the free surface to form reflected P and SV waves, respectively. As a result, the displacement and stress of incident wave field at the artificial boundaries become the superposition of those of the wave field generated by the incident P waves, reflected P, and SV waves, and so it is with obliquely incident SV waves as well, as shown in Figures 1 and 2. The second is that the seismic waves can be obliquely incident into the finite element model of deep-buried underground powerhouses through multiple artificial boundaries. As the distances between the nodes on the artificial boundaries and the epicenter differ from each other, there is a phase difference between different incident seismic waves received by the nodes at the same time. The last one is that the top arch and high sidewalls of an underground powerhouse structure are of larger free surfaces, where the complicated reflection and diffraction phenomena of obliquely incident seismic waves will occur, resulting in the dynamic amplification effect. Considering the above factors, it is assumed that the infinite domain outside the model is a homogeneous elastic medium with a horizontal ground free surface, and the seismic wave is the obliquely incident elastic plane wave. Therefore, the displacement field of the nodes on the artificial boundaries can be deduced as follows:
For obliquely incident P waves:
For obliquely incident SV waves:where , , and are the displacement time histories of incident P and SV waves, reflected P and SV waves, respectively.
There is an assumption that and are the displacement time histories of the incident P and SV waves at zero time, and α is the incident angle between the incident wave front at zero time and the free surface. Accordingly, the incident wave front at zero time is also assumed to be parallel to the longitudinal axis of the underground cavern. The size of the finite element model taken in the z direction is L, and the depth from the top of the model to the free surface is H. Point l (x0, y0, z0) is a node on the artificial boundary of the model. According to the geometrical relationship between the incident wave front and the artificial boundary node l, the displacement time history of the internal wave field at each artificial boundary of the model can be obtained as follows:
For obliquely incident P waves:
For obliquely incident SV waves:where β1 and β2 are the reflection angles of the incident P and SV waves at the free surface, and β1 = arcsin (cssinα/cp), β2 = arcsin (cpsinα/cs); cp and cs are the wave speeds of plane P and SV waves; A1 and A2 are the amplitude amplification ratios of reflected P and SV waves to incident P waves; A3 and A4 are the amplitude amplification ratios of reflected SV and P waves to incident SV waves, and their values can be determined by [23]; Δt denotes the time interval of the incident seismic waves propagating from the wave front at zero time to the boundary node l, so Δt can be described as
Given the displacement time history of the internal wave field, the speed field can be solved by derivation or difference. The corresponding stress field can be confirmed by the internal wave field at artificial boundaries based on the generalized Hook’s law.
According to Equations (1)–(5), the equivalent load fli at direction i of the boundary node l can be obtained by substituting the displacement , the velocity , and the corresponding stress into (6), to realize the input method of obliquely incident seismic waves in deep-buried underground powerhouses:where Al is the equivalent area for the boundary node l; Kli and Cli are the spring and damping coefficients of node l, and the subscript i represents the direction of the component, including the normal and tangential directions. Kli and Cli can be determined by reference [22].
2.3. Numerical Verification of the Input Method
A numerical model is established to analyze the dynamic response of the obliquely incident seismic waves in semi-infinite space to verify the input method. The size of the 3D finite element model is 800 m × 800 m × 800 m. A total of 64000 hexahedral elements are meshed with the maximum mesh size of 20 m, which meets the mesh size requirement, as shown in Figure 3. The elastic modulus of the medium is 10 GPa. Poisson’s ratio is 0.3. The medium mass density is 2000 kg/m3. The 3D viscoelastic artificial boundary is applied for the bottom and sides of the model, while the top is the free surface. The central point A (400, 400, 800) of the free surface is selected as the observation point. Figure 4 is the displacement time history of the input P waves, whose peak displacement is 1.0 m.
Figure 5 exhibits the vertical and horizontal direction displacement time histories at point A under P waves with incident angles of 15° and 30°, respectively. It can be observed clearly that the numerical results of the vertical and horizontal displacements at point A are in good agreement with the theoretical results. The above comparison demonstrates that the developed input method can simulate the free field motion of the obliquely incident seismic waves in semi-infinite space.
(a)
(b)
3. Vibration Deterioration Law of Interface in Layered Rock Mass under Seismic Loading
Because the seismic load is both cyclic load and dynamic load [19], the contact surface between interlayers in layered rock mass under seismic action is simultaneously affected by the cyclic shearing action of cyclic load and the deformation rate of dynamic load, namely, the influence factors of vibration wear and relative shear velocity. The vibration deterioration effect of the seismic action on interlaminar interface is a complex dynamic change process. To quantitatively describe the deterioration degree of interface, the vibration deterioration coefficient D (t) change with time is introduced. Then, the shear strength of the interlaminar interface at time t can be calculated bywhere and are the peak shear strength of interface at the initial time and time t; σ(t) is the normal stress of interface at time t; φ0 and c0 represent the initial internal friction angle and initial cohesion of interface.
3.1. Mathematical Expression of Vibration Wear Influence Coefficient
During the seismic duration, the cyclic shearing action on the interlaminar interface could cause the fatigue wear and passivation of interface, leading to the shear strength of interface decreasing accordingly. It is assumed that the deterioration degree of interface caused by the fatigue wear and passivation decays with a negative exponential function over time under earthquakes [24]. Then, the vibration wear influence coefficient η(t) of interface due to the cyclic shearing action of seismic loads can be written aswhere a and b are dimensionless parameters, δ0 is the convergence value, and all of those can be fitted by experiment; J(t) is the cyclic shear amplitude; K(t) is the cyclic shear number.
In engineering practice, the dynamic shear process of the interlaminar interface in layered rocks is relatively complex and strongly nonlinear. Therefore, it is difficult to obtain the accurate values of the cyclic shear amplitude J(t) and the cyclic shear number K(t), which makes (8) greatly limited in engineering application. In consequence, based on the comprehensive consideration of J(t) and K(t), the dynamic shear process of interface is believed to be a process where the fatigue damage of interface accumulates continuously and the shear strength decreases gradually, which can be assumed to be in connection with the current time and total duration of seismic loading [25]. Thus, the vibration wear influence coefficient η(t) can be obtained aswhere t and ts are the current time and the total duration of seismic loading, respectively; λ is the dimensionless parameter, and R0 is the convergence value of the weakening degree of interface after an earthquake, all of which are fitted by experiment.
3.2. Mathematical Expression of Relative Shear Velocity Influence Coefficient
Li et al. [21] found that the relationship between the peak shear strength of rock joints and the relative shear velocity presents an obvious negative exponential decay form. Then, the influence coefficient γ(t) of the vibration deterioration effect of interface caused by the relative shear velocity can be expressed aswhere is the relative shear velocity between the contact node pairs at time t; m, n and p are dimensionless parameters, all of which are fitted by experiment.
3.3. Vibration Deterioration Coefficient of Interface
Since the mechanism of these two factors on the strength deterioration of interface is different, it is also assumed that the deterioration effect of the cyclic shearing action and relative shear velocity on the shear strength of interface is independent of each other during the earthquake process. So, the vibration deterioration coefficient D (t) of the interlaminar interface at time t is obtained aswhere the parameters R0, λ, m, n, p should be confirmed by the shear test of rock mass material in engineering practice. In this paper, from the perspective of engineering safety, their values are advised as R0 = 0.75, λ = 5, m = 0.883, n = 0.02, p = 0.032, according to the experimental results in references [21, 24].
The relative shear velocity between the contact node pairs in (11) can be solved by the explicit central difference method. At each time step of the dynamic calculation, the vibration deterioration coefficient D(t) of interface is also solved explicitly, and the shear strength parameters get updated in time. Therefore, the dynamic degradation process of the interlaminar interface during the seismic duration can be obtained and applied to the explicit dynamic contact analysis method in Section 4.
4. Dynamic Analysis Model for Layered Rock Mass Contact System
Under seismic action, the soft and hard layered rock mass transfer loads through the contact surface, and the interaction between them belongs to an intricate dynamic contact problem. The evident dislocation displacement occurs easily between interlayers of layered rocks under the combined action of tectonic stress and seismic load, leading to the severe damage of underground powerhouse structure. Aimed at the complex dynamic contact features of the interface between soft and hard rocks under seismic cyclic loading, a dynamic contact force analysis model considering the seismic deterioration effect and bond-slip characteristics of the interface is put forward in this section.
After the finite element discretization of the interlaminar contact system of layered rock mass, the kinematic equation of nodes including dynamic contact force under seismic action has the form ofwhere M, C, and K represent the mass, damping, and stiffness matrix of the contact nodes, respectively; , and are the acceleration, velocity, and displacement vectors, separately; F is the external load vector; R is the dynamic contact force vector, and R = N + T, N and T represent the normal and tangential components of R.
For any contact node pair i and i′, as plotted in Figure 6, given the displacement , velocity and dynamic contact force at time t, the explicit integral equation of contact nodes at time t + Δt can be obtained by central difference method:where is the displacement vector of the contact nodes without considering the dynamic contact force at time t + Δt; is the additional displacement vector caused by the dynamic contact force Rt at time t + Δt.
According to Equations (12)–(15), it is evident that the motion state of contact nodes at time t + Δt can be determined by the motion state and dynamic contact force of contact nodes at time t. The displacement and velocity of contact nodes are known at time t, whereas the dynamic contact force Rt is unknown and depends on the movement state at time t as well as t + Δt. Therefore, the dynamic contact force Rt can be solved based on the contact conditions from time t to t + Δt.
4.1. Dynamic Contact Force Algorithm considering Complex Shear Strength of the Interface
Supposing that the interface between interlayers of layered rock mass is well cemented before an earthquake, the contact node pairs keeps in a bond contact state considering the interface cohesion, as shown in Figure 6. The relative slip may easily occur between interlayers in case of strong seismic events, resulting in the contact node coming into contact with the surface of the corresponding element. At this time, the interface enters into the sliding contact state without considering the interface cohesion.
It is assumed that the contact node pairs between interlayers are still in bond contact state at time t + Δt. Therefore, they should satisfy the contact boundary conditions, including the deformation coordination conditions and the contact force boundary conditions:where ni is the unit normal vector of the contact node pair, from node i′ to node i; ti is the unit tangent vector.
Substituting Equations (13)–(15) into (16),where Mi and are the lumped masses of node i and i′, respectively.
Obviously, the dynamic contact forces and in (17) are solved based on the assumption that the contact node pairs are in bond contact state. Actually, when the dynamic contact force exceeds the shear strength of interface, there are still various contact states such as sliding and separation between soft and hard rocks during strong earthquakes. Hence, it is necessary to check and correct the dynamic contact force according to the ultimate bearing strength of interface. The normal ultimate tensile and tangential ultimate shear force of interface are defined aswhere Ai is the control area by node i; c and μs are the cohesion and static friction coefficient of interface, respectively. D(t) is the vibration deterioration coefficient of interface at time t obtained in Section 3. Before time t + Δt, if the contact node pair is in a bond contact state, the interfacial cohesion should be considered and c > 0. Otherwise, if the contact node pair produces relative slip under seismic action and the interface enters into a sliding contact state, the interfacial cohesion should be no longer considered and c = 0.
The failure modes of interface mainly include the shear slip along the tangential direction and tensile fracture along the normal direction.(1)If and , it indicates that the interface is in tension state and tensile fracture occurs on the interface, so the dynamic contact force should be corrected as(2)If and , or and , it demonstrates that the contact node pairs are in the bond contact state, and the dynamic contact force should not be amended.(3)If and , it is indicated that the contact node pairs are in the sliding contact state, and there is shear slip on the interface. Namely, the interlayers of layered rocks develop relative movement under the shear force action. So, the tangential contact force should be revised: where μd is the dynamic friction coefficient of interface. The dynamic contact force at time t can be solved by Equations (17)–(20) and substituted into Equations (13)–(15), so the motion state of the contact nodes at time t + Δt can be updated in result.
5. Engineering Application and Analysis
5.1. Project Profile and Calculation Model
As illustrated in Figure 7(a), the Azad Pattan hydropower station is situated on the Jhelum river in Pakistan with priority given to power generation. The underground powerhouse is located in the mountain of the left bank. The bedrock where the underground powerhouse is buried shows interlayered distribution between sandstone and mudstone with complex geological conditions [1]. The strata are mainly monoclinal structures with steep dip, whose inclination angles range from 61° to 75°. The interlaminar shear and compressive rupture zone is reasonably common in the project area and has a great influence on the stability of high sidewalls of underground powerhouses.
(a)
(b)
(c)
The Azad Pattan engineering area is a tectonically active area, which is mainly affected by the subduction movement of the Indian plate towards the Eurasian plate. Northern Pakistan and Azad Jammu and Kashmir are areas of intensive earthquake activities, greatly influenced by multiple seismic plate tectonics, as plotted in Figure 7(b). The region is characterized by strong present seismic activity, in which the latest major earthquake was a magnitude 7.6 earthquake on October 8, 2005, with the epicenter near the upper reaches of the Azad project area. In the light of the research of the China Earthquake Administration, the peak acceleration of the 50-year exceedance probability 10% of the project area is 0.32 g, and the corresponding seismic basic intensity is eight degrees. Generally, the seismic and geological background is comparatively complicated. The underground powerhouse adopts a single cavern scheme and is buried in 300–350 m below the ground surface. The size of the main powerhouse cavern is 157.0 m × 24.9 m × 60.15 m (length × width × height), which has the characteristics of a large span and high sidewalls. The thickness of lining structure in the main powerhouse is set as 1.25 m for top arch and 1.0 m for sidewalls. The lining structure is made of C25 concrete. The lower structure of the main powerhouse mainly consisted of mass concrete structure, such as pier base.
A three-dimensional finite element model of the underground powerhouse is set up by selecting the cavern area with interbedded distribution of mudstone and sandstone as shown in Figures 7(c) and 8. The strike is nearly perpendicular to the cavern axis, of which the dip angle is 65°. The computational model is discretized by hexahedron element with 8 nodes, including hard rock, soft rock, lining, and internal mass concrete structure. There are 140524 elements and 156828 nodes divided totally, 23668 of which are concrete structural elements. The model ranges along x-, y-, and z-axes are 190.0 m, 200.0 m, and 170.0 m, respectively. The x-axis of the model is perpendicular to the cavern axis along the horizontal direction, and the y-axis is parallel to the longitudinal axis of the cavern. The z-axis coincides with the geodetic coordinate.
The in situ stress field is confirmed by inversion analysis of the measured data. The lateral pressure coefficients are taken as kx = 1.1, ky = 0.85 and kz = 1.0. The mechanical parameters of the rock mass, interface, and lining material are shown in Table 1. Before the dynamic response analysis of deep-buried underground powerhouses in layered rock mass, the static excavation and support simulation of the caverns are firstly implemented by 3D elastoplastic damage finite element method [26]. Then, the corresponding static calculation results are taken as the initial state of dynamic analysis.
5.2. Calculation Conditions
The calculation program adopts our research group’s in-house dynamic finite element program by a self-developed FORTRAN program [27], in which the input method of obliquely incident seismic waves and the dynamic contact force method are implemented. Before inputting seismic waves, the contact surface is firstly set at the interlaminar interface between hard and soft rocks. The node separation of hard and soft rock elements is accomplished by adding common contact node pairs on both sides of the interface. In numerical calculation, the dynamic elastic-plastic damage constitutive model based on the Mohr-Coulomb yield criterion is used for surrounding rock [28], while the plastic damage constitutive model is used for lining and concrete structure [29]. The scalar damage variable d for lining is defined aswhere dt and dc are the tension and compression damage variables, respectively. s is the stiffness recovery coefficient of the element from tension state to compression state during load reversal. At, Bt and Ac, Bc are the corresponding tension and shear damage parameters, and their values can be obtained with uniaxial compression tests for Ac, Bc and flexion tests for At, Bt. is the equivalent plastic strain. and are the tension and shear damage thresholds. It should be noted that the tensile criteria should be given higher priority in numerical calculation. When an element meets the maximum tensile strain criterion, the lining element is considered as tension damage. And it is unnecessary to judge whether the element satisfies the Mohr–Coulomb yield criterion. Otherwise, the Mohr–Coulomb yield criterion is used to identify the compressive strain state of the element.
Given the large buried depth of the underground cavern, a total of 580,000 elements needs to be added if the model is built to the ground surface. Thus, in order to improve the computational efficiency, the top of the model is built to double the height of the cavern. Accordingly, the 3D viscoelastic artificial boundary is imposed on the bottom, four sides, and top of the model, to absorb the obliquely incident seismic waves and the corresponding reflected waves. The El-Centro wave is employed as the input seismic wave. According to the seismic fortification intensity of the Azad Pattan Hydropower Station, the peak acceleration is adjusted to 3.15 m/s2 and the first 20-second acceleration time history with strong seismic motion is intercepted, which are processed by filtering and baseline correction. The influence of obliquely and vertically incident seismic waves on surrounding rock and lining is considered in this study, where the input angle is 30° under oblique incidence. Besides, the effects of P waves and SV waves on the underground powerhouse are also taken into account for dynamic calculation. Meanwhile, SV waves adopt the incident waves exhibited in Figure 9, and the acceleration of P waves is 2/3 that of SV waves. Since the explicit central difference method is employed to solve the near-field wave propagation problem, it is necessary to determine an appropriate time step Δt to meet the accuracy and stability of the solution. Based on the wave motion principle, Δt can be understood as the shortest time during which the propagation distance of the seismic wave is not greater than any element size. It can be written aswhere α is an empirical coefficient and is generally 0.8 ≤ α ≤ 0.98; Ce is the wave velocity. For the hexahedron element with 8 nodes, le is the ratio of the volume to the maximum area of the element.
(a)
(b)
In order to study the seismic response time-history characteristics of the underground powerhouse structure, the middle section of soft rock is selected as the monitoring section. Four points on the section are picked out as monitoring points to monitor the displacement and stress response of lining structure. The other monitoring points are set in the surrounding rock on both sides of the contact surface between soft and hard rock to monitor the relative displacement and velocity characteristics. The specific arrangement of the monitoring scheme is shown in Figure 10. The seismic calculation is divided into three different cases: ① the vertically incident seismic waves for the model without dynamic contact force; ② the obliquely incident seismic waves for the model without dynamic contact force; ③ the obliquely incident seismic waves for the model with dynamic contact force. The obliquely incident direction vector of seismic waves is (0.5, 0, 0.866), and the incident reference point is the origin of coordinates.
(a)
(b)
5.3. Calculation Results and Analysis
5.3.1. Seismic Deterioration Analysis of the Interlaminar Interface in Layered Rock Mass
Under strong earthquake loading, cyclic and reciprocating interactions occur easily between interlayers of layered rock mass, and the obliquely incident seismic waves aggravate the complex dynamic behavior, leading to vibration deterioration effect on the contact surface. The relative shear velocity time-history curves of contact node pairs near the interface between soft and hard rock at the sidewall and top arch in the case ③ are shown in Figure 11. Thus, the vibration deterioration coefficients of the shear strength are calculated based on the above (11), and the corresponding time-history curves of the seismic deterioration coefficient of interface are exhibited in Figure 12.
It can be found from Figure 11 that the relative shear velocity of the interlaminar interface presents a fluctuating change related to the input seismic waveform over time under cyclic loading and unloading of surrounding rock. In the first 5 s when the input seismic motions are relatively severe, the relative shear velocity of contact node pairs at the sidewall and top arch of the main powerhouse fluctuates violently, ranging from −0.83 to 0.99 mm/s. Moreover, the corresponding amplitudes of the relative shear velocity are larger for the period of more severe vibration. According to Figure 12, the seismic deterioration coefficients of interface between soft and hard rock at the sidewall and top arch generally show an obvious law of negative exponential attenuation over time. In 0 to 5 s, the seismic deterioration coefficient D(t) attenuates rapidly from 1.0 to 0.8, with the evident vibration deterioration effect. In 5 to 10 s, the coefficient D(t) decreases slowly over time and gradually deteriorates to the convergence value 0.75. After 10 s, the coefficient no longer continues to decay and basically tends to be stable, only fluctuating up and down at the convergence value under the influence of the relative shear velocity. Overall, the seismic deterioration effect of the high sidewall in the main powerhouse is much more severe than that of the top arch, which is mainly reflected in the amplitudes of the coefficient fluctuation. Especially in the period of 10∼20 s, the fluctuation amplitude of the seismic deterioration coefficient of the high sidewall is much larger than that of the top arch under the effect of the relative velocity. The former mainly fluctuates in the range of 0.71∼0.78, while the latter fluctuates slightly around 0.75. Consequently, the time-history curves of the seismic deterioration coefficient can intuitively reflect the dynamic deterioration process of the interface between soft and hard rock during a strong earthquake.
5.3.2. Displacement Time-History Analysis of Underground Lining
Figure 13 provides displacement time-history curves of 4 monitoring points for lining structure of the main powerhouse under three various cases. The same seismic response laws under three cases can be summarized: (1) the displacement time-history curves of the top arch, left and right sidewalls, and the bottom of the main powerhouse almost resemble each other in the forms and patterns, indicating that various parts of underground powerhouses keep in a synchronous vibration state. (2) In the first 5 s, the displacement time-history curves of 4 monitoring points fluctuate sharply, and the peak displacements of the sidewalls are significantly larger than those of the top arch and bottom. The motion state of the underground powerhouses is mainly manifested as the spatial fluctuation with seismic waves during earthquake, so the overall deformation may occur on the underground structure. However, the large overall deformation does not mean that such large damage occurs to the underground powerhouse. Therefore, it is necessary to introduce the relative displacement further to characterize the structural deformation. In this paper, the peak displacement difference between monitoring points (A, B, C) and the point D at the bottom of the cavern is used to represent the peak relative displacement, as plotted in Figure 14.
(a)
(b)
(c)
The maximum displacement of 4 monitoring points reaches 6.74 cm, and the peak relative displacements between left and right sidewalls, the top arch and the bottom are 3.1 cm, 2.7 cm and 2.1 cm, which indicates that the relative deformation of the main powerhouse is relatively small in the case ①. When considering the spatial oblique incident characteristics of seismic waves, the impact of the seismic motion on the displacement response of the lining structure under case ② mainly lies in the amplitude of the displacement fluctuation. The maximum displacement of 4 monitoring points is 7.90 cm, and the peak relative displacements between monitoring points (B, C, A) and the point D are 3.7 cm, 3.2 cm, and 2.4 cm, respectively. By contrast, the displacement response of the lining structure is greatly affected by the incident angle of seismic waves, making it extremely easy to cause the inconsistent deformation of underground structure. It is mainly due to the fact that the input methods of obliquely incident seismic waves consider the incident direction and the inconsistency of seismic waves, as well as the influence of the reflected seismic waves on the wave field of deep-buried underground powerhouses. The incident P waves and SV waves will generate a complex wave transformation at the free surface, forming reflected P and SV waves, respectively. Therefore, the seismic wave field at the artificial boundaries becomes the superposition of different incident and reflected waves, and nodes on the boundaries also have different vibration waveforms and vibration directions. The above wave-field characteristics do not exist in the case of vertical incidence; hence, the displacement response of the lining structure under obliquely incident seismic waves is larger than that under vertical incidence.
After considering the obliquely incident seismic waves with dynamic contact force, the maximum displacement of monitoring points is 9.56 cm, and the peak relative displacements of the three monitoring points (B, C, A) are 4.9 cm, 4.2 cm, and 3.0 cm in the case ③. It can be concluded that when considering the vibration deterioration effect of interface between soft and hard rocks, the shear strength parameters of interface tend to deteriorate gradually under earthquakes, and the interface cannot provide the sufficient antisliding force. Furthermore, the obliquely incident seismic waves exacerbate the seismic deterioration effect of interface and shear slip failure occurs between interlayers of layered rocks, leading to a large relative deformation of the high sidewall in the main powerhouse.
5.3.3. Relative Displacement Analysis of the Interlaminar Interface in Layered Rock Mass
If the maximum negative displacement is defined as the peak displacement, the x-direction peak displacement distribution of the high sidewall along the cavern axis is depicted in Figure 15. The x-peak displacement response of surrounding rock in soft rocks is obviously larger than that of other parts. It indicates that earthquake has a significant amplification effect on the rock deformation of soft and weak rocks in layered rocks. When considering the dynamic contact, the discontinuous deformation occurs near the interfaces and both sides of soft rocks in layered rocks, contributing to large relative shear deformation on the lining structure.
In order to further illustrate the influence of the obliquely incident seismic motion on dynamic contact response of the interlaminar interface in layered rocks, Figure 16 offers relative displacement time-history curves between soft and hard rock under different cases. The interlaminar relative displacement between soft and hard rock fluctuates around the 0 line in the case ①, with a fluctuation range from −0.75 to 0.80 cm. Due to the inconsistency of the obliquely incident seismic input, the interlaminar relative displacement in the case ② experiences a sharper fluctuation than that in the case ①, ranging from −1.20 to 1.30 cm. The maximum relative displacement is 2.15 cm at around 6.0 s, and the value gradually decreases to around 0 after earthquake. Considering the seismic deterioration effect and dynamic contact force in the case ③, the evident dislocation displacement occurs between soft and hard rock in 0∼7 s, with the maximum dislocation displacement of −5.01 cm. With seismic loading increasing, the dislocation displacement basically varies in the range of −4.0∼−5.0 cm in 7∼20 s and finally reaches −4.93 cm after the earthquake. It shows that the vibrations of soft and hard rock strata in layered rock mass are not synchronized under strong earthquakes, and the apparent shear slip failure occurs between the strata.
5.3.4. Stress Time-History Analysis of Underground Lining
Since the compressive strength of concrete is far more than its tensile strength, the tensile failure is the main damage and destruction mode of lining structure for underground powerhouse under seismic cyclic loading [24]. Therefore, this paper mainly analyzes the variation laws of maximum principal stress for lining structure under earthquake motion, as displayed in Figure 17. The maximum principal stress time-history curves of the top arch, left and right sidewalls, and the bottom of the main powerhouse under three different cases are basically consistent, which are similar to the input seismic waveform. In 0 to 5 s, the maximum principal stress changes sharply and increases rapidly. Despite several crests, the peak value and fluctuation amplitude reduce remarkably after 5 s, which reveals that the fluctuation law of the maximum principal stress is mainly affected by the input seismic wave.
(a)
(b)
(c)
The magnitude of the maximum principal stress, which can be believed to be the damage indicator of lining structure, attracts the most attention because the concrete material is well known to be of weak tensile capacity. As demonstrated in Figure 17, the maximum tensile stresses of the top arch (point A), sidewalls (points B and C), and the bottom (point D) of lining structure are 1.19 MPa, 1.32 MPa, and 0.76 MPa in the case ①, respectively. Under obliquely incident seismic waves, the maximum tensile stresses of the three parts are separately for 1.31 MPa, 1.54 MPa, and 0.97 MPa in the case ②. When considering the dynamic contact, the maximum tensile stress of the three parts is 1.50 MPa, 1.89 MPa, and 1.27 MPa in the case ③. In contrast with the three cases, the maximum tensile stresses at the top arch and sidewalls are much greater than those of the bottom, and their values even clearly exceed the concrete tensile strength. This result indicates that the top arch and sidewalls of the main powerhouse are damaged by tension cracking, which are the weak parts of lining structure under seismic action. Besides, the maximum tensile stress values of sidewalls and the top arch under case ② increase by 16.7% and 10.1% in comparison with the case ①. It is concluded that the nonuniformity of obliquely incident seismic waves leads to a larger stress response than that under vertical incidence due to the superposition effect of the incident and reflected wave field. When the seismic deterioration effect and the bond-slip characteristics of interface are taken into account, the magnitudes of sidewalls and the top arch in the case ③ are approximately 22.7% and 14.5% greater than those in the case ②. Obviously, it can be indicated that the vibration deterioration and the relative slip between interlayers of layered rock mass have a great significance on the stress response for lining structure, which aggravate the tensile stress fluctuation, resulting in the serious tension damage occurring at the top arch and sidewalls.
5.3.5. Seismic Damage Characteristics of Underground Structure
The damage coefficient d of the lining and concrete structure, which can be acquired by (21), is employed to quantitatively demonstrate how much damage is induced by the earthquake and where it is seriously damaged. The turbine unit where the soft rock passes through is selected for seismic damage analysis, of which the damage coefficient distributions under different cases are depicted in Figures 18 and 19. It can be seen from Figure 18 that the intensity, region, and extent of the damage gradually expand from case ① to case ③. The seismic damage area is mainly distributed in the top arch, high sidewalls, the generator floor, and turbine floor in the case ① of vertical incidence. The damage coefficient in most areas is relatively small and varies from 0.1 to 0.3, with the maximum damage coefficient less than 0.4. After considering oblique incidence, the distribution and magnitude of the damage increase on a large scale, the main display of which is that the damage zone extends further along the longitudinal and vertical direction of the underground cavern separately. More seriously, the damage area at the top arch tends to be connected with that at the upper structure of high sidewalls. Under the conditions of case ③, the lining and concrete structure suffers the most severe damage, and most areas of the main powerhouse are nearly covered by the damage zone. There are rather serious damage and failure occurring at local areas of lining structure, such as the junctions of the top arch and sidewalls, and the generator windshield, with the maximum damage coefficient almost close to 0.9. As a result, the seismic damage degree of underground caverns is not only affected by the properties of seismic waves, such as the incident angle and amplitude, but also related to its spatial distribution characteristics.
(a)
(b)
(c)
The tension and shear damage are the major damage types of the lining and concrete structure, mainly caused by the stresses exceeding the strength limit of concrete under a strong earthquake. The damage type distributions of the underground structure under different cases are exhibited in Figures 20 and 21. Looking horizontally in Figure 20, the upper structure of the main powerhouse is more easily prone to damage and failure than the lower. The tension damage is mainly distributed in the top of the roof arch, upper sidewalls, and floors of the main powerhouse, while the shear damage mostly occurs in the roof arch and the entire sidewall in all three cases. After comparing Figures 18 and 20, it can be found that the damage coefficient of the tension damage zone is larger than that of the shear damage zone, and the damage situation is much more serious as well. Looking longitudinally on Figure 21, the tension damage largely occurred in the center part of soft rock strata in layered rock mass, whereas the shear damage is mainly distributed on the interlaminar dislocation interface and presents a tendency of expanding distribution from the interface to both sides, forming the interlaminar shear and compressive rupture zone.
(a)
(b)
(c)
5.3.6. Damage Mechanism of Underground Powerhouse in Layered Rock Mass
According to the numerical simulation results, two major damage and failure modes of the underground powerhouse in layered rock mass under strong earthquake can be summarized, which are caused by structural deformation and interlaminar dislocation. Figures 22 and 23 provide comparisons between the sketches of damage mechanisms, depicted in Figures 22(a) and 23(a), and the predicted displacement results, plotted in Figures 22(b) and 23(b). As seen in Figure 22(a), the underground powerhouse can be mainly divided into the upper structure composed of the top arch, upper sidewalls, and floors and the lower structure made up of the mass concrete such as the turbine pier. The upper structure is relatively thin and weakly constrained, while the lower structure is poured into the whole concrete with a strong deformation constraint. The obliquely incident seismic waves will produce the amplified horizontal forces acting on the structure in the propagation process, causing the large relative deformation between the upper and lower structures. Figure 22(b) shows that the x-direction displacement of the underground powerhouse tends to decrease gradually from the top arch to the bottom, and the upper structure suffers the largest displacement deformation. The horizontal relative displacement of the upper and lower structures ranges from 3.1 to 4.2 cm, and the underground cavern suffers a large structural deformation, which eventually leads to serious damage of the top arch and high sidewalls.
The comprehensive response of the underground powerhouse under interlayer dislocation is shown in Figure 23(a). Under the seismic action, the displacement deformation of lining structure where the soft rock strata pass through is significantly different from that on both sides of the interlayer dislocation surface (see Figure 23(b)). Therefore, the spatial distribution characteristics are likely to contribute to shear failure of the lining structure.
6. Conclusions
Based on the input method of the obliquely incident seismic waves and the dynamic contact force method considering the vibration deterioration effect and bond-slip characteristics of interface, a dynamic response analysis method for deep-buried underground powerhouses in layered rock mass is constructed. The method is employed to study the seismic response and damage mechanism of the underground powerhouse in Azad Pattan Hydropower Station. The following conclusions can be drawn from the research.(1)The displacement and stress time-history curves of different parts of the underground powerhouse resemble each other under seismic loading. The displacement and stress responses of lining structure are greatly affected by incident angles of seismic waves, which mainly lie in the fluctuation amplitudes. The inconsistency of seismic motion input intensifies the dynamic response of the underground powerhouse. After considering the vibration deterioration effect, the displacement and stress of lining structure around interlayers develop further, and the maximum tensile stress of the top arch and upper sidewalls exceeds the concrete tensile strength.(2)The vibrations of soft and hard rock strata in layered rock mass are not synchronized under strong earthquake, leading to the obvious dislocation displacement between interlayers, with a maximum up to 4.93 cm. It is indicated that the evident shear slip failure and seismic deterioration effect occurring between interlayers and the upper sidewalls suffer a more serious deterioration degree than the top arch. The seismic deterioration coefficient of interface presents a law of negative exponential decay over time, whose time-history curve can directly reflect the dynamic deterioration process of interface at different times.(3)Looking horizontally, the damage area is mainly distributed in the places that have a large free surface, such as the top arch, high sidewalls, and floors. Looking longitudinally, the tension damage is mostly distributed in the place where the soft rock passes through, while the shear damage is largely distributed on the interlaminar dislocation interface and extends from the interface to both sides.(4)Two main damage and failure modes of the underground powerhouse in layered rock mass under seismic action are analyzed in detail. The amplified horizontal force generated by seismic loading on the underground cavern leads to structural deformation and contributes to the distortion and cracks on the top arch and sidewalls. The interlayer dislocation results in the shear deformation of lining structure between soft and hard rock strata.
Data Availability
The proposed numerical model and data used during the current study are available from the corresponding author.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by the Major Program of the National Natural Science Foundation of China (Grant nos. 52079097 and 51579191), the National Key Basic Research Program of China (973 program) (Grant no. 2015CB057904), and China water resources Beifang investigation, design and research co. ltd (BIDR). The supports are greatly acknowledged and appreciated.