An IBEM solution to the scattering of plane SH-waves by a lined tunnel in elastic wedge space

The indirect boundary element method (IBEM) is developed to solve the scattering of plane SH-waves by a lined tunnel in elastic wedge space. According to the theory of single-layer potential, the scattered-wave field can be constructed by applying virtual uniform loads on the surface of lined tunnel and the nearby wedge surface. The densities of virtual loads can be solved by establishing equations through the continuity conditions on the interface and zero-traction conditions on free surfaces. The total wave field is obtained by the superposition of free field and scattered-wave field in elastic wedge space. Numerical results indicate that the IBEM can solve the diffraction of elastic wave in elastic wedge space accurately and efficiently. The wave motion feature strongly depends on the wedge angle, the angle of incidence, incident frequency, the location of lined tunnel, and material parameters. The waves interference and amplification effect around the tunnel in wedge space is more significant, causing the dynamic stress concentration factor on rigid tunnel and the displacement amplitude of flexible tunnel up to 50.0 and 17.0, respectively, more than double that of the case of half-space. Hence, considerable attention should be paid to seismic resistant or anti-explosion design of the tunnel built on a slope or hillside.


Introduction
The scattering of elastic waves by the underground structure and the phenomenon of the dynamic stress concentration is an interesting and important topic in many fields, e.g., in earthquake engineering, non-destructive detection, etc. In general,the solution methods can be divided into the analytical method and numerical method. The analytical methods include the wave function expansion (Lee and Trifunac 1979;Liang et al. 2010;), earthquake coefficient method and response displacement method, etc. The numerical methods include the finite element method (Yang and Liu 1994), boundary element method (boundary integral equation method) (Du et al. 1993;Stamos and Beskos 1996;Liang et al. 2013;Chen et al. 2011) and the hybrid method (Datta et al. 1984), etc.
Note that above studies are mainly restricted to the fullspace or half-space model at present. However, in practical engineering, many tunnels or underground pipes are built on sloping topography or cliffs, and then wedge space model would be more appropriate for preliminary quantitative analysis in such sites. Compared with the half-space model, more difficulties will arise for exactly satisfying the free boundary conditions of the wedge space. Achenbach (Achenbach 1970) studied the transient wave propagation problem in wedge space, considering spatially uniform shear tractions applied to one or both faces of the wedge. Knopoff (1969), Budaev and Bogy (1995) and Gautesen (2002) have investigated the wave reflection and transmission coefficients of the free field in elastic wedge space by numerical or experimental method. Moreover, Li and Gong studied the reflection and transmission of obliquely incident Rayleigh wave in two adjacent rectangular space (Li and Gong 1998).
As for the wave scattered by obstacle in wedge space, available results are rarely seen due to the complicated characteristics of wave propagation and scattering. Lee and Sherif studied the scattering of SH-waves by a canyon in wedge-shape space (Lee and Sherif 1996); Shi et al. studied the scattering solutions with a fixed circular inclusion or a hole in Cartesian space by the method of complex variable function (Shi et al. 2006(Shi et al. , 2007. However, up to date, there is no results published for the wave scattering around a tunnel in wedge space of arbitrary angle, only Liu et al. (2009) presented some results around a cavity in wedge space. Zhang et al. (2013) studied the scattering of SHwave by circular cavity in a right-angle plane.
This paper aims to study the scattering of SH-waves by a lined tunnel in elastic wedge space of any angle by the indirect boundary element method. It is illustrated that this method has several advantages such as reducing dimensions of problems, automatic satisfaction of radiation condition, and high calculation precision. Moreover, in this method, the virtual loads can directly act on the boundary surface, which can be recognized as a direct implementation of the Huygens' principle.
This paper will be arranged as follows. Firstly, the numerical procedure for IBEM solution to SH-waves diffraction in wedge space is presented. Then, the accuracy of this method is verified by the comparison between the degenerated solutions and available solutions. Finally, the effects of key parameters, such as the wedge angle, excitation frequency and the incident angle on dynamic response of tunnel are investigated in detail through numerical examples, and some important conclusions have been obtained. Figure 1 shows a lined tunnel of infinite length and constant cross-section located at arbitrary position of the wedge space. Define the vertex of the wedge as o, the angle between the inclined and horizontal surface as mp, the geometric center of the lined tunnel as o 0 . H and D are the vertically and horizontal distance between o and o 0 , respectively. Assume that the material in the wedge space and lined tunnel is homogeneous and linearly elastic. Define the outer and inner surface of tunnel as C 1 and C 2 , the region of the wedge space and tunnel as D 1 and D 2 , the ground and wedge surface as S and B, respectively. The shear wave velocity of the wedge space is defined as

Model
, here l and q are the shear modulus and the mass density, respectively. Accordingly, b 1 , l 1 , q 1 are corresponding parameters in the tunnel lining. Considering plane SH-waves incidence with angle a to the horizontal surface, the two dimensional anti-plane scattering problem needs to be solved. For simplicity, only the cylindrical tunnel is considered in this paper. Note that the tunnel of arbitrary shape can be treated by the IBEM.

Method and solution
The wave scattering problem can be solved by IBEM (Sanchez-Sesma and Campillo 1991) in the following way: based on the theory of single-layer potential, the scattered field can be formed by applying virtual uniform loads on the surface of the scatterers. Then the density of virtual loads can be solved by boundary conditions. The total wave field can be obtained by the superstition of free-wave field and scattered-wave field.
The displacement wave field u (t) in wedge-shaped space using the polar coordinate system o-h satisfies the wave equation as follows: The displacement at any location of the elastic solid, subject to a time harmonic excitation, can be written by the Somigliana integral representation: In which, u(x) and t(x) are the z-direction component of displacement and stress at arbitrary point x, respectively. f(y) is the body force, G(x,n) and T(x,n) are the corresponding displacement and stress Green's functions, respectively.
Similarly, the stresses can be derived according to Hooke's law: Equations (3) and (4) form the basis of the boundary element formulation. Numerical computation requires a discretization of the boundary surface.

Wave field construction
According to linear elastic theory, the total wave field can be decomposed into the free field and the scattered-wave field. Herein the free field is the solution of SH-waves in elastic wedge space without the lined tunnel. The scatteredwave field, based on the theory of single-layer potential, can be constructed by applying virtual uniform loads on the surface of lined tunnel and near ground surface. Then, the free field and the scattered-wave field add up to the total wave field: 3.1.1 Free field Sanchez presented the expression for the displacement generated by plane SH-waves incidence in wedge space with wedge angle mp (Sanchez-Sesma 1985). The resulting formulation is: where u 0 is the displacement amplitude, e n is the Neumann factor (n = 0, e n = 1; n C 1, e n = 2), a is the incident angle with the horizontal direction, and J n/m is the Bessel function of the first kind of order n/m.
In the polar coordinate system o-h and the Cartesian coordinate system x-y, the stress function can be expressed by a simple derivation as follows: Note that the tangential stress by the free field on the boundary can be written as: where k = x/b is the wave number of the SH-waves in wedge-shaped space (b is the shear wave velocity), and the time factor exp(ixt) has been omitted.

Scattered-wave field
According to the previous discussion, the diffracted field is given by Eqs. (3) and (4), which, in the absence of body forces, can be written as: / j denotes the densities of virtual uniform loads on the boundary. The anti-plane line source Green's function in full-space can be expressed in the following form: 1 ðkrÞðc x n x þ c y n y Þ; ð15Þ , c x = (xx 0 )/r, c y = (yy 0 )/r, (x, y) and (x 0 , y 0 ) are the coordinate of the field point and source point, respectively. H ð2Þ n ðÞ is the Hankel function of the second kind of integer order nÁ(n x , n y ) denotes unit normal vector on boundary surface.

Boundary conditions and solution
The boundary conditions of this problem include the zerotraction condition on the surface of wedge space and inner surface of the tunnel, the continuity of displacements and stresses on the interface between the tunnel and the wedge space.
The zero-traction conditions are expressed as follows: The continuity conditions of displacements and stresses can be written as: According to the free boundary and continuity condition, integral equations can be expressed as: It's a singular Fredholm integral equation of the second kind for the boundary sources. We need to discretize the inner and outer surface of tunnel and the nearby wedge surface, then apply virtual uniform loads on each element. Due to the attenuation characteristics of the scatteredwave, the computational accuracy can easily reach to 10 -3 as the discretization range of the surface of wedge space reaches to 8 times the wavelength near the tunnel. / j (n) is assumed to be constant on each boundary element. Define the discretization numbers of the wedge surface, the outer and inner surface of tunnel are N 1 , N 2 , and N 3 , respectively, then the linear Eqs. (21-24) can be rewritten as: Namely, constructing matrix equation . / (densities of virtual loads) can be solved by the equation, In which, dynamic influence functions can be expressed: Equations (12) and (13) can be calculated directly by using two or three-point Gauss quadrature rules when x 6 ¼ n. Analytical expressions can be obtained by the series expansion of Green's functions when x is in the neighborhood of n. It can be expressed as:  where c is the Euler constant (0.5772), lg is the signs for logarithms, and DS is the length of element.
In summary, the introduction of Green's function for fullspace leads to extra discretization of wedge space surface, but has an advantage of analytically treating the singular integration on each element, which is fairly beneficial for improving the calculation accuracy. Then, the density of virtual loads on each element can be solved through the Eqs. (25)-(28). The total wave field is obtained by the superposition of free field and scattered-wave field. Besides, above calculations are performed in frequency domain, and the time domain solution can be obtained by Fourier transform.

Accuracy verification
Firstly, define non-dimensional frequency as the ratio of the equivalent diameter of scatterer to the wavelength of the incident waves: The degenerated solution in elastic wedge space can be calculated using this method. The result of displacement amplitudes by this method compared with the references (Lee and Sherif 1996) and (Lee and Trifunac 1979) are shown in Fig. 2. In Fig. 2(a), for the model of a canyon in wedge space, the following parameters are set: incident frequency g = 2.0 and wedge angle m = 1/2. In Fig. 2(b), for the model of a tunnel in half-space, q 1 /q = 1/3, l 1 / Fig. 4 Dynamic stress concentration factor at the tunnel surface in wedge space (the outer surface of the tunnel) (g = 0.5, 1.0, 2.0; rigid tunnel) l = 0.35, g = 0.5, 1. It shows that our results by IBEM are in good agreement with the analytical method (Lee and Trifunac 1979). Note that, the scattering of elastic wave by tunnel of arbitrary shape in elastic wedge space can be solved by present method.

Numerical examples
In this part, detailed parameters analysis will be presented considering various wedge angles, the variety of tunnel stiffness, incident frequency and angle of incidence. For the rigid tunnel, the material parameters are: q 1 /q = 5/4, b 1 /b = 5/1. For the flexible tunnel, the material parameters are: q 1 /q = 4/5, b 1 /b = 1/3. The ratio of the inner and outer radius of tunnel is set to be r 1 /r 2 = 10/11. Figure 3 illustrates the surface displacement amplitudes in wedge space for different incident frequencies and angles. Set the angles of wedge to be 90°,120°,150°,180°, the incident frequencies g = 0.5, 1.0, 2.0, and the incident angles a = 0, p/6, p/3, p/2, respectively. The location of the lined tunnel is D = 2a,H = 2a (a is the inner radius of the tunnel). In the figure, the x-axis represents the ratio of the distance between the point on the ground surface and the vertex of the wedge to the inner radius of the tunnel (the negative axis represents the horizontal plane; the positive axis represents the inclined plane). The y-axis represents the displacement amplitude |u (t) |. Obviously, there is large difference between the results of half-space (v = 1) and that of the wedge space. Due to the multiple scattering and interface effect of SH-waves between the surfaces of the wedge and the tunnel, the response characteristic is more complicated and the amplification effect is more significant. As the wedge angle decreases, the displacement amplitude increases significantly. For an example, for m = 1/2 (90°wedge space), the displacement amplitude can reach up to about 6.0 for g = 0.5 (which is about 5.0 in 120°wedge space, and about 3.7 in halfspace). As the frequency increases, the displacement amplitude oscillates more quickly in space and the amplification effect appears to be more obvious, up to about 8.0 for g = 2.0 in 90°wedge space. Note that for horizontally incident waves, the max displacement response usually appears just above the tunnel.
Figures 4 and 5 illustrate the dynamic stress concentration factor (DSCF), defined as the ratio of the total shear stress to the stress of incident waves s h 0 Z =s 0 , at the inner and outer surface of the tunnel for different incident frequencies. The calculation parameters are the same as Fig. 3. It can be seen that as the wedge angle decreases, the wave energy centralization becomes more significant and the peak value of DSCF increases gradually. There is large difference between the 90°wedge space (v = 1/2) and the half-space (v = 1). For an example, for g = 0.5 and vertical incident waves, the DSCF at the outer surface can reach up to about 34.0 for m = 1/2 (90°wedge space), but that is about 10.0 in half-space. Therefore, the phenomenon of the dynamic stress concentration near the tunnel in wedge space cannot be analyzed quantitatively using the model of a half-space. It can also be found that spatial characteristics of DSCF strongly depend on the incident frequency, and the stress oscillates more rapidly for high frequency. In addition, comparing Figs. 4 and 5, the peak of the DSCF at inner surface is a bit bigger than that at outer surface, but the spatial distribution characteristics are similar. From the application perspective, we should pay more attention to seismic design for the tunnel located at steep slope. Figure 6 illustrates the effect of the location of the lined tunnel in wedge space on the displacement response, with incident frequency g = 0.5, wedge angle m = 1/2, incident angle a = 0, and the buried depth of tunnel H = 2a. The horizontal distance from tunnel center to the vertex of wedge space takes: D = 2a, 5a, 10a and 30a, respectively. Obviously, when D = 2a, the wave interference effect between the tunnel and the wedge space is more significant. In general, as the distance D increases, the peak of displacement amplitude decrease gradually, which clearly indicates the great influence of vertical plane on wave scattering in a 90°wedge space. Figures 7, 8 and 9 illustrate the displacement amplitude spectrum around the surface of the tunnel and the dynamic stress concentration factor (DSCF) spectrum on the surface of the tunnel to the incident SH-waves for different wedge angles. The incident frequency takes g[(0,4.0), and observation points for displacement amplitude are located at x/a = -1, -2 (right above the tunnel), -4, and 4 on ground surface, respectively; for DSCF, the points take h 0 = 0°, 90°,180°, and 270°on the surface of the tunnel, respectively. It shows that as the wedge angle decreases, the peak value of the frequency spectrum increases, and the displacement amplitude can reach to 4.3 in the 90°wedge space but that is 1.8 in the half-space. In addition, compared with the half-space, the spectrum curve oscillates more rapidly in the wedge space due to the complex interference effect near the tunnel. As for DSCF spectrum curve, it shows that the dynamic stress concentration effect seems more significant for low frequency waves, and the wedge angle has large influence on the peak value of DSCF spectrum. For an example, in the 90°wedge space, the peak of the stress spectrum inner the tunnel can reach to 55 (g = 0.3) at the point h 0 = 0°(the point at the right of the tunnel). Accordingly, that are just 26 (g = 0.25) in halfspace.
Considering different material properties of the lined tunnel, Figs. 7, 8, and 9 also show the displacement amplitude and DSCF spectrum of the flexible tunnel case. Except for the material parameters, other parameters remain the same as those of the rigid tunnel. Compared with the rigid tunnel, the displacement amplitude on wedge space surface near the tunnel shows more significant amplification effect, which can reach to 7.5 at x/a = -2 (just above the tunnel), but that is 4.3 for the rigid tunnel in 90°wedge space. Additionally, the DSCF spectrums show little amplification effect inside the tunnel for the flexible case. Figures 10, 11, and 12 show the contour pictures of the displacement amplitude both around the rigid and flexible lined tunnel in wedge space. Consider the wedge angle v = 1, 3/4, 1/2, the incident angles a = 0, p/2, and the frequencies g = 0.5, 1.0, 2.0, respectively.
As for the rigid case, in wedge space, the displacement amplification effect becomes more significant, e.g., the displacement amplitude becomes 7.0 in 90°wedge space, and that is just 3.0 in half-space, 4.0 in 135°wedge space. It can be seen that as the frequency increases, the interference effect between the tunnel and the wedge space becomes more notable. There are more ''focusing points'' of wave energy for high incident frequencies. In addition, for the horizontal incident waves in half-space, due to the existence of the lined tunnel, the shielding effect on the displacement response can be seen clearly behind the tunnel, while this phenomenon is not so obvious for the 135°and 90°wedge space due to the existence of the inclined surface of wedge space. As for the flexible tunnel, the displacement amplification effect inner the tunnel should be paid more attention. For example, when g = 1.0 and a = 0 the peak of the displacement amplitude inner the flexible tunnel can reach up to 12, while that of the rigid tunnel is about 2.5. Besides, for the flexible tunnel, in the 90°wedge space, the peak of the displacement amplitude inner the tunnel can reach up to 17.0 for g = 2.0 and a = p/2 , while that is about 6.4 for the half-space case. In general, the peak of displacement appears in the side facing the incoming waves or most close to the wedge space surfaces. Additionally, as the incident frequency increases, the amplification effect seems more significant. Hence, the seismic resistant design of the x/a 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 2ϕ=60 0 2ϕ=75 0 2ϕ=90 0 2ϕ=120 0 L=5aη=0.5 flexible tunnel should try to control the large displacement response of the tunnels. Figure 13 shows the contour picture of dynamic stress concentration factor (DSCF) of the rigid tunnel in wedge space. It can be seen that spatial characteristics of DSCF is more complicated for high frequency, but the peak value increases gradually. Besides, as the wedge angle decreases, the peak value of DSCF increases clearly, e.g., for g = 2.0, a = p/2, the peak value can reach up to about 21.0 for m = 1/2 (90°wedge space), but that is about 8.6 in halfspace. In general, the DSCF decreases from the inner surface to the outer surface gradually.
To study the seismic response of tunnel built in cliffy mountains, Fig. 14 shows a sharply angular wedge space. Define the apex angle of the wedge as 2u, the distance between o and o 0 as L. oo 0 is the bisector of the wedge space, and the incident SH-wave comes along the bisector line oo 0 . Figure 15 illustrates the surface displacement amplitudes both for acute and obtuse angle wedge space with different incident frequencies, apex angles and the location of the lined tunnel. Set the incident frequencies to be g = 0.5, 1.0, 2.0, the apex angle of wedge 2u = 60°, 75°, 90°, 120°, the location of the lined tunnel L = 5a, 10a, 15a, the ratio of mass density and shear wave velocity q 1 /q = 5/4, b 1 /b = 5/1, respectively. Obviously, when the apex angle of wedge is an acute angle, the seismic response is more significant. For example, when g = 0.5, L = 5a, the displacement amplitude can reach up to about 6.0 for 2u = 60°, but that is about 1.4 for 2u = 120°at the apex of wedge. Besides, as the distance L increases, the peak of displacement amplitude decreases gradually.

Conclusions
This paper presents an indirect boundary element method (IBEM) for the scattering of plane SH-waves by a lined tunnel in elastic wedge space based on the theory of singlelayer potential. Compared with the exact analytical solution, the accuracy of this method has been verified. Through detailed parameter analysis, several important conclusions can be drawn as follows: The reflection and diffraction of elastic waves are changed substantially by the inclined surface of wedge space. The wave scattering features become more complicated, and strongly depend on the wedge angle, the angle of incidence, the incident frequency, the location of lined tunnel and material parameters. As the wedge angle decreases, the displacement amplification and dynamic stress concentration effects around the tunnel become more significant. As for the rigid tunnel in 90°wedge space, compared with the half-space case, the peak values of the displacement amplitude on ground surface and the dynamic stress concentration factor inside the tunnel increase more than 100 %. As for the flexible tunnel in wedge space, the dynamic stress concentration effect is not so significant, but the displacement amplitude of the tunnel can reach up to 17 times that of the incident waves. Therefore, the seismic or anti-explosion design of the tunnel or pipes located in slope or hillside should adopt the wedge space model to improve the accuracy.
It is worth mentioning that present method is applicable to the tunnel of arbitrary shape and in arbitrary-angle wedge space. In addition, the solution technique can also be expanded to solve the scattering of P and SV waves in elastic wedge space. The specific solution procedure will be given in another paper.