Parameter Analysis on Seismic Response of Long Lined Tunnel by 2.5D Substructure Method

Abstract

When the numerical analysis of a long lined tunnel is carried out, the calculation amount of the finite element model becomes restricted large-scale parameter analysis. In this paper, an efficient and high-precision 2.5-dimensional (2.5D) frequency-domain finite element method is used to simulate the three-dimensional response of tunnels under the action of oblique incident plane seismic waves. This method can save calculations and avoid the boundary effect caused by the longitudinal truncation of the tunnel. The 2.5D zigzag-paraxial boundary is developed. The artificial boundary is attached to the structure’s surface. The substructure method for oblique plane seismic waves is established. Comparing the substructure method with the analytical solution, the correctness of the site response is verified first. The accuracy of the 2.5D finite element substructure method is further verified. The parameter analysis of different incident angles and conversion angles shows that the underground tunnel does not reach the maximum of structural seismic response when the seismic wave is vertically incident. The location of the soil–rock interface on the tunnel is further discussed. The results show that when the underground tunnel crosses the location of the soil–rock interface, the seismic response of the tunnel will be amplified.

1. Introduction

With the increasing scale of urban underground space construction, the safety performance of underground structures has been of wide concern by researchers [1,2,3]. The observation results of the 1995 Kobe earthquake show that subway stations and interstation tunnels will suffer serious damage under strong earthquakes [4].

At present, researchers mainly use model tests [5], analytical methods [6,7] and numerical methods [8] to study the dynamic interaction of underground structures. Due to the lack of prototype observation data, there is not enough data to reveal the earthquake damage mechanism and failure mode. The model experiment can simulate the damage process during an earthquake. It is a widely used seismic research method at present. However, it is limited by the geometric scale of the model, the boundary effect and other conditions, so the results will have a certain impact. The analytical method established to analyze the soil–structure dynamic interaction under seismic waves has the characteristics of a clear concept and physical meaning, but the research work is only applicable to elastic half-space, simple site conditions or ideal geometric conditions. The rapid development of computer and numerical analysis technology has significantly improved the ability to use numerical methods to obtain quantitative results. Therefore, the numerical simulation method has become one of the most efficient methods in the engineering field [9,10]. When the numerical method is used to solve the soil–structure dynamic interaction problem, the whole analysis method or substructure method based on the finite element method is often used. The whole analysis method refers to the establishment of the near-field whole model including the structure and a part of surrounding media. At the cut-off boundary of the near field model, dynamic artificial boundary conditions are applied to simulate the radiation damping effect of truncated far field media. At the same time, the earthquake input at the artificial boundary of the near field model is realized based on the free field response. Different from the whole analysis method, the substructure analysis method only establishes the structural finite element model. The truncated infinite domain is simulated through the foundation impedance function or the frequency domain dynamic stiffness. The impedance function or the frequency domain dynamic stiffness are integrated into the structural finite element model through the soil–structure interface condition. The earthquake input is realized at the soil–structure interface based on the free field response.

A suitable method for the seismic design of a tunnel isometric underground structure is to use a three-dimensional (3D) finite element model to conduct the whole dynamic time-history analysis of the tunnel–soil interaction system. This method can effectively simulate the dynamic interaction between tunnel and soil [11]. The deformation and internal force at each moment can be obtained in the process of structural seismic response. In order to improve the computational efficiency, Yu et al. [12] proposed a multi-scale finite element method that combines coarse mesh division with local fine mesh division. Based on the multi-scale finite element method, the seismic response of a long water transmission pipeline is analyzed. Huang et al. [13] used three-dimensional dynamic finite element time-history analysis to study the influence of the incident P-wave of seismic wave on the internal force and deformation of a long lined tunnel. At the same time, the truncation effect at the end of the tunnel is discussed in detail. The large size of the three-dimensional finite element model in practical engineering applications, the large number of divided units and the high cost of calculation put forward high requirements of the computing ability, thus limiting the application of the three-dimensional dynamic finite element analysis method [14].

