Zero-stress, cylindrical wave functions around a circular underground tunnel in a flat, elastic half-space: Incident P-waves

Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

Contents lists available at ScienceDirect

Soil Dynamics and Earthquake Engineering journal homepage: www.elsevier.com/locate/soildyn

Zero-stress, cylindrical wave functions around a circular underground tunnel in a flat, elastic half-space: Incident P-waves Chi Hsin Lin a, Vincent W. Lee b, Maria I. Todorovska b, Mihailo D. Trifunac b, a b

Converse Consultants, Monrovia, CA 91016, USA Department of Civil Engineering, University of Southern California, Los Angeles, CA 90089-2531, USA

a r t i c l e in fo

abstract

Article history: Received 11 April 2009 Accepted 23 January 2010

The zero-stress boundary conditions at the surface of the half-space in the presence of surface and subsurface cavities for in-plane, incident cylindrical P- and SV-waves have always posed challenging problems. The outgoing cylindrical P- and SV-waves can be represented by Hankel functions of radial distance coupled with the sine and cosine functions of angle. Together, at the half-space surface the P- and SV-wave functions are not orthogonal over the semi-infinite radial distance from 0 to infinity. Thus, to simultaneously satisfy the zero in-plane, normal, and shear stresses, an approximation of the geometry is often made. This paper presents an analytical formulation of the boundary-valued problem, where the Hankel wave functions are expressed in integral form, changing the representation from cylindrical to rectangular coordinates, so that the zero-stress boundary conditions at the half-space surface can be applied in a more straightforward way. It is sometimes desirable to use alternate, simpler, or approximate solutions to the boundary-valued problems, without going to great lengths to have all of the boundary conditions satisfied. The free-stress boundary conditions at the half-space surface to be solved here are among the most complicated boundary conditions to be satisfied in the present class of problems. It is thus of interest to illustrate what the solution would be like if the half-space boundary conditions are not imposed—in other words, if they are ‘‘relaxed.’’ A section of this paper is devoted to this approach, and the results of the solutions with the half-spaced, stress-free boundary conditions ‘‘imposed’’ and ‘‘relaxed’’ are presented and compared. The method introduced here may also be applied to more complicated wave-propagation problems, like the diffraction problems involving surface and sub-surface inhomogeneities in a poroelastic half-space medium. & 2010 Elsevier Ltd. All rights reserved.

Keywords: In-plane P-wave diffraction Circular tunnel ‘‘Exact’’ & ‘‘Relaxed’’ half-space boundary conditions

1. Introduction Large, underground structures such as cavities, pipes, subways, and tunnels are always present in metropolitan areas. When excited by incoming seismic waves, their presence will generate additional waves by diffraction and scattering, resulting in amplifications and de-amplification of the input motions, which will affect the deformations, distributions, and concentrations of stresses near the ground surfaces and in the structures above. It is thus important to understand the theoretical and engineering aspects of the effects of such scattering and diffraction. The setting for our problem can be illustrated by the January 17, 1994 Northridge earthquake in the Los Angeles metropolitan area of Southern California. The earthquake was centered on a blind thrust fault underneath the San Fernando Valley, and it resulted in significant damage in many parts of Los Angeles County, not just

 Corresponding author.

E-mail address: [email protected] (M.D. Trifunac). 0267-7261/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.soildyn.2010.01.010

in the densely populated San Fernando Valley [53,54,56–60]. It also produced many strong-motion recordings, giving a valuable set of earthquake acceleration data [46,47]. Among the many sites that experienced damage, we mention the Van Norman Complex (VNC) of the Los Angeles Department of Water and Power (LADWP), which supplies water and power to the City of Los Angeles. During the Northridge earthquake, the damage to the power distribution facilities caused power shortages that affected almost 2.5 million customers in Southern California [39]. The Rinaldi Receiving Station was one of the four most severely shaken substations of the LADWP (close to the Van Norman Complex), with the strong-motion recordings there producing a peak acceleration of 0.8 g, and a peak velocity of 170 cm/s [60]. There are many buried corrugated metal pipes (CMP) below the surface in the Van Norman Complex (VNC), and the Northridge earthquake resulted in significant damage to some of them. One such 2.4-m large-diameter pipe suffered complete lateral collapse at the lower San Fernando Dam (LSFD), and Davis et al. [5] studied this case in detail. They noted that underground pipes and tunnel linings are rarely designed to withstand the transverse loads caused by

880

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

earthquake ground motions. This is so because there have been, up until the 1990s, very few documented case histories of transverse failures of underground structures from past earthquakes. The Northridge earthquake demonstrated the need for improved analysis and understanding of the seismic responses of underground buried structures. To this end, it is necessary to develop methods for study of in-plane responses of underground cavities, pipes, and tunnels to the input of seismic waves. For such wave-propagation, boundary-valued problems, it is necessary to select suitable coordinate systems that allow the vector-wave equations to be separable into scalar-wave equations for the P- and SV-wave potentials. However, because of the complexity of the problem, simple analytic solutions to diffraction problems in the vicinity of such underground cavities are limited at present. For such bodies, which occupy a finite volume, the waveboundary-valued problems can be solved by finite-element methods. Many software versions of such formulations are available in industry, but for problems involving infinite or semi-infinite medium, such as wave-propagation problems in a half-space, the finite-element methods are not directly applicable. A typical ‘‘large’’ finite-element mesh may have to be defined, along with adequate transmitting boundaries, in order to simulate the semi-infinite medium. This, however, cannot always be done in a satisfactory manner. Solving for the scattered wave potentials that need to satisfy the zero-stress boundary conditions at the half-space surface is a classical boundary-valued problem, which is of interest to scientists, seismologists, and engineers—and in particular earthquake engineers. Lamb [12] presented a paper of lasting significance, titled ‘‘On the propagation of tremors over the surface of an elastic solid.’’ In this paper, he posed the problem of wave motion generated at the surface of an elastic half-space subjected to a concentrated load at the surface, or inside the halfspace. Both time-harmonic and impulsive loads were considered. Some 30 years ago, we started to develop analytical solutions for an underground, circular, unlined tunnel (cavity) subjected to incident, anti-plane SH waves [13]. Lee and Trifunac [21] extended this solution to the case of an underground circular elastic tunnel. Both solutions used the method of images, which created a line of symmetry on the half-space surface that resulted in zero anti-plane shear stresses there. Such a method of images, however, failed for inplane P- and SV-wave models, although approximations for the reflections from the flat surface were proposed. Lee and Cao [15] and Cao and Lee [2,3] used a large, circular, almost-flat surface to approximate the half-space surface, and they presented numerical solutions to diffraction problems of surface circular canyons of various depths for incident plane P- and SV-waves. Todorovska and Lee [41–43] used the same method for anti-plane SH waves and for incident Rayleigh waves on circular canyons. Later, Lee and Karl [17,18] extended the method to diffraction of P-, and SV-waves by circular underground cavities. Lee and Wu [22,23] used the same method for arbitrary-shaped, two-dimensional canyons. Davis et al. [5] used such an approximate method in their studies of the failure of underground pipes involving transverse SV wave incidence. Liang et al. [25–33] continued the analyses by the same method, for problems involving circular-arc canyons and valleys, as well as underground pipes in elastic and poro-elastic half-space. Recently, the above approximate method has met with criticism. It was noted that the large, circular approximation does not converge quickly to that of the flat half-space when the radius, R, of the circular surface approaches infinity because, as R-N, the Bessel functions used in the transformation would approach zero. It was suggested that without a method to satisfactorily satisfy the free-stress boundary conditions at the half-space surface, it might just be better to relax the free-stress conditions than to use the large circle approximations [24,49]. In the meantime, many computational and numerical methods have also been developed for elastic wave diffraction problems in

the half-space. There are currently two main approaches to deal with such boundary-valued problems: (1) Numerical methods such as finite difference, finite element and boundary element methods. (2) Analytical methods such as the method of wave function expansion, perturbation method, the method of integral equations and integral transforms. In dealing with most practical problems, which involve arbitrary-shaped obstacles, topography with irregular geometry and material nonlinearities, numerical methods are better. To deal with the semi-infinite medium in a finite domain, numerical procedures such as absorbing boundaries [4,8,10], and infinite elements using mapping functions [67] have to be implemented. However, such procedures must be further developed as they are known to produce errors in the solutions [9,11]. In contrast to the numerical methods, analytical methods provide closed-form solutions with better accuracy and relatively simpler numerical implementation. However, the limitation of analytical methods is that those can only treat linearly-elastic or visco-elastic media with simple geometries. Therefore only a limited set of problems has been solved analytically. In spite of this, analytic solutions offer the best benchmarks to test and to verify other approximate solutions. We will use an analytic solution scheme in this study. The method employed in this study is to leave the half-space surface flat, and to represent the diffracted cylindrical waves using polar coordinates, but then to transform their representation into integral form, which uses rectangular coordinates. This approach uses an integral method along the lines proposed by Lamb [12], who represented the cylindrical wave function (namely, the Hankel function of order 0) in an integral form in rectangular coordinates. This study extends such representation to Hankel functions of all integer orders for general cylindrical waves. Such integral representations have been used by Thiruvenkataher and Viswananthan [40], for a boundary-valued problem of an underground cavity with applied traction force, and by Sanchez-Sesma et al. [38]. The method described in this work may also serve as a general method for studies of more complicated wave-propagation problems such as the diffraction problems involving surface and sub-surface inhomogeneities in both elastic and poroelastic half-space media. The method should be of particular value for further developments of analytical studies of complex soil–structure interaction problems in metropolitan areas with dense concentrations of tall buildings [65] with both shallow and deep foundations [62,63], and for buildings with large horizontal dimensions [44,45,55] in the presence of underground tunnels. The method should further be useful in the development of advanced, site-specific procedures when there is, for example, a requirement for detailed probabilistic description of the response spectra [61], for description of the spatial variations in the prediction of synthetic strong ground motion [50,66], or in the prediction of liquefaction hazard [48,52].

2. The model: incident plane P-waves Consider (Fig. 1) a flat, half-space medium with a circular cavity of outer radius ‘‘a’’ and at a depth ‘‘h’’ below the surface. It is lined by a tunnel of outer radius ‘‘a’’ and inner radius ‘‘b’’ made up of another elastic medium. The half-space is assumed to be elastic, isotropic, and homogeneous, with mass density r and Lame´ elastic constants l, m. The tunnel is supposed to be elastic with mass density r1 and Lame´ elastic constants l1, m1. A rectangular coordinate system x1–y1 is defined with the origin O1 at the center of the circular tunnel, and with the y1-axis pointing

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

where

y

sin gb =cb ¼ sin ga =cb ;

x

K1 ¼

O ψr

ψ2 s

φr

φ 2s

ψ1 s

y1

φ1s

x1

φ2t, ψ2t

¼

inga

½e

cr ¼

arn Jn ðka r1 Þeiny1

n ¼ 1 1 X

brn Jn ðkb r1 Þeiny1

ð4aÞ

n ¼ 1

where arn ¼ ð1Þn K1 ei2ka h cos ga einga brn ¼ ð1Þn K2 eiðka h cos ga þ kb h cos gb Þ eingb

