# Analytical Solution for Circular Tunnel under Obliquely Incident P Waves considering Different Contact Conditions

Authors: Yawei Duan; Mi Zhao; Jingqi Huang; Huifang Li; Xiuli Du
Journal: Shock and Vibration
 (2021)

---

## Abstract

An analytical solution for the seismic-induced thrust and moment of the circular tunnel in half-space under obliquely incident P waves is developed in this study, which is the superposition of the solution for deep tunnels under incident and reflected P waves and the reflected SV waves. To consider tangential contact stiffness at the ground-tunnel interface, a spring-type stiffness coefficient is introduced into the force-displacement relationship. Moreover, the tunnel lining is treated as the thick-wall cylinder, providing more precise forecasts than beam or shell models used in previous analytical solution, especially for tunnels with thick lining. The reliability of the proposed analytical solution is assessed by comparing with the dynamic numerical results. Based on the proposed analytical solution, parametrical studies are conducted to investigate the effect of some critical factors on the tunnel’s seismic response, including the incident angles, the tangential contact stiffness at the ground-tunnel interface, and the relative stiffness between the ground and the tunnel. The results demonstrate that the proposed analytical solution performs well and can be adopted to predict the internal forces of circular tunnels under obliquely incident P waves in seismic design.


---

## Body

## 1. Introduction

As the essential part of infrastructures, tunnel structures play an important role in public transportation facilities, sanitation frameworks, irrigation utilities, and so on. Besides thousands of land tunnels, many undersea tunnels are constructed or under the plan for construction with the rapid development of coastal transportation in recent years. However, seismic damage cases in recent years indicate that tunnel structures are vulnerable to earthquake loading [1–4]. Sometimes, a low level of tunnel damage may affect structural serviceability and bring about substantial life and property lost, especially for undersea tunnels. Therefore, the tunnel linings should be carefully designed using an effective design method to withstand the additional loads arising from earthquakes.In the last several decades, plenty of seismic analysis and design methods have been developed for tunnels and other underground structures, which can be divided into analytical solution, quasi-static analysis method, and dynamic time history analysis method. Compared with the quasi-static analysis method and dynamic time history analysis method, the analytical solution can quickly and easily estimate the force and displacement of the tunnel under seismic loadings, which can identify the most relevant variables of the problem. Therefore, analytical solution has been widely used in the preliminary seismic design of tunnels. Many researchers have been conducted to develop analytical solution of circular tunnels for seismic analysis. Wang [5] proposed the analytical solution for thrust, moment, and deformation of the circular tunnel under vertically propagating shear waves. Penzien and Wu [6] put forward a simple formula for calculating the deformation and corresponding stress of bored tunnels resulting from earthquake loading. Then, Penzien [7] supplemented the previously published by Penzien and Wu [6] and proposed a method to calculate the deformation of rectangular and circular tunnel lining under seismic action. The above analytical solution assumes full-slip condition or no-slip condition between the ground and the tunnel lining. However, the contact state at the ground-tunnel interface is between full-slip and no-slip. In this regard, Park et al. [8] established a new analytical solution by using a spring-type flexibility coefficient D to consider the slippage effect of the ground-tunnel interface. Fang et al. [9] introduced the elastic-slip interface model to simulate the ground-tunnel interface properties and analyzed the dynamic response of twin tunnels subjected to blast waves. Then, Li et al.[10] used the same elastic-slip interface model and analyzed the influence of the elastic-slip interface in saturated poroelastic soil on the dynamic response of the underwater tunnel subjected to plane waves. Another important assumption adopted by most analytical solution is dry soil conditions. Bobet [11, 12] considered the influence of pore water pressure and presented an analytical solution for the circular tunnel subjected to seismic loads for drained and undrained loading conditions. All the above analytical solution are derived by assuming the tunnel lining as a beam or shell model, and there is a certain error for thick lining tunnels [13, 14]. Therefore, it is necessary to use a thick-wall cylinder to simulate the tunnel lining. Another important assumption is that the tunnel is assumed to be located in the full-space, and the seismic wave is considered as vertical incidence. In fact, the buried depth of the tunnel is usually limited, and the seismic incident direction has a significant impact on the underground structure in half-space [15–18]. Huang et al. [19] derived an analytical solution of the circular tunnel under obliquely incident SV waves based on a thick-wall cylinder model, considering the no-slip interface condition. Many studies show that the obliquely incident P waves has great influence on the seismic response of structures [20–24]. However, the analytical solution of the circular tunnel under obliquely incident P waves based on thick-wall cylinder model considering different contact conditions of the ground-tunnel interface is limited.In this paper, a new analytical solution for the internal force of the circular tunnel in half-space under obliquely incident P waves is developed by using a thick-walled cylinder to simulate the tunnel lining and introducing a spring-type stiffness coefficientk to consider the different tangential contact stiffness at the ground-tunnel interface. Then, the accuracy and applicability of the new analytical solution are verified by comparison with the dynamic time history analysis method. Finally, the effects of P waves’ incident angle, stiffness coefficient, and flexibility ratio are discussed.


## 2. Formulation of the Problem

Most of works [5, 7, 8, 12] assumed the earthquake motions as the plane body waves propagating along a vertical direction. In fact, the earthquake waves always have an arbitrary propagating direction [21]. Consider a circular tunnel located below the ground and subjected to obliquely incident P waves, as shown in Figure 1. The outer and inner radii are R and r, respectively, with the center O and the depth h under the half-space surface. The seismic excitation is a P wave with an incident angle of α in the x-y plane, and the distance between the wavefront at the zero time of the P waves and the center of the tunnel O is L. For the plane strain problem shown in Figure 1, the tunnel lining and its surrounding ground are assumed to be homogeneous isotropic linear elastic medium. The soft soil has nonlinear behavior under earthquake; linear elastic assumption may limit the practical applicability of the results obtained. However, the analytical solution develops a simple formulation for the preliminary seismic design of the tunnel, which can be used for further analysis. While reducing the calculation effort, some parameters that can reflect the real mechanical behavior of soil are discarded. At the same time, the analytical solution contains the most critical variables in the problem. Moreover, this study aims to the quasi-static analytical solution on the seismic-induced internal forces of the tunnel lining. Therefore, the geostatic stress is not considered. As we all know, the seismic response of underground structure is dominated by the surrounding ground response but not the inertial properties of the structure itself [19, 25]. Therefore, the dynamic inertial effect is ignored in the following derivation.Figure 1
Ground-tunnel interaction system in half-space under obliquely incident P waves.


## 3. New Analytical Solution

### 3.1. Wave Motions in Half-Space Ground and Pseudostatic Analysis Model

The free filed wave motions in the half-space ground under obliquely incident P waves are shown in Figure2. The incident P waves reach the ground surface with an oblique angle of α. According to the wave motion theories [26], reflected P waves and SV waves are generated at the free surface, with the reflected angles of α and β, respectively. As the tunnel being constructed in the ground, the tunnel will be excited by three waves, i.e., incident P waves, reflected P waves, and reflected SV waves. Therefore, the seismic response of the tunnel is superposition of the response of the tunnel under these three waves. In order to calculate the tunnel’s response conveniently, local coordinate systems (η1, ξ1), (η2, ξ2), and (η3, ξ3) are adopted for the three waves, respectively. In local coordinate systems, the velocity-time histories of the three waves are as follows:(1)v1t=v0t−Δt1,v2t=A1v0t−Δt2,v3t=A2v0t−Δt3,where A1 and A2 are the amplitudes of the reflected P waves and the reflected SV wave, Δt1, Δt2, and Δt3 are the wave propagation times of the incident P waves, the reflected P waves, and the reflected SV waves from the zero-time wavefront to the center location of the tunnel, respectively, and t stands for time. These parameters can be expressed as [27](2)A1=cs2sin2αsin2β−cp2cos22βcs2sin2αsin2β+cp2cos22β,A2=−2cpcssin2αcos2βcs2sin2αsin2β+cp2cos22β,β=arcsincssinαcp,Δt1=Lcp,Δt2=L+2hcosαcp,Δt3=L+hcosα−hsinαtanβcp+hcscosβ,Where cs and cp are the shear wave velocity and compression wave velocity in the ground.Figure 2
Wave motions in half-space under obliquely incident P waves.As the wave length of peak velocities is at least 8 times larger than the width of the opening, the free-field stress gradient across the opening is relatively small and the seismic loading can be considered as a pseudo-static load. In most underground openings the quasi-static conditions are usually satisfied [11, 12]. Therefore, the pseudo-static approach is employed to analysis the tunnel response under the three waves. The pseudo-static approach is an accepted method of analysis in current seismic design practice, which has been widely used in earlier studies for analytical solution [6–8, 11, 19, 28]. In the pseudo-static approach, the incident seismic wave is simulated by applying a static far-field stress or strain at the model boundaries to represent the seismic effects. Figure 3 offers the pseudo-static analysis model of the ground-tunnel system for the incident and reflected P waves and reflected SV waves, respectively.Figure 3
The pseudo-static analysis model of the ground-tunnel system: (a) for the incident P waves and reflected P waves and (b) for the reflected SV waves.
(a)(b)The compressive stressσν1 and σh1 caused by incident P waves, the compressive stress σν2 and σh2 caused by reflected P wave, and the shear stress τ caused by reflected SV wave can be expressed as(3)σv1t=Eg1−υg1+υg1−2υgcpv0t−Δt1,σh1t=Egυg1+υg1−2υgcpv0t−Δt1,(4)σv2t=Eg1−υg1+υg1−2υgcpA1v0t−Δt2,σh2t=Egυg1+υg1−2υgcpA1v0t−Δt2,(5)τt=−Eg21+υgcsA2v0t−Δt3,where t stands for time and Eg and νg represent Young’s modulus and Poisson’s ratio of the ground, respectively.In the following section, the analytical solution of the tunnel’s internal force under the incident P waves, the reflected P waves, and the reflected SV waves will be, respectively, derived. Then, the three solution are linearly superimposed to obtain the final solution of the tunnel.


