# Effects of Oblique Incidence of SV Waves on Nonlinear Seismic Response of a Lined Arched Tunnel

Authors: Dunyu Lyu; Sha Ma; Chu Yu; Congcong Liu; Xiaowei Wang; Buyun Yang; Ming Xiao
Journal: Shock and Vibration
 (2020)

---

## Abstract

The incident direction of earthquake motion is an important factor affecting the seismic response of underground structures. In this study, a three-dimensional (3D) oblique incidence method of SV waves is proposed and the effects of incident angles of SV waves on the seismic response of a lined arched tunnel are evaluated. Based on wave field decomposition principle and equivalent node force method and together with viscous-spring artificial boundary, the oblique incidence method of SV waves is implemented by transforming seismic wave field into the equivalent nodal forces acting on the artificial boundaries. By deriving the distance of the incident waves and the reflected wave on free surface to artificial boundaries, this method can comprehensively consider the phase difference of the seismic wave propagation and the influence of the damping effect of the rock medium on the seismic wave propagation. The method is programed into a dynamic finite element program and its effectiveness is examined by a numerical example. Consequently, the oblique incidence method is applied to evaluate the seismic behaviors of the tunnel. The numerical results reveal that (1) the oblique incidence of the seismic wave results in a larger seismic response; (2) the response amplitudes of the stress and displacement increase with the increase of incident angles and reaching the maximum in the case of 30° incident angle; (3) the damage extent increases with an increase in the incident angles, and the oblique incidence of the seismic wave is believed to increase the spatial difference of damage distribution.


---

## Body

## 1. Introduction

Due to the constraints of surrounding rocks or soil, underground structures are well known to suffer less damage induced by earthquakes than surface structures; thus, underground structures are generally believed to be strongly resistant to earthquakes unless they are crossed by active faults [1, 2]. However, many instances of noticeable seismic damage of underground structures were reported during several strong earthquake events, including the Kobe earthquake in Japan [3, 4], the Kocaeli earthquake in Turkey [5, 6], the Chi-Chi earthquake in Taiwan [7, 8], and the Wenchuan earthquake in China [9–11]. These earthquake disasters provide sufficient evidences to demonstrate that the stability of underground structures located in seismically active areas is still an important and emergent issue.The numerical method is most widely used to assess the seismic behaviors of underground structures because it is of low cost and allows complex nonlinear simulations, including material and geometric nonlinearities [12]. Great efforts have been spent by many scholars to study the seismic performances of underground structures using the numerical method. Liu and Song [13] analyzed the effects of vertical earthquake motions and the buried depth on the seismic responses of underground structures in liquefiable soils. Penzien [14] proposed an analytical procedure to evaluate the racking deformation of rectangular and circular tunnel linings during a seismic event. A finite element method in the time domain was adopted by Hatzigeorgiou and Beskos [15] to investigate the impacts of seismic soil-structure interaction on 3D-lined tunnels. Ding et al. [16] conducted a large-scale 3D nonlinear dynamic numerical simulation for an immersed tunnel, and in this simulation, both the nonlinear material behavior and soil-tunnel interaction were considered. A modified discontinuous deformation analysis (DDA) method was presented by Zhang et al. [17] to study the seismic dynamic response of large underground caverns. Cui et al. [18] adopted the 3DEC (three-dimensional distinct element) method to assess the control effects of a large geological discontinuity on the seismic stability of an underground powerhouse cavern. Dowding et al. [19] developed a block model, which was modified to operate in a parallel mode, to investigate the dynamic stability of heavily jointed rock. Those works can effectively enhance the current seismic design, but most of them assume the seismic waves propagating along a vertical direction. Actually, the seismic waves usually have oblique incident directions when underground structures are close to the epicenters [20], and the statistics data of recent strong earthquake records suggest that the incident angle of an actual earthquake motion at the rock site is approximately 60° [21, 22]. The oblique incidence seismic waves are proved to contribute to the spatial nonuniform deformation of large-scale underground structures [23–25]. This is especially true for the large-scale underground structures [26], including tunnels and underground powerhouse of hydropower station. Many scholars have realized the evident impacts of seismic motion incident angles on underground structures. Heymsfield [27] studied the effects of a sloping bedrock half-space on the amplification of shear waves under different incident angles. The analytical expressions for the strains of cylindrical underground structures under random incident angles of seismic waves are derived by Kouretzis et al. [28] based on the 3D shell theory. The effects of the seismic input angle on the dynamic response of underground gas storage salt cavern were evaluated by Wang et al. [29]. Naggar et al. [30] investigated the effects of the seismic input angle on the moments and thrusts of cylinder tunnel linings by a proposed analytical procedure. An input method of P waves and SV waves was proposed by Huang et al. [24, 31] to analyze the impacts of incident angles on the seismic behaviors of the lined tunnel. More relevant works can be found in [32–37]. Clearly, those works have made several achievements; however, most of them were mainly concentrated on the cylindrical lined tunnel, whereas the study of effects of incident angles on the seismic performance of tunnels having vertical walls and arch roof, to date, remains far from adequate.In this work, the effects of oblique incidence of SV waves on the seismic response process and damage characteristics of a lined arched tunnel are evaluated in detail. The rest of this paper is organized as follows. A brief description of the tunnel, including geological conditions, structural arrangement, and support measures is presented and the principle theories of the adopted finite element model are provided in Section2. In Section 3, a 3D oblique incidence method of SV waves is deduced and its effectiveness is verified by a numerical example. In Section 4, the proposed method is used to assess the effects of the incident angles of SV waves on the nonlinear seismic performance of the tunnel.


## 2. Problem Description and Establishment of the FEM Model

### 2.1. Problem Description

A lined tunnel, which is located in Southwest China, is adopted as an engineering example to assess the effects of incident direction of seismic waves on the seismic behaviors of underground structures. The tunnel is embedded in gray limestone, and according to site survey, there is no large geological discontinuity nearby. The tunnel has vertical walls and an arch roof with dimension of 9.2 m wide and 10.5 m high, and its vertical buried depth is approximately 20 m (see Figure1). The lining is made of reinforced concrete and has a thickness of 0.6 m. What makes this tunnel special is that it located in the strong earthquake activity areas of China. According to the seismic ground motion parameters zonation map of China [38], which is published by the China Earthquake Administration, the seismic intensity of the studied area is expected to be as high as seven degree; thus, it is essential to assess the seismic performance of the tunnel.Figure 1
3D finite element model of the tunnel.


