A 2.5D Finite Element Method Combined with Zigzag-Paraxial Boundary for Long Tunnel under Obliquely Incident Seismic Wave

Abstract

Seismic waves propagation with an oblique angle to the tunnel axis will cause asynchronous tunnel motions and have a significant effect on the axial response. A high-precision 2.5D finite element method is established in the frequency domain to simulate the 3D seismic response of the tunnel. This method avoids the disturbance caused by the truncation of the tunnel in the longitudinal direction. Meanwhile, a 2.5D zigzag-paraxial boundary is derived to further improve the calculation efficiency of the 2.5D finite element model. Moreover, by combining the 2.5D finite element method, 2.5D zigzag boundary condition and seismic motion input methods, an obliquely incident substructure method for plane seismic waves is built by converting the plane seismic wave into equivalent nodal forces. The proposed 2.5D finite element method is verified by comparing with a reference solution. Finally, the 2.5D finite element method is applied to study the seismic response of the long lined tunnel. Parameter analyses illustrate that the wave propagation effect to the tunnel axis has a non-negligible influence on the axil deformation of long tunnels.

1. Introduction

In seismic engineering problems, the study of tunnel dynamic characteristics has become increasingly important. The safety of tunnels under earthquakes has attracted extensive attention [1,2,3,4,5]. Some observation results show that when seismic wave propagates along the longitudinal direction, it will cause tension, compression and bending deformation of the tunnel [6]. As an underground structure with a large span in the longitudinal direction, tunnel deformation in the longitudinal direction is more remarkable than that in the transverse direction, which should be highly considered in seismic design.

The common seismic analysis methods of underground tunnels include analytical solutions [7,8], model experiments [9] and numerical methods [10,11,12,13,14]. Compared with the assumption limitation of the analytical method and the long period of experimental testing, numerical methods can efficiently and accurately analyze complex soil–structure dynamic interaction. For a 3D simplified soil–tunnel model, soil around the tunnel is denoted by using soil springs to consider the soil–tunnel interaction [15]. Based on the concurrent discretization of the entire computation domain with coarse and fine-scale finite element meshes, Yu et al. [10] used a multi-scale method to simulate dynamic responses of a long tunnel. Li and Song [11] established a 3D soil–tunnel structure interaction model to simulate the longitudinal response of the tunnel when subjected to an asynchronous earthquake. Based on the time-domain finite element method with the viscous-spring artificial boundary condition, Huang et al. [13] conducted a study of the seismic responses of the long lined tunnels under the obliquely incident earthquake. Fabozzi et al. [12] addressed the problem of underground tunnels when subjected to asynchronous seismic shaking along the axis by a comparison with the synchronous case that is usually considered in the tunnel design. Both sides of the tunnel model are truncated, resulting in some numerical errors in the seismic response of the tunnel. Therefore, the 3D finite element model needs to take a long distance in the longitudinal direction to ensure that the tunnel seismic response results are accurate.

In recent years, some researchers have become increasingly interested in 2.5D modeling methods for tunnel–soil system dynamic interaction analysis, because in long lined tunnels, 2.5D finite element models are more efficient and accurate than 3D models. Hwang and Lysmer [16] first used a special type of wave transmission element to simplify the 3D system into a finite element plane, where each node of the plane element has three degrees of freedom, while retaining the propagation effects of waves. Yang and Hung [17] introduced the 2.5D method, which applies Fourier transform of the longitudinal coordinate to reduce the problem’s dimensionality and named it the 2.5D finite element/infinite element method. According to the different positions of the excitation source, the dynamic interaction analysis of the tunnel–soil system is mainly divided into two categories. One is the internal source problem, such as the vibration analysis of tunnel–soil system under moving train load [17,18,19,20,21,22,23,24,25]. Sheng et al. [19] coupled the finite element boundary element (FE–BE) method to analyze the vibration problem of the soil–tunnel system caused by moving trains load in the ground or tunnels. François et al. [20] studied the dynamic interaction between tunnels and layered soil under train load by using the general 2.5D FE-BE method. Jin et al. [22] used measurement data of dynamic response of tunnel to verify 2.5D coupled FE–BE model under the passage of a train. Lopes et al. [23] proposed a coupled 2.5D FEM–perfectly matched layer method to study the effect of soil stiffness on soil–tunnel interaction in the frequency domain. Recently, Liravi et al. [25] proposed the 2.5D FEM–singular boundary method (FEM–SBM) to deal with soil–tunnel interaction problems under the moving train load. The other category is the study under external source loads such as seismic motions [26,27,28,29]. Lin et al. [26] used the 2.5D finite element/infinite element method to analyze the seismic response of longitudinal invariant tunnels. Zhu et al. [27] used a 2.5D finite element indirect boundary element (FE–IBE) model to analyze the difference in seismic response of tunnels embedded in single-phase and two-phase media. Meanwhile, by using 2.5D FE–BE method, the dynamic response of segmented tunnels [28] and long lined tunnel in a water-saturated poroviscoelastic half-space [29] under earthquake were analyzed.

