#### Abstract

A three-dimensional (3D) detailed numerical model of an immersed tunnel in a horizontally layered site is established in this study. The 3D seismic response of the immersed tunnel in a horizontally layered site subjected to obliquely incident waves is analyzed based on the precise dynamic stiffness matrix of the soil layer and half-space via combined viscous-spring boundary and equivalent node stress methods. The nonlinear effects of external and internal site conditions on the whole model were determined by equivalent linearization algorithm and Mohr–Coulomb model, respectively. The proposed model was then applied to investigate the nonlinear seismic response of an immersed tunnel in the Haihe River subjected to seismic waves of oblique incidence. The dislocation (opening) of pipe joints in the immersed tunnel were analyzed to determine the response characteristics of the shear keys and overall displacement of the tunnel; the dynamic responses of the immersed tunnel subjected to obliquely incident seismic waves markedly differ from those of vertically incident seismic SV waves. The maximum stress value of shear keys and the maximum dislocation of the pipe joint appear as upon critical angle. The overall displacement of the tunnel increases as incident angle increases. Under severe earthquake conditions, both the pipe corners and midspan section of the roof and floor are likely to produce crack. These areas need careful consideration in the seismic design of immersed tunnel structures.

#### 1. Introduction

Immersed tunnel structures have grown increasingly common in cross-sea and cross-river traffic engineering projects. Many of these structures are located in highly seismic regions. The unique characteristics of the immersed tunnel make any structural damage caused by an earthquake particularly difficult to be repaired, so antiseismic design is particularly important. The complexity of the immersed tunnel system makes it difficult to test by shaking table. Numerical analysis methods are an important means to study the seismic performance of immersed tunnels. Immersed tunnels tend to be lengthy and are usually built on weak ground [1, 2], so it is necessary to consider the spatial variation effects of ground motion and the dynamic interaction between the soil and structure in such numerical analyses [3]. Numerical simulation of the whole dynamic analysis model must include both the immersed tunnel and surrounding rock [4].

Recent researchers have conducted many shaking table tests [5–7] and fine numerical simulations [8–10] for investigating the 3D seismic response of immersed tunnels. Hu et al. [11] investigated the mechanical properties of segmental joints in the immersed tunnel of the Hongkong-Zhuhai-Macao Bridge in China; Yu et al. [6] considered the multipoint input of ground motions via experiment. Chen et al. [5] carried out a large immersed tunnel-joint-site soil shaking-table test to determine the seismic response of immersed tunnels and joints.

It is generally assumed that seismic waves in a deep-focus earthquake strike perpendicularly to the bedrock and transfer it to the structure. For shallow-focus earthquakes with small epicentral distance, the influence of seismic wave incident angle on the seismic response of surrounding structures is an important consideration [12]. Du [13] analyzed the oblique incidence of seismic waves by using the transmission viscoelastic boundary conditions. FCPD Barros and Luco [14] obtained the three-dimensional harmonic response in the vicinity of an infinitely long, cylindrical cavity of circular cross section buried in a layered, viscoelastic half-space. Xiao et al. [15] investigated the shear behavior of an immersion joint subjected to compressive shear loads and found that the capacity of the joint is smaller than the sum of the capacities of all shear keys.

Previous researchers have mainly focused on 2D cross sections. However, it is difficult to accurately describe the nonlinear response state by simply decomposing the immersed tunnel into in-plane and out-of-plane components under the oblique incidence of seismic waves. Unfortunately, it is very challenging to secure an accurate 3D model due to its complexity. Huang et al. [16] proposed a method to realize the SV and SH waves with arbitrary incident angles in ABAQUS and indicated that the nonlinear seismic responses of long lined tunnels are highly affected by the incident angles of S waves. The literature mainly centers around the half-space uniform field rather than the 3D horizontally layered nonlinear seismic response of the site. Actual construction sites are usually layered; the amplification effect of the layered site is significantly different from that of the equivalent uniform site due to the resonance amplification and filtering effect of the soil layer [17], and this is proved in our paper [18].