For the long lined underground structure with a constant longitudinal direction, special numerical calculation methods such as the 2.5D finite element method can be considered. At first, Hwang and Lysmer [15] proposed a special wave transmission element to denote the traveling wave effect of seismic wave propagation. Some researchers used this idea to further analyze the scattering effect of seismic waves in irregular sites [16,17,18,19]. Through Fourier transform, Yang and Hung [20] established a two-dimensional (2D) plane element with three degrees of freedom to solve the soil–structure dynamic interaction, and named it the 2.5D method. Subsequently, Alves Costa et al. [21] used the 2.5D finite element method to study the influence of soil nonlinearity on the dynamic response of a high-speed railway track. Sheng et al. [22] and François et al. [23] solved the dynamic interaction between a longitudinally uniform structure and soil under moving load by using the 2.5D finite element–boundary element coupling method. Lin et al. [24] studied the seismic response of underground structures using the 2.5-dimensional finite element infinite element method. Recently, Zhu et al. [25,26] adopted a 2.5D indirect boundary element method to study the dynamic response of structures under earthquake action.

In order to reduce the size of the finite field in the computational model and effectively improve the computational efficiency, the method of setting artificial boundaries is usually used to simulate the radiation damping effect. The purpose is to ensure that scattered waves generated in soil media pass through artificial boundaries from within the finite computational area without reflection. Lysmer [27] proposed a viscous boundary with dampers on the truncated boundary according to the wave equation. Deeks [28] assumed that the scattered wave was a 2D cylindrical wave, and put forward the definition of a 2D viscous spring artificial boundary; that is, the equivalent viscous damper–spring–concentrated mass system is arranged on the truncated boundary. Liu and Lu [29] assumed that the scattered wave was a 3D spherical wave, and proposed a 3D viscous spring artificial boundary. In addition, the Perfectly Matched Layer is an absorbing boundary model, which absorbs external traveling waves by setting a layer of boundary zone with an energy dissipation effect [30]. An effective seismic input method is the key to the dynamic response analysis of the tunnel foundation system. Joney and Chen [31] proposed to transform the incident ground motion into the equivalent load acting on the viscous artificial boundary to realize the wave input of the one-dimensional model. Zhang et al. [32] proposed a two-dimensional substructure input method for two-phase media.

The rest of this paper is as follows. In Section 2, the interaction between soil and long linear underground tunnel under the oblique incidence of a seismic wave is described. The 2.5D finite element substructure analysis method using a high-precision artificial boundary is given in Section 3. The accuracy of the 2.5D finite element substructure model is verified in Section 4 by comparing the site response with the analytical solution and the substructure model with the analytical solution. Section 5 carries out parameter analysis. According to different incident angles and different locations of soil–rock interface, the seismic response of an urban underground long linear tunnel is discussed. The conclusion is given in Section 6.

2. Problem Description

The soil–tunnel interaction model during an earthquake is shown in Figure 1. The longitudinal axis of the tunnel is parallel to the x-axis. The material parameters of the site and tunnel remain unchanged along the longitudinal direction. It is assumed that the soil–tunnel interface is fixed. Sliding and other conditions are not considered. The seismic wave is incident at any angle from the underlying half-space, where θ_h is the included angle between the incident direction of seismic wave and the z-axis. The conversion angle θ_v is the angle between the projection of the incident direction in the x-o-y plane and the x-axis. The initial stress of tunnel and soil are zero.

3. The 2.5D Finite Element Substructure Method

3.1. Review 2.5D Finite Element Method

In the substructure method, the finite domain only includes long lined tunnel structures. Through the 2.5D finite element method, the plane element with three degrees of freedom at each element node is obtained. The seismic response of the tunnel mentioned in the previous section can be expressed as follows:

