# Effects of Spatially Varying Seismic Ground Motions and Incident Angles on Behavior of Long Tunnels

Authors: Yundong Zhou; Yongxin Wu; Ziheng Shangguan; Zhanbin Wang
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

Seismic behavior of long circle tunnels is significantly influenced by the nature of input motion. This study, based on the 3D finite-element method (FEM), evaluates the effects of spatially varying seismic ground motions and uniform input seismic ground motions and their incident angles on the diameter strain rate and tensive/compressive principal stresses under different strata. It is found that (1) the spatially varying seismic ground motions induced larger diameter strain rate (radially deformation) than the uniform input seismic motion, (2) the spatially varying seismic ground motions had an asymmetric effect on the radial strain rate distributions, and (3) the rising incident angles changed the pure shear stress state into a complex stress state for tunnels under specified input motion.


---

## Body

## 1. Introduction

Tunnels, as an integral part of the infrastructure of modern society, would suffer damages when subjected to dynamic loadings during earthquake activity, for example, in Kobe earthquake (1995) [1], Chi-Chi earthquake (1999) [2], Duzce earthquake (1999) [3], Mid-Niigata Prefecture earthquake (2004) [4], Wenchuan earthquake (2008) [5], and Tohoku earthquake (2011) [6]. Nevertheless, the number of large tunnels and underground structures had grown significantly in recent decades. Therefore, seismic evaluation of tunnels in seismically active areas is critical during engineering design.Seismic behavior of tunnels has been widely studied by many researchers [7–11], and these researches have concentrated mostly on 2D analysis. For the analysis of axial and bending deformations of tunnels, it is most appropriate to utilize 3D models.However, for the 3D analysis of seismic behavior of tunnel, soil-tunnel analyses in the past were typically limited to relatively small regions, which made it difficult to fully consider the complex spatial features involved in such large structures [12]. Tunnels often had significant length and could be built on different strata foundations, which made the seismic analysis of tunnels a complex problem and was usually evaluated under idealized conditions by using numerical methods, such as the finite-element method (FEM).For the seismic analysis of complex stress distribution and deformation of long-distance tunnels, it is often more reliable to adopt three-dimensional (3D) methods. With the rapid development of science and technology, it is now possible to use high-performance computers to conduct large-scale 3D FEM seismic analysis for tunnels. There are three kinds of deformation, such as axial compression/extension, longitudinal bending, and ovalling/racking, occurring in tunnels during earthquake [10]. Particularly, the cross-sectional distortion of the tunnel can be related to seismic waves propagating along the tunnel.Yu et al. [13, 14] presented a multiscale 3D FEM analysis of long tunnels under seismic loads where the mechanical characteristics of the tunnel segments and joints under artificial or recorded earthquake loads were evaluated in detail. The model in this paper would take into account not only the motion distribution with time but also the spatial variability (incoherency effect), the wave-passage effect, and the site-response effect.This paper, based on an earlier report by the authors [15], attempts to develop a new model for seismic analysis of long tunnels with multisupport excitations, which properly accounts for the spatial variability. The main focus of this paper is to assess the influence of multisupport input earthquake waves and uniform input earthquake waves as well as their incident angles on the diameter strain rate and tensive/compressive principal stresses under different strata. A full-scale 3D finite-element model is built comprising geological data, tunnel geometry, and so on to farthest simulate actual situation.


## 2. Simulation of Seismic Ground Motions

### 2.1. Simulation of Uniform Seismic Ground Motion

Consider a zero mean Gaussian stationary seismic ground motion with a target autospectral densitySω. The seismic ground motions can be generated through the following expression:(1)xt=2∑l=1NSωlΔω⁡cosωlt+φl,where N is the number of frequency intervals, Δω=ωu/N is frequency increment with ωu as the cutoff frequency, ωl=lΔω, and the φl’s are statistically independent random phase angles uniformly distributed between (0, 2π]. Equation (1) is valid if there is an upper cutoff frequency ωu above which the contribution of the power spectral density (PSD) to the simulations is negligible for practical purposes.


### 2.2. Simulation of Spatially Varying Seismic Ground Motions