In the present study, we established a 3D refined model of horizontal layered site immersed tunnel. We used the dynamic stiffness matrix of the soil layer and half-space combined with viscous-spring boundary and equivalent node stress methods to conduct accurate simulations of 3D seismic responses in a horizontally layered immersed tunnel structure. We applied the equivalent linearization algorithm (self-programmed) to determine the nonlinear effects of the external site of the overall model. We also used the Mohr–Coulomb model to determine the nonlinear effects of soil inside the overall model. The nonlinear effects of concrete were determined by a plasticity damage model of concrete in ABAQUS. The nonlinear seismic response of an immersed tunnel in the Haihe River and the seismic behavior of shear keys and pipe joints under obliquely incident SV waves were then determined based on the proposed model.

#### 2. Obliquely Incident Seismic Waves in ABAQUS

##### 2.1. Viscous-Spring Boundary

We used the 3D viscous-spring boundary set provided by Huang et al. [19] to simulate the stiffness and radiation damping effects of semi-infinite strata in this paper. The viscous-spring boundary reflects the propagation of an advancing wave by superimposing the plane wave and the far-field scattered wave. The characteristics of wave attenuation and multiangle incidence are reflected by two empirical parameters; the wave field of the infinite field at the artificial boundary is transformed to stress on artificial boundaries. When the finite-element integral is stable, the viscous-spring boundary remains stable. The spring stiffness and damping coefficient of the 3D viscous-spring artificial boundary are shown in equations (1) and (2):where *K*_{N} and *C*_{N} are the normal spring stiffness and damping coefficient of the nodes on the artificial boundary, respectively; *K*_{T} and *C*_{T} are the tangential spring stiffness and damping coefficient of the nodes on the artificial boundary, respectively; *A* and *B* are empirical parameters typically set to 0.8 and 1.1, respectively; *G* is the shear modulus; is the first Lame constant; *r* is the distance from the scattered wave source to the artificial boundary, which can be simply taken as the vertical distance from the geometric center of the structure to the artificial boundary; *ρ* is the soil density of the foundation; and and are S wave velocity and P wave velocity, respectively.

##### 2.2. Earthquake Input Method Based on Viscous-Spring Boundary

Here, we combine the finite-element method with the time-domain method of viscous-spring boundaries. The displacement and stress of the artificial boundary caused by equivalent load input on the artificial boundary must be consistent with that of the original free field. In this way, the earthquake input is translated into a free-field motion problem on artificial boundary nodes, whereas free-field motion can be transformed into an equivalent nodal load on artificial boundary nodes [20]. The equivalent nodal load of artificial boundary node *l* in direction *i* is calculated as follows:where is the influence area of node *l* on the artificial boundary surface; is the free-field stress tensor of node *l* in direction *i*; is the free-field displacement tensor of node *l* in direction *i*; is the free-field velocity vector of node *l* in direction *i*; is the viscous-spring boundary spring stiffness of node *l* in direction *i*; and is the viscous-spring boundary damping coefficient of node *l* in direction *i*.

##### 2.3. Plane Wave Input under Oblique Incidence in a 3D Layered Site

The schematic diagram of the input wave field when the SV wave is obliquely incident is shown in Figure 1 [20]. The plane SV wave with the incident angle is reflected on the surface of the bedrock, forming a plane SV wave with the reflection angle and a plane P wave with the reflection angle , where and the critical value of the incident angle of the SV wave . When the problem is solved, the input schedule is first transformed to the frequency domain by Fourier transform, and then all the effective frequency points are calculated in the frequency domain. Finally, the time-domain results are obtained by inverse Fourier transformation. The direct stiffness method can be applied to the free-field solution in the frequency domain.

Wolf [21] gives an accurate dynamic stiffness matrix for layered soil and half-space. The dynamic equilibrium equations of *i* layer and half-space in the plane motion (P-SV wave) are calculated by equations (4) and (5), respectively:where *P*_{i}, *R*_{i}, *P*_{b}, and *R*_{b} are horizontal and vertical load amplitudes of the upper and lower soil layers, respectively; *P*_{o} and *R*_{o} are horizontal and vertical load amplitudes of the free surface in half-space; *u*_{i}, , *u*_{b}, and are horizontal and vertical displacements of the upper and lower surface of the soil layer; and *u*_{o} and are horizontal and vertical displacements of the bedrock surface. The stiffness matrix is a symmetric matrix. The expressions of *k*_{11}–*k*_{44} and *r*_{11}–*r*_{22} are available in our references [21].