### 3.2. Internal Force of Tunnel under Reflected SV Waves

This section gives the analytical solution of a circular tunnel under reflected SV waves, and the pseudo-static ground-tunnel interaction system is shown in Figure3(b). The shear stress can be obtained from Equation (5). Different from previous researchers such as Hoeg [29], Penzien [7], Park et al. [8], and Bobet [11, 12], who regarded the tunnel lining as a beam or a shell, this paper regards the tunnel lining as a thick-wall cylinder, which reduces the errors brought by the thick lining to the analytical solution. Moreover, we introduce a spring-type stiffness coefficient k to consider the tangential contact stiffness at the ground-tunnel interface.The model of the ground-tunnel system in local coordinate system (η, ξ), as shown in Figure 3(b), is transformed to the polar coordinates (ρ, φ), as shown in Figure 4(a). Using coordinate transformation, the far-field shear stress can be expressed as(6)σρ|ρ=∞=τsin2φ,τρφ|ρ=∞=τcos2φ.Figure 4
Decomposition of the ground-tunnel system under S waves: (a) ground-tunnel system in polar coordinates; (b) ground system; (c) tunnel system.For the convenience in analysis, as shown in Figure4(a), the ground-tunnel system was decomposed into two cases: (1) a circular cylindrical cavity subjected to far-field shear stress and contact stress at the ground-tunnel interface (Figure 4(b)); (2) thick-walled cylinder subjected to contact stress at the ground-tunnel interface (Figure 4(c)). The ground-tunnel interface’s normal and tangential contact stresses can be expressed as psin2φ and qcos2φ [6, 28]. In the following formula, subscripts g and l represent the ground and the tunnel, respectively, and subscripts ρ and φ represent the radial and circumferential directions.To consider the relationship between displacement and interaction forces at the ground-tunnel interface, following the work of Park et al. [8], a spring-type stiffness coefficient k is introduced. k ⟶ 0 and k ⟶ ∞ are two extreme states, which represent the full-slip interface condition and the no-slip interface condition, respectively. The displacement and stress at the interface (ρ = R) should be satisfied, namely,(7)uρg|ρ=R=uρl|ρ=R,uφg|ρ=R=uφl|ρ=R+qcos2φk,(8)τρφg|ρ=R=τρφl|ρ=R,σρg|ρ=R=σρl|ρ=R.The displacement and stress components of the cylinder in polar coordinates are obtained from the theory of elasticity [30]:(9)uρ=−m2Bρ−4Cρ−2Dρ3sin2φ+n4Aρ3+2Bρ−2Dρ3sin2φ+Isinφ+Kcosφ,uφ=−m212Aρ3+4Bρ−4Cρ+4Dρ3cos2φ+n24Aρ3+4Bρ+4Cρ+4Dρ3cos2φ+Icosφ−Ksinφ,(10)σρ=−2B+4Cρ2+6Dρ4sin2φ,σφ=12Aρ2+2B+6Dρ4sin2φ,τρφ=−6Aρ2+2B−2Cρ2−6Dρ4cos2φ,where A, B, C, D, I, and K are undetermined constants and I and K represent the rigid body displacement components. m and n are expressed by Young’s modulus E and Poisson’s ratio ν, which are (1 − υ2)/E and  (υ + υ2)/E, respectively.Substituting equation (10) into (6), we can obtain(11)τρφg|ρ=∞=−6Agρ2+2Bg−2Cgρ2−6Dgρ4cos2φ=τcos2φ,σρg|ρ=∞=−2Bg+4Cgρ2+6Dgρ4sin2φ=τsin2φ.Whenρ ⟶ ∞, 1/ρ2 and 1/ρ4 tend to 0. Therefore, it can be concluded that(12)Ag=0,Bg=−τ2.Because the inner surface of the tunnel is free, the following boundary conditions are satisfied:(13)τρφl|ρ=r=−6Alr2+2Bl−2Clr2−6Dlr4cos2φ=0,σρl|ρ=r=−2Bl+4Clr2+6Dlr4sin2φ=0.Substitute equation (10) into (8) to obtain(14)6AlR2+2Bl−2ClR2−6DlR4−6AgR2−2Bg+2CgR2+6DgR4=0,2Bl+4ClR2+6DlR4−2Bg−4CgR2−6DgR4=0.Substitute equation (9) into Equation (7), and it becomes(15)−mg2BgR−4CgR−2DgR3sin2φ+ng4AgR3+2BgR−2DgR3sin2φ+Igsinφ+Kgcosφ=−ml2BlR−4ClR−2DlR3sin2φ+nl4AlR3+2BlR−2DlR3sin2φ+Ilsinφ+Klcosφ−mg212AgR3+4BgR−4CgR+4DgR3cos2φ+ng24AgR3+4BgR+4CgR+4DgR3cos2φ+Igcosφ−Kgsinφ+qcos2φk=−ml212AlR3+4BlR−4ClR+4DlR3cos2φ+nl24AlR3+4BlR+4ClR+4DlR3cos2φ+Ilcosφ−Klsinφ.Equation (15) should be established when φ takes any value. Substituting φ = 0 and φ = 3π/4 into Equation (15), one obtains(16)Ig=Il,Kg=Kl..The tangential contact stress can be expressed as(17)qcos2φ=−τρφl=6AlR2+2Bl−2ClR2−6DlR4cos2φ.Substitute equations (16) and (17) into (15) to simplify(18)4nlR3Al−2ml−nlRBl+4mlRCl+2ml−nlR3Dl−4nlR3Ag+2mg−ngRBg−4mgRCg−2mg−ngR3Dg=0,−6ml−2nlR3+6R2kAl+−2ml−nlR+2kBl+2ml+nlR−2kR2Cl+−2ml−nlR3−6kR4Dl+6mg−2ngR3Ag+2mg−ngRBg−2mg+ngRCg+2mg−ngR3Dg=0.Solving equations (12), (13), (14), and (18) simultaneously, the constants Al, Bl, Cl, and Dl can be solved as(19)A1=−η4qn−qt+3η2qn+qt6η2−13r02τ,B1=η6+η4+2η2qn+2η2qt2η2−13τ,C1=−r122η4qn+η2+1qn+qt2η2−13τ,D1=r043η6qn+η4qn+2η4qt6η2−13τ,where qn, qt, and η are expressed as(20)qn=−4mgmg−ng/3R+δ12+δ22+1/kR4mg+2ng/3R−δ12mg−ng/3R+δ12+δ22+1/k+mg−ng/3R−δ11+δ125mg+ng/3R+δ22+1/k,qt=−4mgmg−ng/3R−δ11+δ12R4mg+2ng/3R−δ12mg−ng/3R+δ12+δ22+1/k+mg−ng/3R−δ11+δ125mg+ng/3R+δ22+1/k,δ11=2rnl−3mlη7+3η3−nl+5ml3η5+η6η2−13,δ12=4rnlη7−3nlη5+3nl−2mlη3−nl+2mlη6η2−13,δ22=2r3ml−nlη7−3η5+3η3+nl+5mlη6η2−13,η=Rr0.The thrust and the moment of the tunnel can be obtained by integration:(21)Tsv=∫rRσφdρ,Msv=∫rRσφρ−R+r2dρ.Substitute equation (19) into (10) and then into equation (21). Finally, the analytical solution for the thrust Tsv and the moment Msv of the circular tunnel under far-field shear stress in the local coordinate system (η, ξ) is obtained:(22)Tsv=R3qn+2qtτsin2φ,Msv=3R2qt−R2+Rr0qn+2qt6τsin2φ.Since the propagation angle of the reflected SV waves in the (X, Y) coordinate system is π − β, the tunnel thrust Tsvr and moment Msvr under the shearing stress τ caused by the reflected SV waves can be obtained from Equation (22) by simply rotating the angle π − β:(23)Tsvr=R3qn+2qtτsin2φ+π−β,Msvr=3R2qt−R2+Rr0qn+2qt6τsin2φ+π−β.


### 3.3. Internal Force of Tunnel under Incident P Waves and Reflected P Waves

