Abstract
The seismic response of hydraulic tunnels is a complex nonlinear process. What makes the case even more interesting is that the large amount of water in hydraulic tunnels which is likely to induce considerable hydrodynamic pressure acted on tunnel structures during earthquakes. In this work, a full three-dimensional (3D) dynamic finite element model is adopted to conduct a comprehensive assessment of the seismic behaviors of a hydraulic tunnel system. In this analysis, the plastic-damage model is employed to reflect the nonlinear mechanical behaviors of the concrete lining, and a fluid-solid coupling method based on an explicit weighted residual approach is proposed to consider the effects of the hydrodynamic pressures on the seismic response of the hydraulic tunnel. The numerical results indicate that the hydrodynamic pressure contributes to a greater seismic response of the tunnel structure. When the hydrodynamic pressure is considered, the magnitudes of the maximum principal stresses are likely to increase by 50% and the displacement amplitudes are approximately 2 cm more than that of without hydrodynamic pressure. The hydrodynamic pressure exacerbates the damage degree of the tunnel structure, and the waist suffers the most severe damage.
1. Introduction
Due to the population growth and urbanization of China in the past dozen years, water resources in some regions are insufficient to meet daily domestic and industrial water use; thus, many hydraulic tunnels have been built to solve the water shortage problem. Generally, tunnels in seismically active areas cannot avoid the earthquake-induced risks. Several strong earthquake events, such as the 1999 Chi-Chi Earthquake [1], the 1999 Turkey Koceali Earthquake [2], the 2004 Mid-Niigata Prefecture Earthquake [3, 4], and the 2008 Wenchuan Earthquake [5–7], demonstrated that those tunnels may suffer different levels of seismic damage, ranging from lightly cracking to severely collapse. Therefore, the seismic stability of the tunnels in seismically active areas becomes an important engineering problem and arouses much attention [8].
Great efforts have been made by many researchers to investigate the seismic behaviors of tunnels using various methods (e.g., the field investigation, scaled model tests, theoretical analysis, and numerical simulation). A systematic damage investigation of 57 tunnels was conducted and their damage patterns were summarized by Wang et al. [1] after the Chi-Chi Earthquake. Tao et al. [9] adopted 3D elastic shell theory to analyze the dynamic response of a tunnel portal and verified the analytical results by shaking table tests. Naggar et al. [10] built a closed form solution for in-plant moments and thrusts in tunnel linings and studied the impacts of the incident angles on the seismic responses. Yashiro et al. [4] proposed a simple classification method for the seismic damage forms of tunnels and summarized their damage mechanism. A 3D input method of shear waves was employed by Huang et al. [11] to investigate the impacts of the incident angles on the seismic response of a long lined tunnel. More relevant works on the seismic behaviors of tunnels can be found in [12–15].
The previous work [15] classified the seismic damage suffered by underground structures into three categories: damage due to shaking; damage due to permanent ground deformations; ground failure, including portal failure. The deformation modes of tunnels under seismic loads are generally considered to be particularly relevant to the seismic damage [16]. Hashash et al. [13] summarized several deformation forms, namely, compression and extension along the tunnel axis; longitudinal bending caused by the components of seismic waves perpendicular to the longitudinal axis; ovaling or racking shear deformation of tunnel sections; compression deformation of tunnel sections. In addition, a series of parameters, such as geological conditions, distance to the epicenter, seismic intensity, maximum ground acceleration/velocity/displacement, etc., were demonstrated to be particularly relevant to the seismic behavior of tunnels [12]. Those works effectively enhance our understanding of the seismic behaviors of tunnels and provide a guidance for the seismic design of tunnels. Most of those works are focused on traffic tunnels; however, hydraulic tunnels are of more complicated stress states because they are subjected to considerable hydrodynamic pressure. Many researchers have realized the significant effect of hydrodynamic pressure on the seismic behavior of hydraulic tunnels. A fluid-solid coupling method in frequency domain was proposed by Chen et al. [17] to study the seismic response of a water-conveyance tunnel. Yu et al. [18] and Chen et al. [19] adopted the added mass method proposed by Housner et al. [20] to simulate the hydrodynamic pressure. Clearly, the added mass method is a simplified method and cannot properly reflect the shaking of the large mass of water in hydraulic tunnels. In fact, the shaking of the inner water is a complicated nonlinear process and has a significant impact on the seismic response; thus, systematic and sufficient studies on the response performance are undoubtedly required.
Due to the complicated nonlinear response of hydraulic tunnels under earthquake loads, which can only be solved accurately and economically with the aid of numerical methods, a full three-dimensional (3D) FEM model for the tunnel system was conducted here to assess the seismic response of a hydraulic tunnel; in this model the fluid-solid coupling was considered and its impacts on the seismic performance were evaluated. Finally, the nonlinear response process was presented and discussed in detail.
2. Problem Layout and Establishing of FEM Model
2.1. Problem Layout
A hydraulic tunnel, which is located in the seismically active areas of China, is taken as a practical example to analyze the nonlinear seismic response process of hydraulic tunnel. The tunnel has a buried depth of approximately 240 m and is designed in circular section with a diameter of 9 m. The surrounding rock of the tunnel is mainly composed of coarse metamorphic volcanic breccia and is not crossed by active faults. Due to the fact that the tunnel is subjected to considerable internal water pressure, whose value is approximately 0.78 MPa, a primary and a secondary support are adopted for the tunnel. The primary support consists of 0.2 m of shotcrete and steel rock bolts. The rock bolts, with length of 4.5 m and diameter of 25 mm, are distributed on a grid of 2.0 m×2.0 m. The secondary support is composed of reinforced concrete with a thickness of 0.6 m. The concrete used for the tunnel support is C25-grade [21]. Two reasons contribute to the evaluation the seismic performance of the tunnel: one is the location of high seismic intensity, the other is the large mass of water in the tunnel, which is expected to induce high hydrodynamic pressure. According to the seismic ground motion parameters zonation map of China [22] prepared by the China Earthquake Administration, the tunnel studied in this work is located in the area with expected basic intensity of seven-degree; thus, the earthquake-induced risks cannot be ignored.
2.2. Geometry of the FEM Model and Input Seismic Motion
In order to simulate and analyze the problem depicted in Section 2.1, a 3D time-history analysis of the hydraulic tunnel has been carried out based on an in-house dynamic finite element numerical simulation platform [23], which is designed for simulating the seismic response and evaluating the seismic damage of large underground structures. A limited computing model is intercepted to establish 3D finite element model based on the actual terrain and size of the tunnel (see Figure 1). The directions of the tunnel longitudinal and vertical axes, as shown in Figure 1, are defined as the x-axis and z-axis in the global coordinate system of the computing model, and the y-axis is determined by the right-hand rule. The length of the model is 183 m, the width 180 m, and the maximum height 274 m. The entire finite element mesh, in which the surrounding rock, the shotcrete, the secondary lining, and the water are included, has been constructed and discretized with eight-node hexahedral elements. A total of 60536 solid elements are meshed, among which 2688 elements represent the secondary lining, and 5264 elements represent the fluid in the tunnel. Due to the lack of measured seismic waves, the acceleration time-history of the 1940 El Centro earthquake is adopted in this paper and its first 20s acceleration time-history is intercepted as the input seismic wave for the dynamic calculation (see Figure 2). The x-axis direction motion is considered in the calculation, and the maximum amplitude is adjusted to of 1.96 m/s2 (0.2 g). The motion is applied to the bottom boundary of the model and propagates upward. The input of the motion is conducted by transforming the motion into the equivalent nodal forces applying to the bottom boundary nodes. More details about the input method can be found in the works of Liu at al. [24]. It should be noted that the tunnel and the surrounding rock are assumed to be tied together on the tunnel-rock interface and the relative displacement on the interface is not taken into consideration [25].
2.3. FEM Model for Reinforced Concrete Lining
2.3.1. Material Model for Concrete
The constitutive model for concrete is the key to describe its mechanical properties and to understand its failure mechanism. The destruction of concrete in earthquakes is a gradual process; from initiation, expansion, and accumulation of microcracks to appearance of macrofractures, the destruction of concrete can be regarded as a process of damage accumulation. Additionally, the mechanical properties of concrete materials are significantly different when subjected to tensile and compression stresses. The nonlinear behavior under compression state is mainly due to stiffness degradation caused by the damage, while under tensile state it is induced not only by the damage but also by the plastic softening; thus, it is necessary to fully consider these nonlinear mechanical characteristics of concrete materials. In this paper, the plastic-damage model proposed by Lubliner at al. [26] is adopted because it can reflect the mechanical characteristics described above (see Figure 3). In the model, the constitutive relationship can be written aswhere is the strain tensor and can be divided to elastic strain tensor and plastic strain tensor ; is the effective stress; is the initial elasticity modulus; means a scalar damage parameter and can be evaluated by interpolating between the tensile damage variable and the shear damage variable where and are functions of the stress introduced to represent stiffness recovery effects associated with stress reversals. The tensile damage variable and the shear damage variable can be obtained by the function proposed by Birtel et al. [27]
The yield function modified by Lee at al. [28] is used to account for different evolution behaviors of the strength under tension and compression. In terms of effective stress, the yield function can be expressed aswhere and are dimensionless material constants; is the Mises equivalent effective stress; is the equivalent hydrostatic pressure; is the algebraically maximum eigenvalue of the effective stress; is the equivalent plastic strain; ; ; is the parameter of the yield curve on the plane of deviatoric stresses.
The plastic strain rate can be calculated according to the nonassociated flow rule, and the plastic potential function is expressed in the form of effective stress space. The plastic strain rate can be described aswhere is the plastic multiplier and can be obtained by the consistency condition; the Drucker-Prager hyperbolic function is adopted as the flow potential where is the parameter related to the properties of concrete materials; is the first invariant of the principal stress; is the second invariant of the deviatoric stress.
2.3.2. Finite Element Simulation of Lining Element
Two assumptions are introduced here to simulate the steel bars in concrete lining: one is that the steel bars are uniformly distributed; the other is that the steels and concrete are well bonded before cracks occur and are under same strain. Based on the assumptions, the stiffness of the lining can be considered to be the superposition of the stiffness of the concrete elements and the steelswhere is the stress transformation matrix; is the elastic matrix of concrete; is the elastic matrix of steel and can be obtained according to the method in [29]; is the reinforcement ratio.
2.4. Analysis Procedure
In order to approximate the actual construction process, the analysis is divided into 2 steps. In the first step, the static excavation analysis is conducted. In this process, the gravity field is adopted as the initial geostress field for the static calculation due to the lack of measured geostress data. The effects of the supports, i.e., lining and rock bolts, are determined by the method of Chen et al. [30]. In the next step, the results of the static calculation are used as the initial conditions for subsequent seismic response calculation. The viscous spring boundary by Liu et al. [31] is applied to the lateral and bottom boundaries of the model, and the Rayleigh damping is adopted with the critical damping ratio of 5%. Moreover, the mechanical behavior of the rock mass is simulated by the Mohr-Coulomb failure criterion with a tension cut-off limit, while that of the concrete is simulated by the plastic-damage model described in Section 2.3.1. Meanwhile, the fluid-solid coupling process of the inner water and the lining is considered by the method in Section 3. The geomechanical properties for the calculations are listed in Table 1. Furthermore, since the explicit method is used in the seismic response calculation to solve the seismic wave field, an appropriate time step is required to ensure the accuracy of the calculation. can be given by [32]where is an empirical coefficient; is the wave velocity; is the minimum size of element.
3. Fluid-Solid Coupling Mechanism
The seismic responses of hydraulic tunnels are highly related to the fluid-solid coupling. Two basic methods are widely used to analysis to the fluid-solid coupling effects: the added mass method proposed by Westergaard et al. [33] and modified by Housner et al. [20] and the multimaterial Arbitrary Lagrangian Eulerian (ALE) method [34]. The added mass method is a simplified method as it does not consider the sloshing of the inner water. The ALE method can effectively simulate the movement and sloshing of the inner water, but its modeling is of much difficulty, and the high computational costs are unacceptable in engineering practice. In this paper, a dynamic contact model is proposed to solve the problems in the above two methods. Firstly, the weighted residual approach is adopted in this model to solve the motion equations of the fluid and the structure separately. And then, the explicit solutions of the motion equations are derived according to the constraint conditions of the fluid-solid coupling interface.
3.1. Explicit Finite Element Time Integration of Fluid Medium
Since the main purpose of this paper is to study the response of the hydraulic tunnel under the horizontal lateral excitation (x-axis direction in Figure 1), it can be generally assumed that the fluid in the tunnel is a steady uniform flow and there is no tangential interaction between the fluid and the structure. Based on the above assumptions, the displacement-based fluid element [35] can be adopted here to simulate the fluid in the tunnel and the dynamic equilibrium equation of the fluid element after finite element discretization can be obtained [36]where the subscript refers to fluid medium; , , and are the vectors for acceleration, velocity, and displacement of the fluid elements, respectively; , , and are the matrixes for mass, damping, and stiffness of the fluid elements; is the contact force vectors on the fluid-solid coupling interface.
Solving the motion equation is the key of the fluid-solid coupling analysis. In general, the motion equation can be solved by the explicit or implicit algorithms. Due to the high efficiency of explicit algorithms in solving complex nonlinear systems, an explicit algorithm, that is weighted residual approach in time domain, is adopted here. The explicit integral equations of fluid medium at time can be expressed aswhere is the time step; the subscript refers to the coordinate direction.
According to (11), it is clear that the motion state of the fluid medium at time can be calculated by the displacement, velocity, and contact force at time . The displacement and velocity at time can be obtained by the recursive method, while the contact force depends on the constraint conditions of the fluid-solid coupling interface.
3.2. Explicit Finite Element Time Integration of Solid Medium
After the finite element discretization, the dynamic equilibrium equation of the solid medium under seismic action can be also defined aswhere the subscript refers to solid medium; is the known external forces vector.
Similarly, the motion equation (12) can be also solved by the weighted residual approach in time domain. The explicit integral format of the motion states (i.e., displacement and velocity) and the contact force can be given by
The stress results of the static excavation calculation are taken as the initial conditions of (11) and (13), and the initial displacement and velocity of the fluid and solid mediums are assumed to be 0; then there are
Obviously, it can be seen from (11) and (13) to (16) that the motion state of the fluid and solid mediums at any time can be derived from the motion state of the previous moment and the contact force at the corresponding moments. Therefore, the key to acquire the movement state is to calculate contact forces and based on the coupling interaction conditions. In addition, it should be noted that the normal component of the contact force at the initial moment is the internal water pressure under static conditions, while the tangential component is ignored.
3.3. Coupling Interaction Conditions
According to the contact characteristics of the fluid and the structure of the pressure tunnel, the normal stress and displacement on the contact surface of the fluid and the structure can be considered to be continuous; thus, several coupling interaction conditions can be introduced.
If the normal stress on the contact surface is continuous and the tangential stress is ignored, here is
If the normal displacement and velocity on the contact surface are continuous, here iswhere the subscripts 1 and 2 are the two tangential directions of the contact surface and the subscript 3 is the normal direction.
Substituting (17), (18), (19), and (20) into (11) and (13), we can have
Substituting (21) and (22) into (11) and (13), the contact force and at time can be obtained accordingly. It should be noted that the dynamic contact model in this paper simplifies the fluid-solid coupling interaction because the effects of tangential component of the contact force and liquid shaking are not fully considered. Nevertheless, since the tangential contact force is not the dominant force affecting the structure, and the liquid shaking is not significant in pressure tunnels, this simplification has little impact on the result of our research and the conclusions that we tried to obtain are still valid.
4. Results and Discussions
The work in this paper aims at the nonlinear seismic response characteristics of the tunnel and the effects of the hydrodynamic pressure on the response process. The detailed analysis process is as follows. First, the stress and displacement response process of the tunnel structure are conducted; in this process, two cases, one considers the hydrodynamic pressure and the other does not, are contrasted to study the impacts of the hydrodynamic pressure. Finally, the damage characteristics of the tunnel structure are analyzed in detail. Specifically, in order to illustrate the results, a typical critical cross section is selected and three monitoring points in this section are chosen; their positions are given in Figure 4. These points are representatives of the following locations: crown (A), waist (B), and inverted arch (C).
4.1. Stress Response Characteristics
The hydrodynamic pressure can be calculated by the fluid-solid coupling method introduced in Section 3. The maximum principle stresses are adopted here to assess the tension states of the tunnel structure in the below two cases: the hydrodynamic pressure is not considered in case 1, while the hydrodynamic pressure is considered in case 2. Figure 5 shows the time histories of the maximum principle stresses of the monitoring points. As plotted in Figure 5, the maximum principal stress of the crown, waist, and inverted arch before the earthquake (t = 0 s) is about 0.20 MPa. At this time, the stress state of the lining structure is mainly controlled by the internal water pressure. After the input of seismic waves, the time histories curves of the stresses of the three monitoring points begin to increase rapidly, and the principal stresses experience a violent fluctuation in the period of 0-5 s. The magnitude of the maximum principal stress, which is usually regarded as the indicator of the damage degree, is of primary concern because the concrete material is well known to be of weakly tensile capacity. The maximum principal stress at the waist (point B) is larger than that of the crown (point A) and the inverted arch (C) under the horizontal excitation (x-axis direction in Figure 4), and their magnitudes even exceed the maximum tensile strength (1.25MPa) when the hydrodynamic pressure is considered (Figure 5(b)). In addition, it can be seen from Figures 5(a) and 5(b) that the magnitudes, in case 2, are nearly 50% larger than that in case 1. Therefore, we can conclude that the waist experiences a larger tension stress than other positions, and the inner water contributes a larger stress response in the tunnel.
(a) Without hydrodynamic pressure
(b) With hydrodynamic pressure
4.2. Displacement Response Characteristics
Figure 6 displays the time-history curves of the displacement of the monitoring points in the two cases: the hydrodynamic pressure is considered or not considered. It can be seen that the waveforms of the curves are roughly similar to each other, and their peaks appear almost at the same time, indicating that the displacement fluctuation is mainly controlled by the input seismic waves. The displacement responses of different positions are evidently varied, and the waist experiences a larger displacement amplitude than other positions. This result is consistent with the stress response in Section 4.1. The different displacement responses of the different positions also suggest that the tunnel suffers structural deformation. Comparing case 1 with case 2, the major difference of the two cases lies in their displacement amplitudes. When the inner water is considered (case 2), the displacement amplitudes of the monitoring points are approximately 2 cm more than that in case 1. The different results of the two cases are mainly due to the hydrodynamic pressure induced by the inner water. Ultimately, it can be concluded that the hydrodynamic pressure induces a larger displacement response of the tunnel and it cannot be neglected.
(a) Without hydrodynamic pressure
(b) With hydrodynamic pressure
4.3. Seismic Damage Characteristics of the Tunnel Structure
The damage factor of the tunnel structure, which can be obtained according to formula (3), is adopted here to quantitatively show the damage degree of the tunnel structure and where suffers the most severe damage. In this paper, parts of the tunnel structure, which are adjacent to the monitoring section, are selected for analysis. The results shown in Figure 7 are obtained with the hydrodynamic pressure considered and are plotted in every 20 seconds during this numerical simulation. Figure 7 offers the damage evolution process of the tunnel structure in the two cases: the hydrodynamic pressure is considered (Figure 7(b)) or not considered (Figure 7(a)). It is clear that the damage intensity increases with time. When t = 5 s, the damage coefficients of most positions are approximately 0.2 in the two cases, and the maximum damage coefficient is close to 0.5 or 0.6. Moreover, it is interesting to note that the waist suffers more serious damage than other positions. As has been expected, the results can be explained that the waist experiences a larger tension stress induced by the horizontal seismic excitation. As time increases, the damage continues to extend from the waist to the crown and the inverted arch. At the end of the numerical simulation, namely, at t = 20 s, the tunnel structure suffers the most severe damage, and the maximum damage coefficient is approximately 0.7 in case 1 and 0.9 in case 2. In addition, the major difference in the results of the two cases lies in the damage degree. When the hydrodynamic pressure is considered, the damage coefficient is obviously larger than that in case 1 at the same time. This result is consistent with the stress and displacement responses. The comparison between the two cases shows that the hydrodynamic pressure has a significant effect on the seismic response of the tunnel structure.
(a) Case 1
(b) Case 2
5. Summary and Conclusions
In order to analyze the nonlinear seismic response process of a hydraulic tunnel, a full 3D finite element model is established. For approximating actual field conditions, the plastic-damage model is adopted to simulate the nonlinear mechanical behavior of the concrete lining under seismic loads, and an explicit fluid-solid coupling method is proposed to consider the effects of the hydrodynamic pressure. When we apply the model to the seismic response analysis of a diversion tunnel, the conclusions of the research are as follows:
The waist of the tunnel lining experiences a larger tension stress than other positions due to the actions of the horizontal lateral seismic loads. The hydrodynamic pressure contributes a larger stress response in the tunnel. And compared with the case of without hydrodynamic pressure, the magnitudes of the maximum principal stresses in the case with hydrodynamic pressure are likely to increase by 35%.
The displacement time-history curves of the monitoring points oscillate in synchrony, and their peaks appear almost at the same time, indicating that the displacement fluctuations are mainly dominated by the input seismic waves. The displacement responses of different positions are evidently varied, and the maximum displacement occurs on the waist of the tunnel lining. This different displacement response suggests the tunnel suffers structural deformation. In the case of with hydrodynamic pressure, the displacement amplitudes of the monitoring points are approximately 2 cm more than that of without hydrodynamic pressure.
The hydrodynamic pressure has an evident influence on the seismic response and the damage distribution of the tunnel structure. When the hydrodynamic pressure is considered, the damage coefficient is obviously larger than that in the case without hydrodynamic pressure. The waist suffers the most severe damage, and the maximum damage coefficient is, respectively, 0.7 and 0.9 in cases without and with hydrodynamic pressure. In addition, although the amplitude of input seismic wave tends to be gentle after t = 10 s, the damage degree of the tunnel lining continues to increase with time, suggesting that the damage degree of the tunnel structure is related not only to the amplitude of the input seismic wave, but also to the seismic wave duration.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study is supported by the National Natural Science Foundation of China (Grant no. 41602273) and the Geological Survey Project “Zhengzhou Urban Geological Survey” of China Geological Survey (Grant no. DD20189262). The supports are acknowledged and greatly appreciated.