The stiffness matrix of each soil layer can be set and adjusted to obtain the overall space stiffness matrix *S*_{P-SV}, which is a (2*n* + 2) order square matrix where *n* is the number of soil layers. The overall dynamic equilibrium equation of the site can be written as follows:where *P* is the load vector, which is only related to the field characteristics, the incident wave amplitude, and the incident angle; the displacement vector *u* at the layer stratification can be solved accordingly.

The horizontal displacement and vertical displacement at any point within each soil layer can be expressed as equations (7a) and (7b). The velocity at any point can be derived from the displacement to the time as expressed by equations (8a) and (8b). Normal stress amplitude and as well as shear stress amplitude , can be derived from the following constitutive equation:where , , , and represent the incident and reflected amplitude of P and SV waves, respectively; is the wave number; is the phase wave velocity; is the circular frequency; and . For the incident SV wave, the direction cosine ; is the P wave propagation direction cosine.

In this study, we investigated the nonlinear seismic response of the site by the equivalent linearization method. The soil characteristics in the matrix are replaced by the equivalent shear modulus and damping ratio which match the equivalent shear strain:where is the maximum shear modulus; is a reference shear strain during the earthquake (the value can be determined experimentally or empirically); is the dynamic shear strain; and is the maximum damping ratio (which also can be obtained empirically or experimentally).

The coordinate system in the direct stiffness method differs from that in ABAQUS. After coordinate transformation, the solution is plugged into equation (3) to obtain the equivalent nodal load for use in ABAQUS.

#### 3. Input Verification

We established a 3D horizontally layered finite-element model (Figure 2(a)) and selected a 240 m × 240 m × 120 m cubic computational domain; 5 m × 5 m × 5 m hexagonal solid elements satisfying the precision requirements of the wave finite-element simulation (equation (11)) were used for discretization, in which A and B are monitoring points. We divided the soil into three layers each 40 m in thickness with the parameters shown in Table 1. The incident plane wave is a pulse wave (equation (12)). The waveform is shown in Figure 3. The SV wave incident angle is . The angle between the incident wave projecting to the *Z*-*X* plane and the *X* axis is , and the angle between the incident wave and the *Y*-axis is (i.e., incident angle in Section 2.3). We used the incident angle of (60°, 30°) for input verification.where is the largest size for the grid, is the minimum shear wave velocity, and is the maximum frequency of fluctuations.where *f* = 2, 0 ≤ *t* ≤ 0.5.

**(a)**

**(b)**

**(a)**

**(b)**

The reflection and refraction waveforms of the SV wave are shown in Figure 4. The viscous-spring boundary absorbs the reflected wave very well. The close agreement between the numerical solution and the analytic solution proves that the oblique-incidence plane wave input method yields accurate simulation results in a 3D layered field with viscous-spring boundaries.

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. Seismic Response of Immersed Tunnel under Seismic Wave Oblique Incidence

##### 4.1. Calculation Model and Parameters

The tunnel immersed below the Tianjin Haihe River is a two-hole, three-pipe gallery structure (Figure 5(a)). The width of the immersed tunnel is 36.6 m, the height of its cross section is 9.65 m, and the thicknesses of the roof, floor, side wall, and partition wall are 1.35 m, 1.4 m, 1.0 m, and 0.6 m, respectively. The span of pipe gallery on both sides is 2.6 m, the span of the middle elements gallery is 2 m, and the span of the dual holes is 12.5 m. The tunnel consists of three pipes (*E*1–*E*3); the length of a single pipe is 85 m. The thickness of the crushed-stone soil cover above the structure is 1.0 m. The tunnel is composed of C40 concrete. The clay layer in the Haihe Tunnel is defined here as the soil layer, which is divided into three sublayers. The material parameters of the soil and structure are shown in Tables 2 and 3.

**(a)**

**(b)**

We established a 3D finite-element model (300 m × 255 m × 55 m) of the Haihe River immersed tunnel using the parameters discussed above. The overall model includes three pipes (*E*1–*E*3), soil layers, and superimposed soil (Figure 5(b)). Around the cut-off boundary and the bottom boundary, a viscous-spring artificial boundary approximates an infinite field. We used the Mohr–Coulomb model for the soil and the plastic damage model in ABAQUS for the concrete.

##### 4.2. Seismic Wave Input

