Abstract

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.

1. 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 [14], 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 fluid-lining coupling interaction belongs to the typical fluid-structure 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 element discretization [9, 10]. In this method, a finite element equation of system can be derived by combining with the boundary and definite conditions [11, 12]. Studies on fluid-structure coupling problem have been quickly developed in hydraulic, marine, and other fields since 1980s, such as the liquid sloshing in the liquid reservoir vessel [1316], dam reservoir water fluctuation [1720], and the structure vibration of deep seabed tunnels [2123]. However, fluid-structure coupling analysis method is rarely applied to the high head and large diameter water conveyance tunnel. In current researches, the effect of water is simply considered with added mass method. However, the interaction law of tunnel structure and fluid cannot be reflected by added mass method, and the seismic response characteristics are not clear as well.

The lining-rock coupling interaction is a typical dynamic contact problem. The current numerical calculation methods are primarily the Lagrange multiplier method [24] and penalty function method [25, 26] and their improved versions [27, 28]. However, these methods tend to either increase the degree of freedom of the system or influence the time step of integration [29, 30], and they can adversely affect the precision and speed of calculation when applied to the dynamic calculation of underground structures, which involve a large number of contact elements and complex contact states. A dynamic contact force method aimed at solving the seismic response problem of the contact crack is put forward by Liu et al. [31, 32]. The convergence and stability of this algorithm are easily satisfied, making it suitable for large and complex contact systems. But it fails to reflect the bond-slip properties and nonlinear deformation characteristics of contact surface.

On the basis of the former researches, a new dynamic explicit finite element analysis method considering the coupling interaction of fluid and lining and the combined bearing interaction of lining and surrounding rocks is presented. This method can directly integrate, solve without simultaneous equations, and greatly simplify the calculation process. Moreover, it can be easily used to analyze 3D wave propagation problems with a variety of complex contact and coupling process and simulate the nonuniform, nonlinear deformation characteristics of contact and structural materials. The research results are of important reference value in solving aseismic safety problem and enhancing the antiearthquake designing level of pipeline structures in fields of interbasin water transfer, hydropower generation, oil-gas transmission, and so forth.

2. Seismic Response Analysis Method of Diversion Tunnel

2.1. 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; and are interaction forces of the fluid-lining coupling interface; and are the dynamic contact forces of the lining-rock coupling interface; and 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.

2.2. 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: no stickiness and no rotation or heat transfer; compressibility or almost incompressibility; small displacement; 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; means the fluid density; 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 aswhere , , and , respectively, refer to mass matrix, damping matrix, and stiffness matrix of fluid medium; is external load matrix of fluid-lining coupling interface; , , and 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 fluid-lining 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.

2.3. 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; and are Lame constants; 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 , , and , respectively, refer to mass matrix, damping matrix, and stiffness matrix of solid medium; for lining medium, , in which and , 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, , in which and , 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 velocity responses of solid medium and node force of fluid-lining coupling interface can be expressed as

With regard to fluid medium and solid medium, the following initial conditions can be introduced:

3. Fluid-Lining-Rock Coupling Interaction Analysis of Diversion Tunnel

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

3.2. Combined Bearing Analysis of Lining and Rock

The stress continuity conditions in the interface of lining and surrounding rocks are satisfied; namely, . It can be seen from (8) and (9) that the displacement of interface node at time step can be solved by displacement and contact force of interface node at time step , but the velocity of interface node at time step should be solved by displacement and contact force of interface node at time step . Therefore, in order to obtain the motion state of interface node at time step , the contact force 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-to-point contact. A certain amount of cohesive force exists between the nodes, which are in cohesive contact state. During dynamic loading process, the stress in some contact nodes would exceed the cohesive force and enter into the state of sliding contact or even separation. 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 and of contact node pair at time step can be obtained by (8). Therefore, the relative normal displacement and tangential displacement 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 aswhere is the initial normal stiffness of contact; and 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 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 into expressions (18) and (19), the normal and tangential stiffness of contact at time step can be solved. Therefore, the dynamic contact force increments of lining-rock coupling interface can be obtained:where and are normal and tangential force increments. On this basis the dynamic contact force of contact node pair at time step can be obtained as follows:where and .

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 , which shows that contact node pair and have the tendency to invade each other and if , then the node pair enter into the state of sliding contact. We havewhere and refer to the static and dynamic friction coefficients, respectively; is the cohesive force between contact surface. Before the calculation of time step , if sliding or separation has not occurred in the contact node pair, ; otherwise, .(2)If , which shows that contact node pair and have the tendency to separate each other, if , then the node pair enter into the state of separation. We havewhere refers to the ultimate tensile strength of contact.

