Numerical SimulationAnalysis ofCombined SeismicResponse for Rock-Lining-Water in Hydraulic Tunnel

In order to explore the influence of internal water on the seismic response of hydraulic tunnel, the combined mechanical analysis models of multimaterial including surrounding rock, lining structure, and internal water are built. Based on the explicit central difference method, the dynamic finite element analysis methods for rock, lining, and water are discussed, respectively. *e dynamic contact force method is used to simulate the rock-lining contact interaction, and the arbitrary Lagrange-Euler (ALE) method is used to simulate the lining-water coupling interaction. *en a numerical simulation analysis method for combined seismic response of rock-lining-water system in hydraulic tunnel is proposed, and the detailed solving steps are given.*ismethod is used to study the seismic stability characteristics of the water diversion tunnel in a hydropower station, and the displacement, stress, and damage failure characteristics of the lining structure under the conditions of no water, static water, and dynamic water are comparatively analyzed. *e results show that the hydrostatic pressure restricts the seismic response of the lining, while the hydrodynamic pressure exacerbates its seismic response and leads to damage, separation, and slip failure appearing on the haunch, which can provide a scientific reference for the seismic design of hydraulic tunnel with high water head and large diameter.


Introduction
In order to ensure the national energy supply and optimize the water resources allocation, a large number of longdistance and cross-basin water diversion projects have been built, under construction and planned in Southwest China, which is rich in water resources. Most of the water transmission lines of these projects have to pass through high mountains, and they mainly rely on hydraulic tunnel for water diversion. Hydraulic tunnel generally inevitably passes through the high seismic intensity area in the construction process due to its long extension, so its seismic safety needs detailed analysis and demonstration [1][2][3][4]. Different from the ordinary underground tunnel, the lining structure of hydraulic tunnel not only interacts with the surrounding rock, but also interacts with the internal water. In particular, the interactions between surrounding rock, lining structure, and internal water are more complex under seismic load. Especially for the hydraulic tunnel with high water head and large diameter, the hydrodynamic pressure formed by the violently shaking internal water has a great influence on the lining stress. erefore, the research on the dynamic response characteristics of hydraulic tunnel structure under the combined actions of seismic load and hydrodynamic pressure is conducive to ensure the long-term safe and stable operation of hydraulic tunnel, which has great theoretical value and social economic benefit.
ere are many research results on the seismic response of ordinary tunnel [5][6][7][8], but the seismic response analysis of hydraulic tunnel is still relatively less. Although the hydraulic tunnel characterized by the internal water is different from the ordinary tunnel, the seismic analysis methods of ordinary tunnels are still suitable for hydraulic tunnels, including analytical method, experimental method, and numerical method. Considering the relatively complex conditions such as the nonlinearity of material parameters, the irregularity of tunnel geometry, and the dynamic interaction of fluid-structure-rock system in hydraulic tunnel, the numerical method is still an efficient and economical way to simulate the seismic response of hydraulic tunnel, which has broad application prospects. In particular, the fluidstructure interaction should be taken into account in the seismic analysis of hydraulic tunnel. Based on the dynamic finite element analysis methods for single-phase fluid and solid mediums, Deng et al. [9] proposed a dynamic explicit finite element integration format for solving the fluid-lining coupling interaction in diversion tunnel. However, this study ignores the rock-lining interaction, which cannot accurately reflect the mechanical characteristics of the lining. Chen et al. [10] built a full three-dimensional (3D) finite element model of a double-line shield tunnel for water supply in Shanghai and studied the seismic performance of the tunnel by defining the nonlinear contact interface of soilstructure interaction. Wang et al. [11] conducted a full 3D seismic analysis of the portal structure of a hydraulic tunnel and investigated the effects of hydrodynamic pressure and rock-structure interaction on the seismic response and damage characteristics of the tunnel portal. Sun et al. [12] built a numerical model considering the fluid-structure-rock interaction in the commercial software ABAQUS and evaluated the nonlinear response effects of hydraulic arched tunnel induced by the mainshock-aftershock excitation. e above studies mainly adopt the simplified added mass method [13] to consider the influence of fluid, which do not reflect the dynamic coupling characteristics of fluid-lining system. erefore, in order to scientifically evaluate the dynamic stability characteristics of hydraulic tunnel under seismic load, it is necessary to propose an effective method for simulating the combined dynamic interaction of rocklining-water system.
On the other hand, the violent shaking of fluid under seismic load belongs to a large deformation problem of mesh.
e methods adopted to describe the large deformation of mesh mainly include Lagrange method and Euler method in the finite element calculation of continuous medium. e pure Lagrange method and the pure Euler method have their own advantages, but also have serious defects. For the Lagrange method, the solution of motion equation is relatively simple, but the large deformation may lead to mesh deformity. For the Euler method, the calculation accuracy will not decrease due to the material distortion, but it needs to introduce a very complex mathematical mapping when dealing with the boundary problem, which may make a large difference between the calculation results and the real solutions. With the continuous development of computer technology, the fluidstructure coupling analysis method based on the ALE description [14,15] has been applied in many fields. e ALE method combines the advantages of the Lagrange method and the Euler method and overcomes the defects of the two methods. For the ALE method, the mesh does not depend on the fluid motion, so it can not only track the free surface of fluid but also ensure the accuracy and stability of solving the large deformation problem. At present, some scholars have tried to apply the ALE method to the seismic analysis of hydraulic tunnel. Lou et al. [16] built a full 3D fluid-tunnelsoil coupling numerical model of a large-diameter parallel diversion tunnel and used the multimaterial ALE method to calculate the seismic response of the tunnel under the action of fluid. Wang et al. [17] used the ALE method to simulate the movement and sloshing of internal water and proposed a parallel calculation method based on coupling load balance to analyze the dynamic response of a water conveyance tunnel under nonuniform seismic excitation. Although the ALE method has the unique advantages in describing the fluid motion and the fluid-structure coupling problems, the research results of its application in the seismic analysis of large-scale hydraulic tunnel are still few, and there is no widely accepted dynamic coupling analysis theory for waterlining system.
In this paper, the Lagrange method is used to describe the motions of surrounding rock and lining structure, and the ALE method is used to describe the motion of internal water. Based on the mechanical analysis models of multimaterial in hydraulic tunnel, the dynamic explicit finite element analysis methods for rock, lining, and water are discussed, respectively. e dynamic contact force method [18] is used to simulate the dynamic interaction process of rock-lining system, and the ALE method is used to simulate the dynamic coupling process of lining-water system. en a combined seismic response analysis method for rock-liningwater system is proposed by combining the above two analysis methods, which is very suitable for simulating the seismic response characteristics of hydraulic tunnel.