$$\tilde{R}(x,y,z,t)=N\overline{R}(y,z)\exp (-i{k}_{x}x)\exp (-i\omega t)$$

(1)

where N is the shape function matrix in 2D finite element analysis; R is the dynamic response amplitude of the structure in the y-o-z plane; the superscript represents the frequency-wavenumber domain; i is an imaginary unit; k_x is the wave number in the x direction and ω is the circular frequency.

Assuming that the urban underground soil layer is a linear elastic solid medium, its motion equation under the global coordinate system (x, y, z) can be expressed as

$$L^{\mathrm{T}}\mathsf{\sigma }=\rho \ddot{u}$$

(2)

where σ is stress; ρ is the density; $\ddot{u}$ is the nodes acceleration; L is a 3D differential operator matrix and the superscript T of the matrix represents transposition.

$$L={\left[\begin{array}{cccccc}\frac{\partial }{\partial x}& 0& 0& \frac{\partial }{\partial y}& 0& \frac{\partial }{\partial z}\\ 0& \frac{\partial }{\partial y}& 0& \frac{\partial }{\partial x}& \frac{\partial }{\partial z}& 0\\ 0& 0& \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]}^{\mathrm{T}}$$

(3)

By replacing the partial derivative of x in the 3D differential operator with—ik_x, the equilibrium equation in the form of 2.5D finite element method is obtained.

$$\left(K-{\omega }^{2}M\right)\overline{u}=\overline{F}$$

(4)

where K and M are the stiffness and mass matrices derived according to the standard finite element process. F is the interaction between the structure and the soil.

So far, the 2.5D equilibrium equation of the finite domain of the tunnel structure has been obtained.

$$\left[\begin{array}{cc}{S}_{nn}& S_{nm}\\ S_{mn}& S_{mm}\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_n\\ {\overline{u}}_m\end{array}\right\}=\left\{\begin{array}{c}0\\ -{\overline{F}}_m\end{array}\right\}$$

(5)

where the subscript n refers to the number of internal nodes of the finite domain except the soil and tunnel interface. The subscript m refers to the number of nodes on the soil and tunnel interface.

The dynamic stiffness matrix is as follow

$$\mathbf{S}=K-{\omega }^{2}M$$

(6)

3.2. 2.5D High-Precision Artificial Boundary Conditions

The artificial boundary diagram of the 2.5D finite element model is shown in Figure 2a. In Figure 2b, a local coordinate system (r, s, t) is established along the outer surface of the tunnel. The red surface in Figure 2a and the red curve in Figure 2b represent artificial boundary conditions. The coordinate conversion relationship of partial derivatives between different coordinate systems is shown below.

$$\left\{\begin{array}{c}\frac{\partial }{\partial y}\\ \frac{\partial }{\partial z}\end{array}\right\}=\left[\begin{array}{cc}1& 0\\ -n_y/n_z& 1/n_z\end{array}\right]\left\{\begin{array}{c}\frac{\partial }{\partial s}\\ \frac{\partial }{\partial t}\end{array}\right\}$$

(7)

The coordinate conversion relationship of the partial derivative is introduced into the wave Equation (2) of solid medium, and the following equation is obtained:

$$L_{st}^{\mathrm{T}}D{L}_{st}\tilde{u}=-{\omega }^2\rho \tilde{u}$$

(8)

where the symbol ~ indicates in the frequency domain wavenumber domain. The differential operator in local coordinate system can be written as

$$L_{st}=-i{k}_x{L}_r+L_s\frac{\partial }{\partial s}+L_t\frac{\partial }{\partial t}$$

(9)

By introducing the differential operator (9) into Equation (8), it is obtained that