### 2.2. FEM Model and Input Seismic Motion

In this paper, the finite element method (FEM) is adopted to conduct the numerical simulation, and a limited finite element model is intercepted to establish computing model based on the actual size of the tunnel (see Figure1). The longitudinal axis of the tunnel and vertical axis are, respectively, defined as the y-axis and z-axis in the global coordinate system of the computing model, and the x-axis is determined by the right-hand rule. The entire finite element mesh is discretized with eight-node hexahedral elements, and a total of 108500 solid elements are meshed. The maximum mesh size is set as 3 m to ensure the accuracy of dynamic calculation [39]. In general, it is appropriate to adopt measured seismic waves or artificial seismic waves synthesized according to site conditions for dynamic calculation; thus, a 20 s acceleration time history of the artificial seismic wave, which is synthesized based on the site condition and seismic design standard of the tunnel, is employed as the input seismic wave (see Figure 2).Figure 2
Time history (a) and Fourier spectrum (b) of input acceleration.
(a)
(b)


### 2.3. Material Model

#### 2.3.1. Material Model for Concrete

An appropriate constitutive model for concrete materials is the key to simulate their seismic behavior. The mechanical response of concrete materials under seismic cyclic loading is a complex multistage process, during which concrete materials may suffer stiffness degradation. In addition, the mechanical behavior of concrete materials under tension and compression is quite different. This difference is mainly reflected in the different strength of concrete materials under tension and compression. Generally, the crack model [40] and plastic damage model [41] are used to describe the mechanical properties of concrete materials. In this paper, the plastic damage model proposed by Lubliner et al. [41] is employed because the model can reasonably simulate the softening process, unloading stiffness degradation, and repeated loading process of concrete materials [42]. Under multiaxial stress conditions, the stress of concrete σ can be expressed as(1)σ=1−dσ¯=1−dE0:ε−εpl,where E0 is the initial elasticity modulus; ε is the strain tensor, which can be divided to elastic part εe and plastic part εpl; and d is a scalar damage coefficient, which can be obtained according to the function of tensile damage variable dt and the shear damage variable dc[43].(2)d=1−1−stdt1−scdc,0≤st,sc≤1,dt=1−σt/E0εtpl1/bt−1+σt/E0,0≤dt≤1,dc=1−σc/E0εcpl1/bc−1+σc/E0,0≤dc≤1,where subscripts t and c represent the tension state and compressive state; εtpl and εcpl are the plastic strain under tension and compressive loadings; bt and bc are dimensionless constants and are suggested to be 0.7 and 0.1 by Birtel and Mark [43], st and sc are functions of effective stress σ¯^ and are introduced to represent stiffness recovery effects(3)st=1−wtrσ¯^,sc=1−wcrσ¯^,(4)rσ¯^=∑i=13σ¯^i∑i=13σ¯^i,i=1,2,3,where wt and wc are the weight factors that controls the stiffness recovery degree; rσ¯^ is the stress weight factor; σ¯^i is the principle value of the effective stress σ¯; ⋅ denotes the McCauley bracket, that is, σ¯^i=0.5σ¯^i+σ¯^i.In order to describe the mechanical behavior after yielding, the yield functionF modified by Lee and Fenves [42] and the plastic flow potential G expressed by the Drucker–Prager hyperbolic function are adopted:(5)Fσ¯,ε˜pl=11−αq¯−3αp¯+βε˜plσ¯^max−γ−σ¯^max−σ¯cε˜cpl≤0,(6)G=αpI1+2J2,where ε˜pl is equivalent plastic strain; q¯ is the Mises equivalent effective stress; p¯ is the equivalent hydrostatic pressure; α is material constant; σ¯^max is the maximum eigenvalue of the effective stress; the function β can be given as βε˜pl=1−ασ¯cε˜cpl/σ¯tε˜tpl, where σ¯c and σ¯t are the effective tensile and compressive stresses; γ=31−Kc/2Kc−1, where Kc is the parameter related to the yield curve; I1 is the first invariant of the principal stress and J2 is the second invariant of the deviatoric stress; and αp is a dimensionless material parameter.According to the nonassociated potential flow, the plastic strain rateε˙p can be given by(7)ε˙p=λ˙∂Gσ¯∂σ¯,where λ is the plastic parameter, which can be obtained based on the consistency condition; λ˙ is the derivative of λ versus time. Obviously, the mechanical behavior of concrete materials can be reasonably described according to equations (1)–(7). In addition, according to the code for design of concrete structures [44] published by China, the tension and compressive strength of concrete with C25-grade can be given as 1.78 MPa and 16.7 MPa, respectively. Consequently, the stress-strain relationship of concrete with C25-grade is shown as Figure 3.Figure 3
Plastic damage evolution curves for C25-grade concrete: (a) Tension state and (b) compressive state.
(a)
(b)


#### 2.3.2. Analysis Procedure

Generally, seismic response analysis is divided into three steps. In the first step, it is necessary to invert the geostress field of the calculated area according to the limited measured ground stress and geological conditions. However, in many cases, the measured geostress data is hard to obtain; thus, the gravity stress field is adopted as an alternative by many scholars. This alternative is acceptable especially for the cases that the geostress field is mainly controlled by gravity stress. In this paper, the gravity stress is employed as the initial geostress field. In the second step, the static excavation calculation is conducted and the computed stress field is adopted as the initial condition for the seismic response calculation of the third step. In this process, the viscous-spring artificial boundary proposed by Liu et al. [44] is applied to the lateral and bottom boundaries, and the Rayleigh damping with a critical damping ratio of 5% is employed. The mechanical behavior of rock mass under dynamic loadings is modeled by the Mohr–Coulomb failure criterion with a tension cut-off limit. The mechanical behavior of concrete is modeled by the method presented in Section 2.3.1. The geomechanical properties of the rock and concrete units are shown in Table 1. Note that, the time step Δt for the dynamic calculation by the explicit method should meet the below condition [10]:(8)Δt≤κmineleCe0.80≤κ≤0.98,where κ is an empirical constant; Ce is the wave velocity; and le is the minimum size of all elements of the model. In general, time step Δt can be valued by the maximum value according to equation (8) to reduce the amount of calculation.Table 1
Geomechanical properties of the materials used in the numerical solutions.
Material
Density (kg/m3)
Elastic modulus (GPa)
Poisson’s ratio
Friction angle (°)
Cohesion (MPa)
Tensile strength (MPa)
Compressive strength (MPa)
Rock
2.75
9
0.27
42
1.05
1.3
35
Concrete
2.5
25
0.167
46
1.30
1.2
16.7