The variations in ground motion are caused by the following four sources: (1) the “incoherence effect,” (2) the “wave-passage effect,” (3) the “site-response effect,” and (4) the “attenuation effect.” The spectral representation method is one of the most widely used methods in simulating the spatially varying seismic ground motions.In practical application, spatially correlated ground motions can be considered as a one-dimensional,n-variety (1D-nV) stochastic vector process Xt with components xjt⋅j=1,2,…,n. Based on the spectral representation method, the jth component of the ground motions can be generated by [13](2)xjt=2∑k=1n∑l=1NHjkωlΔω⁡cosωlt−ϕjkωl+φkl,ωl=lΔωl=0,1,K,…,N−1,Δω=ωuN,where ωu is the upper cutoff frequency beyond which elements of the power spectral can be assumed to be zero for either mathematical or physical reasons and φjk are independent random phase angles uniformly distributed over (0, 2π]. Hjkω and ϕjkω are the modulus and phase parts of Hjkω, respectively, which can be obtained by the root decomposition of power spectral density matrix as follows:(3)Sω=HωHTω,where the superscript T denotes conjugate transpose.The power spectral density matrix is given as(4)Sω=S11ωS12ω⋯S1mωS21ωS22ω⋯S2mω⋮⋮⋱⋮Sm1ωSm2ω⋯Smmω,where Sjjω is the autopower spectral density function and Sjkω,t is the cross-power spectral density function, which can be expressed as(5)Sjkω,t=γjkωSjjωSkkω⁡eiθjkω,where γjkω is the lagged coherence function representing “incoherence” effect and θjkω is composed of wave passage.However, the ground motions simulated by the above method are stationary, while the actual seismic records are nonstationary. Therefore, to obtain nonstationary seismic ground motion, the way of multiplying an envelope function is applied. The envelope function is as follows:(6)ft=t/t12,t<t11,t1≤t<t2e−ct−t2,t≥t2,where t1,t2, and c are three parameters describing the shape of the envelope function. In this study, they are set to be t1=6, t2=10, and c=0.5, respectively.


## 2.1. Simulation of Uniform Seismic Ground Motion

Consider a zero mean Gaussian stationary seismic ground motion with a target autospectral densitySω. The seismic ground motions can be generated through the following expression:(1)xt=2∑l=1NSωlΔω⁡cosωlt+φl,where N is the number of frequency intervals, Δω=ωu/N is frequency increment with ωu as the cutoff frequency, ωl=lΔω, and the φl’s are statistically independent random phase angles uniformly distributed between (0, 2π]. Equation (1) is valid if there is an upper cutoff frequency ωu above which the contribution of the power spectral density (PSD) to the simulations is negligible for practical purposes.


## 2.2. Simulation of Spatially Varying Seismic Ground Motions

The variations in ground motion are caused by the following four sources: (1) the “incoherence effect,” (2) the “wave-passage effect,” (3) the “site-response effect,” and (4) the “attenuation effect.” The spectral representation method is one of the most widely used methods in simulating the spatially varying seismic ground motions.In practical application, spatially correlated ground motions can be considered as a one-dimensional,n-variety (1D-nV) stochastic vector process Xt with components xjt⋅j=1,2,…,n. Based on the spectral representation method, the jth component of the ground motions can be generated by [13](2)xjt=2∑k=1n∑l=1NHjkωlΔω⁡cosωlt−ϕjkωl+φkl,ωl=lΔωl=0,1,K,…,N−1,Δω=ωuN,where ωu is the upper cutoff frequency beyond which elements of the power spectral can be assumed to be zero for either mathematical or physical reasons and φjk are independent random phase angles uniformly distributed over (0, 2π]. Hjkω and ϕjkω are the modulus and phase parts of Hjkω, respectively, which can be obtained by the root decomposition of power spectral density matrix as follows:(3)Sω=HωHTω,where the superscript T denotes conjugate transpose.The power spectral density matrix is given as(4)Sω=S11ωS12ω⋯S1mωS21ωS22ω⋯S2mω⋮⋮⋱⋮Sm1ωSm2ω⋯Smmω,where Sjjω is the autopower spectral density function and Sjkω,t is the cross-power spectral density function, which can be expressed as(5)Sjkω,t=γjkωSjjωSkkω⁡eiθjkω,where γjkω is the lagged coherence function representing “incoherence” effect and θjkω is composed of wave passage.However, the ground motions simulated by the above method are stationary, while the actual seismic records are nonstationary. Therefore, to obtain nonstationary seismic ground motion, the way of multiplying an envelope function is applied. The envelope function is as follows:(6)ft=t/t12,t<t11,t1≤t<t2e−ct−t2,t≥t2,where t1,t2, and c are three parameters describing the shape of the envelope function. In this study, they are set to be t1=6, t2=10, and c=0.5, respectively.


## 3. FE Model

