# Three-Dimensional Nonlinear Seismic Response of Shield Tunnel Spatial End Structure

Authors: Yijing Lu; Yongsheng Song; Yanzhen Wang; Jiale Yuan
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The stiffness mutation of shield tunnel-shaft junction makes the tunnel structure affected by the differential displacement and forms a complex spatial effect. Taking the subsea shield tunnel crossing under the Shantou Gulf, China, as a case study, a three-dimensional finite element global model and a refined local spatial end submodel are established. The nonlinear dynamic behaviors of the seabed soil and concrete, the simulation of the bolt joints between ring segments by using cohesive models and the SMA shape memory alloy flexible joints, and the input ground motions produced by scaling from the high-level earthquake records are considered in detail. The results show that the shield tunnel spatial end structure increases nonlinearly in response to the increase of seismic motion intensity. The opening width and the deformation between ring segments at the vault and the outside spandrel are larger, and serious seismic damage and stress concentration exist at the conjugate 45° directions of shaft. The seismic responses of the tunnel-shaft junction subjected to the seismic motions with rich low-frequency components are much stronger than those of seismic motions with rich high-frequency components. Adding SMA flexible joints, the structural deformation caused by seismic motion propagation can be induced to the preset flexible joint, and the structural damage and stress concentration can be effectively reduced. The seismic response characteristics of shield tunnel spatial end structure calculated by the global model are consistent with those calculated by the submodel, while the seismic response of the submodel is greater than that of the global model.


---

## Body

## 1. Introduction

Maritime trade is the pulse of the contemporary world economy, and the coastal zone is the link to it. This region concentrates on the vast majority of the population of the modern world, carries the dual functions of economy and life, and has a special significance for the world today. The local coastal zone is frequently scattered with islands, and the maritime traffic between the islands will greatly enhance the carrying capacity of this area.The undersea shield tunnel is less affected by tides, currents, waves, and other meteorological and hydrological environments and has strong continuous carrying capacity after completion, which has become the first choice for the construction of cross-sea passage. The shaft usually needs to be set at the shore of the tunnel, for the installation of drainage, power supply, ventilation, and other important facilities inside the tunnel. When attacked by a high-level earthquake, the submarine tunnel-shaft junction is prone to the problem of stress concentration and deformation, which is the weak part of seismic resistance. During the Michoacán earthquake in Mexico in 1985 [1, 2] and the Kobe earthquake in 1995 [3–5], some seismic disasters occurred, such as relative segment dislocation, bolt shearing between the rings, axial cracks in the vault, and water leakage. Therefore, it is vital to conduct special seismic response analysis and seismic design for the submarine shield tunnel-shaft junction.Some scholars have studied and analyzed the seismic response of submarine shield tunnel. At present, the longitudinal seismic response analysis methods of tunnel are mainly divided into response displacement method and soil-structure integral numerical analysis method [6]. Anastasopoulos et al. [7] used the beam-spring model and Zhang et al. [8] derived the seismic analytical solution of the tunnel-shaft junction. Chen et al. [9] and Tang et al. [10] studied the nonlinear longitudinal seismic response characteristics of shield tunnel and utility tunnel based on the generalized response displacement method, respectively. In this method, soil-structure interaction is simulated by soil spring, and the stiffness coefficient of soil spring has a great influence on the results [11, 12]. Numerical analysis method directly divides the strata into units according to the geological profile. In practical engineering, the soil and the structure are mostly simplified as plane strain problems and adopted as two-dimensional models, and the direction of the input seismic motion is limited to the tunnel axial direction [13, 14]. For the seismic direction is perpendicular to the tunnel axial direction [15] and the geological conditions mutation or the tunnel structure changes, it is necessary to establish a three-dimensional model for seismic response analysis. Through the earthquake damage investigation of shield tunnel, Koizumi [16] found that the junction between the shaft (such as the starting shaft, the middle shaft, and the receiving shaft) and the tunnel is prone to damage. Duran et al. [17] found that the torsion of the shaft under earthquake has a significant influence on the deformation and strain of the adjacent shield tunnel rings. Zhao et al. [18], Yang et al. [19], and Song et al. [20] studied the end constraint effect of tunnel structure and believed that since the seismic performance of underground structure is closely related to its geometric shape and stiffness distribution, the mutation of stiffness at the structural changing position will cause the tunnel structure to be affected by differential displacement and form complex spatial effect. At the same time, the material constitutive model used in the existing research is too simple to reflect the real influence of material nonlinear characteristics on the seismic response of the shield tunnel-shaft junction. It is hard to analyze the distribution of cross section stress and opening width without considering the discontinuous deformation between the pipe rings. However, when calculating the numerical model with refined mesh, the computational precision is high and a lot of computing resources are spent. The submodel technology is a method for local precise analysis of large and complex finite element global models [21]. Wu et al. [22] established a three-dimensional fine model of tunnel entrance by using submodel technology to analyze the process and extent of seismic damage of tunnel lining. Banushi et al. [23] conducted earthquake induced permanent ground deformation (PGD) analysis on the submodel of buried pipeline and obtained its accurate and efficient seismic performance. Chen et al. [24] used the generalized response displacement method based on submodel technology to conduct nonlinear seismic analysis of submarine shield tunnel and studied the opening width between pipe rings under high-level earthquake.Based on the ABAQUS software platform and using the spatial end structure of the Suai Undersea Tunnel-shaft junction in Shantou as the engineering background, this paper establishes a three-dimensional model and uses the submodel technology to establish a fine model for the local spatial end structure considering segment longitudinal seam stitching. Considering seabed dynamic nonlinear characteristics of soil, the shield tunnel longitudinal bolt connection between pipe rings, and the flexible connection between the shield tunnel-shaft junctions, making nonlinear dynamic time-history analysis on the tunnel-shaft junction global model and local spatial end structure submodel, and discussing the influence of the characteristics of the input seismic motions and the shock absorption measures on the spatial end structure, the research results can provide a reference for seismic design of shield tunnel-shaft junction.


## 2. Submodel Technology Based on ABAQUS Software

The submodel technology is based on Saint Venant’s principle, which can refine the mesh of a certain part of the global model and further seek a more accurate solution. As shown in Figure1, the relatively rough initial global model is first analyzed, and then the local model that needs to refine the mesh is cut out. Submodel technology is used to drive the local model through the results of the global model. The driving variables refer to the variables that constrain the submodel and generally match the results of the global model. The universal submodel technique is based on the nodes, and the corresponding response of the global model is interpolated to the nodes of the submodel by using the response of nodes (including displacement, temperature, or pressure degrees of freedom). The whole calculation process does not need to extract the response quantity of the global model node and input it to the submodel boundary again. The finite element software ABAQUS identifies the relationship between the submodel and global model by using the spatial positions of them and takes the displacement of the corresponding driving subboundary nodes in the global model as the additional boundary condition of the submodel. The displacement of the driving subboundary nodes are functions of time in dynamic analysis.Figure 1
Diagram of application of submodel technology.


## 3. Analysis Model of Shield Tunnel Spatial End Structure

### 3.1. Tectonic Setting and Engineering Situation

The Suai submarine tunnel project is located between the Shantou Bay Bridge in the downtown area of Shantou and Queshi Bridge, connecting the north and south sides of Shantou. It is the first extra-long and super-large diameter undersea shield tunnel located in the earthquake intensity zone of VIII degree in China. The Suai submarine tunnel site is characterized by variable geological structure, multiple northeast and northwest faults interwoven, and being close to the Taiwan Strait with strong seismic activity. The possible seismic impact will mainly come from the ocean area in the future. Two high-level earthquakes of magnitude 7 or above occurred about 54 km away from the project site in the nearby field. The engineering site has suffered many earthquakes with the intensity of V∼VIII in history, among which there were three earthquakes with the intensity of VIII. Figure2 shows the historical earthquake and fracture distribution in the tunnel site. The starting shaft of the submarine tunnel is located in the cofferdam of the south bank, and the outer diameter of the single pipe of the shield tunnel is 14.5 m and the inner diameter is 13.3 m. The shaft of the Suai submarine shield tunnel is 50 m long and 25 m wide, with 30 m buried depth and 1.2 m wall thickness. The section size of the open cut tunnel is approximately 40 × 20 m, and the outer wall thickness is 1.3 m. Figure 3 shows the geological profile of seabed soil along the longitudinal axis of the tunnel.Figure 2
Regional tectonic setting and historical earthquake data of the Suai submarine tunnel site.Figure 3
Geological sketch of section along tunnel longitudinal axis.


### 3.2. Three-Dimensional Numerical Modeling