ð4bÞ

2.3. The total input: free-field waves

With respect to the {x1, y1} and {r1, y1} coordinates, with origin at O1 (the center of the tunnel), the incident plane P-waves are assumed to have an angle of incidence ga with respect to the y-axis, frequency o and wave number ka ¼ o/ca, where ca is the P-wave speed. The waves are described by the potential function fi of unit amplitude: 1 X

1 X

fr ¼

2.1. Incident plane P-waves

f ¼ f ðx1 ; y1

ð3Þ

which can again be expressed in terms of a Fourier–Bessel series [1] as

vertically upwards into the half-space. A corresponding polar coordinate system r1  y1 is also defined, with the same origin as the rectangular system and r1 being the radial distance and y1 the angle measured counter-clockwise from the x-axis. A corresponding rectangular coordinate system is also defined: x y with origin O at distance h above O1. There is also a polar coordinate r y defined at the same origin. The half-space region is thus  p r y r0, with the half-space surface at y¼ 0 (y1 ¼h), NoxoN or rZ0, y ¼0, p. The excitation consists of a plane incident P-wave.

Þ ¼ eika ðx1 sin ga þ y1 cos ga Þ

ð2bÞ

¼ K2 eiðka h cos ga þ kb h cos gb Þ eikb r sinðy1 gb Þ

φ1t, ψ1t

a

Fig. 1. The model.

i

2 sin 2ga cos 2gb sin 2ga sin 2gb þ k2 cos2 2gb

fr ¼ K1 eika h cos ga eika ðx1 sin ga ðy1 hÞcos ga Þ ¼ K1 ei2ka h cos ga eika r1 sinðy1 ga Þ cr ¼ K2 eika h cos ga eikb ðx1 sin gb ðy1 hÞcos gb Þ

b

φi

i

sin 2ga sin 2gb þ k2 cos2 2gb

Back to the x1, y1 and r1, y1 coordinates, with origin at O1, the reflected plane P- and SV-waves take the form

θ1

O1

K2 ¼

ka =kb ¼ cb =ca

sin 2ga sin 2gb k2 cos2 2gb

k ¼ kb =ka ¼ ca =cb

r1

γα

881

¼

The total free-field P- and SV-wave potentials are the sum of incident and reflected plane waves that satisfy the zero-stress surface boundary conditions. The total free-field potentials are, respectively, fff and cff: 1 X

fff ¼ fi þ fr ¼

eika r1 sinðga þ y1 Þ

an Jn ðka r1 Þeiny1

n ¼ 1

Jn ðka r1 Þeiny1

ð1aÞ

1 X

cff ¼ cr ¼

bn Jn ðkb r1 Þeiny1

ð5aÞ

n ¼ 1

n ¼ 1

where

using the identity [1]: expðikz sin yÞ ¼ expð1=2kzðt1=tÞÞ ¼

1 X

n

t Jn ðkzÞ

an ¼ einga þ arn ¼ einga þ ð1Þn K1 ei2ka h cos ga einga

iy

for t ¼ e

bn ¼ brn ¼ ð1Þn K2 eiðka h cos ga þ kb h cos gb Þ eingb

n ¼ 1

ð5bÞ

With respect to the x, y and r, y coordinates with origin at O, the potential takes the form 2.4. The waves scattered from the underground tunnel

fi ¼ fi ðx; yÞ ¼ eika ðx sin ga þ ðy þ hÞcos ga Þ ¼ eika h cos ga eika ðx sin ga þ y cos ga Þ ¼ eika h cos ga eika r sinðga þ yÞ

ð1bÞ

2.2. Reflected plane waves by the flat ground surface The presence of the half-space surface will result in reflected plane P- and SV-waves. The reflected P-wave will have the same angle of reflection, ga, with respect to the vertical axis, and the reflected SV-wave will have a different (smaller) angle, gb, with respect to the vertical; wave number kb ¼ o/cb, where cb is the SV-wave speed. The reflected waves are represented by their potentials: r

r

ika h cos ga ika ðx sin ga y cos ga Þ

f ¼ f ðx; yÞ ¼ K1 e

e

ika h cos ga ika r sinðga yÞ

¼ K1 e

e

cr ¼ cr ðx; yÞ ¼ K2 eika h cos ga eikb ðx sin gb y cos gb Þ ¼ K2 eika h cos ga eikb r sinðgb yÞ ð2aÞ

The presence of the underground tunnel will result in scattered P- and SV-waves in the half-space. These are represented by Hankel functions of the first kind, and they represent outgoing cylindrical waves. The P- and SV-scattered potentials, measured with respect to the r1, y1 coordinates, take the form:

fs1 ¼ cs1 ¼

1 X n ¼ 1 1 X

A1;n Hnð1Þ ðka r1 Þeiny1 B1;n Hnð1Þ ðkb r1 Þeiny1

ð6Þ

n ¼ 1

The presence of the half-space surface will result in additional P- and SV-waves reflected into the half-space that, with respect to

882

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894 t

the {r1, y1} coordinates, will be finite everywhere 1 X

fs2 ¼ cs2 ¼

A2;n Jn ðka r1 Þeiny1

n ¼ 1 1 X

t

in terms of the transmitted P- and SV-potentials f1 þ f2 and ct1 þ ct2 inside the tunnel.

B2;n Jn ðkb r1 Þeiny1

ð7Þ

3.2. Stress and displacement continuity at the outer surface of the tunnel wall

n ¼ 1

2.5. Total scattered and reflected waves in the half-space The total scattered and reflected waves in the half-space are given by 1 X

f ¼ fff þ fs1 þ fs2 ¼ ff

c¼c

s s þ c1 þ c2

¼

sr ¼ sffr þ ssr;1 þ ssr;2 ¼ str;1 þ str;2

½ðan þ A2;n ÞJn ðka r1 Þþ A1;n Hnð1Þ ðka r1 Þeiny1

n ¼ 1 1 X

½ðbn þ B2;n ÞJn ðkb r1 Þþ B1;n Hnð1Þ ðkb r1 Þeiny1

With respect to the {r1, y1} coordinates, the outer surface of the tunnel at r1 ¼a is the interface between the total waves in the half-space outside the tunnel and the transmitted waves inside the tunnel. The continuity of the stress and displacement boundary conditions give

try ¼ tffry þ tsry;1 þ tsry;2 ¼ ttry;1 þ ttry;2 ð8Þ

n ¼ 1

ur ¼ uffr þ usr;1 þ usr;2 ¼ utr;1 þutr;2 uy ¼ uffy þ usy;1 þusy;2 ¼ uty;1 þ uty;2 ;

0 r y r 2p

at r1 ¼ a

ð11Þ

substituting from Eqs. (5a), (6), and (7). 2.6. The transmitted waves within the tunnel The region within the tunnel wall is b rrra, 0 r y r2p, with the outer circular surface at r ¼a in contact with the half-space medium and the inner circular surface at r¼ b. The presence of two circular surfaces will result in both outgoing and incoming cylindrical waves, which are best represented with respect to the {r1, y1} coordinates as follows: (i) Outgoing waves:

ft1 ¼ ct1

¼

1 X n ¼ 1 1 X

D1;n Hnð1Þ ðkb r1 Þeiny1

ct2 ¼

n ¼ 1 1 X

With respect to the x, y coordinates with origin at O at the halfspace surface right above the underground tunnel, the zero-stress boundary conditions at the half-space surface are given by

1 ox o 1

at y ¼ 0

ð12aÞ

ð9aÞ in terms of the free-field incident, reflected, and scattered waves. Because the free-field waves fff and cff already satisfy the zerostress conditions at the surface, namely, sffy ¼ tffyx ¼ 0, only the s s s s scattered P- and SV-waves, f1 þ f2 and c1 þ c2 , must together satisfy the zero-stress conditions:

and (ii) incoming waves: 1 X

3.3. Zero-stress at the half-space surface

sy ¼ sffy þ ssy;1 þ ssy;2 ¼ 0 tyx ¼ tffyx þ tsyx;1 þ tsyx;2 ¼ 0;

C1;n Hnð1Þ ðka r1 Þeiny1

n ¼ 1

ft2 ¼

Note that for the case of an underground circular cavity (in the absence of a tunnel wall), these two sets of boundary conditions will be reduced to just one set of zero-stresses at the circular surface of the cavity in contact with the half-space medium [5,17,18] and no transmitted waves will be generated inside the cavity.

C2;n Hnð2Þ ðka r1 Þeiny1 D2;n Hnð2Þ ðkb r1 Þeiny1

ð9bÞ

ssy;1 þ ssy;2 ¼ 0 tsyx;1 þ tsyx;2 ¼ 0;

1 o x o1

at y ¼ 0

ð12bÞ

n ¼ 1

4. The solution

with the total transmitted waves inside the tunnel as 1 X

ft ¼ ft1 þ ft2 ¼ ct ¼ ct1 þ ct2 ¼

ðC1;n Hnð1Þ ðka r1 Þ þ C2;n Hnð2Þ ðka r1 ÞÞeiny1

n ¼ 1 1 X

ðD1;n Hnð1Þ ðkb r1 Þ þ D2;n Hnð2Þ ðkb r1 ÞÞeiny1

The solution will be developed in the following subsections. ð9cÞ

4.1. Zero-stress at the tunnel inner surface

n ¼ 1

3. The boundary conditions The boundary conditions are given at the three surfaces: the inner and outer surfaces of the tunnel wall, and the flat surface of the half-space.

At the inner surface of the tunnel, r1 ¼b, 0r y r2p 9 8 ! ! > > @2 f 2m @2 c @c > 2 > > > þ l r f þ 2 m r  ( ) > > = < @r@y @y @r2 r2 sr ! ! ¼ 2 2 2 tr y > 2m @ f @f 1 @ c 1 @c @ c > > > r1 ¼ a > > > ; : r 2 r @r@y  @y þ m r 2 @y2 þ r @r  @r2 >

  ¼

0 0

r1 ¼ b

ð13Þ t

3.1. Zero-stress at the tunnel inner surface With respect to the r1, y1 coordinates, with origin at O1 (the center of the tunnel), the zero-stress boundary conditions at the inner surface r1 ¼b are given by:

sr ¼ str;1 þ str;2 ¼ 0 try ¼ ttry;1 þ ttry;2 ¼ 0;

0 r y r 2p

at r1 ¼ b

ð10Þ

t

t

t

t

t

for f ¼ f ¼ f1 þ f2 and c ¼ c ¼ c1 þ c2 , the transmitted P- and SV-wave potentials inside the tunnel. Substituting their representations from Eq. (9c) into Eq. (13), the following can be obtained, using the orthogonality of the exponential functions in 0r y r2p. For n ¼ N to N: ( ) ( )   C1;n C2;n 0 ð4Þ ½Eð3Þ j j þ½E ¼ ð14Þ r1 ¼ b r1 ¼ b n n D1;n D2;n 0

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