$$A_{ss}\frac{{\partial }^2\tilde{u}}{\partial s^2}+\left(A_{st}+A_{st}^T\right)\frac{{\partial }^2\tilde{u}}{\partial s\partial t}+A_{tt}\frac{{\partial }^2\tilde{u}}{\partial t^2}-i{k}_x\left[(A_{rs}+A_{rs}^T)\frac{\partial \tilde{u}}{\partial s}+(A_{rt}+A_{rt}^T)\frac{\partial \tilde{u}}{\partial t}\right]-k_x^2{A}_{rr}\tilde{u}+{\omega }^2\rho \tilde{u}=0$$

(10)

The coefficient matrices A_ss, A_st, A_tt, A_rs, A_rt and A_rr are as follows:

$$\begin{array}{l}{A}_{ss}=L_s^{T}D{L}_s A_{st}=L_s^{T}D{L}_t A_{tt}=L_t^{T}D{L}_t\\ A_{rs}=L_r^{T}D{L}_s A_{rt}=L_r^{T}D{L}_t A_{rr}=L_r^{T}D{L}_r\end{array}$$

(11)

By discretizing Equation (10) along the s direction, the displacement vector of any thin layer element can be expressed

$$\tilde{u}=N_e\overline{u}$$

(12)

The nodal shape function matrix N_e as

$$N_e=\left[\begin{array}{cccccc}{N}_1& 0& 0& N_2& 0& 0\\ 0& N_1& 0& 0& N_2& 0\\ 0& 0& N_1& 0& 0& N_2\end{array}\right]$$

(13)

where N₁ and N₂ are one-dimensional finite element shape functions.

$$-B_{ss}\frac{{\partial }^2\overline{u}}{\partial s^2}-\left(B_{st}-B_{st}^T\right)\frac{\partial \overline{u}}{\partial s}-(B_{rs}+B_{rs}^T)\frac{\partial \tilde{u}}{\partial s}-(B_{rt}-B_{rt}^T+B_{rr}-B_{tt}+{\omega }^{2}M)\overline{u}={\mathsf{\sigma }}_z$$

(14)

The coefficient matrices B_ss, B_st, B_tt, B_rs, B_rt, B_rr and mass matrix M are as follows:

$$\begin{array}{l}{B}_{ss}=n_z{\displaystyle \int N^T{A}_{ss}Ndt} B_{st}=n_z{\displaystyle \int N^T{A}_{st}\frac{\partial N}{\partial t}dt} B_{tt}=n_z{\displaystyle \int \frac{\partial N^T}{\partial t}{A}_{tt}\frac{\partial N}{\partial t}dt}\\ B_{rs}=-i{k}_x{n}_z{\displaystyle \int N^T{A}_{rs}Ndt} B_{rt}=-i{k}_x{n}_z{\displaystyle \int N^T{A}_{rt}\frac{\partial N}{\partial t}dt}\\ B_{rr}=-k_x^2{n}_z{\displaystyle \int N^T{A}_{rr}Ndt} M=\rho n_z{\displaystyle \int N^{T}Ndt}\end{array}$$

(15)

In the substructure method, it is necessary to establish a high-precision boundary condition to simulate the radiation damping in the infinite domain. The zigzag boundary is converted to the 2.5D form above. At the same time, the paraxial boundary condition is used in the underlying half-space to make up for the deficiency that the zigzag boundary can only be applied to a rigid base.

$${\tilde{\mathsf{\sigma }}}_h={\left\{\begin{array}{ccc}{\tau }_{xz}& {\tau }_{yz}& {\sigma }_z\end{array}\right\}}^T=-A_h\frac{{\partial }^2{\tilde{u}}_h}{\partial y^2}-B_h\frac{\partial {\tilde{u}}_h}{\partial y}+C_h{\tilde{u}}_h$$

(16)

The coefficient of 2.5D paraxial boundary condition is as follows:

$$-A\frac{{\partial }^2\overline{u}}{\partial y^2}-B\frac{\partial \overline{u}}{\partial y}+C\overline{u}=0$$

(17)

where the expressions of matrix A, B and C are as follows:

$$A=B_{ss}\oplus A_h$$

(18)