Figure4(a) shows the global model of the shield tunnel spatial end structure, the longitudinal length of which is 225 m and the transverse width is 180 m. The interface with shear wave velocity ≥500 m/s and no lower wave velocity below it is taken as the seismic bedrock surface, and the depth is 100 m. According to the concept of “equivalent continuum model,” the shield tunnel is equivalent to a homogeneous ring without considering segment splicing in the global model. Each ring of shield tunnel is 2 m long, and 50 rings are taken from a single hole. The ring is equivalent to a uniform ring without considering segment splicing. The soil and structure are simulated by Lagrange solid element (C3D8R). The mesh size of the shaft-tunnel junction structure is 1 m. Meanwhile, the mesh size of the site is refined: the tunnel-shaft junction surrounding soil is 1m, the vertical mesh size (1∼3 m) increases with the depth, and the horizontal mesh size (1∼5 m) increases with the distance from the shaft. The number of elements in the global model is 527616, while the degrees of freedom are 2194357. The contact surface is defined between the structure and the soil, holding “hard contact” in the normal direction, while the tangential friction coefficient is 0.3 [25]. The beam element is used to simulate the steel and embedded into the concrete. Viscoelastic artificial boundaries are set at the side and bottom of the model [26]. The spatial end structure of the one-sided shaft-tunnel junction in the global model is intercepted, and the splitting of longitudinal segments is considered by remodeling the pipe ring. The single pipe ring is composed of 10 segments, as shown in Figure 4(b). In order to retain the influence of soil-structure interaction on the submodel, the 10 m soil layer around the structure is also intercepted to compose the submodel, which should be consistent with the spatial location relationship of the global model. The number of elements in the submodel is 10879, while the degrees of freedom are 44239. All the outer surfaces of the submodel which contact the global model except the upper surface are set as the displacement driving subboundary.Figure 4
Three-dimensional finite element model of the shield tunnel spatial end structure. (a) Global model. (b) Submodel.


### 3.3. Parameters of Analysis Model

#### 3.3.1. Soil Nonlinear Cyclic Behaviors

A non-Masing constitutive model [27, 28] (DCZ model) based on the Davidenkov backbone curve under irregular loading-unloading criteria is used to describe the nonlinear hysteretic behavior of soil under seismic motion. The mathematical description of the DCZ model is as follows.The initial backbone curve uses the Davidenkov backbone curve:(1)τ=G⋅γ=Gmax⋅γ⋅1−Hγ,in which(2)Hγ=γ/γr2B1+γ/γr2BA.If the reversed strain occurs, the subsequent strain-stress curve follows a path along the assigned orientation curve from the stress reversal point to the last extreme value point of the stress-strain space in the previous cycles, as described by the three following equations:(3)τ−τc=Gmax⋅γ−γc⋅1−Hγ−γc2n,(4)2ncγr2B=γex±γc2B⋅1−RR,(5)R=1−τex±τcGmax⋅γex±γc1/A.If the unloading-reloading curve intersects the backbone curve, the stress-strain curve follows the backbone curve until the next stress reversal.τ and γ are shear stress and shear strain, respectively; Gmax is the initial shear modulus; γr is the reference shear strain; both A and B are dimensionless exponents; nc is the proportion coefficient of hysteresis loop; τc and γc are the shear stress and strain at the last stress reversal point, respectively; and τex and γex are the shear stress and strain at the last extreme value point, respectively. During irregular load-unloading period, “±” is positive when the hysteresis loop is bent upward, while it is negative when it is bent downward. Figure 5 is a schematic diagram of the stress-strain curve described by the DCZ model. The shear modulus and the damping ratio curves obtained from the resonant column test for typical undisturbed soil samples are shown in Figure 6, and the parameter values of the adopted soil constitutive model are listed in Table 1.Figure 5
Stress-strain curves of the DCZ model under irregular loading-unloading-reloading.Figure 6
TheG/Gmax∼γ and λ∼γ curves of typical soils for the submarine tunnel site.Table 1
Parameters of the constitutive model of the subsea soils.
SoilABγ0 (×10−4)Silt1.030.405.1Mucky soil1.010.396.7Sandy soil1.200.377.4Silty clay1.120.419.0Weathered rock1.300.4010.0


#### 3.3.2. Tunnel Segments and Ring Joint Bolts

The discontinuous structure of “segment ring + joint” leads to the state nonlinear characteristics of longitudinal deformation of shield tunnel. Segment longitudinal joint is composed of a ring joint surface and a connecting bolt. It is assumed that the connecting bolts bear all the tension in the tensile zone of the ring joint surface, while the segment concrete bears all the pressure in the compression zone. The pretightening force of bolts is not considered temporarily. As shown in Figure7, contact surfaces are set between pipe rings to simulate the state nonlinear behavior of ring joints and connecting bolts. Normal hard contact is used in the compression zone and cohesive model is used in the tension zone. According to the principle that adjacent pipe rings receive the same resultant force with the same opening width amount under unidirectional tensile condition, the bolt stiffness is equivalent to a continuous and uniform distribution along the segment contact area. Physical and mechanical parameters of connecting bolts are listed in Table 2.Figure 7
Deformation characteristics and relation between the forces and deformations between adjacent pipe rings.Table 2
Ring joint bolt parameters of shield tunnel segments.
Diametersd0 (mm)Lengthl0 (mm)Elastic modulusEs (MPa)Tensile yield stressfy (MPa)Tensile ultimate stressfm (MPa)Elastoplastic stiffness ratioα367502.06E56408000.01As the opening width of the ring joints increases, the connecting bolt will yield and the stiffness of the ring joint will decay. Maximum separation criterion is applied for describing the stiffness damage of the cohesive material; that is, the stiffness begins to decay when the normal or tangential deformation reaches the corresponding limit. Meanwhile, the stiffness decay is realized by the damage indexD:(6)K2=1−DK1,where K1 is the stiffness of bolts in the elastic deformation level; K2 is the stiffness of the bolt after yield; and D is the damage index, which changes from “0” to “1.” In ABAQUS, the damage index D is defined by the indirect method of elastic ultimate deformation δy and ultimate deformation δu.In the cohesive model, the elastic ultimate deformationδy is obtained when the bolt reaches the yield stress. The ultimate deformation δu indicates that the bolt stress begins to decrease when the opening width between pipe rings reaches the ultimate displacement:(7)δy=NyK1=nAsfynEsAs/l0=fyl0Es,δu=δy+Nm−NyK2=δy+nAsfm−nAsfynαEsAs/l0=fyEs+fm−fyαEs⋅l0,where Ny is the elastic ultimate tensile force of bolts; Nu is the ultimate tensile force of bolts, and l0 is the length of bolts. For the calculation of cohesive model parameters, the reader is referred to [24], and the relevant parameters are listed in Table 3. 42 bolts are set between rings, and 2 bolts are set between the longitudinal segments. The waterproof limit is 15 mm.Table 3
Parameters of surface-based cohesive behavior in finite element model.
Knn (E8N/m3)Kss (E9N/m3)Ktt (E9N/m3)Elastic deformation limitδy (mm)Ultimate displacementδu (mm)4.4818(4.657)4.4818(4.657)4.4818(4.657)2.360.6


#### 3.3.3. SMA Shape Memory Alloy Flexible Joint

Shape memory alloy is a good shock absorber because of its superelastic properties. The commonly used memory alloy is composed of 55.9% Ni and 44.1% Ti, with a recoverable strain of 8% and an ultimate strain of 17%. In addition, memory alloys can dissipate the energy by forming hysteresis curves within recoverable strains. In the past two decades, shape memory alloys as a new material have attracted extensive attention from scholars at home and abroad, and many scholars [29, 30] have conducted in-depth studies on the numerical simulation of their constitutive relations. In this paper, a flag-SMA piecewise linear simplified constitutive model is adopted [31]. According to the material properties of Ni-Ti shape memory alloy obtained by Liu et al. [32] and Zhou et al. [33], the superelasticity (the material constitutive is shown in Figure 8) is used to simulate shape memory alloy in ABAQUS, 2017 version. The flexible antiseismic joint is shown in Figure 9. A total of 40 sets of energy dissipation flexible joint with a diameter of 25 mm are set in a single ring. The coupling connection of the flexible joint elements with the shaft end wall and the shield tunnel is adopted (i.e., the degree of freedom of the node is consistent), and the interaction setting of the surrounding rock is consistent with that of the shield tunnel. According to the principle that adjacent pipe rings receive the same force with the same opening width, the stiffness of memory alloy material is equivalent to a continuous and uniform distribution along the segment contact area, as described by the following equation:(8)K′=nAsAcK,where K′ is the material stiffness after equivalence; As is the cross-sectional area of a single energy dissipation-damping joint; Ac is the cross-sectional area of pipe rings; n is the number of energy dissipation-damping joints; and K is the material stiffness before equivalence. The parameters of superelasticity are listed in Table 4.EM is Young’s modulus of martensite; γM is Poisson’s ratio of martensite; εL is the uniaxial transformation strain; σtLs is the stress at which the transformation begins during loading in tension; σtLE is the stress at which the transformation ends during loading in tension; σtUS is the stress at which the reverse transformation begins during loading in tension; σtUE is the stress at which the reverse transformation ends during loading in tension; σcLS is the stress at which the transformation begins during loading in compression, as a positive value; T0 is the reference temperature; (δσ/δT)L is the slope of the stress versus temperature curve for loading; and (δσ/δT)U is the slope of the stress versus temperature curve for unloading.Figure 8
Lag-SMA simplified constitutive model.Figure 9
Simulation of the SMA shape memory alloy flexible joint.Table 4
The parameters of superelasticity (unit: MPa).
EMγMεLσtLsσtLEσtUSσtUEσcLST0(δσ/δT)L(δσ/δT)U19.70.330.0670.151.551.40.0370.15150.010.006