where [36], for j ¼1, 2, 3, 4: 2 ðjÞ 3 ðjÞ E11;n E12;n ðjÞ ðjÞ 5 ½En  ¼ ½En ðka r; kb rÞ ¼ 4 ðjÞ ðjÞ E21;n E22;n ðjÞ E11;n ðjÞ E12;n ðjÞ E21;n ðjÞ E22;n

2

( ¼ ½Gð3Þ n ð 1 ; ka1 a; kb1 aÞ

m

C1;n D1;n

883

)

( þ ½Gð4Þ n ð 1 ; ka1 a; kb1 aÞ

m

C2;n D2;n

) ð18bÞ

with

2 2

¼ ðn þ nkb r =2ÞCnðjÞ ðka rÞka rCðjÞ n1 ðka rÞ ðjÞ ¼ i½ðn2 þ nÞCnðjÞ ðkb rÞ þ nkb rCn1 ðkb rÞ ðjÞ ¼ i½ðn2 þ nÞCnðjÞ ðka rÞ þ nka rCn1 ðka rÞ ðjÞ ðjÞ 2 2 2 ¼ ðn þ nkb r =2ÞCn ðkb rÞ þ kb rCn1 ðkb rÞ

ð15Þ

where the EnðjÞ matrix is defined in Eq. (15) above, and for j¼1, 2, 3, 4: 2 3 m½EnðjÞ ðka a; kb aÞ ðjÞ 4 5 and ½Gn ðm; ka a; kb aÞ ¼ ½DnðjÞ ðka a; kb aÞ 2 ðjÞ 3 ðjÞ D11;n D12;n ðjÞ ðjÞ 4 5 where ð19aÞ ½Dn  ¼ ½Dn ðka r; kb rÞ ¼ ðjÞ ðjÞ D21;n D22;n

are the stress functions computed from the wave potentials. Here, for j ¼1, 2, 3, 4, the CnðjÞ for j¼ 1, 2, 3, 4, are the Bessel or ð1Þ ð1Þ ð2Þ ð3Þ ð4Þ Hankel functions, Cð1Þ n ¼ Jn ; Cn ¼ Yn ; Cn ¼ Hn ; and Cn ¼ Hn ; respectively. Eq. (14) can be rearranged to obtain, for each n ¼  N to N: ( ) ( ) ( ) C2;n C1;n C1;n 1 ½Eð3Þ jr1 ¼ b ¼ ½Eð4Þ ¼ ½Fn  ð16Þ n n D2;n D1;n D1;n

with 1 ð3Þ ½Fn  ¼ ½Fn ðka1 b; kb1 bÞ ¼ ½Eð4Þ n  ½En jr1 ¼ b

a matrix defined on the inner radius of the tunnel, so that one set of the coefficients in the tunnel medium can be eliminated.

4.2. Stress and displacement continuity at the tunnel outer surface The boundary conditions at the tunnel outside surface will involve the stress and displacement continuity between the halfspace and the tunnel medium. 9 8 ! ! > > @2 f 2m @2 c @c > > 2 > > l r f þ2 m r  þ > > > > 2 2 > > @r@ y @ y @r r > > 9 > 8 > > > ! ! > > s > > > > r 2 2 2 > > 2m > > > > > > @ f @ f 1 @ c 1 @ c @ c > = < = > < try >  2 > r  þm 2 2 þ 2 r @r@ y @ y @r r r @r @ y ¼ ur > > > > > > > > > > > @f 1 @c > > > > > > : þ > > uy ; > > > > > r @r @ y > > > > > > > > 1 @f @c > > > >  ; r1 ¼ a : r @y @r Outside 9 8 trr > > > > > > > > < tr y = ð17Þ ¼ u > > > > > r > > > : u ; r1 ¼ a y Inside Substituting the waves from Eqs. (8) and (9) into Eq. (17), and again using orthogonality of the trigonometric functions, for each n ¼ 1 to 1 one obtains: 2 3( 3( ) 2 ) m0 ½Eð3Þ m0 ½Eð1Þ A1;n A2;n þan n ðka a; kb aÞ n ðka a; kb aÞ 4 5 4 5 þ B1;n B2;n þbn ½Dð3Þ ½Dð1Þ n ðka a; kb aÞ n ðka a; kb aÞ 2 3( 2 3( ) ) m1 ½Eð3Þ m1 ½Eð4Þ C1;n C2;n n ðka1 a; kb1 aÞ n ðka1 a; kb1 aÞ 5 4 5 ¼4 þ D1;n D2;n ½Dð3Þ ½Dð4Þ n ðka1 a; kb1 aÞ n ðka1 a; kb1 aÞ ð18aÞ or in a simpler matrix form ( ) ( ) A1;n A2;n þ an ð1Þ ½Gð3Þ ð m ; k a; k aÞ ð m ; k a; k aÞ þ ½G b b 0 a 0 a n n B1;n B2;n þ bn

ðjÞ ðjÞ D11;n ¼ nCnðjÞ ðka rÞ þ ka rCn1 ðka rÞ ðjÞ ¼ i½nCnðjÞ ðkb rÞ D12;n ðjÞ ¼ i½nCnðjÞ ðka rÞ D21;n ðjÞ ðjÞ ¼ nCnðjÞ ðkb rÞkb rCn1 ðkb rÞ D22;n

ð19bÞ

are the displacement functions computed from the wave potentials. Note that the left-hand side of Eqs. (18a) and (18b) uses wave functions of the half-space, while the right-hand side uses wave functions inside the elastic medium of the tunnel. 4.3. Zero-stress at the half-space surface (y ¼0) With respect to the x, y coordinates, with origin at O at the half-space surface, the zero-stress boundary conditions at the half-space surface are: !3 2 @2 f @2 c 2 l r f þ 2 m  6 7 ( ) @y2 @x@y 7  0  6 sy 6 ! 7 jy ¼ 0 ¼ 6 ð20Þ 7¼ 2 2 2 txy 6 7 0 4 m 2@ f þ @ c  @ c 5 @x@y @y2 @x2 which are to be applied only to the scattered P- and SV-wave s s s s potentials, f1 þ f2 and c1 þ c2 . They are represented by cylindrical (polar) coordinates, and the above zero-stress conditions are to be applied at the flat surface y ¼ 0 given in rectangular coordinates. As stated in the introduction, applying such equations at the flat half-space surface to cylindrical waves is not a trivial matter, and previous attempts have been made to approximate the flat boundary surface by a large circular surface. The approach in this work is to leave the half-space surface flat and to follow the approach of Lamb’s problem [12] problem to represent each mode or term of the Hankel wave series as an integral in rectangular coordinates. In other words, the Fourier transform of the wave functions is taken. Take the case of the nth mode of the P- and SV-wave scattered potentials:  n  n  R1 i kna Hnð1Þ ðka r1 Þeiny1 ¼ 1 eikx1 na jy1 j dk ipna ka " #  n ð21Þ R1 knb in Hnð1Þ ðkb r1 Þeiny1 ¼ 1 eikx1 nb jy1 j dk ipnb kb qffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where na ¼ k2 k2a , and nb ¼ k2 k2b . The case of n ¼0, in the above integral expressions was first used by Lamb [12] in his paper on the tremors: R1

1 ikx1 na jy1 j e dk ipna R1 1 ikx1 nb jy1 j e dk H0ð1Þ ðkb r1 Þ ¼ 1 ipnb H0ð1Þ ðka r1 Þ ¼

1

ð22Þ

In Eq. (21), the Hankel wave functions in the r1, y1 cylindrical coordinates are transformed to the x1, y1 rectangular coordinates

884

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

suitable for applying the zero-stress boundary conditions at the flat half-space surface. Each integrand in the above integrals of Eq. (21) is the Fourier transform of the nth mode of the wave function. Using the {x, y} coordinate system with origin at O, on the half-space surface the integral takes the following form, with y1 ¼y+ h, in the region 0ry1 rh, or  hryr0:  n  n  R1 i kna ena h eikxna y dk Hnð1Þ ðka r1 Þeiny1 ¼ 1 ipna ka "  n # : ð23Þ R1 knb in Hnð1Þ ðkb r1 Þeiny1 ¼ 1 enb h eikxnb y dk ipnb kb Next, we apply the transform to the whole set of terms for the P- and SV- potentials: 1 X

fs1 ¼ fs1 ðr1 ; y1 Þ ¼ cs1 ¼ cs1 ðr1 ; y1 Þ ¼

n ¼ 1 1 X

s

A1;n Hnð1Þ ðka r1 Þeiny1 ¼ f1 ðx; yÞ ¼ s

B1;n Hnð1Þ ðkb r1 Þeiny1 ¼ c1 ðx; yÞ ¼

Z Z

1

½a1 ðkÞeikxna jyj dk

1 1

½b1 ðkÞeikxnb jyj dk

1

n ¼ 1

ð24Þ such that ( ) a1 ðkÞ b1 ðkÞ

¼

" 1 X in za;n ðhÞ=na ip 0 n ¼ 1

#(

0

zb;n ðhÞ=nb

A1;n

)

B1;n

ð25Þ

with

za;n ðhÞ ¼



kna ka

n

ena h

and

zb;n ðhÞ ¼



knb kb

n

enb h

The presence of both the underground tunnel and of the halfspace flat surface will result in additional scattered and reflected waves (Eq. (7)). Those will now be expressed in integral form, with respect to the rectangular coordinate system (x, y) with origin O at the half-space surface, where for yr0: R ikx þ na y fs2 ¼ fs2 ðx; yÞ ¼ 1 dk 1 ½a2 ðkÞe R1 ð26Þ s s c2 ¼ c2 ðx; yÞ ¼ 1 ½b2 ðkÞeikx þ nb y dk In terms of the rectangular coordinate system (x1, y1) with origin at O1, at the center of the tunnel, for y1 rh in the half-space, they take the form: R ikx1 þ na y1 na h fs2 ¼ fs2 ðx1 ; y1 Þ ¼ 1 e dk 1 ½a2 ðkÞe R1 ð27Þ s s ikx1 þ nb y1 nb h c2 ¼ c2 ðx1 ; y1 Þ ¼ 1 ½b2 ðkÞe e dk

where Z A2;n ¼

  n  kna in ena h a2 ðkÞ dk ka 1 #  Z 1 "  k n b n n b h B2;n ¼ in e b2 ðkÞ dk kb 1 or ( ) A2;n B2;n

1

Z

1

¼ 1

" in

#(

za;n ðhÞ

0

0

zb;n ðhÞ

ð31Þ

a2 ðkÞ b2 ðkÞ

) dk

ð32Þ

with za;n ðhÞ; zb;n ðhÞ defined as in Eq. (25) above. Note that in the half-space: s

s