This section gives the analytical solution of a circular tunnel under incident P waves and reflected P waves, and the pseudo-static ground-tunnel system at the local coordinate system is shown in Figure3(a). The far-field stresses σν and σh induced by incident P waves and reflected P waves are expressed by equations (3) and (4), respectively. The far-field stress of tunnel lining can be decomposed into hydrostatic pressure component and pure shear stress component, as shown in Figure 5; the hydrostatic pressure component p and pure shear stress τ component are expressed in equation (24). Note that the analytical solution for the ground-tunnel system under the pure shear stress τ shown in Figure 5(c) can be obtained by equation (22), which only rotates at an angle of π/4. Then, an analytical solution for tunnel internal force under hydrostatic pressure p is derived. The model of the ground-tunnel system in the local coordinate system (η, ξ), as shown in Figure 5(b), is transformed to the polar coordinates (ρ, φ), as shown in Figure 6(a).Figure 5
Force decomposition of the pseudo-static analysis model under P waves: (a) vertical stressσν and horizontal stress σh; (b) hydrostatic pressure component p; (c) pure shear stress component τ.
(a)(b)(c)Figure 6
Decomposition of the ground-tunnel system under hydrostatic pressure component (p): (a) ground-tunnel system in polar coordinates; (b) ground system; (c) tunnel system.(24)p=σv+σh2,τ=σv−σh2.For the convenience of analysis, the tunnel subjected to hydrostatic pressure (Figure6(a)) can be further decomposed into a cylindrical cavity subjected to an outward internal pressure pn and inward external pressure p (Figure 6(b)); the tunnel lining is subjected to the inward external pressure pn and the inner surface is free (Figure 6(c)). At the ground and the tunnel contact interface, the stress and displacement should be satisfied:(25)σρg|ρ=R=σρl|ρ=R,uρg|ρ=R=uρl|ρ=R.According to the theory of elasticity, the general solution of stress component, and displacement component under axisymmetric stress state [30], the radial and circumferential stress components are expressed as(26)σρ=Aρ2+2C,σφ=−Aρ2+2C.The corresponding radial displacement component is(27)uρ=n−mAρ+2m+nCρ+Icosφ+Ksinφ.For the ground-tunnel system, the stress boundary condition is(28)σρg|ρ=∞=Agρ2+2Cg=−p,σρl|ρ=r=Alr2+2Cl=0.In the contact interface between the ground and the tunnel, the continuity conditions of stress and displacement should be satisfied. Equations (26) and (27) are substituted into equation (5); one obtains(29)ng−mgAgR+2mg+ngCgR+Igsinφ+Kgcosφ=nl−mlAlR+2ml+nlClR+Ilsinφ+Klcosφ,AgR2+2Cg=AlR2+2Cl.Equation (29) should hold for any value of φ. Substituting φ equals 0, π/2, π, and 3π/2 into equation (29), we obtain(30)Ig=Il,Kg=Kl.Equation (29) is simplified as(31)ng−mgAgR+2mg+ngCgR+ml−nlAlR−2ml+nlClR=0,AgR2+2Cg−AlR2−2Cl=0.Combining equations (28) and (31) to solve the equation, the constants Al and Cl can be solved as(32)Al=2A′R2r2σv+σh2,Cl=−A′R2σv+σh2,where A′=mg/−mg+ng+ml−nlr2+mg−ng+ml+nlR2The thrustTh and the moment Mh of the tunnel lining under the hydrostatic pressure p can be obtained by integration:(33)Th=σv+σhA′Rr2−R2,Mh=σv+σhA′R2r2lnrR+R2−r22Rr.The analytical solution for the thrustTp and the moment Mp of the tunnel under the P waves in the local coordinate system (η, ξ) can be obtained by superposing the solution for the internal forces of the tunnel subjected to hydrostatic pressure p and pure shear stress τ:(34)Tpφ=Th+Tsvφ+π4,Mpφ=Mh+Msvφ+π4.Namely,(35)Tp=σv+σhA′Rr2−R2+R6qn+2qtσv−σhsin2φ+π4,Mp=σv+σhA′R2r2lnrR+R2−r22Rr+3R2qt−R2+Rrqn+2qt12σv−σhsin2φ+π4.Since the incident angle of P waves in the (X, Y) coordinate system is α, the analytical solution of tunnel thrust Tpi and moment Mpi under incident P waves can be obtained from equation (35) by transforming φ into φ + α:(36)Tpi=σv1+σh1A′Rr2−R2+R6qn+2qtσv1−σh1cos2φ+α,Mpi=σv1+σh1A′R2r2lnrR+R2−r22Rr+3R2qt−R2+Rrqn+2qt12σv1−σh1cos2φ+α.Similar to the solution of incident P waves, the analytic solution for the thrustTpr and the moment Mpr of tunnel under reflected P waves can be obtained by rotating π−α of φ angle in equation (35):(37)Tpr=σv2+σh2A′Rr2−R2+R6qn+2qtσv2−σh2cos2φ−α,Mpr=σv2+σh2A′R2r2lnrR+R2−r22Rr+3R2qt−R2+Rrqn+2qt12σv2−σh2cos2φ−α.


### 3.4. Internal Force of Tunnel under Obliquely Incident P Waves

By superimposing equations (23), (36), and (37), the analytical solution for the total thrust T and the moment M of the tunnel under obliquely incident P waves can be obtained:(38)T=Tsvr+Tpi+Tpr,M=Msvr+Mpi+Mpr.


## 3.1. Wave Motions in Half-Space Ground and Pseudostatic Analysis Model

The free filed wave motions in the half-space ground under obliquely incident P waves are shown in Figure2. The incident P waves reach the ground surface with an oblique angle of α. According to the wave motion theories [26], reflected P waves and SV waves are generated at the free surface, with the reflected angles of α and β, respectively. As the tunnel being constructed in the ground, the tunnel will be excited by three waves, i.e., incident P waves, reflected P waves, and reflected SV waves. Therefore, the seismic response of the tunnel is superposition of the response of the tunnel under these three waves. In order to calculate the tunnel’s response conveniently, local coordinate systems (η1, ξ1), (η2, ξ2), and (η3, ξ3) are adopted for the three waves, respectively. In local coordinate systems, the velocity-time histories of the three waves are as follows:(1)v1t=v0t−Δt1,v2t=A1v0t−Δt2,v3t=A2v0t−Δt3,where A1 and A2 are the amplitudes of the reflected P waves and the reflected SV wave, Δt1, Δt2, and Δt3 are the wave propagation times of the incident P waves, the reflected P waves, and the reflected SV waves from the zero-time wavefront to the center location of the tunnel, respectively, and t stands for time. These parameters can be expressed as [27](2)A1=cs2sin2αsin2β−cp2cos22βcs2sin2αsin2β+cp2cos22β,A2=−2cpcssin2αcos2βcs2sin2αsin2β+cp2cos22β,β=arcsincssinαcp,Δt1=Lcp,Δt2=L+2hcosαcp,Δt3=L+hcosα−hsinαtanβcp+hcscosβ,Where cs and cp are the shear wave velocity and compression wave velocity in the ground.Figure 2
Wave motions in half-space under obliquely incident P waves.As the wave length of peak velocities is at least 8 times larger than the width of the opening, the free-field stress gradient across the opening is relatively small and the seismic loading can be considered as a pseudo-static load. In most underground openings the quasi-static conditions are usually satisfied [11, 12]. Therefore, the pseudo-static approach is employed to analysis the tunnel response under the three waves. The pseudo-static approach is an accepted method of analysis in current seismic design practice, which has been widely used in earlier studies for analytical solution [6–8, 11, 19, 28]. In the pseudo-static approach, the incident seismic wave is simulated by applying a static far-field stress or strain at the model boundaries to represent the seismic effects. Figure 3 offers the pseudo-static analysis model of the ground-tunnel system for the incident and reflected P waves and reflected SV waves, respectively.Figure 3
The pseudo-static analysis model of the ground-tunnel system: (a) for the incident P waves and reflected P waves and (b) for the reflected SV waves.
(a)(b)The compressive stressσν1 and σh1 caused by incident P waves, the compressive stress σν2 and σh2 caused by reflected P wave, and the shear stress τ caused by reflected SV wave can be expressed as(3)σv1t=Eg1−υg1+υg1−2υgcpv0t−Δt1,σh1t=Egυg1+υg1−2υgcpv0t−Δt1,(4)σv2t=Eg1−υg1+υg1−2υgcpA1v0t−Δt2,σh2t=Egυg1+υg1−2υgcpA1v0t−Δt2,(5)τt=−Eg21+υgcsA2v0t−Δt3,where t stands for time and Eg and νg represent Young’s modulus and Poisson’s ratio of the ground, respectively.In the following section, the analytical solution of the tunnel’s internal force under the incident P waves, the reflected P waves, and the reflected SV waves will be, respectively, derived. Then, the three solution are linearly superimposed to obtain the final solution of the tunnel.


## 3.2. Internal Force of Tunnel under Reflected SV Waves

