Study on Fluid-Lining-Rock Coupling Interaction of Diversion Tunnel under Seismic Load

Fluid-lining-rock coupling interaction of diversion tunnel under seismic load is a critical problem in seismic research which should be solved urgently. Based on the explicit finite element method for dynamic analysis of single-phase fluid and solid medium and combining with the boundary conditions of coupling interface, a dynamic explicit finite element solving format of diversion tunnel considering fluid-lining coupling interaction is established. In light of the basic theory of dynamic contact force method and applying the nonlinear hyperbolic constitutive model of contact surface, a dynamic explicit finite element time-domain integral equation of combined bearing of lining and surrounding rocks, which takes the bond-slip behavior of the contact surface into account, is put forward. Meanwhile, considering the dynamic interaction process of inner water and lining, lining and surrounding rocks, an explicit finite element numerical simulation analysis method of fluid-lining-rock coupling interaction of diversion tunnel under seismic load is presented. The calculation results of case study reasonably reflect the seismic response characteristics of diversion tunnel, and an effective analysis method is provided for the aseismic design of hydraulic tunnel.


Introduction
The spatial distribution of water resources in China is extremely uneven, which highlights the contradiction between water supply and demand.The development and utilization of water resources is gradually turning to the upstream of each watershed.So the long distance basin water diversion and power generation become the important developing directions of water conservancy and hydropower projects in the long run.The long distance diversion tunnels usually cross high earthquake intensity region, facing serious seismic stability problem, especially located in the southwest China [1][2][3][4], like Dianzhong water diversion project, Liyuan diversion tunnel in Yunnan province, Motuo diversion tunnel in Tibet, and the west line of South-to-North water transfer project.Take, for example, Dianzhong water diversion project, almost across the whole Yunnan province, stretching for 800 km, and the earthquake intensity of project area of water diversion routes is in the degree of VII and above.The peak acceleration of the horizontal earthquake exceeding 0.15 g in this area accounts for about 78%.Therefore, researching the antiseismic stability of diversion tunnel is of great significance for guaranteeing the safe and steady operation of the project.
The key to analyze the seismic response of the diversion tunnel under earthquake excitation lies in the research of the fluid-lining-rock coupling interaction process.The fluidlining coupling interaction belongs to the typical fluidstructure interface coupling problem [5,6].Fluid-lining coupling numerical analysis methods of complex systems can be summed up into two categories: one kind is semianalytical method [7,8], which adopts finite element discretization for structure and approximate analytical formulas for fluid or turns fluid into added mass through boundary integral and simplifies the structure by using the assumed mode and aneroid mode method; the other is pure numerical method, which adopts finite element discretization for both structure and fluid or, respectively, adopts finite element and boundary

Seismic Response Analysis Model of Diversion Tunnel.
The seismic response process of diversion tunnel is a coupling interaction process of complex structures which involve three kinds of mediums and two interfaces, including the seismic response process of inner water, lining, and rocks, the FSI coupling process of inner water and lining, and the SSI coupling process of lining and surrounding rocks [33].
According to the mechanical characteristics of structures, the seismic response of diversion tunnel can be regarded as the superposition of conditional seismic response process of three mediums, as shown in Figure 1, in which Figure 1(a) is the sketch of seismic response of diversion tunnel; Figures 1(b), 1(c), and 1(d) are, respectively, sketches of seismic response of surrounding rocks, lining, and inner water;  1 and  2 are interaction forces of the fluid-lining coupling interface;  2 and  3 are the dynamic contact forces of the lining-rock coupling interface;  and  3 refer to seismic input load.
According to the finite element discretization equations of three mediums and combining with the stress and displacement boundary conditions of the coupling interface, the time domain finite element integration scheme of each medium in diversion tunnel is derived.