(1) Waves given by P- and SV-potentials f1 and c1 (Eq. (24)), with Hankel functions of the first kind, represent outgoing waves and satisfy Sommerfeld’s radiation condition and s s (2) Waves given by P- and SV-potentials f2 and c2 (Eq. (30)), with Bessel functions, which are finite everywhere, represent waves reflected from the half-space surface. Together they form two sets of independent potentials for both P- and SV-waves and thus constitute a complete set of solutions of the P- and SV-wave equations. We next apply the zero-stress boundary conditions at the half-space surface. Substituting the integral representations of the waves in Eq. (24) and Eq. (30) into the zero-stress boundary conditions of Eq. (12), 2 3( ) Z 1 2k2 k2 2iknb a1 ðkÞ b 4 5 dk b1 ðkÞ k2 þ n2b 1 2ikna 2 3   Z 1 2k2 k2 2iknb ( a ðkÞ ) 0 b 2 4 5 þ dk ¼ ð33Þ 2 2 ðkÞ b 2ik n k þ n 0 2 a 1 b After arranging Eq. (33), 2 31 2 ( ) 2k2 k2b 2k2 k2b 2iknb a2 ðkÞ 5 4 ¼ 4 b2 ðkÞ 2ikna 2k2 k2b 2ikna ¼

2iknb 2k2 k2b

3( 5

a1 ðkÞ

)

b1 ðkÞ

( ) a1 ðkÞ 1 ½PðkÞ b1 ðkÞ RðkÞ

ð34Þ

where x1 ¼ x, y1 ¼y + h. Writing k ¼ ka cos wa ¼ kb cos wb pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 k2a ¼ þ i k2a k2 ¼ þika sin wa qffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi nb ¼ k2 k2b ¼ þi k2b k2 ¼ þ ikb sin wb

na ¼

where ð28Þ

" ½PðkÞ ¼ ½Pðk; ka ; kb Þ ¼

P 11 ðkÞ

P 12 ðkÞ

P 21 ðkÞ

P 22 ðkÞ

# ;

ð36aÞ

with P 11 ¼ ½ð2k2 k2b Þ2 þ4k2 na nb  P 12 ¼ 4ikð2k2 k2b Þnb P 21 ¼ 4ikð2k2 k2b Þna

Using the identity 1 X

ð35Þ

and

in terms of some imaginary complex angles wa and wb, these reflected waves become R ika r1 cosðy1 wa Þ na h fs2 ¼ 1 e dk 1 ½a2 ðkÞe R1 ð29Þ s c2 ¼ 1 ½b2 ðkÞeikb r1 cosðy1 wb Þ enb h dk

eika r cosðyÞ ¼

RðkÞ ¼ Rðk; ka ; kb Þ ¼ ð2k2 k2b Þ2 4k2 na nb ¼ Rayleigh function

P 22 ¼ ½ð2k2 k2b Þ2 þ4k2 na nb 

in Jn ðkrÞeiny

ð36bÞ

n ¼ 1

4.4. Transformation from integral back to series

the reflected waves take the form Z 1 1 X fs2 ¼ ½a2 ðkÞeikx þ na y dk ¼ A2;n Jn ðka r1 Þeiny1

cs2 ¼

Z

1 1 1

½b2 ðkÞeikx þ nb y dk ¼

n ¼ 1 1 X n ¼ 1

B2;n Jn ðkb r1 Þeiny1