#### 3.3.4. Concrete

Concrete damage model (CDP) is used to describe the irreversible cumulative plastic deformation and stiffness loss of concrete under cyclic loading [34]. The CDP model needs to define plastic parameters and damage parameters. Plastic parameters describe the change form of the material yield surface, while damage parameters describe the loading and unloading characteristics of the material from tensile and compressive stress-strain relations and damage index. The concrete strength grade of the tunnel and shaft is C60, and the parameters of the CDP are listed in Table 5. According to the proportional strain method, the damage index of tension and compression is obtained [35, 36].Table 5
Parameters of the concrete from the plasticity test.
Elastic modulusEc (GPa)Poisson’s ratioμInitial yield stressσco (MPa)Compression strengthσcu (MPa)Tensile strengthσtu (MPa)Dilatation angle Ψ (°)Eccentricityδfbo/fcoKcViscosity coefficientu360.233.2747.573.8838.00.11.162/30.0005fbo/fco is the ratio of biaxial to uniaxial ultimate compressive strength; Kc is the ratio of biaxial to uniaxial ultimate compressive strength; and u is the ratio of the second stress invariant on the tensile and compressive meridian planes.


### 3.4. Selection of Bedrock Input Motion

According to the seismic safety assessment results of Shantou Bay Tunnel site, the exceeded probabilities of 63%, 10%, and 2% of Peak bedrock acceleration are approximately 0.1g, 0.2 g, and 0.4 g, respectively, where “g” is the acceleration of gravity. Due to the lack of strong earthquake records on this site, on the basis of historical seismic data, the bedrock acceleration records of two strong earthquakes in the foreign seismic network are selected as the bedrock input ground motion. The detailed information is listed in Table 6. The acceleration time history of seismic records and its corresponding Fourier spectrum is shown in Figure 10. The Darfield record focuses on low frequency, while the Iwate record develops high-frequency components. There are significant differences in the spectral characteristics of the acceleration records between the two records. In order to consider the influence of the input ground motion intensity, the PGA (peak ground acceleration) of the horizontal direction (vertical to the longitudinal axis of the tunnel) recorded in the original is adjusted to three earthquake levels: 0.1 g (low-level earthquake (LLE)), 0.2 g (moderate-level earthquake (MLE)), and 0.4 g (high-level earthquake (HLE)), and the ratio of horizontal PGA to vertical PGA is set as 1 : 0.65.Table 6
Details of original ground motion records used for site response analysis.
Earthquake eventStation site conditionMagnitudeDistance from the epicenter (km)Ground motion componentPGA (g)Duration D5-95 (s)Mean frequencyfp (Hz)DarfieldPage Road Pumping (PRP) StationMs 7.153.53NS0.23022.290.34UD0.31815.0811.22IwateIWTH27Ms 7.256.85NS0.34720.448.54UD0.12822.414.21Figure 10
Acceleration time history and Fourier spectra of the input bedrock motions.


### 3.5. Monitoring Point Positions

The displacement time history of the submodel should be extracted from the subboundary of the global model, choosing three monitoring points at the subboundary: Point 1 to Point 2 (as shown in Figure4(b)), which are selected to observe its displacement time-history response. The opening width and dislocation are important indexes to evaluate the seismic performance of shield tunnel. Excessive deformation will lead to bolts failure, and water leakage will occur in serious cases. In this paper, the ring joint connecting the first ring of shield tunnel and the wall of shaft in the global model and submodel is referred to as tunnel-shaft joint for short. Based on symmetry, channel “A” is taken as the research object, and the opening width at the tunnel-shaft joint is analyzed from the positions of 8 measuring points shown in Figure 11(a). The peak opening width is defined as the maximum absolute value of the difference between the displacement time history of the right end of the left pipe ring and the corresponding point displacement time history of the left end of the right pipe ring. The first ring of shield tunnel in the submodel is referred to as a tunnel-shaft head ring. The longitudinal seam deformation analysis is conducted on the positions of 10 measuring points between the first ring joints shown in Figure 11(b). The peak deformation is defined as the maximum absolute value of the contact deformation (Copen) time history of adjacent measuring points. The shield tunnel spatial end structure is referred to as a tunnel-shaft junction for short.Figure 11
Monitoring positions at the joint of tunnel-shaft junction.
(a)(b)


## 3.1. Tectonic Setting and Engineering Situation

The Suai submarine tunnel project is located between the Shantou Bay Bridge in the downtown area of Shantou and Queshi Bridge, connecting the north and south sides of Shantou. It is the first extra-long and super-large diameter undersea shield tunnel located in the earthquake intensity zone of VIII degree in China. The Suai submarine tunnel site is characterized by variable geological structure, multiple northeast and northwest faults interwoven, and being close to the Taiwan Strait with strong seismic activity. The possible seismic impact will mainly come from the ocean area in the future. Two high-level earthquakes of magnitude 7 or above occurred about 54 km away from the project site in the nearby field. The engineering site has suffered many earthquakes with the intensity of V∼VIII in history, among which there were three earthquakes with the intensity of VIII. Figure2 shows the historical earthquake and fracture distribution in the tunnel site. The starting shaft of the submarine tunnel is located in the cofferdam of the south bank, and the outer diameter of the single pipe of the shield tunnel is 14.5 m and the inner diameter is 13.3 m. The shaft of the Suai submarine shield tunnel is 50 m long and 25 m wide, with 30 m buried depth and 1.2 m wall thickness. The section size of the open cut tunnel is approximately 40 × 20 m, and the outer wall thickness is 1.3 m. Figure 3 shows the geological profile of seabed soil along the longitudinal axis of the tunnel.Figure 2
Regional tectonic setting and historical earthquake data of the Suai submarine tunnel site.Figure 3
Geological sketch of section along tunnel longitudinal axis.


## 3.2. Three-Dimensional Numerical Modeling

Figure4(a) shows the global model of the shield tunnel spatial end structure, the longitudinal length of which is 225 m and the transverse width is 180 m. The interface with shear wave velocity ≥500 m/s and no lower wave velocity below it is taken as the seismic bedrock surface, and the depth is 100 m. According to the concept of “equivalent continuum model,” the shield tunnel is equivalent to a homogeneous ring without considering segment splicing in the global model. Each ring of shield tunnel is 2 m long, and 50 rings are taken from a single hole. The ring is equivalent to a uniform ring without considering segment splicing. The soil and structure are simulated by Lagrange solid element (C3D8R). The mesh size of the shaft-tunnel junction structure is 1 m. Meanwhile, the mesh size of the site is refined: the tunnel-shaft junction surrounding soil is 1m, the vertical mesh size (1∼3 m) increases with the depth, and the horizontal mesh size (1∼5 m) increases with the distance from the shaft. The number of elements in the global model is 527616, while the degrees of freedom are 2194357. The contact surface is defined between the structure and the soil, holding “hard contact” in the normal direction, while the tangential friction coefficient is 0.3 [25]. The beam element is used to simulate the steel and embedded into the concrete. Viscoelastic artificial boundaries are set at the side and bottom of the model [26]. The spatial end structure of the one-sided shaft-tunnel junction in the global model is intercepted, and the splitting of longitudinal segments is considered by remodeling the pipe ring. The single pipe ring is composed of 10 segments, as shown in Figure 4(b). In order to retain the influence of soil-structure interaction on the submodel, the 10 m soil layer around the structure is also intercepted to compose the submodel, which should be consistent with the spatial location relationship of the global model. The number of elements in the submodel is 10879, while the degrees of freedom are 44239. All the outer surfaces of the submodel which contact the global model except the upper surface are set as the displacement driving subboundary.Figure 4
Three-dimensional finite element model of the shield tunnel spatial end structure. (a) Global model. (b) Submodel.


## 3.3. Parameters of Analysis Model

### 3.3.1. Soil Nonlinear Cyclic Behaviors