## 2.1. Problem Description

A lined tunnel, which is located in Southwest China, is adopted as an engineering example to assess the effects of incident direction of seismic waves on the seismic behaviors of underground structures. The tunnel is embedded in gray limestone, and according to site survey, there is no large geological discontinuity nearby. The tunnel has vertical walls and an arch roof with dimension of 9.2 m wide and 10.5 m high, and its vertical buried depth is approximately 20 m (see Figure1). The lining is made of reinforced concrete and has a thickness of 0.6 m. What makes this tunnel special is that it located in the strong earthquake activity areas of China. According to the seismic ground motion parameters zonation map of China [38], which is published by the China Earthquake Administration, the seismic intensity of the studied area is expected to be as high as seven degree; thus, it is essential to assess the seismic performance of the tunnel.Figure 1
3D finite element model of the tunnel.


## 2.2. FEM Model and Input Seismic Motion

In this paper, the finite element method (FEM) is adopted to conduct the numerical simulation, and a limited finite element model is intercepted to establish computing model based on the actual size of the tunnel (see Figure1). The longitudinal axis of the tunnel and vertical axis are, respectively, defined as the y-axis and z-axis in the global coordinate system of the computing model, and the x-axis is determined by the right-hand rule. The entire finite element mesh is discretized with eight-node hexahedral elements, and a total of 108500 solid elements are meshed. The maximum mesh size is set as 3 m to ensure the accuracy of dynamic calculation [39]. In general, it is appropriate to adopt measured seismic waves or artificial seismic waves synthesized according to site conditions for dynamic calculation; thus, a 20 s acceleration time history of the artificial seismic wave, which is synthesized based on the site condition and seismic design standard of the tunnel, is employed as the input seismic wave (see Figure 2).Figure 2
Time history (a) and Fourier spectrum (b) of input acceleration.
(a)
(b)


## 2.3. Material Model

### 2.3.1. Material Model for Concrete

An appropriate constitutive model for concrete materials is the key to simulate their seismic behavior. The mechanical response of concrete materials under seismic cyclic loading is a complex multistage process, during which concrete materials may suffer stiffness degradation. In addition, the mechanical behavior of concrete materials under tension and compression is quite different. This difference is mainly reflected in the different strength of concrete materials under tension and compression. Generally, the crack model [40] and plastic damage model [41] are used to describe the mechanical properties of concrete materials. In this paper, the plastic damage model proposed by Lubliner et al. [41] is employed because the model can reasonably simulate the softening process, unloading stiffness degradation, and repeated loading process of concrete materials [42]. Under multiaxial stress conditions, the stress of concrete σ can be expressed as(1)σ=1−dσ¯=1−dE0:ε−εpl,where E0 is the initial elasticity modulus; ε is the strain tensor, which can be divided to elastic part εe and plastic part εpl; and d is a scalar damage coefficient, which can be obtained according to the function of tensile damage variable dt and the shear damage variable dc[43].(2)d=1−1−stdt1−scdc,0≤st,sc≤1,dt=1−σt/E0εtpl1/bt−1+σt/E0,0≤dt≤1,dc=1−σc/E0εcpl1/bc−1+σc/E0,0≤dc≤1,where subscripts t and c represent the tension state and compressive state; εtpl and εcpl are the plastic strain under tension and compressive loadings; bt and bc are dimensionless constants and are suggested to be 0.7 and 0.1 by Birtel and Mark [43], st and sc are functions of effective stress σ¯^ and are introduced to represent stiffness recovery effects(3)st=1−wtrσ¯^,sc=1−wcrσ¯^,(4)rσ¯^=∑i=13σ¯^i∑i=13σ¯^i,i=1,2,3,where wt and wc are the weight factors that controls the stiffness recovery degree; rσ¯^ is the stress weight factor; σ¯^i is the principle value of the effective stress σ¯; ⋅ denotes the McCauley bracket, that is, σ¯^i=0.5σ¯^i+σ¯^i.In order to describe the mechanical behavior after yielding, the yield functionF modified by Lee and Fenves [42] and the plastic flow potential G expressed by the Drucker–Prager hyperbolic function are adopted:(5)Fσ¯,ε˜pl=11−αq¯−3αp¯+βε˜plσ¯^max−γ−σ¯^max−σ¯cε˜cpl≤0,(6)G=αpI1+2J2,where ε˜pl is equivalent plastic strain; q¯ is the Mises equivalent effective stress; p¯ is the equivalent hydrostatic pressure; α is material constant; σ¯^max is the maximum eigenvalue of the effective stress; the function β can be given as βε˜pl=1−ασ¯cε˜cpl/σ¯tε˜tpl, where σ¯c and σ¯t are the effective tensile and compressive stresses; γ=31−Kc/2Kc−1, where Kc is the parameter related to the yield curve; I1 is the first invariant of the principal stress and J2 is the second invariant of the deviatoric stress; and αp is a dimensionless material parameter.According to the nonassociated potential flow, the plastic strain rateε˙p can be given by(7)ε˙p=λ˙∂Gσ¯∂σ¯,where λ is the plastic parameter, which can be obtained based on the consistency condition; λ˙ is the derivative of λ versus time. Obviously, the mechanical behavior of concrete materials can be reasonably described according to equations (1)–(7). In addition, according to the code for design of concrete structures [44] published by China, the tension and compressive strength of concrete with C25-grade can be given as 1.78 MPa and 16.7 MPa, respectively. Consequently, the stress-strain relationship of concrete with C25-grade is shown as Figure 3.Figure 3
Plastic damage evolution curves for C25-grade concrete: (a) Tension state and (b) compressive state.
(a)
(b)