$$B=(B_{st}-B_{st}^T+B_{rs}+B_{rs}^T)\oplus B_h$$

(19)

$$C=(-B_{rt}+B_{rt}^T-B_{rr}+B_{tt}-{\omega }^{2}M)\oplus C_h$$

(20)

Based on the partial derivative relationship, a quadratic eigenvalue equation for the wavenumber k_y is obtained in the frequency domain wavenumber domain.

$$\left(k_y^{2}A+i{k}_{y}B+C\right)\overline{\overline{u}}=0$$

(21)

Equation (22) gives the force–displacement relationship at the right artificial boundary in the frequency-wavenumber domain. Similarly, we can also obtain the force–displacement relationship at the left artificial boundary, which will not be repeated here.

As shown in Figure 2c, through the combination of the left and right absorption boundary conditions, a zigzag-paraxial high-precision absorption boundary condition can be obtained in the frequency domain, which can be written as the following division matrix form:

$$\left\{\begin{array}{c}{\overline{F}}_{l,L}\\ {\overline{F}}_{c,L}+{\overline{F}}_{c,R}\\ {\overline{F}}_{r,R}\end{array}\right\}=\left[\begin{array}{ccc}{K}_{ll,L}^{\infty }& K_{lc,L}^{\infty }& 0\\ K_{cl,L}^{\infty }& K_{cc,L}^{\infty }+K_{cc,R}^{\infty }& K_{cr,R}^{\infty }\\ 0& K_{rc,R}^{\infty }& K_{rr,R}^{\infty }\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_l\\ {\overline{u}}_c\\ {\overline{u}}_r\end{array}\right\}$$

(22)

where the subscripts L and R denote the left and right artificial boundaries, respectively; the subscripts c, l and r represent the common nodes of the left and right boundaries, the left boundary except the common nodes and the right boundary except the common nodes, respectively.

So far, the 2.5D zigzag-paraxial boundary condition in the frequency domain is established. Nodes on the artificial boundary are sorted and merged. Equation (22) is rewritten as the following expression:

$$\left\{\begin{array}{c}{\overline{F}}_m^s\\ 0\end{array}\right\}=\left[\begin{array}{cc}{K}_{mm}^{\infty }& K_{mr}^{\infty }\\ K_{rm}^{\infty }& K_{rr}^{\infty }\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_m^s\\ {\overline{u}}_r^s\end{array}\right\}$$

(23)

where K∞ is the stiffness coefficient of the 2.5D zigzag-paraxial boundary. The r represents the remaining artificial boundary nodes other than those at the interface between soil and structure.

3.3. The Equivalent Load Method for 2.5D Earthquake Input

In the substructure analysis method, only the finite element model of the structure is established. The truncated infinite domain of the surrounding medium is simulated by the fundamental impedance function or dynamic stiffness in frequency domain. The impedance function or dynamic stiffness in frequency domain is integrated into the structural finite element model through the structure–soil medium interface condition. The earthquake input is realized at the soil–tunnel interface based on the free field response. In the substructure model, because the artificial boundary is close to the surface of the tunnel, the precision of the artificial boundary condition is required to be higher and the high-precision artificial boundary condition is usually selected. Zhang et al. [30] proposed a high precision artificial boundary condition in the frequency domain based on the paraxial boundary [33] and zigzag boundary [34]. It is used to analyze vector wave propagation in layered half-space. In this section, the artificial boundary condition is further developed for the 2.5D substructure model.

The equivalent nodal force at the free field is expressed as follows:

$${\overline{\mathbf{f}}}^f=L_i{\overline{\mathsf{\sigma }}}^f\mathbf{d}$$

(24)

3.4. System Equation of Substructure Method

The difficulty of the substructure method lies in the calculation of the frequency domain dynamic stiffness of the infinite medium surrounding the structure. To complete the establishment of the above model, the open soil–structure system needs to be decomposed into a closed finite area in the near field and an infinite area in the far field. Secondly, the closed finite region and the far-field infinite region are analyzed using different analysis methods, and the coupling is carried out by the interface conditions. Finally, the ground motion input is realized at the artificial boundary of the closed finite area based on the free field response.