A non-Masing constitutive model [27, 28] (DCZ model) based on the Davidenkov backbone curve under irregular loading-unloading criteria is used to describe the nonlinear hysteretic behavior of soil under seismic motion. The mathematical description of the DCZ model is as follows.The initial backbone curve uses the Davidenkov backbone curve:(1)τ=G⋅γ=Gmax⋅γ⋅1−Hγ,in which(2)Hγ=γ/γr2B1+γ/γr2BA.If the reversed strain occurs, the subsequent strain-stress curve follows a path along the assigned orientation curve from the stress reversal point to the last extreme value point of the stress-strain space in the previous cycles, as described by the three following equations:(3)τ−τc=Gmax⋅γ−γc⋅1−Hγ−γc2n,(4)2ncγr2B=γex±γc2B⋅1−RR,(5)R=1−τex±τcGmax⋅γex±γc1/A.If the unloading-reloading curve intersects the backbone curve, the stress-strain curve follows the backbone curve until the next stress reversal.τ and γ are shear stress and shear strain, respectively; Gmax is the initial shear modulus; γr is the reference shear strain; both A and B are dimensionless exponents; nc is the proportion coefficient of hysteresis loop; τc and γc are the shear stress and strain at the last stress reversal point, respectively; and τex and γex are the shear stress and strain at the last extreme value point, respectively. During irregular load-unloading period, “±” is positive when the hysteresis loop is bent upward, while it is negative when it is bent downward. Figure 5 is a schematic diagram of the stress-strain curve described by the DCZ model. The shear modulus and the damping ratio curves obtained from the resonant column test for typical undisturbed soil samples are shown in Figure 6, and the parameter values of the adopted soil constitutive model are listed in Table 1.Figure 5
Stress-strain curves of the DCZ model under irregular loading-unloading-reloading.Figure 6
TheG/Gmax∼γ and λ∼γ curves of typical soils for the submarine tunnel site.Table 1
Parameters of the constitutive model of the subsea soils.
SoilABγ0 (×10−4)Silt1.030.405.1Mucky soil1.010.396.7Sandy soil1.200.377.4Silty clay1.120.419.0Weathered rock1.300.4010.0


### 3.3.2. Tunnel Segments and Ring Joint Bolts

The discontinuous structure of “segment ring + joint” leads to the state nonlinear characteristics of longitudinal deformation of shield tunnel. Segment longitudinal joint is composed of a ring joint surface and a connecting bolt. It is assumed that the connecting bolts bear all the tension in the tensile zone of the ring joint surface, while the segment concrete bears all the pressure in the compression zone. The pretightening force of bolts is not considered temporarily. As shown in Figure7, contact surfaces are set between pipe rings to simulate the state nonlinear behavior of ring joints and connecting bolts. Normal hard contact is used in the compression zone and cohesive model is used in the tension zone. According to the principle that adjacent pipe rings receive the same resultant force with the same opening width amount under unidirectional tensile condition, the bolt stiffness is equivalent to a continuous and uniform distribution along the segment contact area. Physical and mechanical parameters of connecting bolts are listed in Table 2.Figure 7
Deformation characteristics and relation between the forces and deformations between adjacent pipe rings.Table 2
Ring joint bolt parameters of shield tunnel segments.
Diametersd0 (mm)Lengthl0 (mm)Elastic modulusEs (MPa)Tensile yield stressfy (MPa)Tensile ultimate stressfm (MPa)Elastoplastic stiffness ratioα367502.06E56408000.01As the opening width of the ring joints increases, the connecting bolt will yield and the stiffness of the ring joint will decay. Maximum separation criterion is applied for describing the stiffness damage of the cohesive material; that is, the stiffness begins to decay when the normal or tangential deformation reaches the corresponding limit. Meanwhile, the stiffness decay is realized by the damage indexD:(6)K2=1−DK1,where K1 is the stiffness of bolts in the elastic deformation level; K2 is the stiffness of the bolt after yield; and D is the damage index, which changes from “0” to “1.” In ABAQUS, the damage index D is defined by the indirect method of elastic ultimate deformation δy and ultimate deformation δu.In the cohesive model, the elastic ultimate deformationδy is obtained when the bolt reaches the yield stress. The ultimate deformation δu indicates that the bolt stress begins to decrease when the opening width between pipe rings reaches the ultimate displacement:(7)δy=NyK1=nAsfynEsAs/l0=fyl0Es,δu=δy+Nm−NyK2=δy+nAsfm−nAsfynαEsAs/l0=fyEs+fm−fyαEs⋅l0,where Ny is the elastic ultimate tensile force of bolts; Nu is the ultimate tensile force of bolts, and l0 is the length of bolts. For the calculation of cohesive model parameters, the reader is referred to [24], and the relevant parameters are listed in Table 3. 42 bolts are set between rings, and 2 bolts are set between the longitudinal segments. The waterproof limit is 15 mm.Table 3
Parameters of surface-based cohesive behavior in finite element model.
Knn (E8N/m3)Kss (E9N/m3)Ktt (E9N/m3)Elastic deformation limitδy (mm)Ultimate displacementδu (mm)4.4818(4.657)4.4818(4.657)4.4818(4.657)2.360.6


### 3.3.3. SMA Shape Memory Alloy Flexible Joint

Shape memory alloy is a good shock absorber because of its superelastic properties. The commonly used memory alloy is composed of 55.9% Ni and 44.1% Ti, with a recoverable strain of 8% and an ultimate strain of 17%. In addition, memory alloys can dissipate the energy by forming hysteresis curves within recoverable strains. In the past two decades, shape memory alloys as a new material have attracted extensive attention from scholars at home and abroad, and many scholars [29, 30] have conducted in-depth studies on the numerical simulation of their constitutive relations. In this paper, a flag-SMA piecewise linear simplified constitutive model is adopted [31]. According to the material properties of Ni-Ti shape memory alloy obtained by Liu et al. [32] and Zhou et al. [33], the superelasticity (the material constitutive is shown in Figure 8) is used to simulate shape memory alloy in ABAQUS, 2017 version. The flexible antiseismic joint is shown in Figure 9. A total of 40 sets of energy dissipation flexible joint with a diameter of 25 mm are set in a single ring. The coupling connection of the flexible joint elements with the shaft end wall and the shield tunnel is adopted (i.e., the degree of freedom of the node is consistent), and the interaction setting of the surrounding rock is consistent with that of the shield tunnel. According to the principle that adjacent pipe rings receive the same force with the same opening width, the stiffness of memory alloy material is equivalent to a continuous and uniform distribution along the segment contact area, as described by the following equation:(8)K′=nAsAcK,where K′ is the material stiffness after equivalence; As is the cross-sectional area of a single energy dissipation-damping joint; Ac is the cross-sectional area of pipe rings; n is the number of energy dissipation-damping joints; and K is the material stiffness before equivalence. The parameters of superelasticity are listed in Table 4.EM is Young’s modulus of martensite; γM is Poisson’s ratio of martensite; εL is the uniaxial transformation strain; σtLs is the stress at which the transformation begins during loading in tension; σtLE is the stress at which the transformation ends during loading in tension; σtUS is the stress at which the reverse transformation begins during loading in tension; σtUE is the stress at which the reverse transformation ends during loading in tension; σcLS is the stress at which the transformation begins during loading in compression, as a positive value; T0 is the reference temperature; (δσ/δT)L is the slope of the stress versus temperature curve for loading; and (δσ/δT)U is the slope of the stress versus temperature curve for unloading.Figure 8
Lag-SMA simplified constitutive model.Figure 9
Simulation of the SMA shape memory alloy flexible joint.Table 4
The parameters of superelasticity (unit: MPa).
EMγMεLσtLsσtLEσtUSσtUEσcLST0(δσ/δT)L(δσ/δT)U19.70.330.0670.151.551.40.0370.15150.010.006


### 3.3.4. Concrete

Concrete damage model (CDP) is used to describe the irreversible cumulative plastic deformation and stiffness loss of concrete under cyclic loading [34]. The CDP model needs to define plastic parameters and damage parameters. Plastic parameters describe the change form of the material yield surface, while damage parameters describe the loading and unloading characteristics of the material from tensile and compressive stress-strain relations and damage index. The concrete strength grade of the tunnel and shaft is C60, and the parameters of the CDP are listed in Table 5. According to the proportional strain method, the damage index of tension and compression is obtained [35, 36].Table 5
Parameters of the concrete from the plasticity test.
Elastic modulusEc (GPa)Poisson’s ratioμInitial yield stressσco (MPa)Compression strengthσcu (MPa)Tensile strengthσtu (MPa)Dilatation angle Ψ (°)Eccentricityδfbo/fcoKcViscosity coefficientu360.233.2747.573.8838.00.11.162/30.0005fbo/fco is the ratio of biaxial to uniaxial ultimate compressive strength; Kc is the ratio of biaxial to uniaxial ultimate compressive strength; and u is the ratio of the second stress invariant on the tensile and compressive meridian planes.


## 3.3.1. Soil Nonlinear Cyclic Behaviors