### 2.3.2. Analysis Procedure

Generally, seismic response analysis is divided into three steps. In the first step, it is necessary to invert the geostress field of the calculated area according to the limited measured ground stress and geological conditions. However, in many cases, the measured geostress data is hard to obtain; thus, the gravity stress field is adopted as an alternative by many scholars. This alternative is acceptable especially for the cases that the geostress field is mainly controlled by gravity stress. In this paper, the gravity stress is employed as the initial geostress field. In the second step, the static excavation calculation is conducted and the computed stress field is adopted as the initial condition for the seismic response calculation of the third step. In this process, the viscous-spring artificial boundary proposed by Liu et al. [44] is applied to the lateral and bottom boundaries, and the Rayleigh damping with a critical damping ratio of 5% is employed. The mechanical behavior of rock mass under dynamic loadings is modeled by the Mohr–Coulomb failure criterion with a tension cut-off limit. The mechanical behavior of concrete is modeled by the method presented in Section 2.3.1. The geomechanical properties of the rock and concrete units are shown in Table 1. Note that, the time step Δt for the dynamic calculation by the explicit method should meet the below condition [10]:(8)Δt≤κmineleCe0.80≤κ≤0.98,where κ is an empirical constant; Ce is the wave velocity; and le is the minimum size of all elements of the model. In general, time step Δt can be valued by the maximum value according to equation (8) to reduce the amount of calculation.Table 1
Geomechanical properties of the materials used in the numerical solutions.
Material
Density (kg/m3)
Elastic modulus (GPa)
Poisson’s ratio
Friction angle (°)
Cohesion (MPa)
Tensile strength (MPa)
Compressive strength (MPa)
Rock
2.75
9
0.27
42
1.05
1.3
35
Concrete
2.5
25
0.167
46
1.30
1.2
16.7


## 2.3.1. Material Model for Concrete

An appropriate constitutive model for concrete materials is the key to simulate their seismic behavior. The mechanical response of concrete materials under seismic cyclic loading is a complex multistage process, during which concrete materials may suffer stiffness degradation. In addition, the mechanical behavior of concrete materials under tension and compression is quite different. This difference is mainly reflected in the different strength of concrete materials under tension and compression. Generally, the crack model [40] and plastic damage model [41] are used to describe the mechanical properties of concrete materials. In this paper, the plastic damage model proposed by Lubliner et al. [41] is employed because the model can reasonably simulate the softening process, unloading stiffness degradation, and repeated loading process of concrete materials [42]. Under multiaxial stress conditions, the stress of concrete σ can be expressed as(1)σ=1−dσ¯=1−dE0:ε−εpl,where E0 is the initial elasticity modulus; ε is the strain tensor, which can be divided to elastic part εe and plastic part εpl; and d is a scalar damage coefficient, which can be obtained according to the function of tensile damage variable dt and the shear damage variable dc[43].(2)d=1−1−stdt1−scdc,0≤st,sc≤1,dt=1−σt/E0εtpl1/bt−1+σt/E0,0≤dt≤1,dc=1−σc/E0εcpl1/bc−1+σc/E0,0≤dc≤1,where subscripts t and c represent the tension state and compressive state; εtpl and εcpl are the plastic strain under tension and compressive loadings; bt and bc are dimensionless constants and are suggested to be 0.7 and 0.1 by Birtel and Mark [43], st and sc are functions of effective stress σ¯^ and are introduced to represent stiffness recovery effects(3)st=1−wtrσ¯^,sc=1−wcrσ¯^,(4)rσ¯^=∑i=13σ¯^i∑i=13σ¯^i,i=1,2,3,where wt and wc are the weight factors that controls the stiffness recovery degree; rσ¯^ is the stress weight factor; σ¯^i is the principle value of the effective stress σ¯; ⋅ denotes the McCauley bracket, that is, σ¯^i=0.5σ¯^i+σ¯^i.In order to describe the mechanical behavior after yielding, the yield functionF modified by Lee and Fenves [42] and the plastic flow potential G expressed by the Drucker–Prager hyperbolic function are adopted:(5)Fσ¯,ε˜pl=11−αq¯−3αp¯+βε˜plσ¯^max−γ−σ¯^max−σ¯cε˜cpl≤0,(6)G=αpI1+2J2,where ε˜pl is equivalent plastic strain; q¯ is the Mises equivalent effective stress; p¯ is the equivalent hydrostatic pressure; α is material constant; σ¯^max is the maximum eigenvalue of the effective stress; the function β can be given as βε˜pl=1−ασ¯cε˜cpl/σ¯tε˜tpl, where σ¯c and σ¯t are the effective tensile and compressive stresses; γ=31−Kc/2Kc−1, where Kc is the parameter related to the yield curve; I1 is the first invariant of the principal stress and J2 is the second invariant of the deviatoric stress; and αp is a dimensionless material parameter.According to the nonassociated potential flow, the plastic strain rateε˙p can be given by(7)ε˙p=λ˙∂Gσ¯∂σ¯,where λ is the plastic parameter, which can be obtained based on the consistency condition; λ˙ is the derivative of λ versus time. Obviously, the mechanical behavior of concrete materials can be reasonably described according to equations (1)–(7). In addition, according to the code for design of concrete structures [44] published by China, the tension and compressive strength of concrete with C25-grade can be given as 1.78 MPa and 16.7 MPa, respectively. Consequently, the stress-strain relationship of concrete with C25-grade is shown as Figure 3.Figure 3
Plastic damage evolution curves for C25-grade concrete: (a) Tension state and (b) compressive state.
(a)
(b)


## 2.3.2. Analysis Procedure