The tunnel model adopted in this paper was built with circular appearance, whose outside diameter is 10 m, inside diameter is 9 m, and length is 1000 m. It had a buried depth of 30 m from the tunnel center to soil surface. The simulation setup is shown in Figure1. It consists of a 1000 × 300 × 100 m box with the 1000 m length circular tunnel. The space coordinates were built by taking the length direction as the Z-axis, the width direction as the X-axis, and the height direction as the Y-direction. The soil profile was modeled as four layers of Mohr–Coulomb materials and tunnel as the elastic material. All the property parameters used can be found in Tables 1 and 2. An infinite domain by using the artificial boundary was adopted [12], where both borders of the tunnel were fixed in the Z-direction.Figure 1
Finite-element mesh model adopted in the 3D FE analyses for soil with long tunnels.Table 1
Properties of different soil layers.
Layer number
Young’s modulus (MPa)
Poisson’s ratio
Density (kg/m3)
Cohesion (kPa)
Internal friction angle (°)
Thickness (m)
1
3.029
0.45
1760
10.8
14.2
20.6
2
6.135
0.36
1810
17.7
18.9
23.9
3
6.319
0.47
1870
22
22.9
24.2
4
6.251
0.36
1950
24.9
18.2
31.5Table 2
Properties of the tunnel.
Elastic modulus (kg/m3)
Poisson’s ratio
Unit weight (kN/m3)
3.6e7
0.2
25.7In this work, the constitutive model of the tunnel is an elastic model, which is given as follows:(7)σ=λI1δ+2Gε,in which σ is the stress tensor, ε is the strain tensor, I1 is the first strain invariant, λ is the Lame constant, and G is the shear modulus.The Mohr–Coulomb yield criterion for soil is expressed as follows:(8)m+12maxσ1−σ2+Kσ1+σ2,σ1−σ3+Kσ1+σ3,σ2−σ3+Kσ2+σ3=Syc,where(9)m=SycSyt,K=m−1m+1,and the parameters Syc and Syt are the yield stresses of the material in uniaxial compression and tension, respectively.On the interface, the deformation of the tunnel and soil is in concert and harmony.


## 4. Computation and Analysis

In computation, a regular 3D finite-element model was built based on geotechnical data, tunnel geometry, and so on. The traveling waves velocity motivated from bedrock surface was 500 m/s in the 3D spatially varying seismic ground motion field. The schematic diagram of different incident angles used is shown in Figure2. The simulated spatially varying ground motions are plotted in Figure 3.Figure 2
The schematic diagram at different incident angles.Figure 3
The simulated spatially varying seismic ground motions: (a)x = 0 m, (b) x = 500 m, and (c) x = 1000 m.
(a)
(b)
(c)


### 4.1. Effect on Radial Deformations

The diameter strain rateΔd is defined as(10)Δd=d′−dd,where d and d′ are the tunnel diameters before and after deformation.The positive rate value denotes that the tunnel cross section undergoes tensile deformation in radial direction, and the negative value corresponds to compressive deformation. In this section,dmax is defined as the maximum of Δd in the earthquake time process. dmax in the vertical direction is calculated from displacement time processes at the vertex and nadir of tunnel in the Y-axis direction while the one in the lateral direction is obtained from processes at the right and left sides of the tunnel in the X-axis direction. It was found from the simulation that lateral Δd and vertical Δd were equal and opposite in the direction at any arbitrary time, videlicet, the primary deformation of tunnel under seismic waves was ovalling. In addition, the partial feature results under uniform and spatially varying seismic ground motions, that is, the lateral dmax at Z = 100 m, 300 m, 500 m, 700 m, and 900 m with incident angles of 0°, 30°, 45°, 60°, and 90° can be found in Tables 3 and 4.Table 3
Lateral maximum diameter strain rate of tunnel under uniform seismic motion.
Incident angle
Z = 100 m
Z = 300 m
Z = 500 m
Z = 700 m
Z = 900 m
0°
−0.036%
−0.004%
0
−0.004%
−0.036%
30°
0.087%
0.113%
0.118%
0.113%
0.087%
45°
0.130%
0.158%
0.165%
0.158%
0.130%
60°
0.171%
0.196%
0.204%
0.196%
0.171%
90°
0.214%
0.226%
0.236%
0.226%
0.214%Table 4
Lateral maximum diameter strain rate of tunnel under multisupport seismic motion.
Incident angle
Z = 100 m
Z = 300 m
Z = 500 m
Z = 700 m
Z = 900 m
0°
0.022%
−0.045%
0.033%
0.035%
0.039%
30°
0.075%
−0.125%
−0.169%
−0.169%
−0.098%
45°
0.111%
−0.163%
−0.225%
−0.225%
−0.123%
60°
0.141%
−0.193%
−0.269%
−0.269%
−0.143%
90°
0.170%
−0.208%
−0.293
−0.293%
−0.155%It could be observed that the lateraldmax increases with increasing incident angles for both seismic input forms, but the maximum positions differ from each other. For uniform seismic input simulations, the lateral dmax position is at the middle point (Z = 500 m) of the tunnel and the diameter strain rates are symmetrically distributed on the left and right sides of the tunnel. This can be attributed to the balanced boundary conditions and uniform seismic waves. For spatially varying seismic ground motions, the lateral dmax is located between Z = 700 m and Z = 900 m and exhibits asymmetric features due to the nonuniform waves. The table also implies that the lateral dmax for the whole tunnel under spatially varying seismic ground motions is greater than that under uniform input ground motions.