A Tianjin wave was input at the bottom of the model with a 20 s duration, 0.02 s interval, and peak acceleration of 3.1 m/s^{2} (Figure 6). Per the incident angle definition above (Figure 2) and the critical incident angle of the SV wave in this model (32.3°), the propagation directions of seismic waves are selected as = 0°, 30°, 60° and = 0°, 15°, 30°; this produced nine angle-of-incidence combinations. We used the dynamic relaxation method for static application. The seismic wave was input after 1 s of the static calculation; the entire calculation time was 21 s.

##### 4.3. Results

The pipe joint and shear key are significant components in the immersed tunnel—their inherent weakness may result in serious, even irreversible damage to the whole tunnel in the case of an earthquake. We first calculated and analyzed the mechanical properties of the pipe joint and the shear key, then the deformation of the tunnel (including three pipes) under typical earthquake conditions.

###### 4.3.1. Analysis of Shear Keys

The three shear keys on the left side of the *E*1-*E*2 section interface are numbered A, B, and C. Each shear key was marked with four corners (numbered 1–4) as monitoring points, as shown in Figure 7. The maximum principal stress, minimum principal stress, and shear stress of each monitoring point under the nine different incident angle combinations were extracted as shown in Tables 4–6. Among them, the maximum principal stress, minimum principal stress, and shear stress value of each monitoring point of shear key are maximal when the incidence angle is (30°, 30°), indicating that the angle combination of (30°, 30°) is the most unfavorable (which is close to the critical angle of incidence of SV waves in the half-space). The second-worst is the angle of incidence (60°, 30°). Figure 8 shows a stress cloud chart of the horizontal shear key A and vertical shear key B (30°, 30°).

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Tables 4–6 and Figure 8 altogether indicate the following.(1)The shear key produces significant compressive stress during an earthquake while the shear stress and tensile stress are relatively low.(2)When the seismic wave incidence transforms from vertical to oblique along the pipe cross section (that is, gradually increases), and the incident angle is small (less than 15°), the stress value of the shear key-monitoring point is fairly stable. With further increase in the incident angle, the stress value increases significantly to a maximum stress of 25.53 MPa, which is 7.1 times greater than the corresponding vertical incidence.(3)When the incidence direction of the seismic wave transforms from the tunnel transverse direction to the axial (that is, gradually increases), the maximum principal stress, minimum principal stress, and shear stress of the horizontal shear key A and vertical shear key B increase at first and then decrease. The maximum principal stress of the vertical shear key C increases at first and then decreases while its minimum principal stress and shear force continue to increase.(4)The horizontal shear stress distribution of compressive stress and shear stress is very uniform, the top corner has stress concentration phenomenon, and the tensile stress is mainly concentrated at the root of the shear key. The stress of the vertical shear keys B and C assumes a trapezoidal distribution, that is, the stress at the lower part of the shear key is greater than the upper stress. Tables 4–6 show that the stress of shear key B on the side wall is greater than that of shear key C on the middle partition wall. Seismic designers, to this effect, are recommended to focus on the root of the shear key which connects the shear key and immersed pipe body as well as the vertical shear key; the choice of reinforced concrete materials for the vertical shear key in the outer partition wall is particularly important.

###### 4.3.2. Dislocation (Opening) of Pipe Joint

As shown in Figure 9, we selected five nodes of the roof and the floor, respectively, at the *E*1–*E*2 pipe joint interface and determined the relationship between the dislocation (opening) and the incident angle of the seismic wave. The results are shown in Figure 10.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

The vertical dislocation value can be obtained without indicating the direction of dislocation, as shown in Figure 10. In fact, the vertical dislocation directions at both ends of the cross section of the model are opposite.

The dislocation volume at the joint is significantly affected by the seismic wave incidence angle. The horizontal and vertical dislocations of the roof and floor of the immersed tunnel joints increase as the local seismic wave incidence changes from vertical to oblique along the pipe cross section ( increases) and that of the maximum value is under (30°, 30°). The maximum value of horizontal dislocation at the roof is 10.64 mm, which is 11.9 times greater than the corresponding vertical incidence. The maximum horizontal dislocation at the floor is 2.45 mm, which is 3.7 times greater than the corresponding vertical incidence. The maximum vertical dislocation at the roof is 4.44 mm, which is 6.5 times greater than the corresponding vertical incidence. The maximum vertical dislocation of the floor is 4.47 mm, which is 6.2 times greater than the corresponding vertical incidence.

When the incidence direction of the seismic wave transforms from the tunnel transverse direction to the axial ( increases), the horizontal dislocation at the position of the roof and the floor of the pipe joint interface first increases and then decreases as vertical dislocation continuously increases. The lengthwise opening at the tunnel joint is very small (0.2 mm maximum at the bottom of the joint).