In the domain-based method, it is necessary to introduce artificial boundary conditions to consider the radiation damping effect of infinite domain [30,31,32,33,34,35,36]; this also applies in the 2.5D finite element method. In this case, the artificial boundary needs to be established by 2.5D forms. A 2.5D form of infinite element boundary is first developed by Yang and Hung [17], in which the choice of wave number is different from that of the previous infinite element method. The 2.5D boundary element condition was deduced by Sheng et al. [19], which was used to simulate the ground vibration caused by train moving loads. However, for some complex infinite domain problems, the fundamental solution is complicated, which limits the application of the boundary element method to a certain extent. Zhu et al. [27] established 2.5D indirect boundary element method, which is easier to simulate the layered site than boundary element method in solving tunnel seismic response. In addition to the boundary element method, other artificial boundaries require finite element spatial domain and mesh discretization. At the beginning of this century, a consistent absorbing boundary for horizontally layered soil was derived by Park and Tassoulas [37]. The thin layer method was popularized to obtain the zigzag boundary of multi-layer soil [38,39]. The zigzag boundary can transform the boundary into a nonvertical line. Since the zigzag boundary can adapt to complex geometric shapes, it is easy to attach to the structural surface. However, the zigzag boundary is limited to a rigid base. At the same time, in order to further improve the computational efficiency, Lin et al. [40] established a substructure technique to simulate the wave scattering problem. Combined with the site response analysis method and the seismic input method, Zhang et al. [41] has established a substructure analysis method for seismic response of the soil–structure system, which is suitable for multi-medium coupling conditions. Kausel and Roesset [42] derived an accurate dynamic stiffness matrix of a homogeneous half space, which constitutes an accurate boundary condition of an infinite underlying half space. However, since this matrix is not a polynomial of wave numbers, it cannot be coupled with the zigzag boundary. Later, Andrade [43] deduced the stiffness matrix of the homogeneous half space, which can be assembled in the overall stiffness matrix to deal with the wave problem in layered half space.

Compared with previous numerical analyses of 3D tunnel models, the 2.5D method has a greater advantage in terms of saving computational costs. In addition, the novelty of this article mainly includes the following two aspects. (1) The paper first derives the 2.5D zigzag boundary condition for 2.5D finite element method. Compared with other boundary conditions, the formula is simpler and can be attached to the structure surface, effectively reducing the calculation domain. (2) The paper proposes a 2.5D substructure method by combining the 2.5D finite element method, 2.5D zigzag boundary condition and seismic motion input method, where the 2.5D calculation model only includes structure domain, resulting in fast and efficient computation.

In this study, combined with 2.5D finite element method, zigzag-paraxial boundary condition and seismic wave input method, a 2.5D substructure method for analyzing wave scattering of tunnel in layered half space is proposed. The rest of this article is as follows. The 2.5D underground soil–structure interaction problem is proposed in Section 2. In Section 3, combining the 2.5D finite element method with zigzag-paraxial boundary condition, an efficient 2.5D substructure method is developed. Based on the 2.5D substructure model, it is verified that the 2.5D finite element method combined with zigzag-paraxial boundary is accurate and effective in Section 4. Section 5 discusses the variation law of dynamic response of tunnel under different incident angles of seismic waves and different buried depths. The conclusions are drawn in Section 6.

2. Wave Scattering Problem of Soil–Tunnel Interaction under Oblique Incidence of P-SV Wave

The soil–tunnel interaction system used for seismic analysis of underground long lined tunnel is shown in Figure 1a. The site is divided into the underlying half space bedrock and overlying soil layer. The interface between the soil and the tunnel is fully bonded. In addition, at the interface nodes, the displacement and stress are continuous. The geometry and material properties remain consistent along the longitudinal direction of the tunnel. The underlying half space bedrock, overlying soil layer and tunnel are assumed to be homogeneous, isotropic and linear elastic material. The surface traction and initial stress of the tunnel is zero.

In this study, the seismic wave is input from underlying half space with any angle in the form of a plane wave. The tunnel is subjected to earthquake loading by compression wave and shear wave. In the Cartesian coordinate system (x, y, z), as shown in Figure 1b, θ_v is defined as the incident angle between the x′ axis and the z axis. θ_t is the conversion angle between the x′ axis and the x axis in the x-o-y view in Figure 1c. The wave propagation direction is x′ axis.

3. 2.5D Finite Element Substructure Method

3.1. The Framework of the Proposed 2.5D Finite Element Substructure Method

The framework of the 2.5D finite element substructure method is presented in Figure 2. As shown in Figure 2b, the zigzag-paraxial boundary (denoted by red surface) attached to the surface of tunnel is employed to divide the 3D model into finite domain and infinite domains. The finite domain with complex geometry and material properties is discretized based on the 2.5D finite element method in Figure 2d. The combined zigzag-paraxial boundary condition, as illustrated in Figure 2e, is presented to provide the dynamic stiffness of the truncated infinite rock or soil media and simulate scattered field. As shown in Figure 2f, wave superposition of free field and scattered field on boundary are stated, in which free field displacement and traction on the boundary are obtained by the stiffness matrix method in reference [42]. Figure 2c shows that the final total 2.5D substructure model is established.