This section gives the analytical solution of a circular tunnel under reflected SV waves, and the pseudo-static ground-tunnel interaction system is shown in Figure3(b). The shear stress can be obtained from Equation (5). Different from previous researchers such as Hoeg [29], Penzien [7], Park et al. [8], and Bobet [11, 12], who regarded the tunnel lining as a beam or a shell, this paper regards the tunnel lining as a thick-wall cylinder, which reduces the errors brought by the thick lining to the analytical solution. Moreover, we introduce a spring-type stiffness coefficient k to consider the tangential contact stiffness at the ground-tunnel interface.The model of the ground-tunnel system in local coordinate system (η, ξ), as shown in Figure 3(b), is transformed to the polar coordinates (ρ, φ), as shown in Figure 4(a). Using coordinate transformation, the far-field shear stress can be expressed as(6)σρ|ρ=∞=τsin2φ,τρφ|ρ=∞=τcos2φ.Figure 4
Decomposition of the ground-tunnel system under S waves: (a) ground-tunnel system in polar coordinates; (b) ground system; (c) tunnel system.For the convenience in analysis, as shown in Figure4(a), the ground-tunnel system was decomposed into two cases: (1) a circular cylindrical cavity subjected to far-field shear stress and contact stress at the ground-tunnel interface (Figure 4(b)); (2) thick-walled cylinder subjected to contact stress at the ground-tunnel interface (Figure 4(c)). The ground-tunnel interface’s normal and tangential contact stresses can be expressed as psin2φ and qcos2φ [6, 28]. In the following formula, subscripts g and l represent the ground and the tunnel, respectively, and subscripts ρ and φ represent the radial and circumferential directions.To consider the relationship between displacement and interaction forces at the ground-tunnel interface, following the work of Park et al. [8], a spring-type stiffness coefficient k is introduced. k ⟶ 0 and k ⟶ ∞ are two extreme states, which represent the full-slip interface condition and the no-slip interface condition, respectively. The displacement and stress at the interface (ρ = R) should be satisfied, namely,(7)uρg|ρ=R=uρl|ρ=R,uφg|ρ=R=uφl|ρ=R+qcos2φk,(8)τρφg|ρ=R=τρφl|ρ=R,σρg|ρ=R=σρl|ρ=R.The displacement and stress components of the cylinder in polar coordinates are obtained from the theory of elasticity [30]:(9)uρ=−m2Bρ−4Cρ−2Dρ3sin2φ+n4Aρ3+2Bρ−2Dρ3sin2φ+Isinφ+Kcosφ,uφ=−m212Aρ3+4Bρ−4Cρ+4Dρ3cos2φ+n24Aρ3+4Bρ+4Cρ+4Dρ3cos2φ+Icosφ−Ksinφ,(10)σρ=−2B+4Cρ2+6Dρ4sin2φ,σφ=12Aρ2+2B+6Dρ4sin2φ,τρφ=−6Aρ2+2B−2Cρ2−6Dρ4cos2φ,where A, B, C, D, I, and K are undetermined constants and I and K represent the rigid body displacement components. m and n are expressed by Young’s modulus E and Poisson’s ratio ν, which are (1 − υ2)/E and  (υ + υ2)/E, respectively.Substituting equation (10) into (6), we can obtain(11)τρφg|ρ=∞=−6Agρ2+2Bg−2Cgρ2−6Dgρ4cos2φ=τcos2φ,σρg|ρ=∞=−2Bg+4Cgρ2+6Dgρ4sin2φ=τsin2φ.Whenρ ⟶ ∞, 1/ρ2 and 1/ρ4 tend to 0. Therefore, it can be concluded that(12)Ag=0,Bg=−τ2.Because the inner surface of the tunnel is free, the following boundary conditions are satisfied:(13)τρφl|ρ=r=−6Alr2+2Bl−2Clr2−6Dlr4cos2φ=0,σρl|ρ=r=−2Bl+4Clr2+6Dlr4sin2φ=0.Substitute equation (10) into (8) to obtain(14)6AlR2+2Bl−2ClR2−6DlR4−6AgR2−2Bg+2CgR2+6DgR4=0,2Bl+4ClR2+6DlR4−2Bg−4CgR2−6DgR4=0.Substitute equation (9) into Equation (7), and it becomes(15)−mg2BgR−4CgR−2DgR3sin2φ+ng4AgR3+2BgR−2DgR3sin2φ+Igsinφ+Kgcosφ=−ml2BlR−4ClR−2DlR3sin2φ+nl4AlR3+2BlR−2DlR3sin2φ+Ilsinφ+Klcosφ−mg212AgR3+4BgR−4CgR+4DgR3cos2φ+ng24AgR3+4BgR+4CgR+4DgR3cos2φ+Igcosφ−Kgsinφ+qcos2φk=−ml212AlR3+4BlR−4ClR+4DlR3cos2φ+nl24AlR3+4BlR+4ClR+4DlR3cos2φ+Ilcosφ−Klsinφ.Equation (15) should be established when φ takes any value. Substituting φ = 0 and φ = 3π/4 into Equation (15), one obtains(16)Ig=Il,Kg=Kl..The tangential contact stress can be expressed as(17)qcos2φ=−τρφl=6AlR2+2Bl−2ClR2−6DlR4cos2φ.Substitute equations (16) and (17) into (15) to simplify(18)4nlR3Al−2ml−nlRBl+4mlRCl+2ml−nlR3Dl−4nlR3Ag+2mg−ngRBg−4mgRCg−2mg−ngR3Dg=0,−6ml−2nlR3+6R2kAl+−2ml−nlR+2kBl+2ml+nlR−2kR2Cl+−2ml−nlR3−6kR4Dl+6mg−2ngR3Ag+2mg−ngRBg−2mg+ngRCg+2mg−ngR3Dg=0.Solving equations (12), (13), (14), and (18) simultaneously, the constants Al, Bl, Cl, and Dl can be solved as(19)A1=−η4qn−qt+3η2qn+qt6η2−13r02τ,B1=η6+η4+2η2qn+2η2qt2η2−13τ,C1=−r122η4qn+η2+1qn+qt2η2−13τ,D1=r043η6qn+η4qn+2η4qt6η2−13τ,where qn, qt, and η are expressed as(20)qn=−4mgmg−ng/3R+δ12+δ22+1/kR4mg+2ng/3R−δ12mg−ng/3R+δ12+δ22+1/k+mg−ng/3R−δ11+δ125mg+ng/3R+δ22+1/k,qt=−4mgmg−ng/3R−δ11+δ12R4mg+2ng/3R−δ12mg−ng/3R+δ12+δ22+1/k+mg−ng/3R−δ11+δ125mg+ng/3R+δ22+1/k,δ11=2rnl−3mlη7+3η3−nl+5ml3η5+η6η2−13,δ12=4rnlη7−3nlη5+3nl−2mlη3−nl+2mlη6η2−13,δ22=2r3ml−nlη7−3η5+3η3+nl+5mlη6η2−13,η=Rr0.The thrust and the moment of the tunnel can be obtained by integration:(21)Tsv=∫rRσφdρ,Msv=∫rRσφρ−R+r2dρ.Substitute equation (19) into (10) and then into equation (21). Finally, the analytical solution for the thrust Tsv and the moment Msv of the circular tunnel under far-field shear stress in the local coordinate system (η, ξ) is obtained:(22)Tsv=R3qn+2qtτsin2φ,Msv=3R2qt−R2+Rr0qn+2qt6τsin2φ.Since the propagation angle of the reflected SV waves in the (X, Y) coordinate system is π − β, the tunnel thrust Tsvr and moment Msvr under the shearing stress τ caused by the reflected SV waves can be obtained from Equation (22) by simply rotating the angle π − β:(23)Tsvr=R3qn+2qtτsin2φ+π−β,Msvr=3R2qt−R2+Rr0qn+2qt6τsin2φ+π−β.


## 3.3. Internal Force of Tunnel under Incident P Waves and Reflected P Waves

This section gives the analytical solution of a circular tunnel under incident P waves and reflected P waves, and the pseudo-static ground-tunnel system at the local coordinate system is shown in Figure3(a). The far-field stresses σν and σh induced by incident P waves and reflected P waves are expressed by equations (3) and (4), respectively. The far-field stress of tunnel lining can be decomposed into hydrostatic pressure component and pure shear stress component, as shown in Figure 5; the hydrostatic pressure component p and pure shear stress τ component are expressed in equation (24). Note that the analytical solution for the ground-tunnel system under the pure shear stress τ shown in Figure 5(c) can be obtained by equation (22), which only rotates at an angle of π/4. Then, an analytical solution for tunnel internal force under hydrostatic pressure p is derived. The model of the ground-tunnel system in the local coordinate system (η, ξ), as shown in Figure 5(b), is transformed to the polar coordinates (ρ, φ), as shown in Figure 6(a).Figure 5
Force decomposition of the pseudo-static analysis model under P waves: (a) vertical stressσν and horizontal stress σh; (b) hydrostatic pressure component p; (c) pure shear stress component τ.
(a)(b)(c)Figure 6
Decomposition of the ground-tunnel system under hydrostatic pressure component (p): (a) ground-tunnel system in polar coordinates; (b) ground system; (c) tunnel system.(24)p=σv+σh2,τ=σv−σh2.For the convenience of analysis, the tunnel subjected to hydrostatic pressure (Figure6(a)) can be further decomposed into a cylindrical cavity subjected to an outward internal pressure pn and inward external pressure p (Figure 6(b)); the tunnel lining is subjected to the inward external pressure pn and the inner surface is free (Figure 6(c)). At the ground and the tunnel contact interface, the stress and displacement should be satisfied:(25)σρg|ρ=R=σρl|ρ=R,uρg|ρ=R=uρl|ρ=R.According to the theory of elasticity, the general solution of stress component, and displacement component under axisymmetric stress state [30], the radial and circumferential stress components are expressed as(26)σρ=Aρ2+2C,σφ=−Aρ2+2C.The corresponding radial displacement component is(27)uρ=n−mAρ+2m+nCρ+Icosφ+Ksinφ.For the ground-tunnel system, the stress boundary condition is(28)σρg|ρ=∞=Agρ2+2Cg=−p,σρl|ρ=r=Alr2+2Cl=0.In the contact interface between the ground and the tunnel, the continuity conditions of stress and displacement should be satisfied. Equations (26) and (27) are substituted into equation (5); one obtains(29)ng−mgAgR+2mg+ngCgR+Igsinφ+Kgcosφ=nl−mlAlR+2ml+nlClR+Ilsinφ+Klcosφ,AgR2+2Cg=AlR2+2Cl.Equation (29) should hold for any value of φ. Substituting φ equals 0, π/2, π, and 3π/2 into equation (29), we obtain(30)Ig=Il,Kg=Kl.Equation (29) is simplified as(31)ng−mgAgR+2mg+ngCgR+ml−nlAlR−2ml+nlClR=0,AgR2+2Cg−AlR2−2Cl=0.Combining equations (28) and (31) to solve the equation, the constants Al and Cl can be solved as(32)Al=2A′R2r2σv+σh2,Cl=−A′R2σv+σh2,where A′=mg/−mg+ng+ml−nlr2+mg−ng+ml+nlR2The thrustTh and the moment Mh of the tunnel lining under the hydrostatic pressure p can be obtained by integration:(33)Th=σv+σhA′Rr2−R2,Mh=σv+σhA′R2r2lnrR+R2−r22Rr.The analytical solution for the thrustTp and the moment Mp of the tunnel under the P waves in the local coordinate system (η, ξ) can be obtained by superposing the solution for the internal forces of the tunnel subjected to hydrostatic pressure p and pure shear stress τ:(34)Tpφ=Th+Tsvφ+π4,Mpφ=Mh+Msvφ+π4.Namely,(35)Tp=σv+σhA′Rr2−R2+R6qn+2qtσv−σhsin2φ+π4,Mp=σv+σhA′R2r2lnrR+R2−r22Rr+3R2qt−R2+Rrqn+2qt12σv−σhsin2φ+π4.Since the incident angle of P waves in the (X, Y) coordinate system is α, the analytical solution of tunnel thrust Tpi and moment Mpi under incident P waves can be obtained from equation (35) by transforming φ into φ + α:(36)Tpi=σv1+σh1A′Rr2−R2+R6qn+2qtσv1−σh1cos2φ+α,Mpi=σv1+σh1A′R2r2lnrR+R2−r22Rr+3R2qt−R2+Rrqn+2qt12σv1−σh1cos2φ+α.Similar to the solution of incident P waves, the analytic solution for the thrustTpr and the moment Mpr of tunnel under reflected P waves can be obtained by rotating π−α of φ angle in equation (35):(37)Tpr=σv2+σh2A′Rr2−R2+R6qn+2qtσv2−σh2cos2φ−α,Mpr=σv2+σh2A′R2r2lnrR+R2−r22Rr+3R2qt−R2+Rrqn+2qt12σv2−σh2cos2φ−α.