Combining the high-precision artificial boundary condition Equation (23) and the finite element Equation (5) above, the finite element equation of the whole underground structure is obtained.

$$\left[\begin{array}{ccc}{S}_{nn}& S_{nm}& 0\\ S_{mn}& S_{mm}+K_{mm}^{\infty }& K_{mr}^{\infty }\\ 0& K_{rm}^{\infty }& K_{rr}^{\infty }\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_n^t\\ {\overline{u}}_m^t\\ {\overline{u}}_r^t\end{array}\right\}=\left\{\begin{array}{c}0\\ -{\overline{f}}_m^f+K_{mm}^{\infty }{\overline{u}}_m^f+K_{mr}^{\infty }{\overline{u}}_r^f\\ K_{rm}^{\infty }{\overline{u}}_m^f+K_{rr}^{\infty }{\overline{u}}_r^f\end{array}\right\}$$

(25)

4. Method Verification

4.1. Free Field Response of Elastic Half Space

According to the stiffness matrix method, the 1D site response is obtained. The half-space site response under the oblique incidence of seismic wave is calculated in the frequency domain. Through time delay, the 1D site response can be expanded to the 2.5D finite element response of each node.

To verify the accuracy of 2.5D finite element free field response, select a uniform half-space with density ρ of 2500 kg/m³, Young’s modulus E of 1.323 Gpa and Poisson’s ratio of 0.25. The reference point A taken is located in the middle of the site surface. The size of the site model is 800 m in length and 400 m in height.

As shown in Figure 3, the velocity–time history expression of the input bidirectional pulse is as follows

$$v(t)=A_{max}\left[1-2{\left(\pi f_p\left(t-t_0\right)\right)}^2\right]e^{-{\left(\pi f_p\left(t-t_0\right)\right)}^2}$$

(26)

where the maximum amplitude A_max = 0.001 m/s; the predominant frequency of the Fourier spectrum f_p = 2.0 Hz; time is expressed in t; the time of pulse displacement reaching peak value is t₀ and the calculation time step is 0.005 s.

According to the numerical simulation requirements, the mesh size of plane finite element needs to be discretized as Δy × Δz = 12.5 m × 12.5 m.

It can be seen from Figure 4 that when the incident pulse is P wave, the site response of 2.5D finite element at observation point A is consistent with the analytical solution under the oblique incidence of seismic wave. When the SV wave is inclined to incident, the 2.5D finite element site response and analytical solution are in good agreement at observation point A in Figure 5. The above numerical examples show that the accuracy of site response in the plane 2.5D finite element is verified when the seismic wave is obliquely incident.

4.2. Response of Lining Structure in Elastic Half-Space

For the soil–tunnel interaction problem, firstly, it is compared with the vertical incident seismic wave in the plane. As shown in Figure 6, the radius of concrete circular lining r = 5 m. The dimensionless parameter is η = ωr₁/πc_s = 0.132. The buried depth of lining H = 8.33 r₁. Lining thickness d = 0.1 r₁. The material parameters of concrete lining and uniform half-space site are listed in

Table 1. Site conditions.

	Density ρ (kg/m3)	The First Lamé Constant λ1 (GPa)	The Second Lamé Constant G1 (GPa)
Soil	2665	18	0.25
Lining	2240	22	0.2

. When the seismic wave is obliquely incident as a P-wave, the incident angle θ = 60°. When the seismic wave is an SV wave with oblique incidence, the incidence angle θ = 15°. For the problem of tunnel wave scattering in a 2D homogeneous bedrock half-space, we use the ABAQUS large model as a reference solution in the study, with a finite element model size of 300 m × 140 m. The dotted line in Figure 6b represents a 2D finite element viscous spring artificial boundary.