3.2. Description of the 2.5D Finite Element Method

Through the decomposition of Figure 2b, the finite domain contains the long tunnel, which can be modeled by the 2.5D finite element method. Different from 2D finite element, each 2.5D element node has three degrees of freedom. The tunnel displacement is expressed as follows:

$$u(x,y,z,t)={\displaystyle {\int }_{-\infty }^{\infty }N\overline{u}(y,z)\exp (-i{k}_{x}x)\exp (-i\omega t)\mathrm{d}\omega }$$

(1)

where N is the shape function matrix in the 2D isoparametric finite element analysis. The u-bar is the nodal displacement with three components in the frequency domain and wavenumber domain. i is imaginary. The circular frequency is ω. The wavenumber is given as k_x = cosθ_vsinθ_tω/c_b in the x direction, with c_b being the wave velocity of the elastic underlying half space.

The dry soil layer is assumed to be an elastic solid medium whose governing equation can be written in global coordinate system (x, y, z) as follows:

$$L^{\mathrm{T}}\sigma =-{\omega }^2\rho \overline{u}$$

(2)

where ρ is density; σ is stress matrix; L is a differential operator; the superscript T represents the transpose of the matrix; specifically, as follows:

$$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)

According to the Hooke’s law, the stress–strain relationship is

$$\sigma =D\mathsf{\epsilon }$$

(4)

where ε is strain. It is well known that D is a 3D linear elastic constitutive matrix with Lame constants λ and G.

$$D=\left[\begin{array}{cccccc}\lambda +2G& \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2G& \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2G& 0& 0& 0\\ 0& 0& 0& G& 0& 0\\ 0& 0& 0& 0& G& 0\\ 0& 0& 0& 0& 0& G\end{array}\right]$$

(5)

The strain–displacement relation is

$$\mathsf{\epsilon }=B\overline{u}$$

(6)

where B is the strain matrix.

By replacing partial derivative ∂/∂x with longitudinal wave number −ik_x, the calculation equation of tunnel finite element model in 2.5D form is obtained.

$$K\overline{u}-{\omega }^{2}M\overline{u}=\overline{F}$$

(7)

where $\overline{F}$ is the interaction force vector between the tunnel and the surrounding soil; K and M are the stiffness matrix and mass matrix of the 2.5D tunnel model, respectively

$$K={\displaystyle \iint B^{T}DBJd\xi d\eta }$$

(8)

$$M={\displaystyle \iint N^T\rho NJd\xi d\eta }$$

(9)

where J is the Jacobian matrix which represents dξdη area mapping to dxdy.

As shown in the Figure 2d, the 2.5D finite element equation of finite domain can be rewritten as follows:

$$\left[\begin{array}{cc}{S}_{pp}& S_{pq}\\ S_{qp}& S_{qq}\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_p\\ {\overline{u}}_q\end{array}\right\}=\left\{\begin{array}{c}0\\ -{\overline{F}}_q\end{array}\right\}$$

(10)

with dynamic stiffness matrix

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

(11)

where subscript q represents the nodes at the interface of the soil–tunnel system and subscript p represents all structural nodes except for the structural boundary.

3.3. A Zigzag-Paraxial Boundary Condition Suitable for 2.5D Finite Element Method

The infinite domain includes multilayer soil and underlying half space bedrock. A zigzag-paraxial boundary condition suitable for 2.5D finite element method is presented, which can be used for efficient analysis seismic responses of long lined tunnel.

The establishment process of the 2.5D zigzag-paraxial boundary condition is introduced, leading into an independent non-orthogonal local coordinate system (r, s, t) with zigzag boundaries attached to the surface of the tunnel. Due to the artificial boundary conditions on the left and right having the same building process, the detailed derivation process of the right boundary is introduced in this subsection.

The conversion relationship between the global coordinate system (x, y, z) and the local coordinate system (r, s, t) is shown in Figure 3. The local coordinate system is used in each soil sublayer element to represent the zigzag shape of the boundary, in which the direction of r-axis and s-axis are the same as x-axis and y-axis of global coordinate, respectively. The t-axis is aligned with the zigzag shape of the surface boundary. The coordinate conversion relationship is as follows:

$$\left\{\begin{array}{c}x\\ y\\ z\end{array}\right\}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& n_y\\ 0& 0& n_z\end{array}\right]\left\{\begin{array}{c}r\\ s\\ t\end{array}\right\}$$

(12)

where n_y and n_z are the y and z components of unit directional vector along t axis, respectively.

According to the Equation (12), the derivatives with respective to x, y and z can be obtained from the corresponding ones with respective to r, s and t as follows:

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

(13)

Substituting Equation (13) into Equation (2), the wave equation of multi-layer soil can be expressed as follows:

$$L_k^{\mathrm{T}}D{L}_k\tilde{u}=-{\omega }^2\rho \tilde{u}$$

(14)

where the ~ on u denotes its value in frequency domain.

Since x and r have the same direction, the partial derivative along the r direction ∂/∂r is also replaced by −ik_x. L_st is the differential operator in the local coordinate system. The L_st expression is as follows:

$$L_k={\left[\begin{array}{cccccc}-i{k}_x& 0& 0& \frac{\partial }{\partial s}& 0& -\frac{n_y}{n_z}\frac{\partial }{\partial s}+\frac{1}{n_z}\frac{\partial }{\partial t}\\ 0& \frac{\partial }{\partial s}& 0& -i{k}_x& -\frac{n_y}{n_z}\frac{\partial }{\partial s}+\frac{1}{n_z}\frac{\partial }{\partial t}& 0\\ 0& 0& -\frac{n_y}{n_z}\frac{\partial }{\partial s}+\frac{1}{n_z}\frac{\partial }{\partial t}& 0& \frac{\partial }{\partial s}& -i{k}_x\end{array}\right]}^T$$

(15)

To conveniently discretize the wave equilibrium equation along the s direction, L_k is rewritten as

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

(16)

with

$$L_r={\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]}^T$$

(17a)

$$L_s=\frac{1}{n_z}{\left[\begin{array}{cccccc}0& 0& 0& n_z& 0& -n_y\\ 0& n_z& 0& 0& -n_y& 0\\ 0& 0& -n_y& 0& n_z& 0\end{array}\right]}^T$$

(17b)

$$L_t=\frac{1}{n_z}{\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right]}^T$$

(17c)

Thus, Equation (14) can be rewritten as

$$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$$

(18)

with the coefficient matrices of the Equation (18) being

$$\{\begin{array}{c}{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}$$

(19)

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

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

(20)

with 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]$$

(21)

with

$$N_1=\frac{s_2-s}{s_2-s_1}$$

(22a)

$$N_2=\frac{s-s_1}{s_2-s_1}$$

(22b)

In the frequency domain, the discrete value of displacement $\overline{u}$ at nodes are prescribed as follows:

$$\overline{u}={\left\{\begin{array}{cc}{\overline{u}}_1^T& {\overline{u}}_2^T\end{array}\right\}}^T={\left\{\begin{array}{cccccc}{u}_{1x}& u_{1y}& u_{1z}& u_{2x}& u_{2y}& u_{2z}\end{array}\right\}}^{\mathrm{T}}$$

(23)

Equation (20) is taken into Equation (18). The weak form of the equivalent integral is obtained by using the definition of the stress on the surface.

$$-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}={\sigma }_z$$

(24)

with coefficient matrices of the equation

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

(25)

The stress vector on the right side of the Equation (24) can be expressed as

$${\sigma }_z={\left\{\begin{array}{cc}-n_z{L}_t^T{\sigma }_1& n_z{L}_t^T{\sigma }_2\end{array}\right\}}^T={\left\{\begin{array}{cccccc}-{\tau }_{xz}& -{\tau }_{yz}& -{\sigma }_z& {\tau }_{xz}& {\tau }_{yz}& {\sigma }_z\end{array}\right\}}^T$$

(26)

As shown in Figure 4, the zigzag boundary changes from a red surface to a curve with three degrees of freedom at each node. A zigzag boundary condition suitable for 2.5D finite element method is derived in the local coordinate system.

Due to the zigzag boundary is limited to rigid base, the 2.5D paraxial boundary condition in the underlying half space can improve the boundary accuracy. The paraxial boundary condition is a paraxial approximation of the exact dynamic stiffness matrix. When the seismic wave is incident obliquely, the relationship between the stress and the displacement at z = z_h (as shown in Figure 5) on the underlying half space surface is as follows