Through the above analyses, the flow chart of fluid-lining-rock 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.

4. 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 that, there are FSI of inner water and lining structure and SSI of lining and surrounding rocks, of which the node numbers are both 696. The model ranges along -, -, and -axes are 180.0 m, 180.0 m, and 321.67 m, respectively. The complete finite element model of diversion tunnel is shown in Figure 3, and the local finite element model and layout of monitoring points are shown in Figure 4. The plane  m is chosen as the monitoring section, and point A on the vault, point B on the haunch, and point C on the inverted arch of the lining structure of monitoring section are selected as monitoring points to acquire the quantitative dynamic responses under seismic load.

It should be noted that static calculation of three-dimensional 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.

4.2. 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 , and maximum acceleration = 2.49 m/s2) 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.

4.3. Calculation Results and Analysis
4.3.1. 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: the hydrodynamic pressure increases sharply from  s to  s and diminishes gradually from  s; the hydrodynamic pressure reaches the maximum at  s, and the maximum is 0.36 MPa; 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  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: 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; hydrodynamic pressure distribution of tunnel section is almost symmetric but slightly different under the influence of vibrating direction.

4.3.2. 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 obtained from Figure 8: the displacement time curve of the vault, haunch, and inverted arch almost coincides; -displacement reaches the maximum 6.56 cm at  s and reaches the minimum −5.67 cm at  s; the displacement reaches the wave peak or trough at the same time, indicating that various locations are in a synchronous vibration state.

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.

4.3.3. 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: the maximum principal stress of lining structure is up to 0.6–0.7 MPa under the initial hydrostatic pressure; 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; 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; 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.

4.3.4. 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 load and hydrodynamic pressure, which cause the formation of tensile failure of contact surface and gradually extension to each side.

5. 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 of fluid-lining-rock coupling interaction of diversion tunnel under seismic load is presented. A diversion tunnel of some hydropower station is chosen as an engineering project to analyze the seismic response. Comparing with the condition regardless of hydrodynamic pressure, we can reach the following conclusions:(1)Under seismic load, the hydropower pressure increases sharply in a short time and decreases gradually with diminishing earthquake wave. Due to the horizontal (direction of -axes) excitation, the peak value of hydrodynamic pressure of diversion tunnel always appears on the haunch and decreases along each side to near-zero. Hydrodynamic pressure distribution along the tunnel section is slightly different under the influence of vibrating direction and liquid sloshing of tunnel structure.(2)The displacement time curves of monitoring points are consistent with the input wave, and various locations of diversion tunnel are in a synchronous vibration state. Comparing with the condition regardless of hydrodynamic pressure, the maximum displacements of monitoring points decrease somewhat but obviously on the haunch, and the locations of peak value transfer from the haunch to vault.(3)The maximum principal stress is evenly distributed around the tunnel and does not exceed the tensile strength of concrete without considering the hydrodynamic pressure. Under the opposite condition, the maximum principal stress quickly increases and gradually exceeds the tensile strength of concrete, in which the haunch of lining structure is remarkably influenced, leading to great stress increment and serious damage.(4)The lining structure mainly enters into the stage of damage and cracking under earthquake and hydrodynamic pressure, in which the haunch is greatly affected, resulting in changes of contact states, formation of local cracking and sliding failure, and gradual extension to each side.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This study is supported by the National Key Basic Research Program of China (973 program) (Grant no. 2015CB057904), the Major Program of the National Natural Science Foundation of China (Grant no. 91215301), and the National Natural Science Foundation of China (Grants nos. 51279136 and 51209164). This support is greatly acknowledged and appreciated.