Generally, seismic response analysis is divided into three steps. In the first step, it is necessary to invert the geostress field of the calculated area according to the limited measured ground stress and geological conditions. However, in many cases, the measured geostress data is hard to obtain; thus, the gravity stress field is adopted as an alternative by many scholars. This alternative is acceptable especially for the cases that the geostress field is mainly controlled by gravity stress. In this paper, the gravity stress is employed as the initial geostress field. In the second step, the static excavation calculation is conducted and the computed stress field is adopted as the initial condition for the seismic response calculation of the third step. In this process, the viscous-spring artificial boundary proposed by Liu et al. [44] is applied to the lateral and bottom boundaries, and the Rayleigh damping with a critical damping ratio of 5% is employed. The mechanical behavior of rock mass under dynamic loadings is modeled by the Mohr–Coulomb failure criterion with a tension cut-off limit. The mechanical behavior of concrete is modeled by the method presented in Section 2.3.1. The geomechanical properties of the rock and concrete units are shown in Table 1. Note that, the time step Δt for the dynamic calculation by the explicit method should meet the below condition [10]:(8)Δt≤κmineleCe0.80≤κ≤0.98,where κ is an empirical constant; Ce is the wave velocity; and le is the minimum size of all elements of the model. In general, time step Δt can be valued by the maximum value according to equation (8) to reduce the amount of calculation.Table 1
Geomechanical properties of the materials used in the numerical solutions.
Material
Density (kg/m3)
Elastic modulus (GPa)
Poisson’s ratio
Friction angle (°)
Cohesion (MPa)
Tensile strength (MPa)
Compressive strength (MPa)
Rock
2.75
9
0.27
42
1.05
1.3
35
Concrete
2.5
25
0.167
46
1.30
1.2
16.7


## 3. Input Method of Oblique Incidence SV Waves

The input of seismic waves into the limited calculation model is highly associated to the boundary conditions of the calculation model. In this paper, a widely used 3D viscous-spring artificial boundary proposed by Liu et al. [44] is employed because it can eliminate drift errors in the low-frequency state under viscous boundary.


### 3.1. 3D Viscous-Spring Artificial Boundary Conditions

A spring and damping system is applied at the nodes of viscoelastic artificial boundaries to absorb the scattered waves from the interior of calculation models (see Figure4). The spring parameter K and damping parameter C, according to [45], can be given by the following.In the normal direction,(9)KN=AlαNGr,CN=Alρcp.In the tangential direction,(10)KT=AlαTGr,CT=Alρcs,where subscripts N and T represent normal and tangential directions; G and ρ are shear modulus and density of infinite domain medium, respectively; r is the distance from the load point to the nodes of artificial boundaries; cp and cs are the velocities of the compression wave and shear wave in infinite domain medium, respectively; Al is the equivalent area for the boundary node l; and αN and αT are modified coefficients, and their values are suggested as 1.333 and 0.667 by [45], respectively.Figure 4
Mechanical model for three-dimensional viscous-spring artificial boundary.


### 3.2. Equivalent Node Force

According to Liu et al. [46], the input of seismic waves can be realized by converting the free field displacement uf into equivalent node force acting on the nodes of artificial boundaries. The equivalent node force Pli can be given by(11)Pli=Cliu˙lif+Kliulif+σlifAl,where subscript i represents the direction of node l; u˙lif is the velocity of the free field; and σlif is the stress of input motions and can be obtained according to the generalized Hooke’s law. It can be found from equation (11) that the equivalent node force Pli can be expressed as the function of the displacement of free field ulif, namely, Pli=Pliulif; thus, the key to solving Pli is to obtain ulif.The free fieldulif at the artificial boundary consists of three components: incident SV waves, reflected P, and SV waves by the free surface (see Figure 5). Therefore, the displacement of free field ulif can be written as(12)ulif=uSV1f+uSV2f+uP2fat the lateral boundary,(13)ulif=uSV1fat the bottom boundary,where subscripts SV1, SV2, and P2 stand for incident SV wave, reflected SV wave, and P wave, respectively.Figure 5
Oblique incidence of plane SV waves in three-dimensional space.The propagation process of seismic waves in inelastic medium has two characteristics. One is the phase difference between the seismic waves at different nodes due to the different distances of these nodes to the epicenter. The other is the amplitude attenuation of seismic waves in the propagation process. These two characteristics are likely to contribute to the inconsistent deformation of underground structures. Therefore, it is required to consider the two characteristics when we solve the displacementulif. Assuming u0t is the displacement of the wave front (see Figure 5), uSV1f, uSV2f,and uP2f can be given by(14)uSV1f=ωs1u0t−Δt1cosα,0,−sinαT,uSV2f=ωs2B1u0t−Δt2−cosα,0,−sinαT,uP2f=ωs3B2u0t−Δt3sinθ,0,−cosθT,where α is the incident angle of seismic waves; θ is reflection angle by free surface, and its value can be given by θ=arcsincPsinα/cS; B1 is the amplitude ratio factor of reflected SV waves, as well as B2 are amplitude factors of reflected P waves; ω is the amplitude attenuation coefficient and can be expressed as the function of the propagation distances si(i = 1, 2, 3) of seismic waves; and Δt is the delay time. Assuming the wave front is parallel to the longitudinal axis of the underground structure and the free surface is horizontal (see Figure 4), the propagation distances si and Δt can be determined according to the geometrical relationship between node l and the wave front:(15)s1=xlsinα+zlcosα,s2=xlsinα+2L−zlcosα,s3=xlsinα+Lcosα−L−zltanβsinα+L−zlcosβ,(16)Δt1=xlsinα+zlcosαcs,Δt2=xlsinα+2L−zlcosαcs,Δt3=xlsinα+Lcosα−L−zltanβsinαcs+L−zl/cosβcp,where L is the height of the calculation model and c is wave velocity.


### 3.3. Determination of Seismic Wave Attenuation Coefficientω

Due to the effects of energy diffusion, the amplitudes of seismic waves are expected to attenuate with the increase of their propagation distance in inelastic medium (e.g., rock and soil). A series of parameters, e.g., the frequency of seismic waves, geological conditions, damping parameters of propagation medium, and propagation distance, were demonstrated to be particularly relevant to the attenuate effects of the seismic wave [47]. Complex inversion calculations are required if we try to consider the impacts of all of these parameters. This is unacceptable for large-scale nonlinear dynamic calculations because it undoubtedly contributes to a large amount of additional computations. In fact, the high frequency components of seismic waves are generally removed before the input of seismic waves, and the remaining low-frequency components have little effects on the attenuation coefficient ω; thus, the impacts of the frequency of seismic waves is not considered here. Assuming a linear relationship between attenuation coefficient ω and propagation distance s, the general expression of the attenuation coefficient ω can be given as(17)ω=ω0s,where ω0 is an empirical constant related to geological condition and damping coefficient. Zhang et al. [48] conducted several numerical experiments on the attenuation coefficient of seismic waves propagating in rock medium, and the value of ω0 is suggested to be determined according to Figure 6. In this paper, the values of elastic modulus and critical damping ratio in the dynamic calculation are 8 GPa and 0.05, respectively; thus, ω0 can be valued as 0.048%. Obviously, ulif can be obtained according to equations (12)–(17). Eventually, Pli can be achieved by substituting ulif into equation (11).Figure 6
Empirical relationships ofω0 and critical damping ratio under different elastic modulus.