## 3.4. Internal Force of Tunnel under Obliquely Incident P Waves

By superimposing equations (23), (36), and (37), the analytical solution for the total thrust T and the moment M of the tunnel under obliquely incident P waves can be obtained:(38)T=Tsvr+Tpi+Tpr,M=Msvr+Mpi+Mpr.


## 4. Numerical Verification of Analytical Solution

### 4.1. Dynamic Time History Methodology and Parameter

In previous studies [31–34], validation works always take the numerical results by the quasi-static FE model as the reference results to assess these analytical solution. The quasi-static methodology keeps consistent with the assumptions of these analytical solution, and it does not suffer from the inherent limitations of the forced-based approach (e.g., sensitivity on the dimensions of the model). In practice, the seismic response of underground tunnels under earthquake waves is a dynamic process. Therefore, in this study, the dynamic time history methodology is employed to verify the proposed analytical solution, which can take account of the radiation damping effect and simulates the seismic wave’s propagation and scattering in complex environment well.A two-dimensional finite element model of the ground-tunnel system was established by using the general finite element software ABAQUS [35]. The finite element domain has a width W of 200 m and a height H of 200 m, as illustrated in Figure 7. The tunnel is embedded in domain with an outer radius R of 3 m. The mechanical parameters of the tunnel include density ρl = 2500 kg/m3, Young’s modulus El = 30 GPa, and Poisson’s ratio Ui = 0.2. In the study, the effect of the lining thickness t, buried depth h of the tunnel, ground conditions, and obliquely incident angles α of P waves on the accuracy of proposed solution will be discussed, which will be introduced, respectively, for each verification case in the following sections.Figure 7
finite element model with viscous-spring artificial boundary.Both the ground and the tunnel are simulated by linear elastic solid elements, and the finite element model is discretized according to the accuracy of numerical simulation. The mesh size of the element satisfies the requirement that the mesh size is less than 1/6∼1/8 minimum wavelength in dynamic simulation. At the ground-tunnel interface, the tangential behavior of the interface is modelled by setting equivalent springs to connect the two points on the contact surface for the arbitrary slip state. In terms of the normal contact behavior, the potential separation between the surrounding ground and the tunnel lining is not allowed.The viscous-spring artificial boundary [36] was set at the truncated boundary of the ground to consider the dissipation effect of energy, which can be realized by establishing two pairs of dashpots and springs at the boundary node, i.e., one in the normal direction of the boundary plane and the other in the tangential directions. The spring and damping element parameters of the viscous-spring artificial boundary have the expression as [36](39)KN=Al11+Aλ+2Gr,CN=AlBρcp,KT=Al11+AGr,CT=AlBρcs,where subscripts N and T represent the normal and tangent directions of the truncated boundary, respectively, r can be simply taken as the distance between the geometric center of the structure and the artificial boundary node l, λ, G, and ρ represent the Lame constant, shear modulus, and density of the ground, respectively, cs is the shear wave velocity, cp is the compression wave velocity, Al is the truncated boundary area of all elements including the boundary node l, A and B are correction coefficients, and the recommended values are 0.8 and 1.1 [36].Kobe seismic record of Takatori Station in 1995 Hanshin earthquake is selected as the incident earthquake motions, whose acceleration-time history is shown in Figure8. The earthquake motions are transformed into nodal forces applied at the boundary nodes. The seismic-induced boundary nodal force can be expressed as(40)FBl=Alσlf+KBlulf+CBlu˙lf,where ulf, u˙lf, and σlf are the free-field displacement, velocity, and stress on the artificial boundary node l (xl, yl), respectively. KB and CB are the spring and damping matrices of the viscous-spring boundary. In this paper, the calculation method for the boundary nodal forces induced by obliquely incident plane S waves and P waves is referred to previous works by Zhao et al. [37] and Huang et al. [15, 22]. The physical parameters of the ground and the tunnel are specified in each of the following sections.Figure 8
Acceleration-time history of the Kobe seismic record of Takatori Station.


### 4.2. Verification of Solutions under SV and P Waves in Full-Space

In this study, the proposed analytical solution for tunnels in half-space under obliquely incident P waves is the superposition of the solution for the deep tunnel in full-space under P waves and SV waves. Therefore, the accuracy of the solution for the deep tunnel in full-space under P waves and SV waves is firstly verified. The mechanical parameters of the ground include the densityρg = 2500 kg/m3, Young’s modulus Eg = 4.5 GPa, and Poisson’s ratio νg = 0.25. The effect of the lining thickness on the accuracy of proposed solution is studied here, with values of 0.3 m, 0.6 m, and 0.9 m, respectively. Moreover, different contact conditions at the ground-tunnel interface are also considered. The stiffness coefficient takes values of 1e16, 1e8, and 0, corresponding to the contact conditions of no-slip, middle-slip (midslip), and full-slip, respectively. The Kobe seismic record shown in Figure 8 is input from the model’s bottom as the shear wave and compression wave, respectively.To comprehensively consider the relative error between the analytical solution results and the numerical results at each data point, the errorE is defined as follows:(41)E=X−YY22=∑i=1nxi−yi2∑i=1nyi2,where xi represents the result of the analytical solution at the ith monitoring point of the lining cross section, yi represents the numerical result at the ith monitoring point, and n represents the number of monitoring points.Figures9 and 10, respectively, show the distributions of the peak thrust Tmax and the peak moment Mmax along the cross section of the tunnel under SV waves and P waves using the new analytical solution, the widely used Bobet’s solution [11, 12], and the dynamic numerical solution. As Bobet’s solution can only consider the no-slip and full-slip conditions, the result of Bobet’s solution at midslip is not offered.Figure 9
Distribution of (a) peak thrustTmax and (b) peak moment Mmax in the full-space tunnel under SV waves considering different contact conditions and lining thickness.
(a)(b)Figure 10
Distribution of (a) peak thrustTmax and (b) peak moment Mmax of the tunnel in full-space under P waves considering different contact conditions and lining thickness.
(a)(b)From Figures9 and 10, it can be found that the thrust and the moment have the same distribution form along the cross section of the tunnel. For SV waves, the thrust and the moment get their maximum at the spandrel and the arch foot. However, while for P waves, these two internal forces get their maximum value at the haunches. Whether for P or SV waves, the new analytical solution results are in good agreement with the numerical results for different contact conditions and different lining thicknesses. Although, in some cases, the proposed analytical solution and Bobet’s solution are very close to each other, in many cases, the proposed analytical solution performs far better than Bobet’s solution, especially for tunnels with very thick lining. This is why the tunnel is considered a thick-wall cylinder rather than a beam or shell, which makes the accuracy of the proposed analytical solution not affected by the lining thickness.


### 4.3. Verification of Solutions under Obliquely Incident P Waves in Half-Space