### 4.2. Effect on Stress Distributions

The maximum principal shear (σtmax) and compressive (σcmax) stresses are commonly used to analyze the effect of seismic input ground motions on tunnels. The peak stresses under different incident angles are illustrated in Table 5.Table 5
Maximum stresses of tunnel body at different incident angles.
Motion form
Principal stress
0°
30°
45°
60°
90°
Uniform
σ
t
max
11.60
10.85
9.24
7.18
7.45
σ
c
max
−11.60
−10.58
−10.02
−8.99
−6.81
Multisupport
σ
t
max
17.74
13.92
11.29
8.94
6.66
σ
c
max
−18.84
−15.68
−13.37
−11.16
−8.49It is indicated that the maximum tensile stress is identical to the maximum compressive stress while the incident angle is 0° for uniform seismic input tunnels. This is owing to the tunnel and soil profile placed in pure shear stress state. When incident angles are equal to 30°, 45°, 60°, and 90°, the values of maximum tensile and compressive stresses are different because of the complex shear and normal stress states induced by seismic load. It can be thus concluded that the earthquake acceleration in the vertical direction changes not only the maximum stress but also the stress characteristics. In addition, the spatially varying seismic ground motion is another factor influencing the stress state according to the different values between the maximum tensile stress and maximum compressive stress at an incident angle of 0° under spatially varying seismic ground motions. Meanwhile, the value of maximum tensile/compressive stresses decreases with the increasing incident angles. The maximum tensile stress nephogram at an incident angle of 0° for both uniform and spatially varying seismic ground motions is, respectively, shown in Figures4 and 5.Figure 4
The maximum tensile stress nephogram at an incident angle of 0° under uniform seismic input motion.Figure 5
The maximum tensile stress nephogram at an incident angle of 0° under multisupport seismic input motion.The red upper portion of the tunnel end (Z = 1000 m) in Figure 4 is the maximum tensile stress, responding to the red lower portion near the tunnel end. It is suggested that the spatially varying seismic ground motion affects the stress distribution.


## 4.1. Effect on Radial Deformations

The diameter strain rateΔd is defined as(10)Δd=d′−dd,where d and d′ are the tunnel diameters before and after deformation.The positive rate value denotes that the tunnel cross section undergoes tensile deformation in radial direction, and the negative value corresponds to compressive deformation. In this section,dmax is defined as the maximum of Δd in the earthquake time process. dmax in the vertical direction is calculated from displacement time processes at the vertex and nadir of tunnel in the Y-axis direction while the one in the lateral direction is obtained from processes at the right and left sides of the tunnel in the X-axis direction. It was found from the simulation that lateral Δd and vertical Δd were equal and opposite in the direction at any arbitrary time, videlicet, the primary deformation of tunnel under seismic waves was ovalling. In addition, the partial feature results under uniform and spatially varying seismic ground motions, that is, the lateral dmax at Z = 100 m, 300 m, 500 m, 700 m, and 900 m with incident angles of 0°, 30°, 45°, 60°, and 90° can be found in Tables 3 and 4.Table 3
Lateral maximum diameter strain rate of tunnel under uniform seismic motion.
Incident angle
Z = 100 m
Z = 300 m
Z = 500 m
Z = 700 m
Z = 900 m
0°
−0.036%
−0.004%
0
−0.004%
−0.036%
30°
0.087%
0.113%
0.118%
0.113%
0.087%
45°
0.130%
0.158%
0.165%
0.158%
0.130%
60°
0.171%
0.196%
0.204%
0.196%
0.171%
90°
0.214%
0.226%
0.236%
0.226%
0.214%Table 4
Lateral maximum diameter strain rate of tunnel under multisupport seismic motion.
Incident angle
Z = 100 m
Z = 300 m
Z = 500 m
Z = 700 m
Z = 900 m
0°
0.022%
−0.045%
0.033%
0.035%
0.039%
30°
0.075%
−0.125%
−0.169%
−0.169%
−0.098%
45°
0.111%
−0.163%
−0.225%
−0.225%
−0.123%
60°
0.141%
−0.193%
−0.269%
−0.269%
−0.143%
90°
0.170%
−0.208%
−0.293
−0.293%
−0.155%It could be observed that the lateraldmax increases with increasing incident angles for both seismic input forms, but the maximum positions differ from each other. For uniform seismic input simulations, the lateral dmax position is at the middle point (Z = 500 m) of the tunnel and the diameter strain rates are symmetrically distributed on the left and right sides of the tunnel. This can be attributed to the balanced boundary conditions and uniform seismic waves. For spatially varying seismic ground motions, the lateral dmax is located between Z = 700 m and Z = 900 m and exhibits asymmetric features due to the nonuniform waves. The table also implies that the lateral dmax for the whole tunnel under spatially varying seismic ground motions is greater than that under uniform input ground motions.