### 3.4. Numerical Verification

The accuracy and reliability of the input method of oblique incidence SV waves is verified by a numerical example. A 3D finite element model, as shown in Figure6, is meshed to simulate the propagation of SV waves in elastic medium. The model sizes are 0≤x≤1250, 0≤y≤1250, and 0≤z≤1250. A total of 64000 hexahedral elements are meshed with the maximum mesh size of 31.25 m. The material properties of the medium are listed in Table 2. The central point A (625, 625, 625) on the top boundary, as marked in Figure 7(a), is employed as the monitoring position. A displacement time history having a peak value of 1 m is adopted as the input SV wave and is plotted in Figure 7(b). The incident angle of the SV wave is 20°, and in this case, the normal vector of the wave front is (0.342, 0, 0.940).Table 2
Mechanical properties of the materials used in the numerical example.
Material
Density (kg/m3)
Elastic modulus (GPa)
Poisson’s ratio
Elastic medium
2000
8
0.28Figure 7
(a) 3D Finite element model and (b) incident waves for numerical example.
(a)
(b)Figure8 offers the z-direction and x-direction displacement time histories of point A under 20° incident angle of SV waves. The numerical results in Figure 8 are obtained according to the proposed method, while the theoretical results are obtained based on the theory of elastic wave propagation. The derivation process of theoretical solution can be found in [49, 50]. It can found in Figure 8 that the displacement curves obtained by the two methods are highly similar to each other, and the displacement peaks of the two are almost the same, indicating that the numerical results match the theoretical results well. Therefore, it can be concluded that the present input method is capable of simulating the seismic wave propagation.Figure 8
Displacement time history of point A under 20° incident angle: (a)z-direction and (b) x-direction.
(a)
(b)


## 3.1. 3D Viscous-Spring Artificial Boundary Conditions

A spring and damping system is applied at the nodes of viscoelastic artificial boundaries to absorb the scattered waves from the interior of calculation models (see Figure4). The spring parameter K and damping parameter C, according to [45], can be given by the following.In the normal direction,(9)KN=AlαNGr,CN=Alρcp.In the tangential direction,(10)KT=AlαTGr,CT=Alρcs,where subscripts N and T represent normal and tangential directions; G and ρ are shear modulus and density of infinite domain medium, respectively; r is the distance from the load point to the nodes of artificial boundaries; cp and cs are the velocities of the compression wave and shear wave in infinite domain medium, respectively; Al is the equivalent area for the boundary node l; and αN and αT are modified coefficients, and their values are suggested as 1.333 and 0.667 by [45], respectively.Figure 4
Mechanical model for three-dimensional viscous-spring artificial boundary.


## 3.2. Equivalent Node Force

According to Liu et al. [46], the input of seismic waves can be realized by converting the free field displacement uf into equivalent node force acting on the nodes of artificial boundaries. The equivalent node force Pli can be given by(11)Pli=Cliu˙lif+Kliulif+σlifAl,where subscript i represents the direction of node l; u˙lif is the velocity of the free field; and σlif is the stress of input motions and can be obtained according to the generalized Hooke’s law. It can be found from equation (11) that the equivalent node force Pli can be expressed as the function of the displacement of free field ulif, namely, Pli=Pliulif; thus, the key to solving Pli is to obtain ulif.The free fieldulif at the artificial boundary consists of three components: incident SV waves, reflected P, and SV waves by the free surface (see Figure 5). Therefore, the displacement of free field ulif can be written as(12)ulif=uSV1f+uSV2f+uP2fat the lateral boundary,(13)ulif=uSV1fat the bottom boundary,where subscripts SV1, SV2, and P2 stand for incident SV wave, reflected SV wave, and P wave, respectively.Figure 5
Oblique incidence of plane SV waves in three-dimensional space.The propagation process of seismic waves in inelastic medium has two characteristics. One is the phase difference between the seismic waves at different nodes due to the different distances of these nodes to the epicenter. The other is the amplitude attenuation of seismic waves in the propagation process. These two characteristics are likely to contribute to the inconsistent deformation of underground structures. Therefore, it is required to consider the two characteristics when we solve the displacementulif. Assuming u0t is the displacement of the wave front (see Figure 5), uSV1f, uSV2f,and uP2f can be given by(14)uSV1f=ωs1u0t−Δt1cosα,0,−sinαT,uSV2f=ωs2B1u0t−Δt2−cosα,0,−sinαT,uP2f=ωs3B2u0t−Δt3sinθ,0,−cosθT,where α is the incident angle of seismic waves; θ is reflection angle by free surface, and its value can be given by θ=arcsincPsinα/cS; B1 is the amplitude ratio factor of reflected SV waves, as well as B2 are amplitude factors of reflected P waves; ω is the amplitude attenuation coefficient and can be expressed as the function of the propagation distances si(i = 1, 2, 3) of seismic waves; and Δt is the delay time. Assuming the wave front is parallel to the longitudinal axis of the underground structure and the free surface is horizontal (see Figure 4), the propagation distances si and Δt can be determined according to the geometrical relationship between node l and the wave front:(15)s1=xlsinα+zlcosα,s2=xlsinα+2L−zlcosα,s3=xlsinα+Lcosα−L−zltanβsinα+L−zlcosβ,(16)Δt1=xlsinα+zlcosαcs,Δt2=xlsinα+2L−zlcosαcs,Δt3=xlsinα+Lcosα−L−zltanβsinαcs+L−zl/cosβcp,where L is the height of the calculation model and c is wave velocity.