Dynamic Explicit Finite Element Method of Fluid Medium.
The pure fluid element based on Navier-Stokes equation or the potential fluid element on the basic of Laplace equation can be adopted as the fluid element in the finite element analysis of the interaction systems of inner water, lining, and rocks.The potential fluid element can be divided into two kinds: fluid element based on displacement and fluid element based on potential.The basic assumptions of fluid element based on displacement can be expressed as follows: (1) no stickiness and no rotation or heat transfer; (2) compressibility or almost incompressibility; (3) small displacement; (4) no actual fluid flow.Assuming that the fluid in the tunnel is a steady uniform flow and mainly considering the normal interaction between fluid and structure almost meet the above assumptions.Therefore, the fluid element based on displacement is selected in this paper.
On the basis of above assumptions, the dynamic equation of nonviscous compressible fluid (ideal fluid) is obtained: where the subscript 1 is the fluid medium;  1 means the fluid density;  1 means bulk modulus; ü and  refer to acceleration and displacement, respectively;  and  are Cartesian coordinate components.When the explicit element method is combined with multitransmitting artificial boundary, very little damping which is proportional to the stiffness matrix is introduced to eliminate high frequency instability problem.Meanwhile, the explicit finite element discretization equation of nonviscous compressible fluid can be written as where  1 ,  1 , and  1 , respectively, refer to mass matrix, damping matrix, and stiffness matrix of fluid medium;  1 is external load matrix of fluid-lining coupling interface; ü 1 , u 1 , and  1 are acceleration matrix, velocity matrix, and displacement matrix of fluid medium.
From the viewpoint of stability and calculation workload, the explicit integration scheme based on time domain weighted residuals is adopted.The time domain explicit integral format [34,35] and the node force expression of fluidlining coupling interface for dynamic analysis of nonviscous compressible fluid are obtained: where the subscript  refers to the node number of fluid medium and  is the coordinate direction of node.

Dynamic Explicit Finite Element Method of Solid Medium.
For lining and rocks, in the case of no volume force, a dynamic equation of linear elastic solid medium is obtained according to the classic linear elastic dynamics theory: where the subscript 2 refers to lining medium or rock medium;  2 and  2 are Lame constants;  2 is the density of lining or rock.Discretizing (6) with the explicit finite element method and introducing the Rayleigh damping terms can lead to the following explicit FEM equation of solid medium: where  2 ,  2 , and  2 , respectively, refer to mass matrix, damping matrix, and stiffness matrix of solid medium; for lining medium,  2 =  2 + 2 , in which  2 and  2 , respectively, refer to the external load matrix of fluid-lining coupling interface and dynamic contact force matrix of lining-rock coupling interface; for rock medium,  2 =  3 + 3 , in which  3 and  3 , respectively, refer to the dynamic contact force matrix of lining-rock coupling interface and seismic load matrix of boundary.
Likewise, the explicit integration scheme based on time domain weighted residuals is adopted.The displacement and Shock and Vibration velocity responses of solid medium and node force of fluidlining coupling interface can be expressed as With regard to fluid medium and solid medium, the following initial conditions can be introduced:

Coupling Interaction Analysis of Fluid and Lining.
The continuity conditions of interface stresses and displacements are satisfied for ideal fluid medium and lining medium.The boundary conditions can be expressed as follows.
(1) The normal stress of interface is continuous: (2) The tangential stress of interface is zero: (3) The normal displacement and velocity of interface are continuous: where the subscript 2 refers to lining medium;  is node number of fluid-lining coupling interface;  and  are normal direction and tangential direction of interface node.
Substituting the boundary condition equations from ( 12) to ( 14) into the displacement and velocity expressions (3), ( 4) and ( 8), (9), the normal displacement and velocity of interface node of fluid medium and lining medium at time step  + 1 can be acquired as follows: Substituting the normal displacement and velocity expressions of interface node into ( 5) and (10), the nodal load of interface at next time step can be solved.Then substituting the nodal load expressions into (4) and ( 9), the internal node velocity of fluid medium and lining medium at time step  + 1 can be obtained.
It should be noted that this paper is just focused on the normal load of fluid-lining coupling interaction rather than the tangential load, and the normal load of interface at the initial time is internal water pressure under the static condition.