( ð30Þ

Substituting Eq. (32) into Eq. (34), gives " # ) Z ( ) 1 za;n ðhÞ 0 A2;n a1 ðkÞ in ½PðkÞ ¼ dk: zb;n ðhÞ 0 B2;n b1 ðkÞ 1 RðkÞ

ð37Þ

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

(

5. The numerical integration step

Substituting Eq. (25) into Eq. (37) gives " # ) Z 1 za;n ðhÞ 0 A2;n in ½RðkÞ ¼ zb;n ðhÞ 0 B2;n 1 RðkÞ " #( ) 1 X 0 A1;m im za;m ðhÞ=na  dk zb;m ðhÞ=nb 0 B1;m ip m ¼ 1 ( ) 1 X A1;m ½Iðm; nÞ ; ¼ B1;m m ¼ 1

The numerical implementation of the above equations for the coefficients rests upon the computation and solution of the complex matrices in Eqs. (16), (18), and (38). In particular, the matrices f½Iðm; nÞ; m; n ¼ 0; 1; 2 . . .gin Eq. (38) involve numerical computations of integrals of the form ð38Þ

where

Z

is a 2  2 matrix defined by

½Iðm; n; h1 ; h2 Þ inm ¼ ip

Z

Z

2

P 12 za;n ðh1 Þzb;m ðh2 Þ

P 11 za;n ðh1 Þza;m ðh2 Þ 6 na 1 6 6 6 1 RðkÞ 4 P 21 zb;n ðh1 Þza;m ðh2 Þ

3

7 7 7 dk P 22 zb;n ðh1 Þzb;m ðh2 Þ 7 5

nb

1

na

nb

ð39Þ This completes the derivation. In summary, Eqs. (16) and (18), together with Eq. (38) give the following:

no RðkÞ

dk

is infinite in range:  NokoN, and it has singular points at qffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 k20 ¼ 0 and k¼ 7kR (where

RðkÞ ¼ 0 is the Rayleigh pole k ¼kR), as shown in Fig. 2. These singular points will have to be ‘‘eliminated’’ in the numerical integration step.

Let F ðkÞ=RðkÞ ¼ GðkÞ, and write the integral in the form

½Gð3Þ n ðm0 ; ka a; kb aÞ A1;m B1;m

A1;n

þ ½Gð1Þ n ðm0 ; ka a; kb aÞ

B1;n

½Iðm; n; h; hÞ

 ½Gð3Þ n ðm1 ; ka1 a; kb1 aÞ

þ ½Gð4Þ n ðm1 ; ka1 a; kb1 aÞ½Fn ðka1 b; kb1 bÞ ( ¼ ½Gð1Þ n ðm0 ; ka a; kb aÞ

an bn

)

(

C1;n

Z

k0

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk þ 1 k2 k20

Z

k0 k0

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk þ k2 k20

Z

1 k0

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk; k2 k20

isolating the middle integral Z

k0 k0

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk k2 k20

as the one with the singular point at the left and right end points k¼ 7k0. If we let x ¼k/k0, so k¼ k0x, dk¼k0dx, the integral path is ( k : k0 -k0

k0 k0

)

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk ¼ k2 k20

Z

Z ¼

pffiffiffiffiffiffiffi Z 1GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk ¼ k0 ðk2 k20 Þ k0

k0 k0

Z Gðk0 xÞ i qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k0 dx ¼ i 2 1 k20 ð1x Þ 1

GðkÞ i qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk k20 k2 1 1

Gðk0 xÞ ffi dx: qffiffiffiffiffiffiffiffiffiffiffi 2 1x

ð42Þ

The integral in the form Z i

1 1

m ¼ 1

)

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk ¼ 1 k2 k20

x : 1-1

Together they give the equations that can be used to solve for ( ) A2;n the coefficients. Using Eqs. (38) and (16), the coefficients B2;n ( ) C2;n and can be eliminated from Eq. (18b), resulting in a D2;n ( ) ( ) C1;n A1;n and : system of equations for B1;n D1;n 1 X

1

ð41Þ

from ð38Þ:

)

1 F ðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk k2 k20 RðkÞ

k¼ 7k0 ¼ka or kb (where

Z

(

1 1

Z

(ii) At the outer radius, a, the common interface of the tunnel and half-space: ( ) ( ) A1;n A2;n þan þ ½Gð1Þ ; ½Gð3Þ n ðm0 ; ka a; kb aÞ B n ðm0 ; ka a; kb aÞ B 1;n 2;n þ bn ( ) C1;n ¼ ½Gð3Þ n ðm1 ; ka1 a; kb1 aÞ D 1;n ( ) C2;n þ ½Gð4Þ ð m ; k a; k aÞ from ð18bÞ a b 1 n 1 1 D2;n and (iii) At the half-space medium: ( ) ( ) 1 X A2;n A1;m ¼ ½Iðm; n; h; hÞ B2;n B1;m m ¼ 1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 k2a and nb ¼ k2 k2b , and RðkÞ

5.1. Removing the apparent singularities at k ¼ 7k0 (ka or kb)

(i) At the inner radius b of the tunnel medium: ( ) ( ) C1;n C2;n ¼ ½Fn ðka1 b; kb1 bÞ from ð16Þ; D2;n D1;n



1 F ðkÞ

1

is the Rayleigh function given in Eq. (35). This integral

h2 ¼ h

(

1

where n0 ¼ na or nb, with na ¼

½Iðm; nÞ ¼ ½Iðm; n; h1 ; h2 Þh1 ¼ h;

885

Gðk0 xÞ ffi dx qffiffiffiffiffiffiffiffiffiffiffi 2 1x

qffiffiffiffiffiffiffiffiffiffiffiffi 2 that has wðxÞ ¼ 1= 1x as the weight function can now be numerically computed, say, by the Gauss–Chebyshev quadrature integration scheme [1,37].

D1;n

-∞

-kR

-k0

0

k0

kR

ð40Þ Fig. 2. Poles along the integration path.

∞

886

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

5.2. The integral to infinity involving singularity at k ¼kR

Z

Take the last integral in Eq. (41), where Z 1 GðkÞ 1 F ðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk: k0 k0 k2 k2 k2 k2 RðkÞ

¼

and using L’Hopital’s rule

1

0

¼

0

Itqffiffiffiffiffiffiffiffiffiffiffiffiffiffi has affi singular point at k¼k qffiffiffiffiffiffiffiffiffiffiffiffiffi ffi R. Let n ¼ k ¼ n2 þ k20 , and dn ¼ k=ð k2 k20 Þ dk. And the integral path is ( k : k0 -1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 k20 , so that

1 k0

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dk ¼ k2 k20

Z

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2 k20

1

Z

00 00 ½2F 01 þ ðnnR ÞF 1 R1 þ ½F 01 þ ðnnR ÞF 01 R01 ½F R01 þ ðnnR ÞðF 01 R0 þ F 1 R Þ  n-nR 2R1 R01

and canceling

¼

00 00  ½2F 01 þ ðnnR ÞF 1 R1 ðnnR ÞF 1 R  0 0 n-nR ¼ 0 : 2R1 R1

Again, using L’Hopital’s rule ¼

v : 0-1 Z

½F 1 ðnÞ þ ðnnR ÞF 01 ðnÞR1 ðnÞðnnR ÞF 1 ðnÞR01 ðnÞ 0 ¼ ; 0 R21 ðnÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 þk20 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dv n2 þ k20

1 Gð

GðkÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn ¼ k 0 k2 k20 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z Z 1 F ð n2 þ k2 Þ 1 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn ¼ ¼ 0 Rð n2 þk2 Þ n2 þ k20 0 0

1 0

00 000 00 00 00 000 ½3F 1 þ ðnnR ÞF R1 þ ½2F 01 þ ðnnR ÞF 1 R01 ½F 1 R þðnnR ÞðF 1 R þ F 1 RR Þ  00 n -n R 2ðR0 R0 þ R1 R1 Þ

With many of the terms canceling or dropping out as n-nR,  00 00 2F 01 R01 F 1 R1  2F 01 ðnR ÞR01 ðnR ÞF 1 ðnR ÞR1 ðnR Þ 0 ¼ : h ðnR Þ ¼  n ¼ n R  2½R01 2 2½R01 ðnR Þ2 ð47Þ

F 1 ðnÞ dn R1 ðnÞ

The same technique and conclusion can be applied to determine the first integral of Eq. (41) involving a singularity at k¼  kR. ð43Þ 6. The ‘‘solution with relaxed boundary conditions’’

where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ð n2 þk20 Þ F 1 ðnÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; n2 þ k20 When k¼kR, n ¼ Z

1 0

F 1 ðnÞ dv ¼ R1 ðnÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R1 ðnÞ ¼ Rð n2 þ k20 Þ:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2R k20 ¼ nR

Z n R d Z F 1 ðnÞ dvþ R1 ðnÞ 0

vR þ d vR d

F 1 ðnÞ dn þ R1 ðnÞ

Z

1 vR þ d

F 1 ðnÞ dn R1 ðnÞ ð44Þ

From Eq. (44), it is clear that the singularity exists only in the second term of the integral. The first term can be numerically computed by Romberg integration, and the last term can be computed by the shifted Gauss–Laguerre method.

Let F 1 ðnÞ R1 ðnÞ

I2 ¼ I2 ¼

vR þ d

Z

v R d vR þ d v R d

F 1 ðnÞ dn; R1 ðnÞ hðnÞ dn ¼ nnR

1 X

fff ¼ fi þ fr ¼

Then, let the 2nd integral Z

6.1. The input free-field waves The total input free-field P- and SV-wave potentials are the same as those introduced in Section 2.3 above, and represent the sum of incident and reflected plane waves. Together they satisfy the zero-stress boundary conditions on the flat ground surface. fff and cff are given by:

5.3. Integrating the singular integral around the Rayleigh pole

hðnÞ ¼ ðnnR Þ

It is sometimes desired to produce an alternate, simpler approximate estimate of the solution to the boundary-valued problem without going to great lengths to have all the boundary conditions satisfied. The free-stress boundary conditions at the halfspace surface are often the most complicated boundary conditions to be satisfied in the present problem. It is thus interesting to see what the solution would be like if the half-space boundary conditions are not imposed and not satisfied—in other words, if they are ‘‘relaxed.’’ Recently, some work on diffraction of elastic waves in poroelastic media was directed along this line, namely relaxing the boundary conditions at the half-space surface [24,49]. This section will consider the solution using such an approach.

an Jn ðka r1 Þeiny1

n ¼ 1 1 X

cff ¼ cr ¼

from ð5aÞ

bn Jn ðkb r1 Þeiny1

n ¼ 1

Z

vR þ d vR d

hðnÞhðnR Þ dn þ hðnR Þ nnR

Z nR þ d

Where the last integral above integrates to zero: Z nR þ d 1 dn ¼ 0: nR d nnR

nR d

1 dn nnR

where an ¼ einga þarn ¼ einga þ ð1Þn K1 ei2ka hcosga einga bn ¼ brn ¼ ð1Þn K2 eiðka hcosga þ kb hcosgb Þ eingb

ð45Þ

from ð5bÞ

6.2. The scattered waves in the ‘‘half-space’’ outside the tunnel

with the second term integrated to zero. Using the mean-value theorem of calculus, Z vR þ d hðnÞhðnR Þ I2 ¼ dn ¼ h0 ðn Þð2dÞ  h0 ðnR Þð2dÞ ð46Þ nnR v R d

The same P- and SV-scattered potentials, measured with respect to the r1, y1 coordinates, as in Section 2.4 above, are used:

for some n A ½uR d; uR þ d. In the limit, when d-0; n -nR , and I2 is approximated using the last term in Eq. (46), using a finite interval of width (2d). Evaluating h0 ðnR Þ:   ðnnR ÞF 1 ðnÞ 0 h0 ðnR Þ ¼ limn-nR h0 ðnÞ ¼ limn-nR R1 ðnÞ

cs1 ¼

fs1 ¼

1 X

A1;n Hnð1Þ ðka r1 Þeiny1

n ¼ 1 1 X

from ð6Þ

B1;n Hnð1Þ ðkb r1 Þeiny1

n ¼ 1

Because the free-stress boundary conditions at the half-space surface are to be relaxed in this section, no additional P- and SV-waves are needed. In other words, it is assumed here that the

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

half-space boundary is not present, so the scattered waves given in Eq. (7) of Section 2.4 above are not present here, and the scattered waves subjected to the free-field input waves are thus not different from those scattered from the tunnel in the full-space.

887

or ( ½Gð3Þ n ðm0 ; ka b; kb bÞ

A1;n B1;n

)

( þ½Gð1Þ n ðm0 ; ka b; kb bÞ

(

¼ ½Gð3Þ n ð 1 ; ka1 b; kb1 bÞ

m

C1;n

)

an bn

) (

þ ½Gð4Þ n ð 1 ; ka1 b; kb1 bÞ

m

D1;n

C2;n D2;n

6.3. The transmitted waves within the tunnel

ct ¼ ct1 þ ct2 ¼

1 X

; ð48bÞ

Because the boundary conditions between the half-space and the tunnel remain the same, and the transmitted waves inside the tunnel are the same as in Eq. (9c) of Section 2.6 above:

ft ¼ ft1 þ ft2 ¼

)

ðC1;n Hnð1Þ ðka r1 Þ þ C2;n Hnð2Þ ðka r1 ÞÞeiny1

n ¼ 1 1 X

from ð9cÞ:

ðD1;n Hnð1Þ ðkb r1 Þ þD2;n Hnð2Þ ðkb r1 ÞÞeiny1

n ¼ 1

The boundary conditions to be satisfied at the inner radius of the tunnel (r ¼b) for the transmitted waves will be the same as those given in Section 3, namely zero normal and tangential stress at r ¼b (Eqs. (10) and (13) above). Application of these boundary conditions results in Eq. (16) above for each n ¼  N to N: ( ) ( ) ( ) C2;n C1;n C1;n 1 ð3Þ  ½E j  ¼ ½Eð4Þ ¼ ½F from ð16Þ; n r¼a n n D2;n D1;n D1;n

where the GnðjÞ matrices, for j ¼1, 2, 3, 4 are the same as GnðjÞ given in Eq. (19). With the half-space, free-stress boundary conditions ( ) C2;n in terms of relaxed, the final equation, after expressing D2;n ( ) C1;n , and using Eq. (16), takes a final form that is much D1;n simpler than Eq. (40) for each n: ( ) A1;n ½Gð3Þ ð m ; k a; k aÞ b 0 a n B1;n ð4Þ ð½Gð3Þ n ðm1 ; ka1 a; kb1 aÞþ ½Gn ðm1 ; ka1 a; kb1 aÞ½Fn ðka1 b; kb1 bÞÞ

( ¼ ½Gð1Þ n ðm0 ; ka a; kb aÞ

an bn

(

)

C1;n D1;n

) :

ð49Þ

with 1 ð3Þ ½Fn  ¼ ½Fn ðka1 a; kb1 aÞ ¼ ½Eð4Þ n  ½En jr ¼ a

a matrix defined on the inner radius of the tunnel wall, as in Section 4.1, so that one set of the coefficients in the tunnel medium can be eliminated.

7. Results and discussions

6.4. The ‘‘relaxed boundary conditions’’ Without imposing the free-stress boundary conditions, the remaining boundary conditions require only the stress and displacement continuity equations at the tunnel outer surface (r ¼a), as given in Section 4.2 above: 9 8 ! ! > > @2 f 2m @2 c @c > > 2 > >  > > lr f þ2m 2 þ 2 r > > > > @r@ y @ y @r r > 9 > 8 > > > > ! ! > > sr > > > > 2 2 2 > > > > > > > > 2 m @ f @ f 1 @ c 1 @ c @ c > > > < > = = <  r  m þ þ try 2 2 2 2 r @r @r@y @y r @y @r ¼ r u > > > > > > > > @f 1 @c > > > r > > > > > ; > :u > þ > > y > > > > @r r @ y > > > > > > > > 1 @ f @ c > > > >  ; : r1 ¼ a r @y @r Outside 9 8 t > > > > > rr > > = < tr y > ¼ from ð17Þ; u > > > > > r > > > :u ; y r1 ¼ a Inside

fs1

cs1 the

and only P- and SV-scattered wave potentials with present in the half-space. As in Eqs. (18a) and (17) now gives, for each n ¼  N to N: 2 3( 3( ) ) 2 m0 ½Eð3Þ m0 ½Eð1Þ A1;n an n ðka b; kb bÞ n ðka b; kb bÞ 4 5 5 þ4 B1;n bn ½Dð3Þ ½Dð1Þ n ðka b; kb bÞ n ðka b; kb bÞ 2 3( 3( ) 2 ) m1 ½Eð3Þ m1 ½Eð4Þ C1;n C2;n n ðka1 b; kb1 bÞ n ðka1 b; kb1 bÞ 4 5 4 5 ¼ þ D1;n D2;n ½Dð3Þ ½Dð4Þ n ðka1 b; kb bÞ n ðka1 b; kb bÞ 1

It should be noted that, unlike Eq. (40), which is a system of matrix equations of infinite order, Eq. (49), for each n, is a 4  4 finite system of equations and can be solved easily. This solution with ‘‘relaxed boundary conditions’’ can be thought of as a ‘‘fullspace solution’’ with half-space, free-field waves as input (Eq. (5)).

1

ð48aÞ

7.1. The x- and y-components of displacements A computer program was written in Visual Fortran to calculate the displacement amplitudes and phases on the half-space surface (y¼0) above the tunnel and on the tunnel surface (r¼ a), for incident plane P-waves with various angles of incidence and for various frequencies. As in all previous work and analyses, the significance of the diffraction of the elastic waves in the presence of the elastic tunnel will depend significantly on the wavelength and frequency of the incident waves relative to the size of the underground tunnel. Following all previous work, it is convenient to express the frequency of the incident—i.e., the reflected and scattered waves in terms of the dimensionless parameter Z defined as follows [51]:

Z¼

kb a

p

¼

oa 2a 2a ¼ ¼ pcb cb T l

ð50Þ

The first term, kba/p, in Eq. (50) represents the dimensionless wave number of the shear SV-waves. With kb ¼ o/cb, the term oa/ pcb then corresponds to the dimensionless frequency. With the period of the waves equal to T¼2p/o, the wavelength of the shear waves is l ¼cbT. Then, Z ¼1 will correspond to shear waves with wavelengths equal to one tunnel diameter, whereas Z ¼5 corresponds to shear waves with wavelengths 1/5 that of the diameter of the tunnel. Small Z thus corresponds to waves with wavelengths that are longer compared with the radius of the tunnel. This parameter is used to illustrate how the resulting motions are not dependent separately (or independently) on the size of the tunnel and the wavelength of the waves, but instead on their ratio, given here by the parameter Z, as discussed in Trifunac [51] and Lee [14]. A Poisson ratio of v¼ 0.25 is assumed in all examples presented in this work.

888

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

The horizontal x- and vertical y-components of displacements are computed from the P(f)- and SV(c)-wave potentials: 9 8 @j @c > > > ( ) > þ > > < ux @x @y = : ð51Þ ¼ uy @j @c > > > > > ; : @y  @x > On the half-space ground surface y¼0, where the Cartesian coordinates of displacements will be used, the free-field displacements will be calculated from their Cartesian representations, namely (Eq. (1a) for the incident plane P-waves, fi and Eq. (2a) for the reflected plane P- and SV-waves, fr, cr). The displacements from the scattered waves will be calculated using the integral representations of the scattered potentials, Eq. (24) for f1, c1, and Eq. (26) for f2, c2. Thus at y¼0, 8 9 " #( ( ) ) Z 1 < uffx = ux a1 ðkÞ ik nb ikx þ e j jy ¼ 0 ¼ dk uy b1 ðkÞ : uffy ; y ¼ 0 na ik 1 " #( ) Z 1 a2 ðkÞ ik nb ikx e dk: ð52Þ þ b2 ðkÞ na ik 1 where the first term of free-field waves is: 9 8 9 8 < uffx = < ð1 þK1 Þika sin ga eika x sin ga K2 ikb cos gb eikb x sin gb = ¼ eika h cos ga : uffy ; : ð1K1 Þika cos ga eika x sin ga K2 ikb sin gb eikb x sin gb ; ( ¼

ð1 þ K1 Þika sin ga K2 ikb cos gb

) eika h cos ga eikx x

ð1K1 Þika cos ga K2 ikb sin gb

kx ¼ ka sin ga ¼ kb sin gb ;

ð53Þ

the 2nd term of the first set of scattered waves is: " #( ) Z 1 a1 ðkÞ ik nb dk eikx b1 ðkÞ na ik 1 " # " #( ) Z 1 1 X A1;n ik nb in za;n =na 0 ikx dk e ¼ zb;n =nb B1;n na ik n ¼ 1 ip 0 1 2

Z in ¼ ip n ¼ 1 1 X

ikza;n 6 na 1 6 eikx 6 4 1 za;n

3

( ) 7 A1;n 7 dk ikzb;n 7 B1;n 5 zb;n

ð54Þ

nb

and finally the 3rd term of the second set of scattered waves is: Z

1

" eikx

1

ik

#(

nb

na ik

a2 ðkÞ b2 ðkÞ #

)

dk

" #( ) a1 ðkÞ P11 P12 1 dk b1 ðkÞ na ik RðkÞ P21 P22 1 # " " #" Z 1 1 ikx n X i za;n =na P11 P12 ik nb e ¼ 0 P21 P22 n ¼ 1 ip ik 1 RðkÞ na Z

¼

"

1

eikx

nb

ik

2

Z 1 X in ¼ ip n ¼ 1

ðikP11 þ nb P21 Þza;n 6 1 na eikx 6 6 6 1 RðkÞ 4 ðna P11 ikP21 Þza;n

na

0 zb;n =nb

ðikP12 þ nb P22 Þzb;n

nb ðna P12 ikP22 Þzb;n

nb

#(

A1;n B1;n

)

dk

3 7 (A ) 7 1;n 7dk : 7 B1;n 5

ð55Þ On the surface of the tunnel, r ¼a, displacements will be computed from the potential functions using cylindrical coordinates. From the scattered and diffracted P(f)- and SV(c)-wave potentials, which are all in cylindrical coordinates, the radial and angular components are available from Eq. (17): ur ¼

@f 1 @c þ ; r @y @r

uy ¼

1 @f @c ;  r @y @r

ð56Þ

from which the horizontal x- and vertical y-components of displacements are given by !  !  ux ur cos y sin y ¼ : ð57Þ uy uy sin y cos y For the longitudinal P- and shear SV-waves present, the amplitudes of the displacement components, and of the corresponding phases are:

1=2 jux j ¼ ½Reðux Þ2 þ ½Imðux Þ2 ð58aÞ

1=2 juy j ¼ ½Reðuy Þ2 þ ½Imðuy Þ2   Imðux Þ Reðux Þ   1 Imðuy Þ wy ¼ tan Reðuy Þ

wx ¼ tan1

ð58bÞ

where Re(  ) and Im(  ) are the real and imaginary parts of complex expressions. Note that the inverse tangent (arctangent) always returns a value of the phase in the range of [  p, p] or [0, 2p], but that the actual phases vary continuously beyond the 2p interval. 7.2. The 3D surface displacement plots Fig. 3 shows the horizontal and vertical components of displacement amplitudes (|ux|, left and |uy|, right) along the surface of the half-space, normalized with respect to the amplitude of input plane P-waves |fi|. These are threedimensional plots of displacement amplitudes plotted against the dimensionless distance x/a and the dimensionless frequency Z ¼ oa/pcb. In Fig. 3, the plane P-wave has an incident angle of g ¼301 with respect to the vertical. The displacement amplitudes along the surfaces are plotted from x/a ¼  4 to x/a¼ + 4, and for the dimensionless frequency Z in the range 0 o Z r6. The circular tunnel of the outer radius a is at a depth of h below the surface, with h/a¼ 1.50. The inner radius of the tunnel is b, with the ratio b/a¼0.85. The ratio of the mass densities of the tunnel to that of the half-space is r1/r ¼1.50, and the ratio of their shear moduli is m1/m ¼3.00, with the same Poisson ratio of u¼ u1 ¼ 0.25 for both the half-space and tunnel medium. This makes the tunnel stiffer than the surrounding medium. In the absence of the underground tunnel, for a uniform half-space the free-field amplitude of the ground displacement would be constant at juffx j ¼ 1:12 and juffy j ¼ 1:69, respectively, for displacements in the x- and y-directions. Lee and Cao [15], Cao and Lee [2,3], Trifunac [51] and Wong and Trifunac [64], in their work on wave diffraction around shallow circular and elliptical canyons, noted that the incident waves at low frequency, with wavelengths much longer than the inhomogeneity, ‘‘do not see’’ the presence of the obstacle as much as the higher-frequency (with shorter wavelengths) incident waves. The same was confirmed by Lee and Karl [17,18] for the same underground model studied here for incident P- and SVwaves. There the underground inhomogeneity was a circular cavity, like the case studied here, but without a tunnel wall. As would be expected, similar conclusions can be drawn from Fig. 3 for both the horizontal and vertical components, where at low frequencies, 0 o Z{1, the displacement amplitudes deviate only slightly from the free-field amplitudes (plotted at Z ¼0). As the frequencies increase, the oscillations of the amplitudes increase, and the amplification of the amplitudes at each frequency becomes higher. For an incident angle of g ¼301, a displacement

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

889

Fig. 3. Displacement amplitudes along the half-space surface vs. dimensionless frequency and distance for g ¼ 301 from ‘‘Exact’’ solutions. Left: juffx j ¼ 1:12, jux jmax ¼ 2:90 at Z ¼ 4.4 and x/a¼  1.0. Right: juffy j ¼ 1:69, juy jmax ¼ 1:69 at Z ¼ 5.3 and x/a¼  1.10.

amplitude approaching 4 was observed at high frequencies ðjux jmax ¼ 3:86 at Z ¼ 5:3Þ, corresponding to an amplification almost as high as 2. Trifunac [51], Wong and Trifunac [64], Lee and Trifunac [21], Cao and Lee [2,3], Todorovska and Lee [41], Lee and Wu [22,23], Lee and Manoogian [19,20], Manoogian and Lee [34,35], Lee et al. [16,20], Liang et al. [28,29], Dermendjian et al. [7] and Dermendjian and Lee [6], in their studies of the diffraction and scattering of anti-plane SH waves by surface and sub-surface inhomogeneities, such as surface canyons, valleys, and foundations, underground cavities, and tunnels of circular or arbitrary shapes, all observed this typical diffraction pattern for the SH waves. This pattern shows specific characteristics both in front of and behind the inhomogeneity that result from the scattering and diffraction of incoming waves. At the half-space surface, above and on the left side of the inhomogeneity (the tunnel in our case), with the tunnel acting as a barrier, significant amount of the energy of the incoming waves from the left is reflected back, in our examples in the direction of the negative x-axis. This results in a standingwave pattern superimposed on the motions that are progressing to the right [51]. On the half-space surface to the right of the tunnel, with the tunnel acting as a barrier, a shadow zone is created, with waves of decreasing amplitudes that are also smoother, relative to the standing-wave pattern on the left side. Such wave patterns have long been observed in all of the SH wave papers cited above and by other authors. The present and subsequent figures in this work show that the same wave pattern can also be observed for the in-plane P- and SV-waves. The in-plane waves have two components of motion, and they are further complicated by the presence of waves of both the longitudinal (P-) and transverse (SV-) type, so one of the two components can be more dominant. For the case of P-waves with an incident angle of g ¼301 with respect to the vertical direction, Fig. 3 shows that the free-field waves have larger amplitudes of the vertical (y-) component than of the horizontal (x-) component, juffx j ¼ 1:12 and juffy j ¼ 1:69, and the pattern of standing, oscillatory waves on the left versus shadowy, smoother waves on the right is seen to be more significant for the vertical (y-) than for the horizontal (x-) component.

Fig. 4 shows waves in the same geometry as those in Fig. 3, but with an incident angle of g ¼601 with respect to the y-direction. The only difference is that this time the motions are ‘‘more horizontal’’ than in the previous example, with the free-field waves having higher amplitudes juffx j ¼ 1:73 and juffy j ¼ 1:00, and the pattern of standing, oscillatory amplitudes on the left versus shadowy and smoother motions on the right is seen to be more significant in the horizontal (x-) than in the vertical (y-) direction. For incident angle of, g ¼ 601 in Fig. 4 a displacement amplitude higher than 4 is observed at a high frequency ðjux jmax ¼ 4:43 at Z ¼ 6:0Þ, corresponding to an amplification slightly higher than 2, Figs. 3 and 4 present the x- and y-displacement amplitudes when the zero-stress boundary condition on the half-space surface and the stress and displacement continuity conditions on the surface of the tunnel are imposed. Next we investigate how a simplified solution will change the resulting displacement amplitudes. Figs. 5 and 6 represent the x- and y-displacement amplitudes when the half-space stress-free boundary conditions are ‘‘relaxed,’’ for incident P-waves with angles of g ¼301 and 601, relative to the vertical. Essentially the same results seen in Figs. 3 and 4 are also found in Figs. 5 and 6. The examples shown suggest that ‘‘relaxing’’ the half-space stress-free boundary conditions will not result in serious departures from the calculations based on the exact treatment and solution of the boundary-valued problem. This is an important observation because by allowing the half-space stress-free boundary conditions to be ‘‘relaxed,’’ the numerical calculations are much simpler (Section 6). A physical explanation for justifying a ‘‘relaxed’’ approach may also be sought in the fact that the additional waves resulting from the presence of the scattered waves reflecting from the half-space surface, for the model studied in this paper, contribute only a small part to the total motions. The waves are scattered and diffracted from the surface of the underground tunnel, and as those reach the surface of the half-space, some of their energy has been diminished already by geometric spreading. The first-order, simplified, approximate solution of the present problem, namely the so-called ‘‘relaxed’’ solution, may give a sufficiently good approximation for many applications.

890

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

Fig. 4. Displacement amplitudes along the half-space surface vs. dimensionless frequency and distance for g ¼601 from ‘‘Exact’’ solutions. Left: juffy j ¼ 1:73, jux jmax ¼ 4:43 at Z ¼ 6.0 and x/a ¼  0.5. Right: juffy j ¼ 1:00, juy jmax ¼ 3:43 at Z ¼ 4.4 and x/a¼  0.55.

Fig. 5. Displacement amplitudes along the half-space surface vs. dimensionless frequency and distance for c ¼ 301 from ‘‘Relaxed’’ solutions. Left: juffy j ¼ 1:12, jux jmax ¼ 2:96 at Z ¼1.2 and x/a ¼  0.75. Right: juffy j ¼ 1:69, juy jmax ¼ 3:59 at Z ¼ 5.6 and x/a ¼  1.00.

7.3. The 2-D surface displacement plots Figs. 7a and b show plots of the displacement amplitudes at the half-space surface above the tunnel, for four angles of incidence of P-waves, and for both the x- and y-components of motion. All plots in Fig. 7a are for the same dimensionless frequency of Z ¼ oa/(pcb)¼0.6, while those in Fig. 7b are for Z ¼6.0. In both figures, the same ratios of shear moduli (m1/ m ¼3.00), mass densities (r1/r ¼1.50), and Poisson ratios (u ¼u1 ¼0.25) are used as in the previous examples. Also, the same ratios of depth-to-radius (h/a¼1.50) and inner-to-outer radius (b/a¼ 0.85) are used. The left-hand column shows the horizontal x-component of motions, while the right-hand column is for the vertical motions. The graphs in each plot show the

amplitudes vs. the dimensionless coordinate x/a along the horizontal axis, in the range [ 4, 4]. In each plot, the horizontal dashed line shows the free-field displacement amplitude at the half-space surface in the absence of the underground tunnel. Fig. 7a is for low frequency, Z ¼0.6, which corresponds to waves with wavelength equal to 10/3 times the outer radius of the tunnel. At such a long period, the amplitudes along the half-space surface are less oscillatory and do not deviate much from the free-field amplitude. Fig. 7b illustrates the motions for the incident waves with a much higher frequency, Z ¼6.0, corresponding to waves with wavelength equal to 1/3 of the outer radius of the tunnel. At such frequencies, the amplitudes along the half-space surface are more oscillatory and deviate more from the free-field amplitude.

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

891

Fig. 6. Displacement amplitudes along the half-space surface vs. dimensionless frequency and distance for c ¼ 601 from ‘‘Relaxed’’ solutions. Left: juffx j ¼ 1:73, jux jmax ¼ 3:77 at Z ¼6.0 and x/a ¼  0.50. Right: juffy j ¼ 1:00, juy jmax ¼ 2:82 at Z ¼ 4.4 and x/a¼  0.15.

Surface Displacement Above Circular Tunnel: Incident P Waves Tunnel: ρ1/ρ = 1.50, μ1/μ = 3.00. ν1 = 0.25, ν = 0.25

Surface Displacement Above Circular Tunnel: Incident P Waves Tunnel: ρ1/ρ = 1.50, μ1/μ = 3.00. ν1 = 0.25, ν = 0.25 depth/radius, h/a = 1.50 inner/outer Radius, b/a = 0.85 -2 0

4

depth/radius, h/a = 1.50 inner/outer Radius, b/a = 0.85 -2 0

2

|ux|

h

|uy|

4

η = 0.60

γ = 85°

|uy|

η = 6.00 Half-Space B.C. Relaxed Imposed

2

4

4 Free-field Ave. Amplitudes

γ = 60°

2

4

γ = 30°

2

γ = 60°

2 Displacement Amplitudes

Displacement Amplitudes

h γ = 85°

Half-Space B.C. Relaxed Imposed

2

2

|ux|

4

4

γ = 30° 2

4 γ = 0°

γ = 0°

2

2

0

0 -4

-2

0

2

-2 x/a

0

2

4

-4

-2

0

2

-2

0

2

4

x/a

Fig. 7. (a) Displacement amplitudes along the half-space surface for 4 angles of incidence and Z ¼ 0.6 and (b) displacement amplitudes along the half-space surface for 4 angles of incidence and Z ¼ 6.0.

892

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

Table 1 summarizes the x- and y-components of free-field amplitudes, juffx jandjuffy j, for the four angles of incidence: Plotted in Figs. 7a and b are the amplitudes of the two types of results, calculated with the half-space stress-free boundary conditions: (i) ‘‘Imposed’’ (Section 4) and (ii) ‘‘Relaxed’’ (Section 6).

Those are represented by the solid lines and dashed lines, respectively, and are referred to as the ‘‘imposed’’ and ‘‘relaxed’’ wave solutions. In Figs. 7a and b, the two bottom plots for the x- and y-components are for g ¼ 01, the case of vertical incidence, which produces symmetric plots for both the imposed (solid) and relaxed (dashed) solutions. Both solutions produce very similar Table 1 Free-field displacement amplitudes.

g (deg)

juffx j

juffy j

0 30 60 85

0 1.12 1.73 0.97

2 1.69 1.00 0.35

wiggles; although the relaxed waves appear to have higher maximum amplitudes in the vertical motions, both maxima occur at almost the same x/a on the half-space surface. The maximum amplitude at the dimensionless frequency of Z ¼6.00 (Fig. 7b) is around 3 for the imposed waves and almost 4 for the relaxed waves. An amplitude of 4 here, for a free-field amplitude of 2, thus corresponds to an amplification of 2. The plots for g ¼301 and 601, for both the imposed and relaxed solutions show similar wave forms, and their maximum amplitudes are close in both location and in amplitude. The top plots are for g ¼851, which is a case of almost horizontal incidence. Figs. 8a and b illustrate the displacement amplitudes at the outer surface of the underground tunnel, again for four angles of incidence of P-waves, and for both the x- and y-components of motions. Fig. 8a shows plots for the dimensionless frequency of Z ¼ 0.6, and Fig. 8b presents the plots for Z ¼6.0. Again, all of the results are for the same ratios of shear moduli (m1/m ¼3.00), mass densities (r1/r ¼1.50) and Poisson ratios (u ¼u1 ¼0.25) as in all previous examples. The same ratios of depth-to-radius (h/a ¼1.50) and inner-to-outer radius (b/a¼0.85) are used as well. The graphs in each plot show the amplitudes along the y-axis vs. the angle y in degrees measured along the tunnel surface, using the polar coordinate system shown in Fig. 1. Thus, y ¼ 01 is the rightmost point on the tunnel surface, y ¼901 is its topmost point, y ¼1801is the leftmost point, and y ¼2701 is the bottommost point of the tunnel surface. Each plot shows the amplitudes for the imposed and relaxed solutions, shown by the solid and dashed lines. At a

Tunnel Displacement Below Ground Surface: Incident P Waves

Tunnel Displacement Below Ground Surface: Incident P Waves Tunnel: ρ1/ρ = 1.50, μ1/μ = 3.00. ν1 = 0.25, ν = 0.25

Tunnel: ρ1/ρ = 1.50, μ1/μ = 3.00. ν1 = 0.25, ν = 0.25

depth/radius, h/a = 1.50 inner/outer Radius, b/a = 0.85

depth/radius, h/a = 1.50 inner/outer Radius, b/a = 0.85 -2

4 |ux|

-2

0

180°

γ = 85°

0°

2

|ux|

0°

η = 6.0

270°

2

270°

|uy|

90°

h 180°

Half-Space B.C. Relaxed Imposed

4

4

Displacement Amplitudes

Displacement Amplitudes

0

4

90°

2 h

2

η = 0.60

|uy|

γ = 60°

2

4

γ = 30°

2

2

4

2

4

4

Half-Space B.C. Relaxed Imposed

γ = 0°

2

2

0

0 0

90

180

270

0

90

Angle - degrees

180

270

360

0

90

180

270

0

90

180

270

360

Angle - degrees

Fig. 8. (a) Displacement amplitudes along the tunnel surface for 4 angles of incidence and Z ¼ 0.6 and. (b) displacement amplitudes along the tunnel surface for 4 angles of incidence and Z ¼ 6.0.

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

dimensionless frequency of Z ¼6.00, corresponding to the shear waves with wavelength 1/6 that of the tunnel outer diameter, the amplitudes of both types of waves are very oscillatory along the tunnel surface. The wave amplitudes for dimensionless frequency Z ¼6.0 in Fig. 8b, and for g ¼301, 601, and 851, all show large differences between the imposed and the relaxed solutions. With the tunnel at a depth of h/a ¼1.5 below the half-space surface, the presence of the half-space surface does contribute to the large differences between the imposed waves and the relaxed solutions. In the relaxed solution, the scattered waves do not ‘‘see’’ the half-spaced surface, and hence they do not come back toward the tunnel. A plot of the tunnel surface displacement amplitudes at other, smaller input frequencies shows the similar trends in the differences between the imposed and relaxed solutions, and, as Figs. 3–6 show, they are in better agreement on the half-space surface.

8. Conclusion In summary, the following observations can be made: 1. The incident P-waves with low frequency ‘‘do not see’’ the presence of the obstacle as much as do the higher-frequency (shorter-wavelength) incident P-waves. 2. As the frequencies increase, the variations of the displacement amplitudes increase, along both the half-space and the tunnel surfaces, and the amplification of displacement amplitudes becomes increased. 3. In this work, two types of solutions were considered: with ‘‘imposed’’ (Section 4) and ‘‘relaxed’’ (Section 6) boundary conditions for the reflected waves from the half-space surface. The examples show that the amplitudes of the two types of solutions give similar results. The two solutions do differ, however, along the tunnel surface.

References [1] Abramowitz M, Stegun IA. Handbook of mathematical functions, with formulas, graphs, and mathematical tables. New York, NY: Dover Publications; 1972. [2] Cao H, Lee VW. Scattering of plane SH waves by circular cylindrical canyons with variable depth-to-width ratio. Eur J Earthquake Eng III 1989(2):29–37. [3] Cao H, Lee VW. Scattering and diffraction of plane P-waves by circular cylindrical canyons with variable depth-to-width ratio. Int J Soil Dyn Earthquake Eng 1989;9(3):141–50. [4] Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equations. Bull Seismol Soc Am 1977;67(6):1529–40. [5] Davis CA, Lee VW, Bardet JP. Transverse response of underground cavities and pipes to incident SV waves. Int J Earthquake Eng Struct Dyn 2001;30(3): 383–410. [6] Dermendjian N, Lee VW. Moment solutions of anti-plane (SH) diffraction around arbitrary-shaped rigid foundations on a wedge-shaped half-space. ISET J Earthquake Technol 2003;40(4):161–72. [7] Dermendjian N, Lee VW, Liang J-W. Anti-plane deformations around arbitrary-shaped canyons on a wedge-shaped half-space: moment method solutions. Int J Earthquake Eng Eng Vib 2003;2(2):281–8. [8] Emerman SH, Stephen RA. Comment on ‘‘Absorbing boundary conditions for acoustic and elastic wave equation,’’ by R. Clayton and B. Engquist. Bull Seismol Soc Am 1983;73(2):661–5. [9] Fu L, Bouchon M. Discrete wavenumber solutions to numerical wave propagation in piecewise heterogeneous media—I. Theory of two-dimensional SH case. Geophys J Int 2004;157(2):481–98. [10] Fuyuki M, Matsumoto Y. Finite difference analysis of Rayleigh wave scattering at a trench. Bull Seismol Soc Am 1980;70(6):2051–69. [11] Kawase H. Time-domain response of a semi-circular canyon for incident SV, P, and Rayleigh waves calculated by the discrete wavenumber boundary element method. Bull Seismol Soc Am 1988;78(4):1415–37. [12] Lamb H. On the propagation of tremors over the surface of elastic solid. Phil Trans R Soc London A 1904;203.

893

[13] Lee VW. On deformations near circular underground cavity subjected to incident plane SH waves. In: Symposium of applications of computer methods in engineering, August 23–26. Los Angeles, USA: Univ of Southern California; 1977. p. 951–61. [14] Lee VW. Three-dimensional diffraction of plane P, SV and SH waves by a hemispherical alluvial valley. Int J Soil Dyn Earthquake Eng 1984;3(3): 133–44. [15] Lee VW, Cao H. Diffraction of SV Waves by circular cylindrical canyons of various depths. J Eng Mech ASCE 1989;115(9):2035–56. [16] Lee VW, Chen S, Hsu LR. Antiplane diffraction from canyon above a subsurface unlined tunnel. J Eng Mech ASCE 1999;125(6):668–75. [17] Lee VW, Karl J. Diffraction of SV waves by underground, circular, cylindrical cavities. Int J Soil Dyn Earthquake Eng 1992;11(8):445–56. [18] Lee VW, Karl J. Diffraction of elastic plane p waves by circular, underground unlined tunnels. Eur J Earthquake Eng 1993;6(1):29–36. [19] Lee VW, Manoogian ME. Surface motion above an arbitrary shape underground cavity for incident SH waves. Eur J Earthquake Eng 1995;7(1):3–11. [20] Lee VW, Manoogian ME, Chen S. Antiplane deformations near a surface rigid foundation above a subsurface rigid circular tunnel. Int J Earthquake Eng Eng Vib 2002;1(1):27–35. [21] Lee VW, Trifunac MD. Response of tunnels to incident SH waves. J Eng Mech ASCE 1979;105:643–59. [22] Lee VW, Wu XY. Application of the weighted residual method to diffraction by 2-D canyons of arbitrary shape, I: Incident SH waves. Int J Soil Dyn Earthquake Eng 1994;13(5):355–64. [23] Lee VW, Wu XY. Application of the weighted residual method to diffraction by 2-D canyons of arbitrary shape, II: incident P, SV & Rayleigh waves. Int J Soil Dyn Earthquake Eng 1994;13(5):365–73. [24] Liang JW, Ba Z, Lee VW. Diffraction of SV waves by a shallow circular-arc canyon in a saturated poroelastic half-space. Int J Soil Dyn Earthquake Eng (Spec Issue Honoring Biot) 2006;26(6–7):582–610. [25] Liang JW, Yan L, Lee VW. Scattering of plane P waves around circular-arc layered alluvial valleys: analytical solution. ACTA Seismol Sinica 2001;14(2):176–95. (Chinese Translation: ACTA Seismol Sinica 23(2): 167–84). [26] Liang JW, Yan L, Lee VW. Effect of a covering layer in a circular-arc canyon on incident plane SV waves. ACTA Seismol Sinicaica 2001;14(6):660–75. [27] Liang JW, Yan L, Lee VW. Scattering of incident plane P-waves by a circular-arc canyon with a covering layer. ACTA Mech Solida Sinica 2002;23(4). [28] Liang JW, Zhang Y, Gu X, Lee VW. Surface motion of circular-arc layered alluvial valley for incident plane SH waves. Chin J Geotech Eng 2000;22(4):396–401. [29] Liang JW, Zhang Y, Gu X, Lee VW. Scattering of plane elastic SH waves by circular-arc layered canyons. J Eng Vib 2003;16(2):158–65. [30] Liang JW, Zhang H, Lee VW. A series solution for surface motion amplification due to underground twin tunnels: incident SV waves. Int J Earthquake Eng Eng Vib 2003;2(2):289–98. [31] Liang JW, Zhang H, Lee VW. A series solution for surface motion amplification due to underground group cavities: incident P waves. ACTA Seismol Sinicaica 2004;17(3):296–307. [32] Liang JW, Zhang H, Lee VW. An analytical solution for dynamic stress concentration of underground twin cavities due to incident SV waves. J Vib Eng 2004;17(2):132–40. [33] Liang JW, Zhang H, Lee VW. An analytical solution for dynamic stress concentration of underground cavities under incident P waves. J Geotech Eng 2004;26(6):815–9. [34] Manoogian ME, Lee VW. Diffraction of SH-waves by subsurface inclusions of arbitrary shape. J Eng Mech ASCE 1996;122(2):123–9. [35] Manoogian ME, Lee VW. Antiplane deformations near arbitrary-shape Alluvial valleys. ISET J Earthquake Technol 1999;36(2):107–20. [36] Pao Y-H, Mow CC. Diffraction of elastic waves and dynamics stress concentrations. New York, N.Y: Crane, Russak & Company Inc.; 1973. [37] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran, the art of scientific computing, 2nd ed. Cambridge, UK: Cambridge University Press; 1992. [38] Sanchez-Sesma FJ, Campillo M, Irikura K. A note on the Rayleigh hypothesis and the Aki-Larner method. Bull Seism Soc Am 1989;79:1995–9. [39] Schiff AJ. Northridge earthquake—lifeline performance and post-earthquake response. In: ASCE Technical council on lifeline earthquake engineering, Monograph no. 8, 1995. p. 339. [40] Thiruvenkataher VR, Viswananthan K. Dynamic response of an elastic half space with cylindrical cavity to time dependent surface tractions over the boundary of the cavity. J Math Mech 1965;12:541–71. [41] Todorovska MI, Lee VW. A note on response of shallow circular valleys to Rayleigh waves: analytical approach. Earthquake Eng Eng Vib 1990;10(1):21–34. [42] Todorovska MI, Lee VW. Surface motion of shallow circular alluvial valleys for incident plane SH waves—analytical solution. Int J Soil Dyn Earthquake Eng 1991;10(4):192–200. [43] Todorovska MI, Lee VW. A note on scattering of Rayleigh waves by shallow circular canyons: analytical approach. ISET J Earthquake Technol 1991;28(2):1–16. [44] Todorovska MI, Trifunac MD. Antiplane earthquake waves in long structures. ASCE, EMD 1989;115(2):2687–708. [45] Todorovska MI, Trifunac MD. Propagation of earthquake waves in buildings with soft first floor. ASCE, EMD 1990;116(4):892–900.

894

C.H. Lin et al. / Soil Dynamics and Earthquake Engineering 30 (2010) 879–894

[46] Todorovska MI, Trifunac MD. Distribution of pseudo spectral velocity during Northridge, California earthquake of 17 January 1994. Soil Dyn Earthquake Eng 1997;16(3):173–92. [47] Todorovska MI, Trifunac MD. Amplitudes, polarity and time of peaks of strong ground motion during the 1994 Northridge, California earthquake. Soil Dyn Earthquake Eng 1997;16(4):235–58. [48] Todorovska MI, Trifunac MD. Liquefaction opportunity mapping via seismic wave energy. J Geotech Geoenviron Eng ASCE 1999;125(12):1032–42. [49] Todorovska MI, Al Rjoub Y. Plain strain soil-structure interaction model for a building supported by a circular foundation embedded in a poroelastic halfspace. Soil Dyn Earthquake Eng 2006;26(6–7):694–707. [50] Trifunac MD. A method for synthesizing realistic strong ground motion. Bull Seism Soc Am 1971;61(6):1739–53. [51] Trifunac MD. A note on scattering of plane SH waves by a semi-cylindrical canyon. Int J Earthquake Eng Struct Dyn 1973;1:267–81. [52] Trifunac MD. Empirical criteria for liquefaction in sands via standard penetration tests and seismic wave energy. Soil Dyn Earthquake Eng 1995;14(6):419–26. [53] Trifunac MD, Todorovska MI. Northridge, California, earthquake of 17 January 1994: density of pipe breaks and surface strains. Soil Dyn Earthquake Eng 1997;16(3):193–207. [54] Trifunac MD, Todorovska MI. Northridge, California, earthquake of 1994: density of red-tagged buildings versus peak horizontal velocity and intensity of shaking. Soil Dyn Earthquake Eng 1997;16(3):209–22. [55] Trifunac MD, Todorovska MI. Response spectra and differential motion of columns. Earthquake Eng Struct Dyn 1997;26(2):251–68. [56] Trifunac MD, Todorovska MI. Nonlinear soil response as a natural passive isolation mechanism—the 1994 Northridge, California Earthquake. Soil Dyn Earthquake Eng 1998;17(1):41–51.

[57] Trifunac MD, Todorovska MI. Damage distribution during the 1994 Northridge, California, earthquake in relation to generalized categories of surface geology. Soil Dyn Earthquake Eng 1998;17(4):239–53. [58] Trifunac MD, Todorovska MI, Ivanovic´ SS. A note on distribution of uncorrected peak ground accelerations during the Northridge, California, Earthquake of 17 January 1994. Soil Dyn Earthquake Eng 1994;13(3):187–96. [59] Trifunac MD, Todorovska MI, Ivanovic SS. Peak velocities and peak surface strains during Northridge, California, earthquake of 17 January 1994. Soil Dyn Earthquake Eng 1996;15(5):301–10. [60] Trifunac MD, Todorovska MI, Lee VW. The Rinaldi strong motion accelerogram of the Northridge, California, Earthquake of 17 January 1994. Earthquake Spectra 1998;14(1):225–39. [61] Udwadia FE, Trifunac MD. Characterization of response spectra through statistics of oscillator response. Bull Seism Soc Am 1974;64(1):205–19. [62] Wong HL, Trifunac MD. Interaction of a shear wall with the soil for incident plane SH-waves: elliptical rigid foundation. Bull Seism Soc Am 1974;64(6): 1825–42. [63] Wong HL, Trifunac MD. Scattering of plane SH-waves by a semi-elliptical canyon. Int J Earthquake Eng Struct Dyn 1974;3(2):157–69. [64] Wong HL, Trifunac MD. Surface motion of a semi-elliptical alluvial valley for incident plane SH-waves. Bull Seism Soc Am 1974;64(5):1389–408. [65] Wong HL, Trifunac MD. Two-dimensional, antiplane, building-soil-building interaction for two or more buildings and for incident plane SH-waves. Bull Seism Soc Am 1975;65:1863–85. [66] Wong HL, Trifunac MD. Generation of artificial strong motion accelerograms. Int J Earthquake Eng Struct Dyn 1979;7(6):509–27. [67] Zienkiewicz OC, Taylor RL. The finite element method, vol. 1, 5th ed. Oxford: Butterworth-Heinemann; 2000.