## 3.3. Determination of Seismic Wave Attenuation Coefficientω

Due to the effects of energy diffusion, the amplitudes of seismic waves are expected to attenuate with the increase of their propagation distance in inelastic medium (e.g., rock and soil). A series of parameters, e.g., the frequency of seismic waves, geological conditions, damping parameters of propagation medium, and propagation distance, were demonstrated to be particularly relevant to the attenuate effects of the seismic wave [47]. Complex inversion calculations are required if we try to consider the impacts of all of these parameters. This is unacceptable for large-scale nonlinear dynamic calculations because it undoubtedly contributes to a large amount of additional computations. In fact, the high frequency components of seismic waves are generally removed before the input of seismic waves, and the remaining low-frequency components have little effects on the attenuation coefficient ω; thus, the impacts of the frequency of seismic waves is not considered here. Assuming a linear relationship between attenuation coefficient ω and propagation distance s, the general expression of the attenuation coefficient ω can be given as(17)ω=ω0s,where ω0 is an empirical constant related to geological condition and damping coefficient. Zhang et al. [48] conducted several numerical experiments on the attenuation coefficient of seismic waves propagating in rock medium, and the value of ω0 is suggested to be determined according to Figure 6. In this paper, the values of elastic modulus and critical damping ratio in the dynamic calculation are 8 GPa and 0.05, respectively; thus, ω0 can be valued as 0.048%. Obviously, ulif can be obtained according to equations (12)–(17). Eventually, Pli can be achieved by substituting ulif into equation (11).Figure 6
Empirical relationships ofω0 and critical damping ratio under different elastic modulus.


## 3.4. Numerical Verification

The accuracy and reliability of the input method of oblique incidence SV waves is verified by a numerical example. A 3D finite element model, as shown in Figure6, is meshed to simulate the propagation of SV waves in elastic medium. The model sizes are 0≤x≤1250, 0≤y≤1250, and 0≤z≤1250. A total of 64000 hexahedral elements are meshed with the maximum mesh size of 31.25 m. The material properties of the medium are listed in Table 2. The central point A (625, 625, 625) on the top boundary, as marked in Figure 7(a), is employed as the monitoring position. A displacement time history having a peak value of 1 m is adopted as the input SV wave and is plotted in Figure 7(b). The incident angle of the SV wave is 20°, and in this case, the normal vector of the wave front is (0.342, 0, 0.940).Table 2
Mechanical properties of the materials used in the numerical example.
Material
Density (kg/m3)
Elastic modulus (GPa)
Poisson’s ratio
Elastic medium
2000
8
0.28Figure 7
(a) 3D Finite element model and (b) incident waves for numerical example.
(a)
(b)Figure8 offers the z-direction and x-direction displacement time histories of point A under 20° incident angle of SV waves. The numerical results in Figure 8 are obtained according to the proposed method, while the theoretical results are obtained based on the theory of elastic wave propagation. The derivation process of theoretical solution can be found in [49, 50]. It can found in Figure 8 that the displacement curves obtained by the two methods are highly similar to each other, and the displacement peaks of the two are almost the same, indicating that the numerical results match the theoretical results well. Therefore, it can be concluded that the present input method is capable of simulating the seismic wave propagation.Figure 8
Displacement time history of point A under 20° incident angle: (a)z-direction and (b) x-direction.
(a)
(b)


## 4. Results and Discussions

The work in this paper aims at the evaluation of the effects of incident angles of SV waves on the nonlinear response characteristics of the tunnel. The analysis process is organized as follows. Firstly, a detailed displacement and stress response analysis of the tunnel is conducted; in this process, the effects of the incident angles on the response process are assessed. Finally, the damage characteristics of the tunnel under different incident angles are studied. Specifically, a typical cross section of the 3D model plotted in Figure1 is chosen and five points on this cross section are selected to illustrate the results; the positions of these points are displayed in Figure 9. In this section, four incident angles α of SV waves are modeled. They are 0°, 10°, 20°, and 30°, respectively.Figure 9
Positions of the monitoring points.


### 4.1. Stress Response Characteristics

The maximum principle stress is employed here to assess the tension state of the concrete structure. Figure10(a) shows the time histories of the maximum principle stresses of point A (at the top arch) under different incident angles. As plotted in Figure 10(a), the stresses experience a violent fluctuation in the period of 2.5–15 s. The fluctuation in this period ranges from −1.16 to 1.638 MPa. The magnitude of the maximum principal stress, which can be regarded as a damage indicator, is what we are extremely concerned about as concrete materials are generally of weak tensile strength. All of the magnitudes of the maximum principle stresses under different incident angles are larger than the tensile strength (1.2 MPa), which indicates the top arch suffers tension damage. In addition, the stress responses of point A under different input angles are clearly different from each other. The primary difference lies on the fluctuation magnitude, which increases as the increasing of the input angles. Figure 10(b) shows the peaks of the maximum principle stresses on the five monitoring points. It can be found that the spandrel (point B) and base corner (point D) of the tunnel experience larger tensile stresses than that of other positions. The peaks of the maximum principle stresses at the same position increase when the input angles increase, and reach the maximum value in the case of 30° input angle. According to these results, we can conclude that the oblique incidence seismic waves contribute to a larger stress response, and this seismic response continues to increase as the increasing input angles of SV waves.Figure 10
Maximum principal stresses under different incident angles: (a) stress time histories of point A and (b) peaks of maximum principal stresses.
(a)
(b)


### 4.2. Displacement Response Characteristics

Figure11(a) displays the time-history curves of the displacement of point A under different input angles. It can be found that the fluctuation laws of the curves are similar to each other, and their peaks appear nearly at the same time. The magnitudes of these displacement curves under different input angles are evidently varied. The magnitudes increase as the input angle increases and reach a maximum of approximately 0.15 m in the case of 30° input angle. This result is consistent with the results of stress response in Section 4.1. The displacement peaks of the five monitoring positions under different input angles are plotted in Figure 11(b). It can be seen that the displacement peaks experience an increase as the increasing of the input angles, and the tunnel suffers the largest displacement in the case of 30° input angle. In addition, the displacement peak of the sidewall (point C) is larger than that of other positions, which is a typical deformation characteristic for underground structures with high sidewalls.Figure 11
Displacement response under different incident angles: (a) displacement time histories of pointA and (b) displacement peaks.
(a)
(b)