Combined Bearing Analysis of Lining and Rock.
The stress continuity conditions in the interface of lining and surrounding rocks are satisfied; namely,  2 +  3 = 0.It can be seen from ( 8) and ( 9) that the displacement of interface node at time step  + 1 can be solved by displacement and contact force of interface node at time step , but the velocity of interface node at time step  + 1 should be solved by displacement and contact force of interface node at time step +1.Therefore, in order to obtain the motion state of interface node at time step +1, the contact force  +1 should be solved primarily by the contact conditions [36].
Before application of the dynamic loading, the contact nodes of lining and surrounding rocks belong to point-topoint 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.Under normal conditions, if relative slide does not occur between contact nodes or the relative slide is small, the contact nodes are in or approximately in point-to-point contact, which are assumptions made in this paper.
If the displacements of contact node pair are     and    after the calculation time step  and the contact force   is known, the displacements  +1   and  +1  of contact node pair at time step +1 can be obtained by (8).Therefore, the relative normal displacement Δ 1 and tangential displacement Δ 2 of contact node pair can be obtained as follows: where   is the unit normal vector of contact node and  and   are node number of lining medium and rock medium on the contact surface.
The stress-strain relationship of contact is usually nonlinear under seismic load.Meanwhile, the failure phenomena of sliding, shearing, crushing, wear, and cracking may arise along the rough surface in the case of shear deformation and the high compressive stress.Based on the nonlinear deformation characteristics of contact, the nonlinear hyperbolic constitutive model [37,38] is used to describe the change law of normal and tangential stiffness of contact.The basic equation can be expressed as where   is the initial normal stiffness of contact;   and  max are the normal relative displacement and maximum permissible normal relative displacement of contact nodes pair;  is Poisson's ratio of concrete material; among them,   and  max are concerned with the opening, strength, and roughness of contact surface, which can be obtained by empirical formula [39].Substituting displacements of contact node pair after the calculation time step  + 1 into expressions ( 18) and ( 19), the normal and tangential stiffness of contact at time step +1 can be solved.Therefore, the dynamic contact force increments of lining-rock coupling interface can be obtained: where Δ +1  and Δ +1  are normal and tangential force increments.On this basis the dynamic contact force of contact node pair at time step +1 can be obtained as follows: where Furthermore, the dynamic contact force should be discriminated and corrected according to stress state of contact node pair after every calculation time step.The specific method is as follows: (1) If ‖Δ 1 ‖ ≤ 0, which shows that contact node pair  and   have the tendency to invade each other and if ‖ +1  ‖ >   ‖ +1  ‖ + , then the node pair enter into the state of sliding contact.We have where   and   refer to the static and dynamic friction coefficients, respectively;  is the cohesive force between contact surface.Before the calculation of time step  + 1, if sliding or separation has not occurred in the contact node pair,  > 0; otherwise,  = 0.
(2) If ‖Δ 1 ‖ > 0, which shows that contact node pair  and   have the tendency to separate each other, if  +1  >   , then the node pair enter into the state of separation.We have where   refers to the ultimate tensile strength of contact.
Through the above analyses, the flow chart of fluid-liningrock coupling interaction analysis of diversion tunnel under seismic load can be shown in Figure 2. It should be noted that the subscript 3 in the flow chart refers to rock medium.