$${\tilde{\mathbf{\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$$

(27)

The matrices A_h, B_h, and C_h as the 3 × 3 matrices are given below

$$A_h=\frac{i{G}_h{c}_s}{2\omega }\left[\begin{array}{ccc}1& 0& 0\\ 0& \left(1-2\alpha \right)/{\alpha }^3& 0\\ 0& 0& \left(\alpha -2\right)/\alpha \end{array}\right]$$

(28a)

$$B_h=\left[\begin{array}{ccc}0& 0& -\frac{G_h{c}_s{k}_x}{2\omega }\left(\frac{{\alpha }^3-2{\alpha }^2}{{\alpha }^3}+1\right)\\ 0& 0& -\frac{G_h\left(1-2\alpha \right)}{\alpha }\\ \frac{G_h{c}_s{k}_x}{2\omega }\left(\frac{{\alpha }^3-2{\alpha }^2}{{\alpha }^3}+1\right)& \frac{G_h\left(1-2\alpha \right)}{\alpha }& 0\end{array}\right]$$

(28b)

$$C_h=\left[\begin{array}{ccc}i\omega \rho C_s+\frac{i{G}_h{c}_s{k}_x^2}{2\omega }(\frac{{\alpha }^3-2{\alpha }^2}{{\alpha }^3})& -\frac{i{G}_h{k}_x(1-2\alpha )}{\alpha }& 0\\ \frac{i{G}_h{k}_x(1-2\alpha )}{\alpha }& i\omega \rho c_p+\frac{i{G}_h{c}_s{k}_x^2}{2\omega }(\frac{1-2\alpha }{{\alpha }^3})& 0\\ 0& 0& i\omega \rho c_s+\frac{i{G}_h{c}_s{k}_x^2}{2\omega }\end{array}\right]$$

(28c)

where the matrices subscript h indicates the half space; the P and SV wave velocities are c_p and c_s, respectively; α is the ratio of c_s to c_p. G_h is the shear modulus of the half-space.

The thin layer element Equation (24) is assembled by displacement and stress continuity conditions in the soil layer. According to the relationship between the global coordinate system and the local coordinate system, ∂/∂s = ∂/∂y can be obtained. Finally, combined with the relationship between the surface stress and displacement of the underlying half space Equation (27), the assembled equation of Equations (24) and (27) can be obtained

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

(29)

with

$$A=B_{ss}\oplus A_h\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B=(B_{st}-B_{st}^T+B_{rs}+B_{rs}^T)\oplus A_h\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}C=(-B_{rt}+B_{rt}^T-B_{rr}+B_{tt}-{\omega }^{2}M)\oplus C_h$$

where the notation $\oplus$ is used to describe 2.5D finite element assembly in the matrix.

Taking the Fourier transform of Equation (29) into the wavenumber domain and transforming into solving the quadratic eigenvalue problem,

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

(30)

where k_y is the wavenumber along the y direction and $\stackrel{=}{\mathbf{u}}$ is nodal displacement vector in the wavenumber-frequency domain.

According to Zhang et al. [44], when the seismic wave propagates in the positive direction of the y-axis, the wave number satisfies the condition that Re(k_y) < 0 and Im(k_y) > 0. Select the corresponding k_y on the right boundary to superimpose the displacement mode as

$$\overline{\mathbf{u}}=\mathsf{\Phi }{E}_y\mathsf{\Gamma }$$

(31)

where E_y is a diagonal matrix and the diagonal elements are exp(−ik_yy); Φ and Г are the eigenvector corresponding to the selected eigenvalue and the undetermined coefficient, respectively.

Taking the artificial boundary on the right as an example, the surface stress of the thin layer element is discretized. The element node stress vector can be expressed as follows

$$\overline{\mathbf{f}}={\displaystyle \int N^T\tilde{T}ds}=-B_{ss}\frac{\partial \overline{\mathbf{u}}}{\partial s}-(B_{rs}^T+B_{st})\overline{\mathbf{u}}$$

(32)

with the surface stress $\tilde{T}$ of artificial boundary given as

$$\tilde{T}=-n_z{L}_s^T\sigma =i{k}_x{n}_z{A}_{rs}^T\tilde{u}-n_z{A}_{ss}\frac{\partial \tilde{u}}{\partial s}-n_z{A}_{st}\frac{\partial \tilde{u}}{\partial t}$$

(33)

The force of the bedrock half space on the bottom of the overlying soil layer is represented by ${\overline{\mathbf{f}}}_h$

$${\overline{\mathbf{f}}}_h=-A_h\frac{\partial {\overline{\mathbf{u}}}_h}{\partial y}-\frac{B_h}{2}{\overline{\mathbf{u}}}_h$$

(34)

Assembling the above $\overline{\mathbf{f}}$ and ${\overline{\mathbf{f}}}_h$, the right artificial boundary node force can be obtained

$$\overline{\mathbf{F}}=-A\frac{\partial \overline{\mathbf{u}}}{\partial y}-\mathbf{H}\overline{\mathbf{u}}$$

(35)

with

$$H=B_{rs}^T+B_{st}+\frac{B_B}{2}$$

(36)

Taking Equation (31) into Equation (35), the expression of the right boundary force-displacement relationship in the frequency-wavenumber domain is as follows

$$\overline{\mathbf{F}}=-A\mathsf{\Phi }\frac{\partial E_y}{\partial y}\mathsf{\Gamma }-H\mathsf{\Phi }{E}_y\mathsf{\Gamma }=\left[iA\mathsf{\Phi }K{\mathsf{\Phi }}^{-1}-H\right]\mathsf{\Phi }{E}_y\Gamma =K^{\infty }\overline{\mathbf{u}}$$

(37)

with

$$K^{\infty }=iA\mathsf{\Phi }K{\mathsf{\Phi }}^{-1}-H$$

(38)

where K is a diagonal matrix consisting of diagonal elements k_y and $K^{\infty }$ is the stiffness coefficient of the 2.5D zigzag-paraxial boundary.

As shown in Figure 2c, a 2.5D zigzag-paraxial boundary condition be obtained, which can be written as the following

$$\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\}$$

(39)

where the subscripts L and R denote the left and right artificial boundaries, respectively; the subscripts c represent the common nodes of the left and right boundaries. The l and r only refer to the left and right boundaries of the structure.

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

(40)

where the r represents the remaining artificial boundary nodes except the interface between soil and structure.

3.4. Seismic Wave Input

In this subsection, the 2.5D seismic wave input method is used to convert the incident seismic wave into the equivalent node load on the artificial boundary surface.

When solving the free field response [41], we use the stiffness matrix method [42] to calculate the free field response of the 1D site. The 3D site seismic response can be obtained through time delay. In the local coordinate system (x′, y′, z) shown in Figure 6, the displacement vector of the i-node on the artificial boundary of the soil structure interaction model can be written as

$${\overline{u}}_i^{f(1)}({x^{\prime }}_i,y_i,{z^{\prime }}_i,t)=\{\begin{array}{cc}0& 0\le t\le t_i\\ {\overline{u}}_j(0,y_i,0,t-t_i)& t>t_i\end{array}$$

(41)

with

$$t_i=\frac{\left|{x^{\prime }}_i\cos \phi +{z^{\prime }}_i\sin \phi \right|}{c_{x^{\prime }}}$$

(42)

where the superscript (1) indicates the local coordinate system (x′, y′, z); u_i and u_j are the displacement of the node i and j on the 1D site response analysis; t_i is the time of along the z direction from j to i.

The coordinate transformational relationship of site response on the truncated boundary is as follows

$${\overline{u}}^f=Q{\overline{u}}^{f(1)}$$

(43)

$${\overline{\sigma }}^f=Q{\overline{\sigma }}^{f(1)}{Q}^T$$

(44)

where ${\overline{u}}^f$ and ${\overline{\sigma }}^f$ are the displacement vector and stress vector in the global coordinates (x, y, z), respectively; the free field response is denoted by the superscript f, which is determined according to the site response. Q is the coordinate transform matrix, which can be written as

$$Q=\left[\begin{array}{ccc}\cos \alpha & -\sin \alpha & 0\\ \sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right]$$

(45)

The node force of the element can be obtained from stress by finite element discretization along the boundary.

$${\overline{F}}^f=L_i{\overline{\sigma }}^{f}d$$

(46)

where d is the element outer normal vector of the artificial boundary surface; L is the element length on the 2.5D artificial boundary surface.

The scattered field response can be simulated by the proposed 2.5D zigzag-paraxial boundary conditions.

3.5. System Equation

The total force and displacement can be decomposed into the following expressions, respectively

$$\overline{F}={\overline{F}}^f+{\overline{F}}^s$$

(47)

$$\overline{u}={\overline{u}}^f+{\overline{u}}^s$$

(48)

where the scattered field, denoted by the superscript s, is unknown.

Equations (47) and (48) are brought into Equation (40). The unknown scattering field response at the artificial boundary can be represented by the known free field response and the total response.

$$\left[\begin{array}{cc}{K}_{qq}^{\infty }& K_{qr}^{\infty }\\ K_{rq}^{\infty }& K_{rr}^{\infty }\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_q\\ {\overline{u}}_r\end{array}\right\}=\left\{\begin{array}{c}-{\overline{f}}_q^f+K_{qq}^{\infty }{\overline{u}}_q^f+K_{qr}^{\infty }{\overline{u}}_r^f\\ K_{rq}^{\infty }{\overline{u}}_q^f+K_{rr}^{\infty }{\overline{u}}_r^f\end{array}\right\}$$

(49)

Considering the continuous conditions of displacement and stress between layers, the 2.5D finite element Equation (10) and the 2.5D zigzag-paraxial boundary condition Equation (49) are assembled, and the equation of the 2.5D substructure method is established as follows

$$\left[\begin{array}{ccc}{S}_{pp}& S_{pq}& 0\\ S_{qp}& S_{qq}+K_{qq}^{\infty }& K_{qr}^{\infty }\\ 0& K_{rq}^{\infty }& K_{rr}^{\infty }\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_p\\ {\overline{u}}_q\\ {\overline{u}}_r\end{array}\right\}=\left\{\begin{array}{c}0\\ -{\overline{f}}_q^f+K_{qq}^{\infty }{\overline{u}}_q^f+K_{qr}^{\infty }{\overline{u}}_r^f\\ K_{rq}^{\infty }{\overline{u}}_q^f+K_{rr}^{\infty }{\overline{u}}_r^f\end{array}\right\}$$

(50)

By inverse Fourier transform of the results, we will obtain the displacement time history of the total response of the tunnel.

4. Method Verification

The substructure models are established by using the 2.5D method combined with zigzag-paraxial boundary condition, which verifies the accuracy of the proposed method under P-SV wave incidence.

4.1. Oblique Incidence of Seismic Waves in the Plane

For the problem of P-SV wave scattering with oblique incidence at any angle in the plane, the geometric characteristics are shown in Figure 7a. The lining depth of the tunnel is H = 8.33r₁, and the thickness is 0.1r₁. The numerical example takes r₁ = 5.0m. The material parameters are shown in

Table 1. Material constants of lining and rock.

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

. The 2.5D finite element substructure model is compared with the 2D finite domain extension model. The size of the 2D finite domain is 300 m × 140 m.

When the P wave is incident, the angle with the vertical direction is θ_v = 60°. For the incident SV wave, the angle with the vertical direction is θ_v = 15°. Due to the verification of the 2D oblique incidence problem, the conversion angle is θ_t = 90°. The dashed line in Figure 7a represents the 2D viscous spring artificial boundary. The radial displacement of the substructure model is compared with that of the 2D finite element model, as shown in Figure 8. It can be seen from the results that the 2.5D finite element substructure method has a high consistency with the reference solution under 2D oblique incidence.

To improve the comprehensiveness of the validation results, the data are further processed and analyzed [45]. The validation results are converted into a scatter plot based on relevant features (maximum normalized displacement amplitude for seismic simulations), as shown in Figure 9. The scatter plot uses the 3D model reference solution as the x-axis, and the 2.5D model solution is used as the y-axis.

Based on validation example 4.1, the R², maximum error, average error and standard deviation of the correlation coefficient k_c are calculated. The results are listed in

Table 2. R², maximum error, mean error and standard deviation of the ratio between maximum seismic response calculated by 2.5D finite element method and 3D reference solution U_2.5D/U_3D.

Verification	U2.5D/U3D
	R2	Maximum Error	Mean Error	Standard Deviation
P wave Ur − 1	0.9998	0.0240	0.0057	0.0075
SV wave Ur − 1	0.9999	0.0317	0.0049	0.0064

. It can be seen from Table 2 that all R² values are close to 1, indicating that the 2.5D substructure method proposed in the paper is consistent with the 3D reference solution.

4.2. Incidence along the Longitudinal Direction

For the accuracy of the longitudinal response of the method proposed, according to the results proposed by Debarros and Luco [46], the P-SV waves incident obliquely along the tunnel axis were analyzed. The substructure model is shown in Figure 10b; the oblique incidence angle θ_v = 30° and the conversion angle θ_t = 0° of the seismic wave are input. The corresponding dimensionless frequency is η = ωr₁/πc_s. Tunnel are buried with a depth of H = 5.0r₁ = 4.545r₀ and a lining thickness of d = 0.1r₁ = 0.0909r₀. The material parameters are shown in

Table 3. Material constants of lining and rock.

Soil	Elastic Modulus E1 (GN/m3)	Poisson’s Ratio ν1	Density ρ1 (kg/m3)
Lining	16	0.2	2240
Rock	7.567	0.333	2664

.

By comparing the normalized radial displacement u_r(r₀) in Figure 11 and the normalized longitudinal displacement u_l(r₀) in Figure 12, the longitudinal seismic response of the 2.5D finite element substructure calculation method is highly in agreement with the reference solution.

Figure 13 and Figure 14 show the comparison between the 2.5D calculated value and the 3D reference solution for each node when seismic waves are obliquely incident along the longitudinal direction. The closer the data points are to the 1:1 line, the more accurate the proposed method out of the plane.

It can be seen from

Table 4. R², maximum error, mean error and standard deviation of the ratio between maximum seismic response calculated by 2.5D finite element method and 3D reference solution U_2.5D/U_3D.

Verification	U2.5D/U3D
	R2	Maximum Error	Mean Error	Standard Deviation
P wave Ur – 2	0.9990	0.1255	0.0610	0.0386
SV wave Ur – 2	0.9986	0.1011	0.0347	0.0359
P wave Ul	0.9954	0.1406	0.0440	0.0360
SV wave Ul	0.9993	0.0730	0.0244	0.0197

that all R² values are close to 1, indicating that the 2.5D substructure method proposed in the paper is consistent with the 3D reference solution. On the whole, the statistical indicators used for evaluation show that the 2.5D model proposed in this article has high accuracy.

4.3. Accuracy of Longitudinal Responses

According to the 3D finite element calculation model established by Huang et al. [13], the length, width and height are 220 m, 60 m and 80 m, respectively. The length of the tunnel is 200 m. In the middle section of the numerical model, the response of the tunnel is accurate and is not affected by the truncation of both sides. The calculated material parameters are the same as Table 3. The 2.5D substructure method proposed in this paper is verified with a numerical example of the 3D finite element method. The substructure model and the incident angle of seismic wave are the same as those in Section 4.2.

Figure 15 shows the peak value of radial displacement u_r along the longitudinal direction of 3D tunnel finite element model and 2.5D tunnel finite element model respectively. In fact, when seismic motion is incident on the underlying half space at any incident angle, according to Snell’s law of wave propagation, the seismic peak response along the longitudinal direction of the tunnel is the same. There is only a delay of occurred time of dynamic response, so the peak value should be constant along the longitudinal direction of the tunnel, as shown in the 2.5D finite element method results presented in Figure 15. In the middle section of the tunnel, the seismic response results of the two methods are basically consistent. The two sides of the 3D long tunnel model are truncated, resulting in the error of the numerical calculation results (as shown by the red dashed box in Figure 15). It is found that the 2.5D finite element model can avoid the disturbance caused by the truncation of the tunnel in the longitudinal direction. At the same time, compared with the calculation time of 3D finite element model, the time cost of using 2.5D substructure model is only 1/6, which reflects the efficiency of 2.5D substructure calculation method.

5. P-SV Wave Scattering of Long Lined Tunnel

5.1. Influence of Seismic Wave Incidence Angle and Conversion Angle on Seismic Response of Tunnels

The computational model is shown in Figure 16. The tunnel and site parameters are shown in Table 1. Kobe wave is selected for the incident seismic wave record. The acceleration time history curve and spectrum curve are shown in Figure 17. For the incident with the P wave and the SV wave, respectively, 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.

At dimensionless frequency η = 0.132, Figure 18 shows the variation of radial displacement when the P wave is incident at different incident angles and conversion angles. With the gradual increase of the incident angle, the position of the maximum radial displacement gradually shifted from the vault to the haunch. Except for vertical incidence, for the same incidence angle, the radial displacement of the tunnel lining increases gradually with the increase of the conversion angle. The radial displacement values of the SV wave at different incident angles are given in Figure 19. At this time, with the increase of the incident angle, the position of the maximum radial displacement is gradually shifted from the haunch to the arch bottom. For the axial displacement, Figure 20 shows that when the P wave is incident and the conversion angle is fixed, the axial displacement increases gradually with the increase of the incident angle. In the case of the same incident angle, the larger the conversion angle, the smaller the radial displacement. When the SV wave is incident and the conversion angle remains the same, as the angle of incidence increases, the radial displacement decreases gradually in Figure 21. Therefore, when the incident angle is large, the longitudinal seismic response of the tunnel should be given enough attention. With the change of the transfer angle, the position of the maximum seismic response of the tunnel in the cross section is changed.

5.2. Influence of Buried Depths of Tunnels

The parameters and dimensions of tunnel material, thickness of soil layer and seismic waves selected for ground motion input are the same as those in Section 5.1, and the parameters of site and bedrock are shown in Table 3. The depth of the tunnel center is 10~40 m. According to the shear wave velocity of the medium and the dominant frequency of the seismic wave, the mesh size of the finite element is determined. The seismic wave is obliquely incident on the bedrock surface. According to the discussion in Section 5.1, the fixed conversion angle is 45° and the oblique incidence inclination angle is changed. The seismic response law of the tunnel lining at different burial depths and the influence of the change of the oblique incident angle on the seismic response of the lining at the same burial depth were investigated. To ensure the generalization, normalized depth is the ratio of tunnel burial depth D to inner diameter R, which is used to describe the influence of the embedded depth.

Figure 22 and Figure 23 show the maximum values of normal strain and shear strain of the tunnel lining under the incidence of P wave and SV wave, respectively. As can be seen from Figure 22a–c, when the P wave is incident, as the distance from the center point of the tunnel to the surface gradually increases, the normal strains in the three directions gradually decrease in most calculated conditions; this illustrates that the seismic response of the shallow buried structure cannot be ignored. The shear strain increases with the burial depth at a small incident angle (less than 45°). When the incident angle is larger, the shear strain decreases with the increase of burial depth. SV wave results are given in Figure 23.

As the burial depth increases, the internal force of the tunnel will gradually decrease. Figure 24 shows the time history curves of the maximum internal force calculated from the burial depth of 10 m when the P wave and SV wave are incident, respectively.

6. Conclusions

In this study, combined with 2.5D finite element method, zigzag-paraxial boundary condition and seismic wave input method, a 2.5D substructure method for analyzing wave scattering of tunnel in half space is proposed. Compared with the 3D time history analysis, the computational efficiency is greatly improved by using the 2.5D finite element substructure method. It provides a high efficiency analysis method for calculating large-scale seismic response of long linear underground structures.

The high precision artificial boundary condition can be placed on the surface of the structure, which greatly saves the computational cost. The dynamic stiffness of the surrounding infinite domain is provided. Under the incident P-SV wave, the proposed method is compared with the numerical model and the analytical solution. The correctness of the uniform half space wave scattering problem is verified by the calculation results. Further comparison shows that the 2.5D substructure method can effectively avoid the disturbance caused by the truncation of the tunnel in the longitudinal direction. In order to verify the capability of the proposed method, the actual seismic records are used in the numerical examples. When the P-SV wave is incident, the influence law of the incident angle of the seismic wave is mainly discussed. The effect of burial depth on the seismic response of lining structures is presented.

At the same time, the proposed method can also adapt to the complex structural shape of the cross section and the multi-structure soil–structure interaction problem. Due to the high precision of this method, discussing the changes of different model factors can provide accurate analysis for engineering design. The dynamic stiffness matrix of the underlying half space is obtained by the second-order paraxial approximation. In the next work, researchers can provide a new way to further improve the accuracy of the paraxial boundary.