## 4.2. Effect on Stress Distributions

The maximum principal shear (σtmax) and compressive (σcmax) stresses are commonly used to analyze the effect of seismic input ground motions on tunnels. The peak stresses under different incident angles are illustrated in Table 5.Table 5
Maximum stresses of tunnel body at different incident angles.
Motion form
Principal stress
0°
30°
45°
60°
90°
Uniform
σ
t
max
11.60
10.85
9.24
7.18
7.45
σ
c
max
−11.60
−10.58
−10.02
−8.99
−6.81
Multisupport
σ
t
max
17.74
13.92
11.29
8.94
6.66
σ
c
max
−18.84
−15.68
−13.37
−11.16
−8.49It is indicated that the maximum tensile stress is identical to the maximum compressive stress while the incident angle is 0° for uniform seismic input tunnels. This is owing to the tunnel and soil profile placed in pure shear stress state. When incident angles are equal to 30°, 45°, 60°, and 90°, the values of maximum tensile and compressive stresses are different because of the complex shear and normal stress states induced by seismic load. It can be thus concluded that the earthquake acceleration in the vertical direction changes not only the maximum stress but also the stress characteristics. In addition, the spatially varying seismic ground motion is another factor influencing the stress state according to the different values between the maximum tensile stress and maximum compressive stress at an incident angle of 0° under spatially varying seismic ground motions. Meanwhile, the value of maximum tensile/compressive stresses decreases with the increasing incident angles. The maximum tensile stress nephogram at an incident angle of 0° for both uniform and spatially varying seismic ground motions is, respectively, shown in Figures4 and 5.Figure 4
The maximum tensile stress nephogram at an incident angle of 0° under uniform seismic input motion.Figure 5
The maximum tensile stress nephogram at an incident angle of 0° under multisupport seismic input motion.The red upper portion of the tunnel end (Z = 1000 m) in Figure 4 is the maximum tensile stress, responding to the red lower portion near the tunnel end. It is suggested that the spatially varying seismic ground motion affects the stress distribution.


## 5. Conclusions

Numerical simulation was an effective method for studying the tunnel responses under seismic event. In this paper, a novel large-scale analytical method was established for estimating the seismic response of long tunnels within soil foundations. Based on the 3D FEM platform, this method could offer a reliable way for investigating the nonlinear seismic behavior of long tunnels. Main findings are summarized as follows:(1)
The spatially varying seismic ground motions induced larger diameter strain rates (radial deformation) than the uniform input seismic motion. The maximum radial strain rate increased with the increasing incident angles for both seismic input forms.(2)
For uniform seismic input simulations, the maximum radial strain rate occurred at the middle point of the tunnel and the radial strain rates were symmetrically distributed on the left and right sides of the tunnel. Moreover, the spatially varying seismic ground motions had an asymmetric effect on the radial strain rate distributions.(3)
The rising incident angles changed the pure shear stress state into a complex stress state for tunnels under specified input motion. The spatially varying seismic ground motions also influenced the stress state relative to the uniform seismic input. Meanwhile, the values of maximum tensile/compressive stresses decreased with the increasing incident angles.



---

*Source: 8195396-2018-06-11.xml*