Engineering Application and Analysis
4.1.Engineering Background.Some hydropower station located in western Yunnan province of China is the fifth cascade hydropower plant of upper reaches of Lancang River, which adopts the underground powerhouse structure with four generating units, and the installed capacity is 1900 MW (4 × 475 MW).The diversion tunnel of the hydropower station spans 7126 m, which goes across different geological tectonic zones.The tunnel section from 2 + 500.00 to 2 + 680.00 located in the middle of water diversion system is chosen for analysis in this paper, of which the seismic intensity is up to VII.The tunnel section is almost horizontal, and it is designed in circular structure, with the excavation diameter 10.2 m, inner diameter 9.0 m after secondary lining, center elevation of tunnel 32.00 m, buried depth of tunnel 241.67 m, and maximum inner head 135.24 m.The reinforced concrete structure is applied in the tunnel with C25 concrete, and the lining thickness is 0.6 m.
Three-dimensional finite element model of diversion tunnel including fluid, lining structure, and surrounding rocks is established, which consists of 64003 nodes and 59584 hexahedron elements (eight nodes), and the lining and fluid element numbers are 2688 and 5264, respectively.Besides

Yes Yes
No Yes Calculate R n+1 of FSI using (10) Calculate f n+1 2 of FSI using (5) Calculate f n+1 1 of FSI using (16) of FSI using ( 15 It should be noted that static calculation of threedimensional finite element should be conducted firstly to simulate the excavation and lining support of diversion tunnel in order to provide the seismic analysis with initial stress and displacement conditions.
The mechanical parameters of the rock mass and lining material are shown in Table 1.The strain rate strengthening effect of dynamic parameters under seismic load is ignored.

Calculation Conditions.
The dynamic analysis of this paper is performed with an in-house FEM numerical simulation platform, SUCED [35], which is designed for estimating seismic damage of large underground caverns group.In the process of numerical calculation, the elastic-plastic-damage constitutive model based on Mohr-Coulomb yield criterion is used for surrounding rocks, while the plastic-damage model is adopted for lining structure, which is proposed by  Lubliner et al. [40], improved by Lee and Fenves [41,42], and generalized for three-dimensional model by Omidi and Lotfi [43].
Free field boundary condition is applied to absorb the reflection waves along the two vertical boundaries which are perpendicular to -axes, and viscoelastic artificial boundary condition is used to absorb the incident waves from the bottom of the model [44].
In this study, the acceleration time history curve of the wave of the 1940 EI Centro earthquake (north-south, magnitude  = 6.7, and maximum acceleration = 2.49 m/s 2 ) is adopted.According to existing codes, the peak value is adjusted to 0.17 g (g is gravity acceleration), which is equivalent to the seven-degree fortification criterion, and duration time is 20 s.The horizontal (direction of -axes) vibration is considered in this paper, and after pretreatment, the -acceleration time history curve which is input from the bottom of the limited area is shown in Figure 5.

Hydrodynamic Pressure Time History Analysis.
The hydrodynamic pressure on the haunch is greater than other parts when considering the horizontal (direction of -axes) excitation.Through calculation, the hydrodynamic pressure time history curve of monitoring point B is shown in Figure 6, which indicates the following: (1) the hydrodynamic pressure increases sharply from  = 2 s to  = 6 s and diminishes gradually from  = 14s; (2) the hydrodynamic pressure reaches the maximum at  = 4.72 s, and the maximum is 0.36 MPa; (3) the hydrodynamic pressure time history curve of monitoring point does not conform with the acceleration time history curve, which indicates that water compressibility has a great influence on hydrodynamic pressure.
In order to reflect the distribution laws of hydrodynamic pressure around the monitoring section in the process of earthquake, three moments  = 4.72 s, 10.06 s, and 14.56 s are chosen, and the pressure distribution of each moment at monitoring section is obtained in Figure 7, in which monitoring points A, B, and C correspond to 90 ∘ , 180 ∘ , and 270 ∘ , respectively.The hydrodynamic pressure distribution around the monitoring section indicates the following: (1) in the case of vertical incidence and horizontal excitation of earthquake wave, the peak value of hydrodynamic pressure always appears on the haunch, decreases along each side, and is close to zero on the vault and inverted arch; (2) hydrodynamic pressure distribution of tunnel section is almost symmetric but slightly different under the influence of vibrating direction.

Displacement Time History Analysis of Lining Structure.
-displacement time curves of monitoring points of lining structure with considering hydrodynamic pressure or not are shown in Figures 8 and 9.The conclusions can be By comparing with the displacement time curves of monitoring points under two conditions, we can reach the following conclusions.The displacement of each monitoring point decreases when considering the hydrodynamic pressure, in which the displacements of vault and inverted arch slightly change, while the haunch remarkably decreases, and the maximum of decreasing displacement is about 10 mm, which seems small, but the relative deformation of 10 mm can make great influence on the mechanical properties of lining structure.These conclusions indicate that the haunch of tunnel structure is significantly affected by the hydrodynamic pressure, which impedes the deformation of lining structure.

Principal Stress-Time History Analysis of Lining Structure.
Since the compressive strength of concrete is far more than the tensile strength, the damage and failure modes of lining structure under seismic load mainly refer to tensile failure.This paper is focused on the change laws of maximum principal stresses of monitoring points.The maximum principal stress-time curves of monitoring points of lining structure with considering hydrodynamic pressure or not are shown in Figures 10 and 11, respectively, which indicate the following: (1) the maximum principal stress of lining structure is up to 0.6-0.7 MPa under the initial hydrostatic pressure; (2) the maximum principal stress quickly increases and gradually exceeds the tensile strength of concrete when considering the hydrodynamic pressure, in which the haunch of lining structure is remarkably influenced, leading to great stress increment and serious damage; (3) the maximum principal stress gradually decreases and is evenly distributed around the tunnel without considering the hydrodynamic pressure, in which the maximum principal stress of haunch changes slightly and is almost under 0.6 MPa; (4) the hydrodynamic pressure has a significant influence on the stress state of lining structure, and it is of great importance in the antiearthquake design of tunnel to consider the coupling interaction of inner water and lining structure under seismic load.

Failure Modes of Contact Surface between Lining and
Surrounding Rocks.The interaction process of lining and surrounding rocks under seismic load is simulated with three-dimensional explicit dynamic finite element methods.Before seismic action, the stress state of contact of lining and rocks is perfect, and the contact lies in cohesive contact state.Under cyclic seismic load, the damage and failure of lining structure gradually occur and develop, leading to the failure of contact surface between lining and rocks.The haunch of tunnel structure is remarkably influenced by earthquake and hydrodynamic pressure, causing serious damage, which finally leads to a range of separating and sliding zones.The failure modes of contact surface under seismic action can be shown in Figure 12.
It can be seen from Figure 12 that under the earthquake and hydrodynamic pressure, damage and cracking of lining structure develop gradually, leading to loads transfer to the deep layer, and make the contact surface of lining and surrounding rocks bear more load, which result in changes of contact states and formation of local cracking, sliding phenomena.The cracking failure of lining structure on the haunch is serious because of the great influence of earthquake

Conclusions
Based on the seismic response mechanism of diversion tunnel, the coupling interaction process of fluid-lining and lining-rock under seismic load is studied, and then an explicit finite element numerical simulation analysis method

Figure 1 :
Figure 1: Seismic response analysis model of diversion tunnel.

1 ,Figure 2 :
Figure 2: The flow chart of fluid-lining-rock coupling interaction analysis of diversion tunnel.

Figure 5 :
Figure 5: -acceleration time history curve of input seismic wave.

Figure 6 :Figure 7 :
Figure 6: Hydrodynamic pressure time history curve of monitoring point B.

Figure 8 :
Figure 8: Displacement time curves of monitoring points considering hydrodynamic pressure.

Figure 9 :
Figure 9: Displacement time curves of monitoring points regardless of hydrodynamic pressure.

Figure 10 :
Figure 10: Maximum principal stress-time curves of monitoring points considering hydrodynamic pressure.

Figure 11 :
Figure 11: Maximum principal stress-time curves of monitoring points regardless of dynamic water pressure.
load and hydrodynamic pressure, which cause the formation of tensile failure of contact surface and gradually extension to each side.

Table 1 :
Mechanical parameters of materials.