Combined Seismic Response Analysis Method for Rock-Lining-Water in Hydraulic Tunnel
Hydraulic tunnel is a complex underground structure composed of multimaterial including surrounding rock, lining structure, and internal water, so the seismic response of hydraulic tunnel can be regarded as the conditional combined seismic response of the three materials [19]. e mechanical analysis models of multimaterial in hydraulic tunnel under seismic load are shown in Figure 1. Figure 1( where ρ s is the density of rock or lining. u is the displacement. t is the time.
x is the spatial coordinate. i and j are the coordinate directions. b is the body force of unit mass. σ is the stress. After the rock or lining is discretized by finite element in space, the dynamic balanced equation of finite element based on Lagrange can be obtained: where M and C are the mass and damping matrixes, respectively. _ U and € U are the velocity and acceleration vectors of nodes, respectively. F int is the internal force vector of nodes, which is obtained by the internal force integration of elements. F ext is the external force vector of nodes. For rock, F ext � F + R 1 , in which F is the load vector of seismic wave and R 1 is the contact force vector exerted on rock by lining. For lining, F ext � R 2 + R 3 , in which R 2 is the contact force vector exerted on lining by rock and R 3 is the coupling force vector exerted on lining by water. e explicit center difference method is used here to solve (2), which does not require iteration and has no convergence problem. erefore, this method is very suitable for dealing with the highly nonlinear problem with complex boundary and contact conditions. When the motion state of rock or lining at time t is known, the timedomain integration formats of displacement and velocity at time t + Δt can be obtained: where Δt is the time step. e initial conditions are given as

Dynamic Finite Element Analysis Method for Water.
For internal water, its mass and momentum conservation equations based on the ALE description can be obtained when it is assumed to be incompressible viscous fluid:

Shock and Vibration
where ρ w is the water density. v is the water velocity. c is the convective velocity, and c � v − v, in which v is the mesh velocity. σ is the water stress, which obeys the Stokes constitutive equation: where μ is the dynamic viscosity of water. p is the water pressure. δ ij is the Kronecker function. Similarly, the dynamic discrete equation of finite element based on ALE can be obtained by discretizing the water in space: where L is the transfer matrix. F ext is the external force vector of nodes, and F ext � R 4 , in which R 4 is the coupling force vector exerted on water by lining. e explicit center difference method is still used to solve (8).
From the above analysis, it can be seen that, for each material in hydraulic tunnel, the motion state of nodes at time t + Δt can be obtained according to that at time t. e top priority is to solve the external load F t ext of each material. F t is known here, while R t 1 , R t 2 , R t 3 , and R t 4 are unknown. erefore, the key to the seismic response analysis of hydraulic tunnel lies in solving the interaction forces between rock, lining, and water at each time.

Dynamic Contact Analysis Method for Rock-Lining.
Many researches [20,21] have shown that the dynamic contact force method is very suitable for simulating the contact behaviors such as bond, separation, and slip between rock and lining, and this method is easy to be combined with the central difference method. Node 1 and node 2 are assumed to be a contact node pair of rock-lining system. U t 1 and U t 2 are recorded as their displacements, and R t 1 and R t 2 are recorded as their contact forces at time t.
en the rock-lining contact interface should satisfy the force balance condition: When F t is known, the predicted displacements U t+Δt 1 and U t+Δt 2 of the contact nodes at time t + Δt without considering the dynamic contact force can be calculated by using (3). en, the normal contact gap Δ 1n and tangential contact gap Δ 1t of the contact nodes can be determined: where n 1 and t 1 are the unit normal and tangential vectors of the contact nodes, respectively. e direction of n 1 is from node 2 to node 1.
In order to solve the dynamic contact force R t 1 , it is necessary to assume first that the contact nodes are in bond contact state at time t + Δt. After R t 1 is calculated, the contact state must be checked according to Δ 1n and R t 1 , and then R t 1 must be corrected accordingly. e calculation formula of R t 1 can be expressed in the form of normal component N t 1 and tangential component T t 1 : where K n and K t are the normal and tangential stiffness of contact interface, respectively.
(1) If Δ 1n < 0 and ‖N t 1 ‖ > f c A 1 , it can be judged that the contact nodes enter into separation state, and then N t where f c is the cohesive force of contact interface. A 1 is the control area of the contact nodes.
it can be judged that the contact nodes enter into slip state, and then T t 1 is corrected as where μ s and μ d are the coefficients of static and dynamic friction for contact interface, respectively. When F t , R t 1 , and R t 2 are known, the modified displacements U t+Δt of the contact nodes considering the dynamic contact force can be calculated by using (3) and (4).

Dynamic Coupling Analysis Method for Lining-Water.
Node 3 and node 4 are assumed to be a coupling node pair of lining-water system. _ U t 3 and _ U t 4 are recorded as their velocities, and R t 3 and R t 4 are recorded as their coupling forces at time t. en the lining-water coupling interface should satisfy the motion continuity condition and the force balance condition: where X is the Lagrange coordinate. χ is the ALE coordinate. e calculation formula of R t 3 can be expressed as [22] 4 Shock and Vibration where N w is the interpolation matrix of water shape function. n w is the unit external normal vector of lining-water coupling interface. σ w is the stress of water element in direct contact with the lining. e lining-water coupling process based on ALE can be briefly described as follows [23]: the water transfers force to the lining on coupling interface, while the lining transfers velocity to the water on coupling interface. erefore, the coupling force R t 4 is reflected as the velocity boundary condition _ U t+Δt 4 exerted on the water by the lining. When the motion state of the coupling nodes at time t is known, the prediction-correction thought based on ALE proposed by Wei et al. [22] is used for a reference to solve the seismic response problem of lining-water coupling system at time t + Δt. is thought is revised reasonably in order to match the dynamic contact analysis of rock-lining system, and its detailed calculation steps are as follows: (1) When F t is known, calculate the predicted displacement U t+Δt 3 and predicted velocity _ U t+Δt 3 of the lining node on lining-water coupling interface without considering the coupling force by using (3) and (4) of the lining node, the updated water mesh, and the corrected velocity of water by repeating the previous steps.

Seismic Response Analysis Steps of Hydraulic Tunnel.
It must be noted that the lining-water coupling force is ignored in the above rock-lining contact analysis, and the rock-lining contact force is ignored in the above lining-water coupling analysis. erefore, in order to correctly simulate the seismic response characteristics of each material in hydraulic tunnel, it is necessary to combine the rock-lining contact analysis and the lining-water coupling analysis. e combined seismic response analysis steps of rock-liningwater system at time t + Δt are given as follows: (1) When F t is known, calculate U

Calculation Model and Parameters.
e installed capacity of the expansion project of a hydropower station is 80 MW, and the water diversion and power generation system is arranged according to a water diversion tunnel, a surge shaft, and two units. e water intake is about 2 km away from the upstream left bank of the dam. e total length of the diversion tunnel from the intake to the surge shaft is 7126 m, whose center elevation is reduced from 80.7 m to 32.0 m, and the total length of the diversion tunnel from the surge shaft to the power house is 148.1 m. e water diversion tunnel is located in the eroded middle and low mountain landform with undulating ground along the route, and the geological conditions are relatively poor.
In this paper, a tunnel section within 80 m range downstream of the surge shaft is selected for the seismic response analysis, whose center elevation is 32.0 m. e excavation section is circular, and the tunnel diameter is 10.5 m. e lining adopts C25 concrete structure, whose thickness is 0.75 m. e maximum piezometric head of internal water is 123.0 m. A 3D finite element model of the tunnel section is built, as shown in Figure 2. e model is discretized completely by the 8-node hexahedrons, which consists of 58080 rock elements, 6336 lining elements, and 14256 water elements. e x-axis of the model is vertical to the tunnel axis along the horizontal direction, and the calculation range is from −60.0 m to 60.0 m. e y-axis is consistent with the tunnel axis, and the calculation range is from 0.0 m to 80.0 m. e z-axis is consistent with the water level elevation, and the calculation range is from −28.0 m to 154.7 m. e contact node pairs are arranged between rock and lining, and the coupling node pairs are arranged between lining and water. e lateral pressure coefficients of initial in-situ stress field are k x � 1.1, k y � 0.8, and k z � 1.0. e values of the physical and mechanical parameters of rock, lining, and water are shown in Table 1. Both the coefficients of static and dynamic friction for rock-lining contact interface are 0.6, and the cohesive force of the contact interface is 0.5 MPa.

Calculation Conditions.
e self-developed 3D dynamic finite element program for underground cavern [24] is used for calculation, in which the proposed combined seismic response analysis method for rock-lining-water system is embedded. Both rock and concrete belong to quasibrittle materials, which show obvious nonlinear deformation characteristics under dynamic cyclic load. ey will produce certain plastic permanent deformation and show obvious fatigue damage characteristics before destruction. erefore, the elastoplastic damage constitutive relation based on the

Shock and Vibration
Mohr-Coulomb yield criterion is used for both surrounding rock and lining structure. e yield criterion can be expressed as where f s � 0 represents the shear yield criterion, and f t � 0 represents the tensile yield criterion. σ 1 and σ 3 are the minimum and maximum principal stresses, respectively. σ t is the tensile strength. c is the cohesive force. N ϕ is a parameter related to the internal friction angle. D is the damage coefficient, which is expressed as where D i is the damage coefficient in the direction of i-th principal stress, i � 1, 2, 3. ε p i is the cumulative plastic deviator strain in the direction of i-th principal strain. R is the damage constant.
Viscoelastic artificial boundary condition is applied at the bottom of the calculation model, and free field artificial boundary condition is applied at the four vertical boundaries. e top of the model is built to the free surface. e representative Imperial-Valley wave is truncated for 20 s as the input load, and its peak acceleration is adjusted to 1.575 m/s 2 to match the seismic intensity of the project area, as shown in Figure 3. e input wave is incident vertically from the bottom of the calculation model. Both x-direction and z-direction seismic excitation are considered in the calculation. e x-direction input load adopts the Imperial-Valley wave, and the z-direction input load adopts 2/3 that of the x-direction.
In order to ensure the stability of dynamic calculation, the mesh size in the direction of wave propagation should be smaller than 1/12-1/8 of the shortest wavelength [25]. Most energy of the Imperial-Valley wave is concentrated within 20 Hz (the high frequency components over 20 Hz are filtered out), so the maximum mesh size of the calculation model is determined to be 6.75 m, while the maximum mesh size of the finite element model built in this paper is within 6.75 m, which can satisfy the requirements of wave propagation in surrounding rock. e calculation is divided into three working conditions: (1) no water condition; (2) static water condition (the hydrostatic pressure is considered, while the hydrodynamic pressure is not considered); and (3) dynamic water condition (both the hydrostatic and hydrodynamic pressures are considered). ree monitoring points are arranged on the top arch, haunch, and bottom arch at the middle section of the lining to monitor the hydrodynamic pressure exerted on the lining, the   Figure 4. It must be noted that the tunnel excavation, lining support, and hydrostatic pressure application should be carried out before the seismic load is applied.

Hydrodynamic Pressure Analysis.
When the dynamic coupling interaction between the internal water and the lining is considered under seismic load, the hydrodynamic pressure produced by the internal water exerted on the lining can be calculated. e hydrodynamic pressure timehistory of the three monitoring points is shown in Figure 5. It can be seen from Figure 5 that the variation laws of the hydrodynamic pressure time-history of the three points are basically consistent. At time 0-7 s, the hydrodynamic pressure gradually increases with time and increases to the maximum value due to the violent fluctuation of seismic wave. At time 7-11 s, it gradually decreases due to the decreasement of seismic strength. After time 11 s, it tends to be stable, while the seismic strength is already very small. e maximum hydrodynamic pressures of the top arch, haunch, and bottom arch are 0.256 MPa, 0.466 MPa, and 0.326 MPa, respectively, which all appear at time 5.2 s. On the whole, the hydrodynamic pressure of the haunch is larger than that of the other parts.

Lining Displacement Analysis.
e peak displacements of the three monitoring points all appear at time 14.9 s both in x-direction and z-direction. Taking the peak displacements at time 14.9 s as the maximum displacements, the maximum displacements of the monitoring points under the three working conditions in x-direction and z-direction can be obtained, as shown in Figure 6.
It can be seen from Figure 6 that in x-direction, the maximum displacements of the top arch, haunch, and bottom arch under condition (1) are 6.00 cm, 6.19 cm, and 5.88 cm, respectively. e maximum displacements of the three parts under condition (2) are 0.33 cm, 0.31 cm, and 0.27 cm smaller than that under condition (1), respectively. While the maximum displacements of the three parts under condition (3) are 0.31 cm, 0.60 cm, and 0.39 cm larger than that under condition (1), respectively. In z-direction, the maximum displacements of the three parts under condition (1) are 3.55 cm, 3.74 cm, and 3.86 cm, respectively. e maximum displacements of the three parts under condition (2) are 0.15 cm, 0.19 cm, and 0.17 cm smaller than that under condition (1), respectively. While the maximum displacements of the three parts under condition (3) are 0.19 cm, 0.24 cm, and 0.38 cm larger than that under condition (1), respectively. ese show that the lining deformation is restricted by the internal water when only the hydrostatic pressure is considered, while the hydrodynamic pressure exacerbates the displacement response of the lining. In addition, the transverse deformation of the haunch and the vertical deformation of the bottom arch are the mostly affected by the hydrodynamic pressure.

Lining Stress Analysis.
In order to explore the influence of hydrodynamic pressure on the mechanical characteristics of the lining, the maximum principal stress is taken as an example to illustrate. e maximum principal     (17) and (18)), the damage coefficient time-history of the three monitoring points can be obtained, as shown in Figure 8. It can be seen from Figure 8 that at time 0∼3 s, the damage coefficient of each part is 0, which shows that the lining is in elastic stress state. At time 3∼7 s, it gradually increases with time due to the violent fluctuation of seismic wave. After time 7 s, it basically remains unchanged. Due to the plastic strain accumulation of concrete, the damage developing process of the lining shows irreversibility. On the whole, the damage of the haunch is more serious than that of the other parts. e damage coefficient distribution of the lining (taking the middle section of the lining) under the hydrodynamic pressure after the earthquake is shown in Figure 9.
It can be seen from Figure 9 that all parts of the lining are damaged, and the damage degree of each part is different. From the view point of radial direction, the lining part in direct contact with the internal water is seriously damaged, and the closer the lining is to the rock, the smaller its damage degree is. From the view point of tangential direction, the haunch is the most seriously damaged, and the damage tends to expand to its both sides. ese show that the hydrodynamic pressure has a great influence on the mechanical behaviors of the lining and can directly affect its seismic safety. Shock and Vibration e rock-lining contact interface is basically in bond contact state before the earthquake. As the lining damage develops gradually under seismic load, the rock-lining contact interface is also gradually destructive. e separation and slip zone distributions of contact interface after the earthquake are shown in Figure 10.
It can be seen from Figure 10 that the separation and slip zones of rock-lining contact interface are mainly distributed in the haunch and its sides, while the top arch and the bottom arch are well bonded. is is mainly because under the combined actions of seismic load and hydrodynamic pressure, the contact interface of the haunch bears more load, and the damage of the haunch is the most serious, which lead to the separation and slip failure appearing on this part. erefore, the influence of hydrodynamic pressure cannot be ignored on the seismic design of hydraulic tunnel, and the haunch is the seismic weak part of lining structure.

Conclusions
(1) Combining the dynamic contact analysis of rocklining interaction and the dynamic coupling analysis of lining-water interaction, a combined seismic response analysis method for rock-lining-water system in hydraulic tunnel is proposed. Taking the water diversion tunnel project in a hydropower station as an example, the dynamic response characteristics and failure mechanisms of the lining structure under seismic load are simulated and analyzed, and the research results can provide a scientific reference for the seismic design of hydraulic tunnel with high water head and large diameter. (2) e hydrodynamic pressure produced by the internal water increases with the increasement of seismic strength, and the hydrodynamic pressure of the haunch is larger than that of the other parts. Compared with the no water condition, both the maximum displacement and maximum tensile stress of the lining decrease under the static water condition, while they both increase obviously under the dynamic water condition. ese show that the hydrostatic pressure restricts the deformation and stress response of the lining, while the hydrodynamic pressure exacerbates the deformation and stress response of the lining.  direct contact with the internal water is seriously damaged, and the haunch is the mostly seriously damaged. e separation and slip zones of rocklining contact interface are mainly distributed in the haunch and its sides under the combined actions of seismic load and hydrodynamic pressure. erefore, the haunch is the seismic weak part of lining structure.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.