This section verifies the accuracy and applicability of the new analytical solution for the internal forces of the tunnel under obliquely incident P waves in a half-space. In order to study the influence of different ground conditions on the accuracy of the new analytical solution, three ground conditions of soft soil, medium stiff soil, and rock were selected by referring to Seismic Design of Urban Rail Transit Structures (GB50909-2014) [38]. The ground physical parameters are shown in Table 1. The outer radius of the tunnel is R = 3 m, the thickness t is 0.3 m, and the physical parameters of the lining are the same as Section 4.1. The contact condition between the ground and the tunnel is set as no-slip condition. The Kobe seismic record, as shown in Figure 8, is used as the compression wave, which is inputted with the angle of α = 0°, 30°, 45°, and 60°, respectively.Table 1
Physical parameters of the ground.
Ground casesDensityρg (kg/m3)Young’s modulusEg (MPa)Poisson’s ratioνgShear velocitycs (m/s)Case 1: soft soil (cs ≤ 150 m/s)2000540.35100Case 2: medium stiff soil (250 m/s <cs ≤ 500 m/s)23009000.3388Case 3: rock (800 m/s <cs)250045000.25848.5Figure11 shows the distribution of peak thrust Tmax and peak moment Mmax along the tunnel’s cross section, with considering different incident angles and ground conditions. It can be seen from Figure 11 that the peak thrust Tmax and peak moment Mmax by the new analytical solution and numerical solution under obliquely incident P waves are very close to each other. According to the above analysis, the new analytical solution can accurately predict the thrust and the moment of the tunnel under obliquely incident P waves and is suitable for different ground and different lining thickness.Figure 11
Distribution of (a) peak thrustTmax and (b) peak moment Mmax of the tunnel in half-space under obliquely incident P waves with considering different incident angles and ground conditions.
(a)(b)The analytical solution for the internal force of the tunnel in the half-space under the obliquely incident P waves is obtained by superimposing the analytical solution for the internal force of the tunnel in the full-space. Moreover, Park et al. [8] and Bobet [11] pointed out that the analytical solution is based on the condition that the buried depth of the tunnel is greater than 2R. Therefore, it is necessary to discuss the applicability of the new analytical solution under different buried depths. The dynamic time history analysis Method described in Section 4.1 is used as the exact solution. The buried depth h of the model is set to 6 m, 9 m, 12 m, 15 m, 20 m, 30 m, 40 m, and 50 m, and the incident angles α of the P waves are set to 0° and 30°. The errors E (%) of the new analytical solution for the thrust and moment under different depth ratios h/R and different ground conditions are shown in Table 2.Table 2
Error indexE (%) for tunnel thrust and moment in half‐space under obliquely incident P waves considering different buried depths and ground conditions.
h/RCase 1 (soft soil)Case 2 (medium stiff soil)Case 3 (rock)0°300°30°0°30°TMTMTMTMTMTM2.031.5936.4915.617.4427.5527.5511.810.7327.5518.567.7613.253.016.5422.6811.457.2515.4015.404.795.5615.4014.735.326.364.010.7413.718.817.0510.6410.643.844.5310.649.804.505.755.08.7711.397.546.778.098.093.364.308.099.694.374.826.76.147.296.386.105.685.682.954.135.687.494.254.7510.03.273.255.125.643.643.642.673.853.645.054.204.6613.31.892.153.995.052.612.612.553.622.613.714.113.7816.71.281.822.584.641.921.922.493.551.923.053.973.48It can be obtained from Table2 that the error of analytical solution decreases with the increase of h/R, which shows that the influence of free surface boundary condition at the top of the ground on the tunnel is weakened. When the buried depth is 4 times greater than the tunnel radius, the error of the analytical solution is less than 15%.Figure12 shows the effect of buried depth on thrust and moment of tunnels. One can clearly observe that, as the buried depth increases, there is a most dangerous buried depth which causes the maximum internal force of the tunnel. Moreover, as the ground changes from soft ground to hard ground (Case 1 to Case 3), the most dangerous buried depth increases. It can be noted that even for soft soil ground, the most dangerous buried depth is more than 50 m, which is 16 times greater than the tunnel radius. Therefore, the analytical solution’s accuracy at the most dangerous buried depth can meet the engineering requirements, as shown in Table 2.Figure 12
Effect of buried depth on thrust and moment of tunnel. (a) P waves’ incident angle is 0°. (b) P waves’ incident angle is 30°.
(a)(b)


## 4.1. Dynamic Time History Methodology and Parameter

In previous studies [31–34], validation works always take the numerical results by the quasi-static FE model as the reference results to assess these analytical solution. The quasi-static methodology keeps consistent with the assumptions of these analytical solution, and it does not suffer from the inherent limitations of the forced-based approach (e.g., sensitivity on the dimensions of the model). In practice, the seismic response of underground tunnels under earthquake waves is a dynamic process. Therefore, in this study, the dynamic time history methodology is employed to verify the proposed analytical solution, which can take account of the radiation damping effect and simulates the seismic wave’s propagation and scattering in complex environment well.A two-dimensional finite element model of the ground-tunnel system was established by using the general finite element software ABAQUS [35]. The finite element domain has a width W of 200 m and a height H of 200 m, as illustrated in Figure 7. The tunnel is embedded in domain with an outer radius R of 3 m. The mechanical parameters of the tunnel include density ρl = 2500 kg/m3, Young’s modulus El = 30 GPa, and Poisson’s ratio Ui = 0.2. In the study, the effect of the lining thickness t, buried depth h of the tunnel, ground conditions, and obliquely incident angles α of P waves on the accuracy of proposed solution will be discussed, which will be introduced, respectively, for each verification case in the following sections.Figure 7
finite element model with viscous-spring artificial boundary.Both the ground and the tunnel are simulated by linear elastic solid elements, and the finite element model is discretized according to the accuracy of numerical simulation. The mesh size of the element satisfies the requirement that the mesh size is less than 1/6∼1/8 minimum wavelength in dynamic simulation. At the ground-tunnel interface, the tangential behavior of the interface is modelled by setting equivalent springs to connect the two points on the contact surface for the arbitrary slip state. In terms of the normal contact behavior, the potential separation between the surrounding ground and the tunnel lining is not allowed.The viscous-spring artificial boundary [36] was set at the truncated boundary of the ground to consider the dissipation effect of energy, which can be realized by establishing two pairs of dashpots and springs at the boundary node, i.e., one in the normal direction of the boundary plane and the other in the tangential directions. The spring and damping element parameters of the viscous-spring artificial boundary have the expression as [36](39)KN=Al11+Aλ+2Gr,CN=AlBρcp,KT=Al11+AGr,CT=AlBρcs,where subscripts N and T represent the normal and tangent directions of the truncated boundary, respectively, r can be simply taken as the distance between the geometric center of the structure and the artificial boundary node l, λ, G, and ρ represent the Lame constant, shear modulus, and density of the ground, respectively, cs is the shear wave velocity, cp is the compression wave velocity, Al is the truncated boundary area of all elements including the boundary node l, A and B are correction coefficients, and the recommended values are 0.8 and 1.1 [36].Kobe seismic record of Takatori Station in 1995 Hanshin earthquake is selected as the incident earthquake motions, whose acceleration-time history is shown in Figure8. The earthquake motions are transformed into nodal forces applied at the boundary nodes. The seismic-induced boundary nodal force can be expressed as(40)FBl=Alσlf+KBlulf+CBlu˙lf,where ulf, u˙lf, and σlf are the free-field displacement, velocity, and stress on the artificial boundary node l (xl, yl), respectively. KB and CB are the spring and damping matrices of the viscous-spring boundary. In this paper, the calculation method for the boundary nodal forces induced by obliquely incident plane S waves and P waves is referred to previous works by Zhao et al. [37] and Huang et al. [15, 22]. The physical parameters of the ground and the tunnel are specified in each of the following sections.Figure 8
Acceleration-time history of the Kobe seismic record of Takatori Station.


## 4.2. Verification of Solutions under SV and P Waves in Full-Space

In this study, the proposed analytical solution for tunnels in half-space under obliquely incident P waves is the superposition of the solution for the deep tunnel in full-space under P waves and SV waves. Therefore, the accuracy of the solution for the deep tunnel in full-space under P waves and SV waves is firstly verified. The mechanical parameters of the ground include the densityρg = 2500 kg/m3, Young’s modulus Eg = 4.5 GPa, and Poisson’s ratio νg = 0.25. The effect of the lining thickness on the accuracy of proposed solution is studied here, with values of 0.3 m, 0.6 m, and 0.9 m, respectively. Moreover, different contact conditions at the ground-tunnel interface are also considered. The stiffness coefficient takes values of 1e16, 1e8, and 0, corresponding to the contact conditions of no-slip, middle-slip (midslip), and full-slip, respectively. The Kobe seismic record shown in Figure 8 is input from the model’s bottom as the shear wave and compression wave, respectively.To comprehensively consider the relative error between the analytical solution results and the numerical results at each data point, the errorE is defined as follows:(41)E=X−YY22=∑i=1nxi−yi2∑i=1nyi2,where xi represents the result of the analytical solution at the ith monitoring point of the lining cross section, yi represents the numerical result at the ith monitoring point, and n represents the number of monitoring points.Figures9 and 10, respectively, show the distributions of the peak thrust Tmax and the peak moment Mmax along the cross section of the tunnel under SV waves and P waves using the new analytical solution, the widely used Bobet’s solution [11, 12], and the dynamic numerical solution. As Bobet’s solution can only consider the no-slip and full-slip conditions, the result of Bobet’s solution at midslip is not offered.Figure 9
Distribution of (a) peak thrustTmax and (b) peak moment Mmax in the full-space tunnel under SV waves considering different contact conditions and lining thickness.
(a)(b)Figure 10
Distribution of (a) peak thrustTmax and (b) peak moment Mmax of the tunnel in full-space under P waves considering different contact conditions and lining thickness.
(a)(b)From Figures9 and 10, it can be found that the thrust and the moment have the same distribution form along the cross section of the tunnel. For SV waves, the thrust and the moment get their maximum at the spandrel and the arch foot. However, while for P waves, these two internal forces get their maximum value at the haunches. Whether for P or SV waves, the new analytical solution results are in good agreement with the numerical results for different contact conditions and different lining thicknesses. Although, in some cases, the proposed analytical solution and Bobet’s solution are very close to each other, in many cases, the proposed analytical solution performs far better than Bobet’s solution, especially for tunnels with very thick lining. This is why the tunnel is considered a thick-wall cylinder rather than a beam or shell, which makes the accuracy of the proposed analytical solution not affected by the lining thickness.