As shown in Figure 7, the radial displacement of the substructure model is compared with the normalized radial displacement of the structure under the soil–structure interaction of the 2D finite element model. From the results, it can be seen that the 2.5D finite element substructure method is highly consistent with the reference scheme under 2D oblique incidence, regardless of whether the P wave or SV wave are obliquely incident. The effectiveness of the proposed finite element substructure method has been verified under the condition of oblique incidence of seismic waves.

So far, the accuracy of the in-plane substructure method has been verified. In order to further verify the accuracy of the tunnel along the axis direction, the oblique incidence angle of seismic wave is taken as θ_v = 30°, as shown in Figure 8a. As in the plane, the radius of concrete circular lining r = 5 m, lining thickness d = 0.1 r₁, and lining buried depth H = 5 r₁. The material parameters of circular lining and single half-space site are listed in

Table 2. Site conditions.

	Density ρ (kg/m3)	Young’s Modulus E (GPa)	Poisson’s Ratio υ
Soil	2664	7.567 × 109	0.333
Lining	2240	16 × 109	0.2

. The dimensionless parameter is η = ωr₁/πc_s = 0.105.

The results are compared with the analytical solutions of Debarros and Luco [35]. By comparing the results, it can be found that the 2.5D substructure calculation method used in this paper is in good agreement with the normalized radial displacement amplitude of the analytical solution under the vertical incidence of seismic wave, whether P waves or SV waves are used, as shown in Figure 9. In the direction of the tunnel axis, Figure 10 shows that the normalized axial displacement amplitude of the substructure model is consistent with that of the analytical solution. The accuracy of the proposed substructure method in the longitudinal direction of the tunnel is verified.

5. P–SV Wave Scattering of Long Lined Tunnel

According to the verification in Section 4, the 2.5D finite substructure method proposed in this paper can accurately calculate the seismic response of the tunnel. In this section, the impact of different seismic wave incident angles and conversion angles, and soil–rock interface position on the seismic response of the tunnel lining is further analyzed.

5.1. Influence of Different Incidence Angle on Structural Seismic Response

Taking a circular section tunnel as an example, the tunnel section form and lining material parameters remain unchanged along the longitudinal direction. The calculation model is the same as Figure 6b. The buried depth H of the tunnel is 15 m. The thickness of the soil layer is 80 m. The outer diameter r₀ = 5.5 m, and the inner diameter r₁ = 5 m. Assuming that the tunnel lining maintains linear elastic deformation, the lining materials and circular tunnel parameters are shown in Table 1. The Kobe wave is selected for the incident seismic wave record, and the acceleration–time curve and spectrum curve are shown in Figure 11. For the incident with the P wave and the SV wave, the incident angle θ_v ranges from 0° to 90°, and the conversion angle θ_t ranges from 0° to 90°. The increment interval between the incident angle and the conversion angle is 15°. The external surface seismic responses of tunnel lining are calculated under different incident angles and conversion angles.

Figure 12 is a schematic diagram of the internal force of the lining; N, τ and M_c represent the axial force, shear force and bending moment of the lining in the cross section, T represents the longitudinal axial force of the lining. M_v and M_h represent the vertical and horizontal bending of the lining in the horizontal plane moment.

By changing the incident angle of the seismic wave, the law of the maximum bending moment values, M_c, M_v and M_h corresponding to the lining changing with the incident angle, are obtained. When the P wave is incident, as shown in Figure 13a, M_c increases with the increase of the incident angle and the conversion angle, and the maximum value is reached when both the incident angle and the conversion angle are 90°. The M_v and M_h results are shown in Figure 13c,e, the bending moment values are mainly affected by the transition angle. Both of them demonstrate peak value at vertical incidence and the transition angle of 90°. When the incident seismic wave is an SV wave, as shown in Figure 13b, M_c is also mainly affected by the conversion angle. The maximum lining bending moment first increases and then decreases with the conversion angle. At the same conversion angle, M_c increases with the increase of the incident angle. Figure 13d,f shows that the maximum value of M_v and M_h increase first and then decrease with the increase of the conversion angle under seismic waves’ vertical incidence. It can be seen that when the P wave is incident, considering the seismic response, the maximum bending moment of the lining is at the conversion angle of 90°. When the SV wave is incident, the unfavorable bending moment increases first and then decreases with the same incident angle. The different seismic waves’ incident angles also affect the internal force of the tunnel, which provides a way for us to further analyze the seismic response of underground structures.

