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

Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences, Shijiazhuang 050061, China North China University of Water Resources and Electric Power, Zhengzhou 450046, China Changjiang Institute of Survey, Planning, Design and Research, 1863 Jiefang Avenue, Jiang’an District, Wuhan, China State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China


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][10][11]. ese 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. e 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 soilstructure 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. ose 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]. e oblique incidence seismic waves are proved to contribute to the spatial nonuniform deformation of large-scale underground structures [23][24][25].
is 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. e 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. e 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][33][34][35][36][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. e 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 Section 2. 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.

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. e tunnel is embedded in gray limestone, and according to site survey, there is no large geological discontinuity nearby. e 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 Figure 1). e 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.

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 Figure 1). e 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 righthand rule. e entire finite element mesh is discretized with eight-node hexahedral elements, and a total of 108500 solid elements are meshed. e 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).

Material Model for
Concrete. An appropriate constitutive model for concrete materials is the key to simulate their seismic behavior. e 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. is 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 where E 0 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 d t and the shear damage variable d c [43].
where subscripts t and c represent the tension state and compressive state; ε pl t and ε pl c are the plastic strain under tension and compressive loadings; b t and b c are dimensionless constants and are suggested to be 0.7 and 0.1 by Birtel and Mark [43], s t and s c are functions of effective stress σ and are introduced to represent stiffness recovery effects where w t and w c 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, In order to describe the mechanical behavior after yielding, the yield function F modified by Lee and Fenves [42] and the plastic flow potential G expressed by the Drucker-Prager hyperbolic function are adopted: 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 , where σ c and σ t are the effective tensile and compressive stresses; where K c is the parameter related to the yield curve; I 1 is the first invariant of the principal stress and J 2 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 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.

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. is 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. e mechanical behavior of rock mass under dynamic loadings is modeled by the Mohr-Coulomb failure criterion with a tension cut-off limit. e mechanical   Table 1. Note that, the time step Δt for the dynamic calculation by the explicit method should meet the below condition [10]: where κ is an empirical constant; C e is the wave velocity; and l e 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.

Input Method of Oblique Incidence SV Waves
e 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 lowfrequency state under viscous boundary.

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 Figure 4). e spring parameter K and damping parameter C, according to [45], can be given by the following.
In the normal direction, In the tangential direction,  Shock and Vibration 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; c p and c s are the velocities of the compression wave and shear wave in infinite domain medium, respectively; A l 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.

Equivalent Node
Force. According to Liu et al. [46], the input of seismic waves can be realized by converting the free field displacement u f into equivalent node force acting on the nodes of artificial boundaries. e equivalent node force P li can be given by where subscript i represents the direction of node l; _ u f li is the velocity of the free field; and σ f li 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 P li can be expressed as the function of the displacement of free field u f li , namely, P li � P li (u f li ); thus, the key to solving P li is to obtain u f li . e free field u f li at the artificial boundary consists of three components: incident SV waves, reflected P, and SV waves by the free surface (see Figure 5). erefore, the displacement of free field u f li can be written as where subscripts SV 1 , SV 2, and P 2 stand for incident SV wave, reflected SV wave, and P wave, respectively. e 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. e other is the amplitude attenuation of seismic waves in the propagation process.
ese two characteristics are likely to contribute to the inconsistent deformation of underground structures. erefore, it is required to consider the two characteristics when we solve the displacement u f li . Assuming u 0 (t) is the displacement of the wave front (see Figure 5), u f SV 1 , u f SV 2 ,and u f P 2 can be given by  Shock and Vibration where α is the incident angle of seismic waves; θ is reflection angle by free surface, and its value can be given by θ � arcsin(c P sin α/c S ); B 1 is the amplitude ratio factor of reflected SV waves, as well as B 2 are amplitude factors of reflected P waves; ω is the amplitude attenuation coefficient and can be expressed as the function of the propagation distances s i (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 s i and Δt can be determined according to the geometrical relationship between node l and the wave front: where L is the height of the calculation model and c is wave velocity.

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. is 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 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, u f li can be obtained according to equations (12)- (17). Eventually, P li can be achieved by substituting u f li into equation (11).

Numerical Verification.
e 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 Figure 6, is meshed to simulate the propagation of SV waves in elastic medium. e 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. e material properties of the medium are listed in Table 2. e 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). e 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). Figure 8 offers the z-direction and x-direction displacement time histories of point A under 20°incident angle of SV waves. e 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. e 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. erefore, it can be concluded that the present input method is capable of simulating the seismic wave propagation.

Results and Discussions
e 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. e 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 Figure 1 is chosen and five points on this cross section are selected to illustrate the results; the positions of these points are 6 Shock and Vibration displayed in Figure 9. In this section, four incident angles α of SV waves are modeled. ey are 0°, 10°, 20°, and 30°, respectively.

Stress Response Characteristics.
e maximum principle stress is employed here to assess the tension state of the concrete structure. Figure 10   Shock and Vibration 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. e 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 11 is result is consistent with the results of stress response in Section 4.1.

Displacement Response Characteristics.
e 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. e damage coefficient d, 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. e 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. e 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. e major difference in the results of the different input angles lies in the damage degree. e overall damage degree of the tunnel increases with the increasing of the input angles of SV waves.

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) e oblique incidence seismic waves contribute to a larger stress response. e incident angles of SV waves affect the stress response of the underground structure by increasing the stress fluctuation amplitudes. e stress fluctuation amplitudes increase with the increasing of the incident angles and reaching the maximum in the case of 30°incident angle.
(2) e incident angles have an evident influence on the displacement of the tunnel. e 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) e 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. e 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.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.