## 4.3. Verification of Solutions under Obliquely Incident P Waves in Half-Space

This section verifies the accuracy and applicability of the new analytical solution for the internal forces of the tunnel under obliquely incident P waves in a half-space. In order to study the influence of different ground conditions on the accuracy of the new analytical solution, three ground conditions of soft soil, medium stiff soil, and rock were selected by referring to Seismic Design of Urban Rail Transit Structures (GB50909-2014) [38]. The ground physical parameters are shown in Table 1. The outer radius of the tunnel is R = 3 m, the thickness t is 0.3 m, and the physical parameters of the lining are the same as Section 4.1. The contact condition between the ground and the tunnel is set as no-slip condition. The Kobe seismic record, as shown in Figure 8, is used as the compression wave, which is inputted with the angle of α = 0°, 30°, 45°, and 60°, respectively.Table 1
Physical parameters of the ground.
Ground casesDensityρg (kg/m3)Young’s modulusEg (MPa)Poisson’s ratioνgShear velocitycs (m/s)Case 1: soft soil (cs ≤ 150 m/s)2000540.35100Case 2: medium stiff soil (250 m/s <cs ≤ 500 m/s)23009000.3388Case 3: rock (800 m/s <cs)250045000.25848.5Figure11 shows the distribution of peak thrust Tmax and peak moment Mmax along the tunnel’s cross section, with considering different incident angles and ground conditions. It can be seen from Figure 11 that the peak thrust Tmax and peak moment Mmax by the new analytical solution and numerical solution under obliquely incident P waves are very close to each other. According to the above analysis, the new analytical solution can accurately predict the thrust and the moment of the tunnel under obliquely incident P waves and is suitable for different ground and different lining thickness.Figure 11
Distribution of (a) peak thrustTmax and (b) peak moment Mmax of the tunnel in half-space under obliquely incident P waves with considering different incident angles and ground conditions.
(a)(b)The analytical solution for the internal force of the tunnel in the half-space under the obliquely incident P waves is obtained by superimposing the analytical solution for the internal force of the tunnel in the full-space. Moreover, Park et al. [8] and Bobet [11] pointed out that the analytical solution is based on the condition that the buried depth of the tunnel is greater than 2R. Therefore, it is necessary to discuss the applicability of the new analytical solution under different buried depths. The dynamic time history analysis Method described in Section 4.1 is used as the exact solution. The buried depth h of the model is set to 6 m, 9 m, 12 m, 15 m, 20 m, 30 m, 40 m, and 50 m, and the incident angles α of the P waves are set to 0° and 30°. The errors E (%) of the new analytical solution for the thrust and moment under different depth ratios h/R and different ground conditions are shown in Table 2.Table 2
Error indexE (%) for tunnel thrust and moment in half‐space under obliquely incident P waves considering different buried depths and ground conditions.
h/RCase 1 (soft soil)Case 2 (medium stiff soil)Case 3 (rock)0°300°30°0°30°TMTMTMTMTMTM2.031.5936.4915.617.4427.5527.5511.810.7327.5518.567.7613.253.016.5422.6811.457.2515.4015.404.795.5615.4014.735.326.364.010.7413.718.817.0510.6410.643.844.5310.649.804.505.755.08.7711.397.546.778.098.093.364.308.099.694.374.826.76.147.296.386.105.685.682.954.135.687.494.254.7510.03.273.255.125.643.643.642.673.853.645.054.204.6613.31.892.153.995.052.612.612.553.622.613.714.113.7816.71.281.822.584.641.921.922.493.551.923.053.973.48It can be obtained from Table2 that the error of analytical solution decreases with the increase of h/R, which shows that the influence of free surface boundary condition at the top of the ground on the tunnel is weakened. When the buried depth is 4 times greater than the tunnel radius, the error of the analytical solution is less than 15%.Figure12 shows the effect of buried depth on thrust and moment of tunnels. One can clearly observe that, as the buried depth increases, there is a most dangerous buried depth which causes the maximum internal force of the tunnel. Moreover, as the ground changes from soft ground to hard ground (Case 1 to Case 3), the most dangerous buried depth increases. It can be noted that even for soft soil ground, the most dangerous buried depth is more than 50 m, which is 16 times greater than the tunnel radius. Therefore, the analytical solution’s accuracy at the most dangerous buried depth can meet the engineering requirements, as shown in Table 2.Figure 12
Effect of buried depth on thrust and moment of tunnel. (a) P waves’ incident angle is 0°. (b) P waves’ incident angle is 30°.
(a)(b)


## 5. Discussion of Critical Factors

### 5.1. Effect of Incident Angle of P Waves

In this section, the effect of the incident angle of P waves on the seismic response of the tunnel is discussed, considering different ground conditions. The incident anglesα takes values of 0°, 15°, 30°, 45°, 60°, 75°, and 85°. Three ground cases are selected, whose physical parameters are listed in Table 1. Figures 9 and 10 show that the thickness of the lining has little influence on the distribution of the tunnel’s internal force under the obliquely incident P waves. Referring to the parameters adopted in previous studies [8, 31, 33], the thickness of lining used in the parameter analysis part is 0.3 m. The tunnel’s buried depth h is 50 m. Moreover, no-slip contact condition at the ground-tunnel interface is used. Other physical parameters of the tunnel are the same as those in Section 4.3.Figure13 shows the distribution of the peak thrust and peak moment along the tunnel’s cross section considering different incident angles. To keep the figure concise, only 0°, 30°, 60°, and 85° are plotted. It can be found that the incident angle of P waves has great influence on the distribution of the thrust and moment along the tunnel’s cross section. In detail, for the vertically incident case with α = 0°, the tunnel haunches suffer the largest thrust and moment. With increasing the incident angle α from 0° to 85°, the critical position for the thrust and moment has the trend of moving from the haunches to the vault and arch bottom. The larger the incident angle is, the closer to the vault and arch bottom the critical position is. Moreover, the ground conditions also have effect on the distribution of the thrust and moment of tunnels under obliquely incident angle. For a given oblique angle, the critical position is getting closer to the vault and arch bottom, as the ground changing from soft ground (ground Case 1) to hard ground (ground Case 3).Figure 13
Distribution of peak thrustTmax and peak moment Mmax considering different incident angles and ground conditions.To evaluate the effect of the incident angle on the amplitude of the thrust and the moment of the tunnel, the maximum of the thrust and moment along the tunnel’s cross section is plotted in Figure14, considering different incident angles and ground conditions. From Figure 14, one can find that the maximum thrust that the tunnel suffers decrease with the increase of the incident angle for soft ground (ground Case 1). However, for hard ground (ground Case 2 and Case 3), the maximum thrust increases first and decreases with the incidence angle subsequently. The maximum moment increases first and then decreases with the incident angle for all ground conditions. From the analysis above, one can find that there is a most unfavorable incident angle where the tunnel suffers the largest thrust or largest moment. Moreover, the most unfavorable incident angle becomes larger as the ground changes from soft ground to hard ground.Figure 14
Extreme values of thrust and moment of tunnel considering different incident angles and ground conditions.


### 5.2. Effect of the Tangential Contact Stiffness at the Ground-Tunnel Interface

In this study, the stiffness coefficientk is introduced to simulate the tangential contact stiffness at the ground-tunnel interface. In this section, the influence of the stiffness coefficient k on the normalized thrust and moment is studied. The physical parameters and dimensions of the ground and tunnel are the same as those of Section 5.1.Figure15 shows the variation of normalized thrust and moment with the stiffness coefficient k under three ground conditions. It can be found that when k changes from full-slip to no-slip, the normalized thrust increases with the increase of k, and the normalized moment decreases with the increase of k. Theoretically, the influence range of k is 0 − ∞. However, it can be seen from Figure 15 that the stiffness coefficient has an apparent effective influence range, but outside the effective influence range, it has little influence on the thrust and moment. Moreover, the influence range of k is different under different ground conditions. As the ground changes from soft to hard, the influence range of the stiffness coefficient gradually expands.Figure 15
Effect ofk on normalized thrust and moment distribution under different ground conditions.


### 5.3. Effect of the Relative Stiffness between the Ground and the Tunnel