The horizontal dislocations at joints are similar at a certain angle of incidence. The end nodes are larger than the intermediate nodes in the vertical dislocation, and the direction of dislocation is opposite; in other words, the immersed tunnel shows torsional deformation. When the incident angle is (30°, 30°), the D-value of the vertical dislocation between the end node and the intermediate node of the floor joint is maximal (0.014°), which is in the scope of security.

###### 4.3.3. Overall Deformation of Immersed Tunnel

The nodes from the immersed tunnel’s central axis including three pipes were used to draw horizontal and vertical dislocation curves at the maximum stress (Figure 11).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

Figure 11 shows that as the seismic wave incidence changes from vertical to oblique along the pipe cross section (that is, increases), the horizontal and vertical displacement of the immersed tunnel increases. When the incidence direction of seismic wave transforms from the tunnel transverse direction to the axial (that is, increases), the overall horizontal and vertical displacement of the tunnel first increases and then decreases.

Furthermore, when the seismic wave is transversely incident (), the horizontal and vertical displacement distributions at intermediate tunnel joints are larger than at both sides; the displacement of both sides of tunnel joints at the pipe ends is basically the same. The horizontal displacement difference and vertical displacement difference of both sides of the tunnel (*E*1 and *E*3) joint increase as *α*_{1} increases, resulting in axial tensile deformation. When the angle of incidence is (60°, 30°), the horizontal dislocation difference of the end tunnel section reaches its maximum value. The axial tensile strain is approximately 1.73*e *−* *7, which is in the scope of security.

###### 4.3.4. Distribution of Cracks in Pipe Joint Cross Section

The distribution of tensile damage factors of the cross section of the immersed tunnel joint under incident angle (30°, 30°) is shown in Figure 12. This is an approximate representation of the distribution and progression of cross section cracks.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

As shown in Figure 12, the end of the side wall produces cracks first and is followed by the midspan of the floor, which is connected to horizontal shear key cracks; the midspan of the roof then cracks followed by the midspan of the roof in the middle pipe gallery and four corner points of the tunnel roof. Key parts that should be strengthened in an actual engineering scenario include the tunnel corner, the midspan of the roof, and the floor—that is to say, the connections between the horizontal shear key and floor.

#### 5. Conclusion

The 3D seismic response of an underwater immersed tunnel under obliquely incident seismic waves was investigated in this study based on simulations of the viscous-spring artificial boundary and seismic load. The influence of seismic SV wave incident angle on the 3D immersed tunnel joint and overall deformation of the tunnel is discussed by modelling the Tianjin Haihe River tunnel structure. Our main conclusions can be summarized as follows.(1)The oblique incidence of seismic SV waves, unlike the vertical incidence, has a marked effect on the seismic response of the immersed tunnel.(2)Among nine seismic wave combinations, the stress of the shear key is the largest under 30°, 30°; this is the worst-case angle to the site discussed in this paper. The tensile stress of the horizontal shear key is mainly concentrated where the shear key is connected to the main body of the immersed tube. The stress in the lower part of the vertical shear key is larger than that in the upper part, which makes it an important design point.(3)The horizontal and vertical dislocation of the roof and the floor of the immersed tunnel joints increase as the seismic wave incidence angle increases; they are maximal at an incident angle of 30°, 30°. The maximum relative dislocation of tunnel joints observed in this study was about 10 mm, while the average residual compression of GINA water stops is 15 mm–25 mm. In other words, the dislocation is within a safe range. The vertical dislocation of the end section is greater than that of the intermediate node and the dislocation directions are opposite, which suggests that the deformation of the immersed tunnel is torsional.(4)The overall horizontal and vertical displacement of the immersed tunnel increases as the seismic wave incidence angle increases. The displacement difference between the two ends of the same tunnel segment also increases, which can cause axial tensile deformation of the tunnel joints.(5)Under the oblique incidence of SV waves, cracks appear in the ends of the side walls (the corners of the tunnel), the midspan of the roof and the floor, and the part of the horizontal shear key connected to the floor. In an actual engineering scenario, construction measures and anticracking treatment should be targeted to these key areas.

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

The research described in this paper was financially supported by the National Natural Science Foundation of China under grant numbers 51678389 and 51678390.