### 4.3. Seismic Damage Characteristics

The damage coefficientd, which can be determined by equation (3), is employed here to show the damage characteristics of the tunnel. Figure 12 shows the evolution process of d under different input angles. The damage evolution characteristic is similar to each other under different input angles. When time is 5 s, the damage mainly occurs at the spandrel and base corner with damage coefficients of less than 0.3. With time increasing, the damage zone continues to expand to the sidewalls, the top arch, and the floor (point E). At the end of the simulation (i.e., t = 20 s), the tunnel suffers the most severe damage. Most areas of the tunnel are covered by the damage zone. The most severely damaged positions occur at the spandrel and the base corner with a maximum damage coefficient of 0.9, which means these positions are more vulnerable to be damaged during earthquakes. Moreover, it can be found in Figure 12 that the seismic damage of the tunnel is affected by the input angle. The major difference in the results of the different input angles lies in the damage degree. The overall damage degree of the tunnel increases with the increasing of the input angles of SV waves.Figure 12
Damage coefficients of the tunnel under different input angles: (a) 0° input angle; (b) 10° input angle; (c) 20° input angle; (d) 30° input angle.
(a)
(b)
(c)
(d)


## 4.1. Stress Response Characteristics

The maximum principle stress is employed here to assess the tension state of the concrete structure. Figure10(a) shows the time histories of the maximum principle stresses of point A (at the top arch) under different incident angles. As plotted in Figure 10(a), the stresses experience a violent fluctuation in the period of 2.5–15 s. The fluctuation in this period ranges from −1.16 to 1.638 MPa. The magnitude of the maximum principal stress, which can be regarded as a damage indicator, is what we are extremely concerned about as concrete materials are generally of weak tensile strength. All of the magnitudes of the maximum principle stresses under different incident angles are larger than the tensile strength (1.2 MPa), which indicates the top arch suffers tension damage. In addition, the stress responses of point A under different input angles are clearly different from each other. The primary difference lies on the fluctuation magnitude, which increases as the increasing of the input angles. Figure 10(b) shows the peaks of the maximum principle stresses on the five monitoring points. It can be found that the spandrel (point B) and base corner (point D) of the tunnel experience larger tensile stresses than that of other positions. The peaks of the maximum principle stresses at the same position increase when the input angles increase, and reach the maximum value in the case of 30° input angle. According to these results, we can conclude that the oblique incidence seismic waves contribute to a larger stress response, and this seismic response continues to increase as the increasing input angles of SV waves.Figure 10
Maximum principal stresses under different incident angles: (a) stress time histories of point A and (b) peaks of maximum principal stresses.
(a)
(b)


## 4.2. Displacement Response Characteristics

Figure11(a) displays the time-history curves of the displacement of point A under different input angles. It can be found that the fluctuation laws of the curves are similar to each other, and their peaks appear nearly at the same time. The magnitudes of these displacement curves under different input angles are evidently varied. The magnitudes increase as the input angle increases and reach a maximum of approximately 0.15 m in the case of 30° input angle. This result is consistent with the results of stress response in Section 4.1. The displacement peaks of the five monitoring positions under different input angles are plotted in Figure 11(b). It can be seen that the displacement peaks experience an increase as the increasing of the input angles, and the tunnel suffers the largest displacement in the case of 30° input angle. In addition, the displacement peak of the sidewall (point C) is larger than that of other positions, which is a typical deformation characteristic for underground structures with high sidewalls.Figure 11
Displacement response under different incident angles: (a) displacement time histories of pointA and (b) displacement peaks.
(a)
(b)


## 4.3. Seismic Damage Characteristics

The damage coefficientd, which can be determined by equation (3), is employed here to show the damage characteristics of the tunnel. Figure 12 shows the evolution process of d under different input angles. The damage evolution characteristic is similar to each other under different input angles. When time is 5 s, the damage mainly occurs at the spandrel and base corner with damage coefficients of less than 0.3. With time increasing, the damage zone continues to expand to the sidewalls, the top arch, and the floor (point E). At the end of the simulation (i.e., t = 20 s), the tunnel suffers the most severe damage. Most areas of the tunnel are covered by the damage zone. The most severely damaged positions occur at the spandrel and the base corner with a maximum damage coefficient of 0.9, which means these positions are more vulnerable to be damaged during earthquakes. Moreover, it can be found in Figure 12 that the seismic damage of the tunnel is affected by the input angle. The major difference in the results of the different input angles lies in the damage degree. The overall damage degree of the tunnel increases with the increasing of the input angles of SV waves.Figure 12
Damage coefficients of the tunnel under different input angles: (a) 0° input angle; (b) 10° input angle; (c) 20° input angle; (d) 30° input angle.
(a)
(b)
(c)
(d)


## 5. Summary and Conclusions

In order to evaluate the effects of incident angles of SV waves on the nonlinear seismic response process of a lined tunnel, a 3D input method of oblique incident SV waves is proposed based on the wave field decomposition principle and the viscous-spring artificial boundary. Both the phase difference and amplitude attenuation during the propagation of seismic waves are considered in this method. When we apply the method to investigate the seismic response of a lined tunnel under different input angles of SV waves, the conclusions of the research are as follows:(1)
The oblique incidence seismic waves contribute to a larger stress response. The incident angles of SV waves affect the stress response of the underground structure by increasing the stress fluctuation amplitudes. The stress fluctuation amplitudes increase with the increasing of the incident angles and reaching the maximum in the case of 30° incident angle.(2)
The incident angles have an evident influence on the displacement of the tunnel. The values of the displacement peak of the tunnel increase at the increase of incident angles, which indicates that the larger the incident angles are, the larger displacement the tunnel suffers.(3)
The damage zone of the tunnel continues to extend with the increase of time and incident angles. It is interesting to note that the damage distribution has an obvious spatial difference as the spandrel and base corner of the tunnel suffer the most severe damage than other positions. The oblique incidence of SV waves is likely to increase the spatial difference of the damage distribution and result in a larger intensity and extent of the damage zone.



---

*Source: 8093804-2020-02-20.xml*