5.2. Influence of Soil–Rock Interface on Structural Seismic Response

When the P–SV wave is inclined to incidence, the incidence angle is 30° and the conversion angle is 45°. The Kobe wave in Section 5.1 is selected for the seismic wave time history record. The tunnel inner diameter R = 5 m, and the lining thickness r = 0.1 R. The parameters of soil and rock materials are given in

Table 3. Site conditions.

	The First Lamé Constant λ (GPa)	The Second Lamé Constant G (GPa)	Density ρ (kg/m3)
Lining	4.44	6.67	2240
Soil	2.14	0.24	2665
Rock	2.592	1.728	2700

. The location of soil–rock interface is shown by the red lines in Figure 14a–c, which respectively represent that the tunnel is above the soil–rock interface, the tunnel is on the soil–rock interface and the tunnel is below the soil–rock interface. The longitudinal seismic response of the tunnel at different locations of the soil–rock interface is focused on.

The method proposed in this paper is applied to analyze the influence of different soil–rock interface positions on the tunnel seismic response. When the P–SV wave is incident, the maximum bending moment values on the inner surface of the lining are shown in Figure 15 and Figure 16, respectively. According to Figure 15a, when the tunnel is located at the soil–rock interface, the in-plane bending moment M_c reaches the maximum value. For the longitudinal bending moments M_v and M_h in Figure 15b,c, when the tunnel passes through the soil–rock interface, the longitudinal bending moments reach the maximum.

However, as shown in Figure 16a,b, when the SV wave is incident, the in-plane bending moment value M_c and the longitudinal bending moment value M_v are close to the maximum value when the soil layer thickness is 5 m and the tunnel passes through the soil–rock interface. When considering the longitudinal bending moment value M_h, as shown in Figure 16c, the bending moment reaches the maximum value when crossing the soil–rock interface. Compared with numerical examples of different soil–rock interfaces, the internal force value of tunnel crossing soil–rock interface is still close to or reaches the maximum when the SV wave is incident.

6. Conclusions

In this paper, a 2.5D substructure method is proposed for the analysis of half-space tunnel wave scattering. On the surface of the structure, the high-precision 2.5D zigzag boundary is used. In the dynamic analysis of soil–structure interaction, the finite domain can only be composed of tunnels, which saves the calculation cost. At the same time, the 2.5D finite element method can effectively avoid the longitudinal truncation effect. Actual seismic records are used in numerical examples. The influence of the incident angle when P–SV waves are incident is mainly discussed. The influence of the position of the soil–rock interface on the seismic response of the lining structure is introduced.

The seismic response of underground tunnels is unique; for example, tunnels are less affected by their own inertia force, so are mainly controlled by the free field deformation of surrounding media. Their seismic response analysis needs to consider the dynamic interaction between the tunnel and surrounding media. When the incidence angle of seismic wave is considered, the seismic response of the vertical incidence is not necessarily the maximum value, and the seismic response may be amplified in the case of oblique incidence. For different site conditions, the same underground tunnel will show different seismic response characteristics. The seismic response of the underground structure near the soil–rock interface of a long linear tunnel may be greater than that of the underground structure in a uniform site.

In addition, in the next work, this method can be used to consider the dynamic interaction between soil and structure of multiple structures. It provides an effective calculation method for studying the seismic response interaction between 3D structures.