As we all know, the seismic response of the underground structures depends on the ground deformation and the relative stiffness between the ground and the structure. In this section, the effect of the relative stiffness between the ground and the tunnel on the seismic response of the tunnel is studied. The flexibility ratioF proposed by Peck et al. [39] is employed to represent the relative stiffness, which was defined as(42)F=Eg1−υl2R36ElI1+υg,where Eg and νg, and El and νl, are Young’s modulus and Poisson’s ratio of the ground and tunnel, respectively, and I and R are the moment of inertia and the outer radius of the tunnel, respectively.Values of the flexibility ratioF with values over 1 suggest that the tunnel lining has the lower stiffness with the surrounding ground. The tunnel lining and incident earthquake intensity are fixed, while Young’s modulus of the ground is varied parametrically for capturing a wide range of flexibility ratios F. In other words, the seismic response of a tunnel under a certain earthquake exciting is studied, which is located under different ground conditions. The ground parameters are shown in Table 3, and the corresponding flexibility ratio ranges from 0.1 to 1000. Other parameters, such as the size, depth, and physical parameters of the tunnel, keep step with Section 5.1. Moreover, different contact conditions are considered, i.e., no-slip, midslip, and full-slip, respectively.Table 3
Physical parameters of the ground and corresponding flexibility ratio.
Density,ρ (kg/m3)Poisson’s ratio,νYoung’s modulus,Eg (MPa)Flexibility ratio,F25000.251.9530.125000.257.8130.425000.2519.531125000.2548.8282.525000.25117.188625000.25195.3131025000.25585.9373025000.251367.1887025000.251953.12510025000.253906.25020025000.259765.62550025000.2519531.251000Figure16 shows the effect of the flexibility ratio on the peak thrust and peak moment of the tunnel considering different contact conditions. With the increase of the flexibility ratio, as shown in Figure 16, the internal forces of the tunnel lining increase first and then decrease. That is to say, there is a most unfavorable flexibility ratio where the tunnel suffers the largest thrust or the largest moment. It should be noticed that the most unfavorable flexibility ratio for the thrust is larger than that for the moment. Moreover, the contact conditions have little effect on the value of the most unfavorable flexibility ratio.Figure 16
Peak thrust and peak moment of tunnel considering different contact conditions and flexibility ratio.


## 5.1. Effect of Incident Angle of P Waves

In this section, the effect of the incident angle of P waves on the seismic response of the tunnel is discussed, considering different ground conditions. The incident anglesα takes values of 0°, 15°, 30°, 45°, 60°, 75°, and 85°. Three ground cases are selected, whose physical parameters are listed in Table 1. Figures 9 and 10 show that the thickness of the lining has little influence on the distribution of the tunnel’s internal force under the obliquely incident P waves. Referring to the parameters adopted in previous studies [8, 31, 33], the thickness of lining used in the parameter analysis part is 0.3 m. The tunnel’s buried depth h is 50 m. Moreover, no-slip contact condition at the ground-tunnel interface is used. Other physical parameters of the tunnel are the same as those in Section 4.3.Figure13 shows the distribution of the peak thrust and peak moment along the tunnel’s cross section considering different incident angles. To keep the figure concise, only 0°, 30°, 60°, and 85° are plotted. It can be found that the incident angle of P waves has great influence on the distribution of the thrust and moment along the tunnel’s cross section. In detail, for the vertically incident case with α = 0°, the tunnel haunches suffer the largest thrust and moment. With increasing the incident angle α from 0° to 85°, the critical position for the thrust and moment has the trend of moving from the haunches to the vault and arch bottom. The larger the incident angle is, the closer to the vault and arch bottom the critical position is. Moreover, the ground conditions also have effect on the distribution of the thrust and moment of tunnels under obliquely incident angle. For a given oblique angle, the critical position is getting closer to the vault and arch bottom, as the ground changing from soft ground (ground Case 1) to hard ground (ground Case 3).Figure 13
Distribution of peak thrustTmax and peak moment Mmax considering different incident angles and ground conditions.To evaluate the effect of the incident angle on the amplitude of the thrust and the moment of the tunnel, the maximum of the thrust and moment along the tunnel’s cross section is plotted in Figure14, considering different incident angles and ground conditions. From Figure 14, one can find that the maximum thrust that the tunnel suffers decrease with the increase of the incident angle for soft ground (ground Case 1). However, for hard ground (ground Case 2 and Case 3), the maximum thrust increases first and decreases with the incidence angle subsequently. The maximum moment increases first and then decreases with the incident angle for all ground conditions. From the analysis above, one can find that there is a most unfavorable incident angle where the tunnel suffers the largest thrust or largest moment. Moreover, the most unfavorable incident angle becomes larger as the ground changes from soft ground to hard ground.Figure 14
Extreme values of thrust and moment of tunnel considering different incident angles and ground conditions.


## 5.2. Effect of the Tangential Contact Stiffness at the Ground-Tunnel Interface

In this study, the stiffness coefficientk is introduced to simulate the tangential contact stiffness at the ground-tunnel interface. In this section, the influence of the stiffness coefficient k on the normalized thrust and moment is studied. The physical parameters and dimensions of the ground and tunnel are the same as those of Section 5.1.Figure15 shows the variation of normalized thrust and moment with the stiffness coefficient k under three ground conditions. It can be found that when k changes from full-slip to no-slip, the normalized thrust increases with the increase of k, and the normalized moment decreases with the increase of k. Theoretically, the influence range of k is 0 − ∞. However, it can be seen from Figure 15 that the stiffness coefficient has an apparent effective influence range, but outside the effective influence range, it has little influence on the thrust and moment. Moreover, the influence range of k is different under different ground conditions. As the ground changes from soft to hard, the influence range of the stiffness coefficient gradually expands.Figure 15
Effect ofk on normalized thrust and moment distribution under different ground conditions.


## 5.3. Effect of the Relative Stiffness between the Ground and the Tunnel

As we all know, the seismic response of the underground structures depends on the ground deformation and the relative stiffness between the ground and the structure. In this section, the effect of the relative stiffness between the ground and the tunnel on the seismic response of the tunnel is studied. The flexibility ratioF proposed by Peck et al. [39] is employed to represent the relative stiffness, which was defined as(42)F=Eg1−υl2R36ElI1+υg,where Eg and νg, and El and νl, are Young’s modulus and Poisson’s ratio of the ground and tunnel, respectively, and I and R are the moment of inertia and the outer radius of the tunnel, respectively.Values of the flexibility ratioF with values over 1 suggest that the tunnel lining has the lower stiffness with the surrounding ground. The tunnel lining and incident earthquake intensity are fixed, while Young’s modulus of the ground is varied parametrically for capturing a wide range of flexibility ratios F. In other words, the seismic response of a tunnel under a certain earthquake exciting is studied, which is located under different ground conditions. The ground parameters are shown in Table 3, and the corresponding flexibility ratio ranges from 0.1 to 1000. Other parameters, such as the size, depth, and physical parameters of the tunnel, keep step with Section 5.1. Moreover, different contact conditions are considered, i.e., no-slip, midslip, and full-slip, respectively.Table 3
Physical parameters of the ground and corresponding flexibility ratio.
Density,ρ (kg/m3)Poisson’s ratio,νYoung’s modulus,Eg (MPa)Flexibility ratio,F25000.251.9530.125000.257.8130.425000.2519.531125000.2548.8282.525000.25117.188625000.25195.3131025000.25585.9373025000.251367.1887025000.251953.12510025000.253906.25020025000.259765.62550025000.2519531.251000Figure16 shows the effect of the flexibility ratio on the peak thrust and peak moment of the tunnel considering different contact conditions. With the increase of the flexibility ratio, as shown in Figure 16, the internal forces of the tunnel lining increase first and then decrease. That is to say, there is a most unfavorable flexibility ratio where the tunnel suffers the largest thrust or the largest moment. It should be noticed that the most unfavorable flexibility ratio for the thrust is larger than that for the moment. Moreover, the contact conditions have little effect on the value of the most unfavorable flexibility ratio.Figure 16
Peak thrust and peak moment of tunnel considering different contact conditions and flexibility ratio.


## 6. Conclusions

In this study, an analytical solution for the seismic-induced internal forces of the circular tunnel in half-space, which is subjected to obliquely incident P waves, is developed. The new analytical solution is different from other previous analytical solution in three respects. (1) Previous analytical solutions were mainly developed for the tunnel’s internal force under vertically incident shear waves, while the new analytical solution is for the tunnel’s internal force under obliquely incident P waves in half-space. (2) Different from previous analytical solutions that regarded the tunnel lining as a beam or a shell, the new solution treats the tunnel lining as a thick-wall cylinder, which can provide more precise predictions than beam or shell models used in previous analytical solutions, especially for tunnels with thick lining. (3) Most of previous solution only consider the full-slip and no-slip contact states at the ground-tunnel interface. In this study, a spring-type stiffness coefficient is introduced to consider the arbitrary contact state of the ground-tunnel interface.The reliability of the proposed analytical solution is assessed by comparing with the dynamic numerical results. The results demonstrate that the proposed analytical solution performs well and can be adopted to predict the internal forces of circular tunnel under obliquely incident P waves in seismic design. Based on the proposed analytical solution, parametrical studies are conducted to investigate the effect of some critical parameters on the tunnel’s response, including the incident angles, the tangential contact stiffness at the ground-tunnel interface, and the relative stiffness between the ground and the tunnel. Some conclusions can be drawn as follows:(1)
With the increasing of the incident angle, the critical position for thrust and moment tend to move from the haunches to the vault and arch bottom. For a given oblique angle, the critical position is getting closer to the vault and arch bottom, as the ground changing from soft to hard. Moreover, there is a most unfavorable incident angle where the tunnel suffers the largest thrust or moment, which becomes being larger with the ground changing from soft to hard.(2)
There is a limited range for the tangential contact stiffness at the ground-tunnel interface, where the tangential contact stiffness has a significantly increasing and decreasing effect on the tunnel’s thrust and moment, respectively. Beyond the range, the tunnel responses almost keep constant. Moreover, as the ground varying from soft to hard, the limited range for the tangential contact stiffness gradually expands.(3)
With the increasing of the relative stiffness between the ground and the tunnel, the thrust and the moment of the tunnel increase first and then decrease. There is a most unfavorable relative stiffness where the tunnel suffers the largest thrust or moment. It should be noticed that the most unfavorable relative stiffness for the thrust is larger than that of the moment.



---

*Source: 1946184-2021-12-22.xml*