A non-Masing constitutive model [27, 28] (DCZ model) based on the Davidenkov backbone curve under irregular loading-unloading criteria is used to describe the nonlinear hysteretic behavior of soil under seismic motion. The mathematical description of the DCZ model is as follows.The initial backbone curve uses the Davidenkov backbone curve:(1)τ=G⋅γ=Gmax⋅γ⋅1−Hγ,in which(2)Hγ=γ/γr2B1+γ/γr2BA.If the reversed strain occurs, the subsequent strain-stress curve follows a path along the assigned orientation curve from the stress reversal point to the last extreme value point of the stress-strain space in the previous cycles, as described by the three following equations:(3)τ−τc=Gmax⋅γ−γc⋅1−Hγ−γc2n,(4)2ncγr2B=γex±γc2B⋅1−RR,(5)R=1−τex±τcGmax⋅γex±γc1/A.If the unloading-reloading curve intersects the backbone curve, the stress-strain curve follows the backbone curve until the next stress reversal.τ and γ are shear stress and shear strain, respectively; Gmax is the initial shear modulus; γr is the reference shear strain; both A and B are dimensionless exponents; nc is the proportion coefficient of hysteresis loop; τc and γc are the shear stress and strain at the last stress reversal point, respectively; and τex and γex are the shear stress and strain at the last extreme value point, respectively. During irregular load-unloading period, “±” is positive when the hysteresis loop is bent upward, while it is negative when it is bent downward. Figure 5 is a schematic diagram of the stress-strain curve described by the DCZ model. The shear modulus and the damping ratio curves obtained from the resonant column test for typical undisturbed soil samples are shown in Figure 6, and the parameter values of the adopted soil constitutive model are listed in Table 1.Figure 5
Stress-strain curves of the DCZ model under irregular loading-unloading-reloading.Figure 6
TheG/Gmax∼γ and λ∼γ curves of typical soils for the submarine tunnel site.Table 1
Parameters of the constitutive model of the subsea soils.
SoilABγ0 (×10−4)Silt1.030.405.1Mucky soil1.010.396.7Sandy soil1.200.377.4Silty clay1.120.419.0Weathered rock1.300.4010.0


## 3.3.2. Tunnel Segments and Ring Joint Bolts

The discontinuous structure of “segment ring + joint” leads to the state nonlinear characteristics of longitudinal deformation of shield tunnel. Segment longitudinal joint is composed of a ring joint surface and a connecting bolt. It is assumed that the connecting bolts bear all the tension in the tensile zone of the ring joint surface, while the segment concrete bears all the pressure in the compression zone. The pretightening force of bolts is not considered temporarily. As shown in Figure7, contact surfaces are set between pipe rings to simulate the state nonlinear behavior of ring joints and connecting bolts. Normal hard contact is used in the compression zone and cohesive model is used in the tension zone. According to the principle that adjacent pipe rings receive the same resultant force with the same opening width amount under unidirectional tensile condition, the bolt stiffness is equivalent to a continuous and uniform distribution along the segment contact area. Physical and mechanical parameters of connecting bolts are listed in Table 2.Figure 7
Deformation characteristics and relation between the forces and deformations between adjacent pipe rings.Table 2
Ring joint bolt parameters of shield tunnel segments.
Diametersd0 (mm)Lengthl0 (mm)Elastic modulusEs (MPa)Tensile yield stressfy (MPa)Tensile ultimate stressfm (MPa)Elastoplastic stiffness ratioα367502.06E56408000.01As the opening width of the ring joints increases, the connecting bolt will yield and the stiffness of the ring joint will decay. Maximum separation criterion is applied for describing the stiffness damage of the cohesive material; that is, the stiffness begins to decay when the normal or tangential deformation reaches the corresponding limit. Meanwhile, the stiffness decay is realized by the damage indexD:(6)K2=1−DK1,where K1 is the stiffness of bolts in the elastic deformation level; K2 is the stiffness of the bolt after yield; and D is the damage index, which changes from “0” to “1.” In ABAQUS, the damage index D is defined by the indirect method of elastic ultimate deformation δy and ultimate deformation δu.In the cohesive model, the elastic ultimate deformationδy is obtained when the bolt reaches the yield stress. The ultimate deformation δu indicates that the bolt stress begins to decrease when the opening width between pipe rings reaches the ultimate displacement:(7)δy=NyK1=nAsfynEsAs/l0=fyl0Es,δu=δy+Nm−NyK2=δy+nAsfm−nAsfynαEsAs/l0=fyEs+fm−fyαEs⋅l0,where Ny is the elastic ultimate tensile force of bolts; Nu is the ultimate tensile force of bolts, and l0 is the length of bolts. For the calculation of cohesive model parameters, the reader is referred to [24], and the relevant parameters are listed in Table 3. 42 bolts are set between rings, and 2 bolts are set between the longitudinal segments. The waterproof limit is 15 mm.Table 3
Parameters of surface-based cohesive behavior in finite element model.
Knn (E8N/m3)Kss (E9N/m3)Ktt (E9N/m3)Elastic deformation limitδy (mm)Ultimate displacementδu (mm)4.4818(4.657)4.4818(4.657)4.4818(4.657)2.360.6


## 3.3.3. SMA Shape Memory Alloy Flexible Joint

Shape memory alloy is a good shock absorber because of its superelastic properties. The commonly used memory alloy is composed of 55.9% Ni and 44.1% Ti, with a recoverable strain of 8% and an ultimate strain of 17%. In addition, memory alloys can dissipate the energy by forming hysteresis curves within recoverable strains. In the past two decades, shape memory alloys as a new material have attracted extensive attention from scholars at home and abroad, and many scholars [29, 30] have conducted in-depth studies on the numerical simulation of their constitutive relations. In this paper, a flag-SMA piecewise linear simplified constitutive model is adopted [31]. According to the material properties of Ni-Ti shape memory alloy obtained by Liu et al. [32] and Zhou et al. [33], the superelasticity (the material constitutive is shown in Figure 8) is used to simulate shape memory alloy in ABAQUS, 2017 version. The flexible antiseismic joint is shown in Figure 9. A total of 40 sets of energy dissipation flexible joint with a diameter of 25 mm are set in a single ring. The coupling connection of the flexible joint elements with the shaft end wall and the shield tunnel is adopted (i.e., the degree of freedom of the node is consistent), and the interaction setting of the surrounding rock is consistent with that of the shield tunnel. According to the principle that adjacent pipe rings receive the same force with the same opening width, the stiffness of memory alloy material is equivalent to a continuous and uniform distribution along the segment contact area, as described by the following equation:(8)K′=nAsAcK,where K′ is the material stiffness after equivalence; As is the cross-sectional area of a single energy dissipation-damping joint; Ac is the cross-sectional area of pipe rings; n is the number of energy dissipation-damping joints; and K is the material stiffness before equivalence. The parameters of superelasticity are listed in Table 4.EM is Young’s modulus of martensite; γM is Poisson’s ratio of martensite; εL is the uniaxial transformation strain; σtLs is the stress at which the transformation begins during loading in tension; σtLE is the stress at which the transformation ends during loading in tension; σtUS is the stress at which the reverse transformation begins during loading in tension; σtUE is the stress at which the reverse transformation ends during loading in tension; σcLS is the stress at which the transformation begins during loading in compression, as a positive value; T0 is the reference temperature; (δσ/δT)L is the slope of the stress versus temperature curve for loading; and (δσ/δT)U is the slope of the stress versus temperature curve for unloading.Figure 8
Lag-SMA simplified constitutive model.Figure 9
Simulation of the SMA shape memory alloy flexible joint.Table 4
The parameters of superelasticity (unit: MPa).
EMγMεLσtLsσtLEσtUSσtUEσcLST0(δσ/δT)L(δσ/δT)U19.70.330.0670.151.551.40.0370.15150.010.006


## 3.3.4. Concrete

Concrete damage model (CDP) is used to describe the irreversible cumulative plastic deformation and stiffness loss of concrete under cyclic loading [34]. The CDP model needs to define plastic parameters and damage parameters. Plastic parameters describe the change form of the material yield surface, while damage parameters describe the loading and unloading characteristics of the material from tensile and compressive stress-strain relations and damage index. The concrete strength grade of the tunnel and shaft is C60, and the parameters of the CDP are listed in Table 5. According to the proportional strain method, the damage index of tension and compression is obtained [35, 36].Table 5
Parameters of the concrete from the plasticity test.
Elastic modulusEc (GPa)Poisson’s ratioμInitial yield stressσco (MPa)Compression strengthσcu (MPa)Tensile strengthσtu (MPa)Dilatation angle Ψ (°)Eccentricityδfbo/fcoKcViscosity coefficientu360.233.2747.573.8838.00.11.162/30.0005fbo/fco is the ratio of biaxial to uniaxial ultimate compressive strength; Kc is the ratio of biaxial to uniaxial ultimate compressive strength; and u is the ratio of the second stress invariant on the tensile and compressive meridian planes.


## 3.4. Selection of Bedrock Input Motion

According to the seismic safety assessment results of Shantou Bay Tunnel site, the exceeded probabilities of 63%, 10%, and 2% of Peak bedrock acceleration are approximately 0.1g, 0.2 g, and 0.4 g, respectively, where “g” is the acceleration of gravity. Due to the lack of strong earthquake records on this site, on the basis of historical seismic data, the bedrock acceleration records of two strong earthquakes in the foreign seismic network are selected as the bedrock input ground motion. The detailed information is listed in Table 6. The acceleration time history of seismic records and its corresponding Fourier spectrum is shown in Figure 10. The Darfield record focuses on low frequency, while the Iwate record develops high-frequency components. There are significant differences in the spectral characteristics of the acceleration records between the two records. In order to consider the influence of the input ground motion intensity, the PGA (peak ground acceleration) of the horizontal direction (vertical to the longitudinal axis of the tunnel) recorded in the original is adjusted to three earthquake levels: 0.1 g (low-level earthquake (LLE)), 0.2 g (moderate-level earthquake (MLE)), and 0.4 g (high-level earthquake (HLE)), and the ratio of horizontal PGA to vertical PGA is set as 1 : 0.65.Table 6
Details of original ground motion records used for site response analysis.
Earthquake eventStation site conditionMagnitudeDistance from the epicenter (km)Ground motion componentPGA (g)Duration D5-95 (s)Mean frequencyfp (Hz)DarfieldPage Road Pumping (PRP) StationMs 7.153.53NS0.23022.290.34UD0.31815.0811.22IwateIWTH27Ms 7.256.85NS0.34720.448.54UD0.12822.414.21Figure 10
Acceleration time history and Fourier spectra of the input bedrock motions.


## 3.5. Monitoring Point Positions

The displacement time history of the submodel should be extracted from the subboundary of the global model, choosing three monitoring points at the subboundary: Point 1 to Point 2 (as shown in Figure4(b)), which are selected to observe its displacement time-history response. The opening width and dislocation are important indexes to evaluate the seismic performance of shield tunnel. Excessive deformation will lead to bolts failure, and water leakage will occur in serious cases. In this paper, the ring joint connecting the first ring of shield tunnel and the wall of shaft in the global model and submodel is referred to as tunnel-shaft joint for short. Based on symmetry, channel “A” is taken as the research object, and the opening width at the tunnel-shaft joint is analyzed from the positions of 8 measuring points shown in Figure 11(a). The peak opening width is defined as the maximum absolute value of the difference between the displacement time history of the right end of the left pipe ring and the corresponding point displacement time history of the left end of the right pipe ring. The first ring of shield tunnel in the submodel is referred to as a tunnel-shaft head ring. The longitudinal seam deformation analysis is conducted on the positions of 10 measuring points between the first ring joints shown in Figure 11(b). The peak deformation is defined as the maximum absolute value of the contact deformation (Copen) time history of adjacent measuring points. The shield tunnel spatial end structure is referred to as a tunnel-shaft junction for short.Figure 11
Monitoring positions at the joint of tunnel-shaft junction.
(a)(b)


## 4. Displacement Time History of the Monitoring Points in Subboundary

Node-based submodel technique is an analysis technique that uses the global model’s node displacement field to interpolate to the boundary nodes of the submodel. Figure12 shows the displacement time history of the monitoring nodes of the global model and the submodel subboundary. It can be seen from the figure that, under the same seismic motions, the trend of displacement time-history curve at the same monitoring point is consistent. This indicates that the seismic displacement response of the global model can be transferred to the corresponding nodes of the submodel structure, which means the global model can drive the submodel to perform seismic response analysis effectively.Figure 12
The displacement time history of the monitoring points in subboundary.


## 5. Joint Seam Deformation Analysis

### 5.1. Peak Opening Width of Tunnel-Shaft Joint

Figure13 shows the peak opening width of the tunnel-shaft joint calculated by the global model and the submodel. Figure 14 shows the time-history curve of the outside spandrel ⑧ joint seam opening width. It can be seen from the figure that the phenomenon of opening and closing between the pipe rings occurred repeatedly under the seismic motion, and the opening width of the ring joints seam is always positive, which indicates that the ring segments are not intruding into each other due to compressional crushing. The peak opening width is symmetric on the cross section of the tunnel-shaft joint under the LLE of different frequency spectrum. With the increase of input ground motion intensity, the peak opening width of tunnel-shaft joint at monitoring points increases, which shows the annular asymmetry characteristics, and the maximum values occur at the outside spandrel. This is owing to the tunnel-shaft junction being affected by the deformation of the surrounding soil; meanwhile, the strong nonlinear characteristics of marine soil increase more strongly along with the seismic motion intensity. When inputting the HLE-Darfield seismic motion, the peak opening width of the outside spandrel is 17.94 mm in the global model, while it is 19.41 mm in the submodel, which has potential safety hazards for exceeding the waterproof limit. On the whole, the peak opening width of ring joints mostly occurs at the vault or the outside spandrel. Two aspects are responsible for the phenomenon. Firstly, when the shear wave is perpendicular or approximately perpendicular to the tunnel axis, the circular tunnel will undergo elliptical deformation [37]. Secondly, the overlying soil of the tunnel is soft, and the stiffness difference between the tunnel and the soil layer is great, which leads the relative slip between the vault and the soil layer to be greater under the seismic motion. It is suggested that special antiseismic measures should be taken to reduce the opening width of the vault and the outside spandrel during the design of the submarine tunnel.Figure 13
Peak opening width of tunnel-shaft joint.Figure 14
Time histories of opening width at the outer spandrel of tunnel-shaft joint.Compared with the Iwate record with abundant high-frequency components, Darfield record with rich low-frequency components causes a larger opening width at the ring joints, which is closely related to the nonlinearity of soil. The phenomenon of high-frequency filtering and low-frequency amplification will occur, when the seismic wave propagates through the soft seabed soil layer. However, the structure depends on the displacement of the surrounding soil layer, and the dynamic characteristics of the tunnel itself have little influence on it. Therefore, the opening width of the ring joints is particularly significant under the low-frequency seismic motion.As clearly indicated in Figures13 and 14, the opening width of the submodel is generally larger than that of the global model. This is due to the damage of underground structures in soft soil, which is mainly controlled by the deformation of surrounding soils; the acceleration and velocity response of the structure itself has no apparent effect on it. Meanwhile, the artificial boundary of the global model has a constraint effect on the structure, while the submodel is surrounded by the driving subboundary, lacking corresponding boundary transmission effect, which leads to greater openings.


### 5.2. Peak Deformation of Tunnel-Shaft Head Ring

Aiming at the deformation concentrated part in the spatial end structure of shield tunnel, a mesh refinement submodel structure was established. At the same time, considering longitudinal seam splice, the homogeneous rings are replaced by segment rings. Figure15 shows the peak deformation of longitudinal joints between all segments in the tunnel-shaft head ring of the submodel. As shown in the figure, the deformation characteristics of the head ring are basically the same as the joint shown in Figure 13. Both of them show the deformation of the outside spandrel and the vault is larger. As the direction of input seismic motion is the two tangential directions along the cross section of the tunnel, the ring segments are mostly squeezed, causing the mutual dislocation to be relatively small. Therefore, the deformation of the tunnel-shaft head ring is less than the opening width of the tunnel-shaft joint.Figure 15
Peak deformation of tunnel-shaft head ring.


## 5.1. Peak Opening Width of Tunnel-Shaft Joint

Figure13 shows the peak opening width of the tunnel-shaft joint calculated by the global model and the submodel. Figure 14 shows the time-history curve of the outside spandrel ⑧ joint seam opening width. It can be seen from the figure that the phenomenon of opening and closing between the pipe rings occurred repeatedly under the seismic motion, and the opening width of the ring joints seam is always positive, which indicates that the ring segments are not intruding into each other due to compressional crushing. The peak opening width is symmetric on the cross section of the tunnel-shaft joint under the LLE of different frequency spectrum. With the increase of input ground motion intensity, the peak opening width of tunnel-shaft joint at monitoring points increases, which shows the annular asymmetry characteristics, and the maximum values occur at the outside spandrel. This is owing to the tunnel-shaft junction being affected by the deformation of the surrounding soil; meanwhile, the strong nonlinear characteristics of marine soil increase more strongly along with the seismic motion intensity. When inputting the HLE-Darfield seismic motion, the peak opening width of the outside spandrel is 17.94 mm in the global model, while it is 19.41 mm in the submodel, which has potential safety hazards for exceeding the waterproof limit. On the whole, the peak opening width of ring joints mostly occurs at the vault or the outside spandrel. Two aspects are responsible for the phenomenon. Firstly, when the shear wave is perpendicular or approximately perpendicular to the tunnel axis, the circular tunnel will undergo elliptical deformation [37]. Secondly, the overlying soil of the tunnel is soft, and the stiffness difference between the tunnel and the soil layer is great, which leads the relative slip between the vault and the soil layer to be greater under the seismic motion. It is suggested that special antiseismic measures should be taken to reduce the opening width of the vault and the outside spandrel during the design of the submarine tunnel.Figure 13
Peak opening width of tunnel-shaft joint.Figure 14
Time histories of opening width at the outer spandrel of tunnel-shaft joint.Compared with the Iwate record with abundant high-frequency components, Darfield record with rich low-frequency components causes a larger opening width at the ring joints, which is closely related to the nonlinearity of soil. The phenomenon of high-frequency filtering and low-frequency amplification will occur, when the seismic wave propagates through the soft seabed soil layer. However, the structure depends on the displacement of the surrounding soil layer, and the dynamic characteristics of the tunnel itself have little influence on it. Therefore, the opening width of the ring joints is particularly significant under the low-frequency seismic motion.As clearly indicated in Figures13 and 14, the opening width of the submodel is generally larger than that of the global model. This is due to the damage of underground structures in soft soil, which is mainly controlled by the deformation of surrounding soils; the acceleration and velocity response of the structure itself has no apparent effect on it. Meanwhile, the artificial boundary of the global model has a constraint effect on the structure, while the submodel is surrounded by the driving subboundary, lacking corresponding boundary transmission effect, which leads to greater openings.


## 5.2. Peak Deformation of Tunnel-Shaft Head Ring

Aiming at the deformation concentrated part in the spatial end structure of shield tunnel, a mesh refinement submodel structure was established. At the same time, considering longitudinal seam splice, the homogeneous rings are replaced by segment rings. Figure15 shows the peak deformation of longitudinal joints between all segments in the tunnel-shaft head ring of the submodel. As shown in the figure, the deformation characteristics of the head ring are basically the same as the joint shown in Figure 13. Both of them show the deformation of the outside spandrel and the vault is larger. As the direction of input seismic motion is the two tangential directions along the cross section of the tunnel, the ring segments are mostly squeezed, causing the mutual dislocation to be relatively small. Therefore, the deformation of the tunnel-shaft head ring is less than the opening width of the tunnel-shaft joint.Figure 15
Peak deformation of tunnel-shaft head ring.


## 6. Stress and Seismic Damage Analysis of Spatial End Structural

### 6.1. Peak Mises Stress of Tunnel-Shaft Joint

The peak Mises stress (Mises stress time-history peak) of monitoring points in tunnel-shaft joint is calculated by the global model and the submodel, as shown in Figure16. As can be seen, the peak Mises stress distribution on the ring cross section is “X” shape, which means the peak Mises stress in the conjugate with the vertical axis of 45° is larger. Furthermore, the max Mises stress is located at the outer arch foot of the cross section mostly and shows symmetry with the max peak opening width. Because of the two-direction seismic motion input, the tunnel-shaft joint will undergo annular torsion. When the spandrel has a larger opening width, the deformation of the pipe ring will be transferred to the other end, resulting in a larger compression at the arch foot. Due to the fact that the compressive stiffness of concrete pipe ring is much larger than the tensile stiffness of connecting bolt, the peak stress of arch foot is larger than the peak stress of spandrel. Figure 17 shows the Mises stress time history of the outer arch foot ⑥ at the tunnel-shaft joint; the monitoring point is in a state of low stress level before 20 s. However, with the enhancement of seismic motion, the stress has a period of sharp increase. After that, the stress decreases slowly due to structural damage and stiffness degrades.Figure 16
The peak Mises stress of the tunnel-shaft joint.Figure 17
Time histories of Mises stress at the outer arch foot of the tunnel-shaft joint.


### 6.2. Earthquake Damage of Tunnel-Shaft Junction

The damage index between “0” (blue: no damage state) and “1” (red: complete damage state) is used to describe the seismic damage degree at tunnel-shaft junction. The distribution cloud diagram of seismic damage index at the spatial end of the shield tunnel is shown in Figure18. The damage intensity and damage range of the global model and the submodel increase with the increase of the seismic motion intensity. Due to the fact that the geometry size and stiffness from shaft structure to shield tunnel change significantly, leading to large longitudinal tension and compression between shaft wall and shield tunnel, when the seismic motion of bedrock with different intensity and spectrum characteristics is excited, the stress concentration is more serious, the damage is more obvious, and the damage evolution of shaft wall is more intense than that of shield tunnel. In the case of LLE, only the shaft wall is partially damaged, and the tunnel-shaft junction is still in elastic working state. In the case of MLE, the damage of the shaft walls at the junction increases slightly under the excitation of the Iwate record with rich high-frequency components. Under the excitation of Darfield record with rich low-frequency components, the damage degree of tunnel within the range of 2∼4 rings and shaft increase significantly, and the junction has entered into the elastic-plastic working state. In the case of HLE, under the excitation of Darfield recorded, the concrete is seriously damaged and the bearing capacity is completely lost within the range of 4∼6 tunnel rings and the shaft. It can be basically determined that the tunnel-shaft junction has overall destruction. According to the damage degree and distribution of the tunnel-shaft junction, the damage first occurs on the wall of the conjugate with the vertical axis of 45°, with the damage distribution on the cross section present as “X” shape, which is consistent with the peak Mises stress distribution on the cross section in Figure 14. Therefore, some measures should be taken in the structure design to improve the antiseismic performance of the junction.Figure 18
Cloud diagram of the damage indexes and time-history curve of the damage factor of tunnel-shaft junction.From the time-history curve of damage index, the degree and speed of structural damage increase sharply with the increase of seismic motion intensity. Taking Darfield record as an example, the damage index gradually increases and tends to be stable between 17 s and 33 s when under the excitation of LLE. The damage index increases rapidly between 19 s and 30 s during MLE. In the case of HLE, the damage index increases rapidly between 19 s and 23 s and reaches the state of complete destruction. The final damage distribution and the evolution rule, calculated by the global model and the submodel, is consistent, which indicates that the submodel technique can be used to analyze the seismic damage of the locally refined structures effectively.


## 6.1. Peak Mises Stress of Tunnel-Shaft Joint

The peak Mises stress (Mises stress time-history peak) of monitoring points in tunnel-shaft joint is calculated by the global model and the submodel, as shown in Figure16. As can be seen, the peak Mises stress distribution on the ring cross section is “X” shape, which means the peak Mises stress in the conjugate with the vertical axis of 45° is larger. Furthermore, the max Mises stress is located at the outer arch foot of the cross section mostly and shows symmetry with the max peak opening width. Because of the two-direction seismic motion input, the tunnel-shaft joint will undergo annular torsion. When the spandrel has a larger opening width, the deformation of the pipe ring will be transferred to the other end, resulting in a larger compression at the arch foot. Due to the fact that the compressive stiffness of concrete pipe ring is much larger than the tensile stiffness of connecting bolt, the peak stress of arch foot is larger than the peak stress of spandrel. Figure 17 shows the Mises stress time history of the outer arch foot ⑥ at the tunnel-shaft joint; the monitoring point is in a state of low stress level before 20 s. However, with the enhancement of seismic motion, the stress has a period of sharp increase. After that, the stress decreases slowly due to structural damage and stiffness degrades.Figure 16
The peak Mises stress of the tunnel-shaft joint.Figure 17
Time histories of Mises stress at the outer arch foot of the tunnel-shaft joint.


## 6.2. Earthquake Damage of Tunnel-Shaft Junction

The damage index between “0” (blue: no damage state) and “1” (red: complete damage state) is used to describe the seismic damage degree at tunnel-shaft junction. The distribution cloud diagram of seismic damage index at the spatial end of the shield tunnel is shown in Figure18. The damage intensity and damage range of the global model and the submodel increase with the increase of the seismic motion intensity. Due to the fact that the geometry size and stiffness from shaft structure to shield tunnel change significantly, leading to large longitudinal tension and compression between shaft wall and shield tunnel, when the seismic motion of bedrock with different intensity and spectrum characteristics is excited, the stress concentration is more serious, the damage is more obvious, and the damage evolution of shaft wall is more intense than that of shield tunnel. In the case of LLE, only the shaft wall is partially damaged, and the tunnel-shaft junction is still in elastic working state. In the case of MLE, the damage of the shaft walls at the junction increases slightly under the excitation of the Iwate record with rich high-frequency components. Under the excitation of Darfield record with rich low-frequency components, the damage degree of tunnel within the range of 2∼4 rings and shaft increase significantly, and the junction has entered into the elastic-plastic working state. In the case of HLE, under the excitation of Darfield recorded, the concrete is seriously damaged and the bearing capacity is completely lost within the range of 4∼6 tunnel rings and the shaft. It can be basically determined that the tunnel-shaft junction has overall destruction. According to the damage degree and distribution of the tunnel-shaft junction, the damage first occurs on the wall of the conjugate with the vertical axis of 45°, with the damage distribution on the cross section present as “X” shape, which is consistent with the peak Mises stress distribution on the cross section in Figure 14. Therefore, some measures should be taken in the structure design to improve the antiseismic performance of the junction.Figure 18
Cloud diagram of the damage indexes and time-history curve of the damage factor of tunnel-shaft junction.From the time-history curve of damage index, the degree and speed of structural damage increase sharply with the increase of seismic motion intensity. Taking Darfield record as an example, the damage index gradually increases and tends to be stable between 17 s and 33 s when under the excitation of LLE. The damage index increases rapidly between 19 s and 30 s during MLE. In the case of HLE, the damage index increases rapidly between 19 s and 23 s and reaches the state of complete destruction. The final damage distribution and the evolution rule, calculated by the global model and the submodel, is consistent, which indicates that the submodel technique can be used to analyze the seismic damage of the locally refined structures effectively.


## 7. Shock Absorption Measure Analysis

### 7.1. Peak Opening Width of Tunnel-Shaft Joint under Shock Absorption Measure

From the above analysis cases, it can be seen that the structure has the most intense seismic response under the Darfield record with excitation of LLE, so this case is selected as the standard case. Figure19 shows the distribution of peak opening width of tunnel-shaft joint along the cross section with SMA flexible shock absorption measure when Darfield record of LLE is input. The max peak opening width of tunnel-shaft joint at monitoring points is greatly increased; when adding SMA shape memory alloy connection, the global model is 124.7% of that under the standard case, while the submodel is 134.5% of that under the standard case. As the input external excitation is shear two directions, the peak opening width of tunnel-shaft joint at the outside spandrel increases sharply; the global model and the submodel reach 1.49 and 1.74 times of the waterproofing limit, respectively. Moreover, the peak opening width at several monitoring points exceeds the waterproofing limit. Because of the additional SMA shape memory alloy connection, the energy dissipation-damping characteristic of flexible joint leads the shaft and the tunnel into relatively independent parts, which can induce structural deformation caused by seismic wave propagation to the preset flexible joint, making the flexible joint deformation obviously and increasing the peak opening width sharply. Therefore, the flexible joint becomes the weak part of the structure in seismic resistance and should be able to meet the expected seismic deformation requirement.Figure 19
Peak opening width of tunnel-shaft joint under shock absorption measure.Figure19 also shows the distribution of residual peak opening width of tunnel-shaft joint; the residual openings of each measuring point are reduced significantly after adding the SMA shape memory alloy connection. Owing to the reset characteristics of the shape memory alloy connection, the residual opening width can be effectively controlled after seismic excitation is over, ensuring the recovery of the joint.


### 7.2. Earthquake Damage and Stress of Tunnel-Shaft Junction under Shock Absorption Measure

Figure20 shows, after adding the SMA flexible damping connection, the seismic damage cloud diagram of the spatial end structure and the peak Mises stress of the tunnel-shaft joint at measuring points on the joint cross section, when Darfield record of LLE is input. Compared with Figure 18, the damage range and degree of shield tunnel are reduced significantly, indicating that flexible connection can enhance the structural differential displacement and stress concentration caused by the stiffness mutation of shaft and shield tunnel effectively. However, for all cases, the seismic damage degree of the shaft end wall is still relatively severe. The existence of shield tunnel weakens the integral seismic performance of the shaft end wall, which is the place that should be paid attention to in the seismic design.Figure 20
Cloud diagram of the earthquake damage and the max Mises stress of the tunnel-shaft junction under shock absorption measure.The distribution of peak Mises stress at measuring points of tunnel-shaft joint cross section is shown in Figure20. The stress concentration at the spandrel and arch foot of tunnel-shaft joint is significantly greater than that at other positions. Combining Figures 11 and 14, the distribution of peak Mises stress is consistent with the distribution of peak opening width; both of them have a strong correlation, which proves the cohesive model in this paper can effectively simulate the spatial deformation characteristics of bolts. The addition of flexible joint can effectively reduce the seismic stress and stress concentration of the junction. Compared with the standard case, the maximum reductions of the peak Mises stress at tunnel-shaft joint of the global model and the sub model are 52.6% and 57.4%, respectively.


## 7.1. Peak Opening Width of Tunnel-Shaft Joint under Shock Absorption Measure

From the above analysis cases, it can be seen that the structure has the most intense seismic response under the Darfield record with excitation of LLE, so this case is selected as the standard case. Figure19 shows the distribution of peak opening width of tunnel-shaft joint along the cross section with SMA flexible shock absorption measure when Darfield record of LLE is input. The max peak opening width of tunnel-shaft joint at monitoring points is greatly increased; when adding SMA shape memory alloy connection, the global model is 124.7% of that under the standard case, while the submodel is 134.5% of that under the standard case. As the input external excitation is shear two directions, the peak opening width of tunnel-shaft joint at the outside spandrel increases sharply; the global model and the submodel reach 1.49 and 1.74 times of the waterproofing limit, respectively. Moreover, the peak opening width at several monitoring points exceeds the waterproofing limit. Because of the additional SMA shape memory alloy connection, the energy dissipation-damping characteristic of flexible joint leads the shaft and the tunnel into relatively independent parts, which can induce structural deformation caused by seismic wave propagation to the preset flexible joint, making the flexible joint deformation obviously and increasing the peak opening width sharply. Therefore, the flexible joint becomes the weak part of the structure in seismic resistance and should be able to meet the expected seismic deformation requirement.Figure 19
Peak opening width of tunnel-shaft joint under shock absorption measure.Figure19 also shows the distribution of residual peak opening width of tunnel-shaft joint; the residual openings of each measuring point are reduced significantly after adding the SMA shape memory alloy connection. Owing to the reset characteristics of the shape memory alloy connection, the residual opening width can be effectively controlled after seismic excitation is over, ensuring the recovery of the joint.


## 7.2. Earthquake Damage and Stress of Tunnel-Shaft Junction under Shock Absorption Measure

Figure20 shows, after adding the SMA flexible damping connection, the seismic damage cloud diagram of the spatial end structure and the peak Mises stress of the tunnel-shaft joint at measuring points on the joint cross section, when Darfield record of LLE is input. Compared with Figure 18, the damage range and degree of shield tunnel are reduced significantly, indicating that flexible connection can enhance the structural differential displacement and stress concentration caused by the stiffness mutation of shaft and shield tunnel effectively. However, for all cases, the seismic damage degree of the shaft end wall is still relatively severe. The existence of shield tunnel weakens the integral seismic performance of the shaft end wall, which is the place that should be paid attention to in the seismic design.Figure 20
Cloud diagram of the earthquake damage and the max Mises stress of the tunnel-shaft junction under shock absorption measure.The distribution of peak Mises stress at measuring points of tunnel-shaft joint cross section is shown in Figure20. The stress concentration at the spandrel and arch foot of tunnel-shaft joint is significantly greater than that at other positions. Combining Figures 11 and 14, the distribution of peak Mises stress is consistent with the distribution of peak opening width; both of them have a strong correlation, which proves the cohesive model in this paper can effectively simulate the spatial deformation characteristics of bolts. The addition of flexible joint can effectively reduce the seismic stress and stress concentration of the junction. Compared with the standard case, the maximum reductions of the peak Mises stress at tunnel-shaft joint of the global model and the sub model are 52.6% and 57.4%, respectively.


## 8. Conclusions

Aiming at the weak seismic part of shield tunnel spatial end structure, this paper takes Suai undersea tunnel project as the background, considering dynamic nonlinear characteristics of seabed soil, bolt connection between pipe rings, and SMA shape memory alloy joint. The three-dimensional finite element global model and local refined submodel of the tunnel-shaft junction are established for exploring the influence of input seismic motion characteristics, as well as shock absorption measures. The following conclusions are reached:(1)
The cohesive model can effectively simulate the three-dimensional deformation characteristics between segments of shield tunnel. The peak opening width of tunnel-shaft joint and the peak deformation of tunnel-shaft head ring are the largest at the vault and the outer spandrel of the tunnel-shaft junction. The seismic stress concentration of the shaft wall and the vertical conjugate direction at 45° are significant, and the stress at the arch foot is slightly larger than that at the spandrel. The seismic damage to the shaft wall of this structure is serious.(2)
Because of the characteristics of high-frequency filtering and low-frequency amplification, the seismic response of the tunnel-shaft junction is stronger when the seismic motion is excited with rich low frequency.(3)
Adding the flexible connection between the shaft and the shield tunnel can induce structural deformation caused by seismic wave propagation to the preset flexible joint, reducing the seismic damage range and degree of the tunnel-shaft junction. However, the seismic safety of the flexible joint is sacrificed for reducing the seismic stress of the shield tunnel. The flexible joint becomes the weak part of the structure, and the design should be able to meet the expected seismic deformation requirement.(4)
The basic characteristics of seismic deformation and stress at the tunnel-shaft junction calculated by the global model are consistent with the submodel, but the calculated value of the submodel is greater than that of the global model.For the selected input seismic motion, the 16-CPU parallel algorithm was used to carry out three-dimensional nonlinear dynamic time-history analysis on the global model, and each case takes about 20 hours, while each refined submodel case takes about 4 hours. If the segment longitudinal seam splicing and refined mesh is considered in the global model, it will take about 70 hours for every single case. Therefore, the submodel technology can effectively improve computational efficiency and obtain the local location refinement results, but there are some safety errors compared with the original global model.



---

*Source: 8441325-2021-09-15.xml*