A 2.5-D dynamic model for a saturated porous medium. Part II: Boundary element method

Available online at www.sciencedirect.com

International Journal of Solids and Structures 45 (2008) 359–377 www.elsevier.com/locate/ijsolstr

A 2.5-D dynamic model for a saturated porous medium. Part II: Boundary element method Jian-Fei Lu a, Dong-Sheng Jeng

b,*

, Sally Williams

c

a

b

Department of Civil Engineering, Jiangsu University, Zhenjiang, Jiangsu 212013, PR China Division of Civil Engineering, School of Engineering, Physics and Mathematics, University of Dundee, Dundee DD1 4HN, Scotland, UK c School of Civil Engineering, The University of Sydney, NSW 2006, Australia Received 15 July 2006; received in revised form 8 May 2007 Available online 14 September 2007

Abstract The 3-D boundary integral equation is derived in terms of the reciprocal work theorem and used along with the 2.5-D Green’s function developed in Part I [Lu, J.F., Jeng, D.S., Williams, S., submitted for publication. A 2.5-D dynamic model for a saturated porous medium: Part I. Green’s function. Int. J. Solids Struct.] to develop the 2.5-D boundary integral equation for a saturated porous medium. The 2.5-D boundary integral equations for the wave scattering problem and the moving load problem are established. The Cauchy type singularity of the 2.5-D boundary integral equation is eliminated through introduction of an auxiliary problem and the treatment of the weakly singular kernel is also addressed. Discretisation of the 2.5-D boundary integral equation is achieved using boundary iso-parametric elements. The discrete wavenumber domain solution is obtained via the 2.5-D boundary element method, and the space domain solution is recovered using the inverse Fourier transform. To validate the new methodology, numerical results of this paper are compared with those obtained using an analytical approach; also, some numerical results and corresponding analysis are presented.  2007 Elsevier Ltd. All rights reserved. Keywords: 2.5-D boundary element method; Discrete wavenumber; The porous medium; Moving loads; Biot’s theory; The Fourier transform

1. Introduction Boundary element method (BEM) has been demonstrated as an effective tool in the numerical solution of complex engineering problems across a wide variety of disciplines [Beskos, 1987, 1997]. In particular, the application of BEM to geotechnical problems has been well researched. A complete boundary integral equation for saturated porous media was formulated by Dominguez (1991, 1992) using three solid displacements and a fourth variable proportional to the pore pressure. By considering the same fundamental variables,

*

Corresponding author. E-mail addresses: [email protected] (J.-F. Lu), [email protected] (D.-S. Jeng), [email protected] (S. Williams). 0020-7683/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2007.07.026

360

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

Cheng et al. (1991) developed a boundary integral equation for porous media in the frequency domain using a 3-D boundary integral equation. The 3-D problem of wave scattering within a porous medium in the frequency domain was later addressed by Zimmerman and Stern (1993). By incorporating the Green’s function into the transformed Laplace domain and using a convolution quadrature method, a time domain boundary integral equation method was developed by Schanz (2001). The BEM for porous media has also been successfully applied in the dynamic analysis of the interaction between structures and their surrounding soils (Kattis et al., 2003; Maeso et al., 2005). In contrast to the more conventional Finite Element Method (FEM), a complete solution can be obtained using BEM from boundary values only. Although this achieves substantial savings in modeling effort, a full 3D analysis using BEM can still require considerable memory availability and computation time. The 2.5-D BEM (Luco et al., 1990; Zhang and Chopra, 1991; Pedersen et al., 1994; Papageorgiou and Pei, 1998) is a valuable technique which can be used to simplify a problem when the applied load and structural response are 3-D, yet the structure itself can be considered as 2-D. The 2.5-D BEM code and discretisation are much simpler than that for the full 3-D BEM and computational requirements are significantly reduced. The 2.5-D BEM has been well developed for problems of elastic media. Tadeu and Kausel (2000) derived a 2.5-D Green’s function for an isotropic elastic medium to facilitate the associated 2.5-D BEM analysis. Based on Green’s function for a moving point load with a constant velocity, Luco et al. (1990) proposed a 2.5-D Indirect Boundary Integral Equation Method to evaluate the 3-D response of an infinitely long canyon when subjected to a wave with an arbitrary incident angle. By using the full space Green’s function, Zhang and Chopra (1991) developed a direct BEM to address the 3-D response of a 2-D structure embedded in a viscoelastic half-space. By means of direct 2.5-D BEM, Stamos and Beskos (1996) calculated the 3-D seismic response of an infinitely long lined tunnel in a viscoelastic half space. Papageorgiou and Pei (1998) presented a discrete wavenumber direct BEM for an infinitely long canyon when subjected to a plane elastic wave with an arbitrary incident angle. More recently, Sheng et al. (2005) developed a coupled discrete wavenumber BEM and FEM to evaluate the ground vibration due to the moving load from a railway. They retrieved the ground response in the 3-D space using an inverse Fourier transform method. Although 2.5-D BEM is well established for isotropic elastic media and has been widely applied, the 2.5-D BEM for a saturated porous medium has yet to be developed. In part I of this work (Lu et al., submitted for publication), the 2.5-D Green’s function for a saturated porous medium was derived and discussed. We now follow by establishing the 2.5-D boundary integral equation using the 2.5-D Green’s function and Biot’s theory. It is noted that in the literature ‘2.5-D BEM’ often refers to the discrete wavenumber BEM in which the Green’s function is evaluated by synthesizing a discrete wavenumber Green’s function in order to reduce the order of the singularity (Bouchon, 2003; Gaffet and Bouchon, 1989). In this paper, the 2.5-D BEM is formulated in the axial wavenumber domain. The numerical characteristics of the new 2.5-D BEM for the saturated porous medium are then investigated by consideration of appropriate examples and comparison with previous analytical studies (Lu and Jeng, 2006). 2. Two-and-a-half-dimensional (2.5-D) boundary integral equation for the porous medium In this section, the 3-D boundary integral equation in the frequency domain for a poroelastic medium is derived. Based on the 3-D boundary integral equation and using the Fourier transform method, a new 2.5D boundary integral equation in the frequency-wavenumber domain will be developed. 2.1. Biot’s theory In this section we outline the fundamental equations of Biot’s theory for a saturated porous medium (Biot, 1956, 1962). A more detailed discussion can be found in Section 2 of Part I of this paper (Lu et al., submitted for publication). The constitutive equations for a homogeneous porous medium (Biot, 1956, 1962) are: rij ¼ 2leij þ kdij e  adij p

ð1aÞ

p ¼ aMe þ Mf

ð1bÞ

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

361

And the equations of motion for the bulk porous medium and the pore fluid can be expressed as (Biot, 1956, 1962): €i lui;jj þ ðk þ a2 M þ lÞuj;ji þ aMwj;ji þ F i ¼ qb € u i þ qf w g wi þ KðtÞ  w_ i ui þ m€ aMuj;ji þ Mwj;ji þ fi ¼ qf € k

ð2aÞ ð2bÞ

The physical meanings of variables in (1)–(2) are the same as those in Part I (Lu et al., submitted for publication). Following the procedure outlined in Part I, performing the Fourier transform with respect to time Z þ1 Z þ1 1 f ðtÞeixt dt; f ðtÞ ¼ ð3Þ f^ ðxÞ ¼ f^ ðxÞeixt dx 2p 1 1 on Eqs. (1) and (2), the governing equations for Biot’s theory can be reduced to (Bonnet, 1987) uj;ji þ qg x2 ^ l^ ui;jj þ ðk þ lÞ^ ui  ag ^ p;i ¼ F^ gi ^ p  b3 ^ uj;j ¼ #^f p;jj þ b2 x2 ^

ð4aÞ ð4bÞ

where F^ gi ¼ F^ i  b4 f^ i , #^f ¼ f^ j;j , qg = qb  b4qf, ag = a  b4, b2 = 1/(b1x2), b3 = qfx2  aM/b1, ^ ^ b4 = qfx2b1/M, b1 ¼ M=½mx2  ixðg=kÞKðxÞ, KðxÞ denotes the Fourier transform of K(t) and a caret represents the Fourier transform with respect to time. The infiltration displacement of the pore fluid relative to the solid skeleton can be expressed in terms of the solid skeleton displacement and the pore pressure as ^i ¼ w

b1 ½^ p;i  qf x2 ^ ui  f^ i  M

ð5Þ

2.2. Three-dimensional (3-D) boundary integral equation for the porous medium The 3-D frequency domain boundary integral equation for a poroelastic medium can be derived from the reciprocal work theorem for a porous medium (Cheng et al., 1991) or alternatively by the weighted residual formulation (Dominguez, 1991, 1992). In this study, we apply the reciprocal work theorem by considering two states, represented by superscripts ð1Þ ð2Þ ð2Þ ð1Þ ^ij ^eij þ ^ ^ij ^eij þ ^ r pð1Þ fð2Þ ¼ r pð2Þ fð1Þ

ð6Þ

Integration of Eq. (6) over an arbitrary domain X gives the reciprocal work theorem for the porous medium as Z Z ð1Þ ð2Þ ð2Þ ð1Þ ð^ rij ^eij þ ^ pð1Þ fð2Þ Þ dX ¼ ð^ rij ^eij þ ^ pð2Þ fð1Þ Þ dX ð7Þ X

X

Using (7) and the governing equations for the porous medium, the generalized Betti’s formula for the porous medium has the form Z Z Z Z ð1Þ ð2Þ ð2Þ ð1Þ ð1Þ ð2Þ ð2Þ ð1Þ ð2Þ ð1Þ ð1Þ ^ ð2Þ ð2Þ ^ ð1Þ ^ ð2Þ ^ ^ ^ V V uj  ^tj ^uj  dS  ½^ w ½^tj ^ pð1Þ w  p  dS ¼ ½^ p  p  dX  ½F^ gj ^uj  F^ gj ^uj  dX f f n n S

S

X

X

ð8Þ ^jk nk , w ^ n ¼ ðb1 =MÞð^p;j  qf x2 u ^j Þnj and V^ f ¼ ðb1 =MÞ#^f . where ^tj ¼ r To derive the Somigliana identity for the porous medium, state 1 is defined as the real state, assuming that there are no body forces acting on the solid skeleton and the pore fluid, while state 2 corresponds to the Green’s function of the porous medium. The displacement of the solid frame is obtained by considering the case of a unit point force acting in the i-direction at a point x for equation (8) Z Z Gs Gs ^ ^ ^ Gs ðx; nÞ^pðnÞ  P^ Gs ðx; nÞ^ ^ ^ ui ðx; xÞ ¼ ½U ij ðx; n; xÞtj ðn; xÞ  T ij ðx; n; xÞ^uj ðn; xÞ dSðnÞ  ½W wn ðnÞ dSðnÞ ni i S

S

ð9Þ

362

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

^ Gnis ðx; nÞ are given by where the Green’s function T^ Gij s ðx; nÞ and W ^ Gij;ks þ U ^ Gik;js  þ kU ^ Gil;ls djk  aP^ Gs djk ; ^Gijks ¼ l½U r i

^Gijks ðx; n; xÞnk ðnÞ T^ Gij s ðx; n; xÞ ¼ r

^ Gik;js ðx; n; xÞnk ðnÞ þ kU ^ Gil;ls ðx; n; xÞnj ðnÞ  aP^ Gs ðx; n; xÞnj ðnÞ; W ^ Gs ðx; n; xÞ ^ Gij;ks ðx; n; xÞ þ U ¼ l½U i ni  Gs  ^ b oP i ðx; n; xÞ ^ Gs ðx; n; xÞ nj ðnÞ  qf x 2 U ¼ 1 ij onj M

ð10Þ

Similarly, considering a dilatation source acting at a point x within the pore fluid for equation (8) gives Z Z = = = G G G ^ j f ðx; n; xÞ^tj ðn; xÞ  T^ j f ðx; n; xÞ^uj ðn; xÞ dSðnÞ  ½W ^ n f ðx; n; xÞ^pðn; xÞ ^ pðx; xÞ ¼ ½U S

S =

wn ðn; xÞ dSðnÞ  P^ Gf ðx; n; xÞ^ =

ð11Þ

=

G T^ j f

=

Gf n

G ^j f, U

=

G T^ j f ,

=

Gf n

^ have similar expressions as in (10). It is noted that the Green’s function ^ and where and W W = P^ Gf in (11) is different from the 3-D Green’s function defined in Lu et al. (submitted for publication), the relationship between the two being =

^ Gf ; ^j f ¼ M U U b1 j G

=

M G T^ j f ¼ T^ j f ; b1 G

=

^ Gf ; ^ nf ¼ M W W b1 n

M P^ Gf ¼ P^ Gf b1 =

G

ð12Þ =

G

^ i f ðx; n; xÞ can be derived Moreover, in terms of (8), the following relation between P^ Gi s ðx; n; xÞ and U =

G

^ i f ðx; n; xÞ P^ Gi s ðx; n; xÞ ¼ U

ð13Þ

which is useful for the calculation of the Green’s function for the porous medium. Introducing the following Green’s functions, pseudo-displacement and pseudo-traction vectors ^u3 u ^4 T 2 G ^ ij s U ^ G  ¼ 4 = ¼ ½U mn G ^j f U

^ u ¼ ½^ u1 ^ G U

^ u2

T ¼ ½^ u1 ^ u2 ^ u3 ^ p  ; ^t ¼ ½ ^t1 3 2 G P^ Gi s T^ ij s  ^ G ¼ ½T^ G  ¼ 4 = 5; T mn = G P^ Gf T^ j f

T ^t2 ^t3 ^t4 T ¼ ½ ^t1 ^t2 ^t3 w ^n  3 ^ Gnis W 5; m; n ¼ 1; 2; 3; 4 = G ^ nf W

and combining the integral boundary Eqs. (9) and (11) yields Z Z ^ G ðx; n; xÞ^tj ðn; xÞ dSðnÞ  T^ G ðx; n; xÞ^uj ðn; xÞ dSðnÞ; ^ U ui ðx; xÞ ¼ ij ij S

i; j ¼ 1; 2; 3; 4

ð14Þ

ð15Þ

S

Performing a limit process as outlined by Zimmerman and Stern (1993), the following boundary integral equation for the porous medium can be obtained from Eq. (15) Z Z ^ G ðx; n; xÞ^tj ðn; xÞ dSðnÞ  T^ G ðx; n; xÞ^uj ðn; xÞ dSðnÞ; i; j ¼ 1; 2; 3; 4 cij ðxÞ^ U uj ðx; xÞ ¼ ð16Þ ij ij S

S

in which the free term cij depends on the local boundary geometry around the point x, with x belonging to a  smooth part of the boundary, cij = 1/2dij. Note that the second integral in (16) involving T^ Gij ðx; n; xÞ should be understood in the sense of the Cauchy principal value. 2.3. Two-and-a-half-dimensional (2.5-D) boundary integral equation for the porous medium In this section, the 2.5-D boundary integral equation for a 2-D porous structure experiencing a 3-D deformation is derived from the 3-D boundary integral Eq. (16) using the Fourier transform method. The Fourier transforms with respect to the x3(z)-coordinate and the wavenumber k3(kz) are required, and defined as follows (Sneddon, 1951) Z þ1 Z þ1 1 f ðzÞeikz z dz; f ðzÞ ¼ ð17Þ f ðk z Þeikz z dk z f ðk z Þ ¼ 2p 1 1 where an overbar denotes the Fourier transform with respect to the x3(z) coordinate.

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

363

It is supposed that the porous medium is infinitely long and has a uniform cross-section along the z-axis. Moreover, letting C denote the intersection of the boundary of the porous medium with the x  y plane and taking into account the translation identity of the Green’s function (Kobayashi, 1987), the boundary integral Eq. (16) can be recast in the following form cij ðx? Þ^ uj ðx? ; z; xÞ ¼

Z Z C



þ1 1

Z Z C

^ G ðn?  x? ; n3  z; xÞ^tj ðn? ; n3 ; xÞ dn3 dCðn? Þ U ij

þ1

1

 T^ Gij ðn?  x? ; n3  z; xÞ^uj ðn? ; n2 ; xÞ dn3 dCðn? Þ;

i; j ¼ 1; 2; 3; 4

ð18Þ

where x? = xix + yiy, n? = n1ix + n2iy, and ix,iy are the base vectors. Carrying out the Fourier transform with respect to the z-coordinate on (18), we obtain   Z Z þ1 Z þ1 G ik z ðzn3 Þ ik n Þ z 3 ^ cij ðx? Þ~ uj ðx? ;k z ;xÞ ¼ dz  ½^tj ðn? ;n3 ;xÞe dn3  dCðn? Þ U ij ðn?  x? ;n3  z;xÞe C 1 1   Z Z þ1 Z þ1 G ik z ðzn3 Þ ik z n3 ^  dz  ½^ uj ðn? ;n3 ;xÞe dn3  dCðn? Þ; T ij ðn?  x? ; n3  z;xÞe 1

C

1

i; j ¼ 1; 2; 3;4

Using the definition of the Fourier transform (17), the above equation can be further reduced to Z ~ G ðn?  x? ; k z ; xÞ~tj ðn? ; k z ; xÞ dCðn? Þ cij ðx? Þ~ U uj ðx? ; k z ; xÞ ¼ ij C Z   T~ Gij ðn?  x? ; k z ; xÞ~uj ðn? ; k z ; xÞ dCðn? Þ; i; j ¼ 1; 2; 3; 4

ð19Þ

ð20Þ

C

~ G , T~ G are the Eq. (20) is the 2.5-D integral equation for the porous medium, in which the Green’s function U ij ij 2.5-D Green’s function derived by Lu et al. (submitted for publication). Since half of the boundary variables ~tj ðn? ; k z ; xÞ and u˜j (n?, kz, x) are given a priori by the boundary condition, the boundary integral equation can be solved by an appropriate discretization scheme. If the boundary condition is given in the space domain, then the wavenumber domain boundary condition can be obtained using the Fourier transformation. Once the boundary variables are determined by solving (20), the frequency-wavenumber domain variables can be evaluated by the following formula Z ~ G ðn?  x? ; k z ; xÞ~tj ðn? ; k z ; xÞ dCðn? Þ ~ ui ðx? ; k z ; xÞ ¼ U ij C Z   T~ Gij ðn?  x? ; k z ; xÞ~ uj ðn? ; k z ; xÞ dCðn? Þ; x? 62 C; i; j ¼ 1; 2; 3; 4 ð21Þ C

Once the variables in the wavenumber domain are obtained by (20) and (21), the space domain solution can be determined by synthesizing the wavenumber domain solution via the inverse Fourier transform, which is achieved by the FFT method in this study. 3. The 2.5-D boundary integral equation for wave scattering and moving loading problems Using the 2.5-D boundary integral equation developed in Section 2, we are now able to formulate the 2.5-D boundary integral equations for the wave scattering problem and the moving load problem. 3.1. The 2.5-D boundary integral equation for the wave scattering problem The boundary integral Eq. (20) can be used in the calculation of the wave scattering problem encountered in earthquake engineering. For the wave scattering problem, the total wave field of the porous medium is the

364

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

superposition of the free wave field due to incident waves, i.e. the wave field in the absence of scatter, and the scattered wave field Þ ~tj ðx? ; k z ; xÞ ¼ ~tðf ~ðsÞ j ðx? ; k z ; xÞ þ t j ðx? ; k z ; xÞ; ðf Þ

ðsÞ

~ uj ðx? ; k z ; xÞ þ ~ uj ðx? ; k z ; xÞ uj ðx? ; k z ; xÞ ¼ ~

ð22Þ

It is emphasized that the scattered wave field described in (22) satisfies the boundary integral Eq. (20). Substituting the scattered wave field (22) into (20) and using the wavenumber domain expression of the Green’s function and the free wave field yields the boundary integral equation for the total wave field Z ~ G ðn?  x? ; k z ; xÞ~tj ðn? ; k z ; xÞ dCðn? Þ cij ðx? Þ~ U uj ðx? ; k z ; xÞ ¼ ij C Z  ðf Þ  T~ Gij ðn?  x? ; k z ; xÞ~uj ðn? ; k z ; xÞ dCðn? Þ þ ~ui ðx? ; k z ; xÞ; i; j ¼ 1; 2; 3; 4 C

ð23Þ ðf Þ ~ uj ðn? ; k z ; xÞ

Since the free wave field is known a priori, and half of the total wave field solution u˜j (n?, kz, x) and ~tj ðn? ; k z ; xÞ is given by the boundary condition, the boundary integral equation can be solved. Based on the solution of Eq. (23), the total wave field for a point x? in the calculation domain is obtained as Z ~ G ðn?  x? ; k z ; xÞ~tj ðn? ; k z ; xÞdCðn? Þ ~ U ui ðx? ; k z ; xÞ ¼ ij C Z  ðf Þ  T~ Gij ðn?  x? ; k z ; xÞ~ uj ðn? ; k z ; xÞdCðn? Þ þ ~ui ðx? ; k z ; xÞ; x? 62 C; i; j ¼ 1; 2; 3; 4 C

ð24Þ For a single-phase 2-D elastic medium subjected to an oblique incident wave, translational invariance holds in a reference frame moving with the apparent velocity along the axis of the 2-D structure (Stamos and Beskos, 1996; Papageorgiou and Pei, 1998). Thus, this kind of scattering problem can be reduced to a quasi 2-D problem. However, due to the attenuation of the body waves of the porous medium, translational invariance is lost. Consequently, even for plane waves of the porous medium with an oblique incident angle, the scattering problem cannot be reduced to a quasi 2-D problem as is the case for a singlephase elastic medium. 3.2. The 2.5-D boundary integral equation for moving loading problems Moving load problems are very common in civil and transportation engineering (Fryba, 1999). To date, moving loads problems for porous media have been mainly addressed analytically (Burke and Kingsbury, 1984; Jin et al., 2004). However, for complex topographies, numerical approaches such as FEM and BEM are necessary. The 2.5-D boundary integral Eq. (20) can also be used to resolve moving load problems if the properties of moving loads are appropriately taken into account. For a moving load with a constant velocity v, all the variables have the following expressions uj ðx; tÞ ¼ uj ðx? ; z; tÞ ¼ uj ðx? ; z  vtÞ;

tj ðx; tÞ ¼ tj ðx? ; z; tÞ ¼ tj ðx? ; z  vtÞ

ð25Þ

Applying the Fourier transform with respect to time and the z-coordinate on (25), the following equations are obtained Z þ1 Z þ1 ~ uj ðx? ; k z ; xÞ ¼ uj ðx? ; z  vtÞeiðkz zxtÞ dz dt 1 1 Z þ1 Z þ1  j ðx? ; k z Þ ½uj ðx? ; z  vtÞeikz ðzvtÞ dzeitðxkz vÞ dt ¼ 2pdðx  k z vÞU ð26aÞ ¼ 1

1

~tj ðx? ; k z ; xÞ ¼ 2pdðx  k z vÞT j ðx? ; k z Þ

ð26bÞ

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

365

where Z þ1 0 uj ðx? ; z  vtÞeikz ðzvtÞ dðz  vtÞ ¼ uj ðx? ; z0 Þeikz z dz0 ; 1 1 Z þ1 Z þ1 0 tj ðx? ; z  vtÞeikz ðzvtÞ dðz  vtÞ ¼ tj ðx? ; z0 Þeikz z dz0 T j ðx? ; k z Þ ¼

 j ðx? ; k z Þ ¼ U

Z

þ1

1

ð27Þ

1

Substituting (26) into (20), the following boundary integral equation is obtained Z ~ G ðn?  x? ; k z ; xÞT j ðn? ; k z Þ dCðn? Þ  cij ðx? ÞU j ðx? ; k z Þdðx  k z vÞ ¼ dðx  k z vÞ U ij C  Z G ~   T ij ðn?  x? ; k z ; xÞU j ðn? ; k z Þ dCðn? Þ ; i; j ¼ 1; 2; 3; 4

ð28Þ

C

Applying the inverse Fourier transform with respect to x and considering the properties of Dirac-d function, one has Z ~ G ðn?  x? ; k z ; k z vÞT j ðn? ; k z Þ dCðn? Þ  cij ðx? ÞU j ðx? ; k z Þ ¼ U ij C  Z   j ðn? ; k z Þ dCðn? Þ ; i; j ¼ 1; 2; 3; 4  T~ Gij ðn?  x? ; k z ; k z vÞU ð29Þ C

Eq. (29) is the 2.5-D integral equation for a saturated porous medium in the wavenumber domain when subjected to moving loads. This can be solved to give the wavenumber domain variables within the calculation domain as Z  Z  j ðx? ; k z Þ ¼ ~ G ðn?  x? ; k z ; k z vÞT j ðn? ; k z Þ dCðn? Þ  T~ G ðn?  x? ; k z ; k z vÞU  j ðn? ; k z Þ dCðn? Þ ; U U ij ij C

x? 62 C;

C

i; j ¼ 1; 2; 3; 4

ð30Þ

 j ðn? ; k z Þ, T j ðn? ; k z Þ in the wavenumber domain. The freEqs. (29) and (30) determine the solutions for U quency-wavenumber domain solutions for u˜j (x?, kz, x) and ~tj ðx? ; k z ; xÞ are given by (26), while the space-time domain solutions for ui (x?, z, t) and ti (x?, z, t) are as follows  Z þ1 Z þ1 1 iðk z zxtÞ ~ ui ðx? ; z; tÞ ¼ u ðx ; k ; xÞe dx dk z i ? z 2 ð2pÞ 1 1 Z þ1 Z þ1   1  i ðx? ; k z Þeiðkz zxtÞ dx dk z ¼ dðx  k z vÞU 2p 1 1 Z þ1 Z þ1 1  i ðx? ; k z Þeikz ðzvtÞ dk z ; ti ðx? ; z; tÞ ¼ 1 U T i ðx? ; k z Þeikz ðzvtÞ dk z ¼ ð31Þ 2p 1 2p 1 Moreover, the frequency domain response of the porous medium subjected to moving loads can also be obtained by changing the sequence of integrations in (31), which yields 1 x  x 1 x  x ^ ui ðx? ; z; xÞ ¼ ei v z U ; ^ti ðx? ; z; xÞ ¼ ei v z T i x? ; ð32Þ i x? ; v v v v Eq. (32) is useful in the analysis of the frequency spectrum of the response of the porous medium due to traveling loads. It follows from (32) that the amplitudes of the frequency domain response are independent of the z-coordinate; the variation of the z-coordinate only results in a phase variation. 4. Desingularisation of the 2.5-D boundary integral equation for the porous medium In this section, the treatment of Cauchy principal value integrals in the boundary integral Eq. (20) is discussed. Also, to improve integration efficiency, weakly singular integrals involving logarithmic function are addressed.

366

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

4.1. Desingularisation of Cauchy principal value integrals in the 2.5-D boundary integral equation for the porous medium The Cauchy principal value integrals in (20) need special treatment to be evaluated numerically. There are two approaches for evaluating the Cauchy principal value integrals involved in the boundary integral equation: the first being the direct calculation of the Cauchy principal value integrals (Hall, 1989); the second is the indirect method, which evaluates the Cauchy principal value integrals via the rigid displacement methodology (Banerjee, 1981). To avoid direct evaluation of the Cauchy integral, in this paper an auxiliary problem is introduced to eliminate the Cauchy singularity from the 2.5-D boundary integral equation of the porous medium. According to the analysis of Lu et al. (submitted for publication), the asymptotic expression for the  2.5-D Green’s function T~ Gij (n?  x?, kz, x) in (20) when n? ! x? is equivalent to the combination of the solutions for a static plane strain problem, a static anti-plane elastic problem and a 2-D potential problem. Thus, to eliminate the Cauchy singularity of the 2.5-D Green’s function in (20), the following Green’s functions are introduced for the auxiliary problem 2 Ga 3 2 Ga 3 a a 0 0 T 12 T G12 U 12 U G12 0 0 a a a a 6 UG UG 6TG TG 0 0 7 0 0 7 a a 6 7 6 21 7 Ga 22 22 Ga UG ¼ ½U Gij  ¼ 6 21 ; T ð33Þ ¼ ½T  ¼ 7 6 7 a a ij G G 4 0 5 4 0 0 T 33 0 5 0 U 33 0 a a 0 0 0 W Gn 0 0 0 PG a

a

a

a

a

a

in which U Gij , T Gij , i, j = 1,2, U G33 , T G33 , P G and W Gn are given by   1 1 1 1 1 1 a a Ga ð3  4mÞ log dij þ r;i r;j ; i; j ¼ 1; 2; log ; PG ¼ log U ij ¼ U G33 ¼ 8plð1  mÞ r 2pl r 2p r   1 or a ½ð1  2mÞdij þ 2r;i r;j   ð1  2mÞðr;i nj  r;j ni Þ ; i; j ¼ 1; 2 T Gij ¼  4pð1  mÞr on 1 1 a a r;k nk ; W Gn ¼  r;k nk ; k ¼ 1; 2 ð34Þ T G33 ¼  2pr 2pr pffiffiffiffiffiffiffiffi where nk is the direction cosine for the outward unit normal, r ¼ xk xk , k = 1,2 and m is Poisson’s ratio of the porous medium. Introducing the following pseudo-displacement and pseudo-traction vectors corresponding to the Green’s function (34) ðaÞ

uðaÞ ¼ ½ui T ¼ ½ u1ðaÞ

ðaÞ

u2

ðaÞ

u3

ðaÞ

T pðaÞ  ;

tðaÞ ¼ ½ti T ¼ ½ tðaÞ 1

ðaÞ

t2

ðaÞ

t3

wnðaÞ 

T

ð35Þ

ðeÞ with wðeÞ n ¼ op =on, then, the boundary integral equation for the auxiliary problem has the following form Z a a ðaÞ ðaÞ ðaÞ ðaÞ ð36Þ cij ðx? Þuj ðx? Þ ¼ ½U Gij ðn?  x? Þtj ðn? Þ  T Gij ðn?  x? Þuj ðn? Þ dCðn? Þ C

Note that the 2-D boundary C in the above equation is the same as in the original problem. However, the boundary unit outward normal of the auxiliary problem is opposite to that of the original problem. Consequently, if the original problem occupies a finite region, then, the auxiliary problem will occupy the outer infinite region surrounding the original region and vice versa. ðaÞ If u˜j (x?, kz, x) at point x? in (20) is used as a constant solution for (36), then, tj ðn? Þ in (36) will vanish and the boundary integral Eq. (36) is reduced to Z a ðaÞ cij ðxÞ~ uj ðx? ; k z ; xÞ  c1 ~ T Gij ðn?  x? Þ~uj ðx? ; k z ; xÞ dCðn? Þ; i; j ¼ 1; 2; 3; 4 ð37Þ ui ðx? ; k z ; xÞ ¼ C

in which c1 = 1, 1/2 or 0 for an auxiliary problem with an infinite, a half-infinite or a finite region, respectively. It is noted that in order to be consistent with the sign conventions of the original poroelastic problem a the sign of the traction T Gij ðn?  x? Þ has been reversed.

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

367

Adding (37) to (20) gives the 2.5-D integral equation for the porous medium Z ~ G ðn?  x? ; k z ; xÞ~tj ðn? ; k z ; xÞ dCðn? Þ ð1  c1 Þ~ U ui ðx? ; k z ; xÞ ¼ ij C Z  a  ½T~ Gij ðn?  x? ; k z ; xÞ~uj ðn? ; k z ; xÞ  T Gij ðn?  x? Þ~uj ðx? ; k z ; xÞ dCðn? Þ ð38Þ C ðaÞ cij ðx? Þ

in which the identity cij ðx? Þ þ ¼ dij has been applied. The desingularised characteristic of (38) can be  demonstrated by adding and subtracting T~ Gij (n?  x?, kz, x) u˜j (x?, kz, x) from the second term in (38), which generates Z ~ G ðn?  x? ; k z ; xÞ~tj ðn? ; k z ; xÞ dCðn? Þ ð1  c1 Þ~ U ui ðx? ; k z ; xÞ ¼ ij C Z   fT~ Gij ðn?  x? ; k z ; xÞ½~uj ðn? ; k z ; xÞ  ~uj ðx? ; k z ; xÞ dCðn? Þ ZC  a  ½T~ Gij ðn?  x? ; k z ; xÞ  T Gij ðn?  x? Þ~uj ðx? ; k z ; xÞ dCðn? Þ ð39Þ C

As it is assumed that ~ uj ðn? ; k z ; xÞ can be expanded into a Taylor series in the neighborhood of the point x? for  the 2.5-D problem, consequently, the singularity of T~ Gij (n?  x?, kz, x) · [u˜j (n?, kz, x)  u˜j (x?, kz, x)] is a  O(r0), which is integrable. Moreover, since the order of the singularity of T~ Gij (n?  x?, kz, x) and T Gij (n?  x?) is equal, the Cauchy singularity has been removed from the third term of (39). Therefore, (38) or (39) only contain integrable singular terms which facilitates numerical integration. 4.2. Treatment of the weakly singular kernel in the boundary integral equation ~ G (n?  x?, kz, x) is weakly singular, having an asymptotic expression equivalent to the It is noted that U ij ~ G (n?  x?, kz, x) is represented by logarithmic function when n? ! x?. However, as the Green’s function U ij the Hankel function, it is difficult to obtain the logarithmic function explicitly. Consequently, to improve the integration efficiency, Eq. (38) is rewritten in the following form Z ~ G ðn?  x? ; k z ; xÞ  U Ga ðn?  x? Þ~tj ðn? ; k z ; xÞ dCðn? Þ ui ðx? ; k z ; xÞ ¼ ½U ð1  c1 Þ~ ij ij C Z Z a  þ U Gij ðn?  x? Þ~tj ðn? ; k z ; xÞ dCðn? Þ  ½T~ Gij ðn?  x? ; k z ; xÞ~uj ðn? ; k z ; xÞ C

C a

 T Gij ðn?  x? Þ~ uj ðx? ; k z ; xÞ dCðn? Þ

ð40Þ

Thus eliminating the logarithmic singularity from the first integral in (40) allows the first integral to be calculated using the standard Gaussian quadrature formula, while the second integral can be evaluated by the a Gauss–Laguerre quadrature formula (Stroud and Secrest, 1966). It is noted that due to Eq. (12), U G44 in a Eq. (40) assumes the form U G44 ¼ M=ð2pb1 Þ logð1=rÞ. 5. Discretization of the 2.5-D boundary integral equation for the porous medium Suppose the boundary C in Eq. (20) is discretized by Ne 1-D iso-parametric elements and each element has M nodes, where the value of M may vary for different elements. For an arbitrary 1-D element, supposing that the point x belongs to the element and its intrinsic (local) coordinate is g, we have the following interpolation formulae xðgÞ ¼

M X l¼1

N l ðgÞxðlÞ ;

~ ui ðgÞ ¼

M X l¼1

ðlÞ

N l ðgÞ~ ui ;

~ti ðgÞ ¼

M X

ðlÞ

N l ðgÞ~ti

ð41Þ

l¼1

where x(l) is the lth node coordinate of the element in the global coordinate system, Nl(g) is the lth shape funcðlÞ ðlÞ tion, ~ ui and ~ti represent the values for the pseudo-displacement and the pseudo-traction at the lth node.

368

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

For a 1-D boundary element, each node can be shared by two elements or located within one element. Thus, each node is shared by r elements, where r = 2 for the end nodes and r is equal to unity for an interior node. As with conventional BEM, the nodes of the boundary elements are taken as the collocation points, thus, a collocation point xm is shared by elements with number e1 , . . ., er(r = 2 or r = 1). Supposing that for the elements sharing the collocation point, the collocation point xm is located at the nk th node of the element ek (k = 1  r), using (41), Eq. (40) is discretized as follows  X Z Z Ne X Ne X M  M X a ðklÞ ðklÞ G Ga ~ ~ ~ ui ðxm Þ ¼ ðU ij  U ij ÞN l ðgÞJ ðgÞdg þ U Gij N l ðgÞJ ðgÞdg tj tj ð1  c1 Þ~ k¼1 l¼1



Ne X

Sk

M  X



k¼e1 ;...;er l¼1

ðklÞ

~ uj

 T~ Gij N l ðgÞJ ðgÞdg

~ uj

k¼1ðk6¼e1 ;...;er Þ l¼1 M  X X

ðklÞ

Z



k¼1 l¼1

Sk

Sk

Z h

 i  a T~ Gij N l ðgÞ  dnk l T Gij J ðgÞdg þ ~uj ðxm Þ

Sk

Ne X k¼1ðk6¼e1 ;...;er Þ

Z

a

T Gij J ðgÞdg

Sk

ð42Þ where the superscripts k, l of u˜j, ~tj denote the element and the node number, respectively, and dnk l is the Kronecker d. It is worth noting that, to achieve high integration efficiency, the second integral in (42) should be evaluated by the Gauss–Laguerre quadrature formula (Stroud and Secrest, 1966), while the other integrals in (42) can be calculated by the conventional Gaussian quadrature formula. Following integration of (42) over all the elements for the collocation points, a set of linear equations for the unknown boundary variables can be constructed. Solution of the obtained linear equations yields the unknown boundary variables. Using the obtained boundary variables, the displacement at an interior point x is obtainable via discretisation of Eq. (21) and has the following discrete form  X  Z Z Ne X Ne X M  M  X ðklÞ ðklÞ G G ~ ~ ^ ^ ^uj ui ðxÞ ¼ U ij N l ðgÞJ ðgÞ dg  T ij N l ðgÞJ ðgÞ dg ð43Þ tj k¼1

l¼1

Sk

k¼1

l¼1

Sk

The stresses of the porous medium can be further obtained by substituting (21) into the constitutive relation (1). Once the wavenumber domain solution is determined, the space domain solution can be recovered by the inverse Fourier transform, which is accomplished using a FFT method. If the number of the wavenumber sample points is taken to be N, and sample spacing for the space and the wavenumber domain is Dz and Dkz, respectively, then the following relation holds Dk z ¼

2p N Dz

ð44Þ

The sample spacing in the space domain should satisfy Dz 6 kmin/2, where kmin is the shortest wavelength for the waves propagating within the calculation domain. 6. Numerical results In this section, two numerical examples are presented. In the first example, the dynamic response of an infinitely long circular tunnel embedded in a full-space porous medium and subjected to a harmonic Bell-shaped distributed load with a fixed location will be calculated by the proposed 2.5-D BEM. The obtained numerical results are compared with those obtained from an analytical approach (Lu and Jeng, 2006). In the second example, the response of an infinitely long semi-circular tunnel embedded in a half-space porous medium and subjected to a Bell-shaped moving load is calculated. Note that in the two examples of this section, Biot’s theory for low frequency waves (Biot, 1956) is used. 6.1. Validation of the proposed 2.5-D boundary element method Fig. 1 shows an infinitely long circular tunnel embedded in a full-space porous medium and subjected to a harmonic Bell-shaped distributed load with a fixed location. This problem can be solved in the wavenumber

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

369

x surrounding porous medium

R

eiwt

z

o ring harmonic bell-shaped load

Fig. 1. An infinitely long circular tunnel embedded in a full space porous medium and subjected to a harmonic Bell-shaped normal ring load.

domain using Eqs. (20) and (21), and the space domain solution can be retrieved by performing the FFT on the discrete wavenumber domain solutions. The radius of the tunnel is R = 3 m. The parameters for the surrounding porous medium are assumed as: k = 4.0 · 107 Pa, l = 2.0 · 107 Pa, qs = 2.5 · 103 kg/m3, qf = 1.0 · 103 kg/m3, a1 = 3, a = 0.95, M = 4.0 · 108 Pa, / = 0.3, g = 1.0 · 103 Pa s, and k = 1.0 · 1012 m2. The distributed Bell-shaped ring normal load is given by ^ðbÞ ðz; xÞ ¼ rn eðzzs Þ r

2

=a2l ixt

ð45Þ

e

where al = 6 m and zs = Dz(N  1)/2, where N is the number of sample points in the wavenumber domain and Dz is the sample spacing in the space domain. The tunnel surface is assumed to be completely permeable. Thus, in the space domain, the axisymetrical boundary condition for the tunnel boundary is ^rr ðR; z; xÞ ¼ r ^ðbÞ ðz; xÞ; r

^zr ðR; z; xÞ ¼ 0; r

^rh ðR; z; xÞ ¼ 0; r

^pðR; z; xÞ ¼ 0

ð46Þ

The discrete Fourier transform accomplished by FFT with respect to the z-coordinate is performed on Eq. (46) to yield the boundary condition for the 2.5-D boundary integral equation in the wavenumber domain. In this example, the number of sample points in the wavenumber domain is set as N = 501, while the sample spacing in the space domain is Dz = 5.0, 1.0 m, for frequencies of 10 and 100 1/s, respectively. The sample spacing in the wavenumber domain is given by (44). As the tunnel is symmetrical with respect to the x-axis and the y-axis, only 1/4 of the tunnel boundary needs to be discretized. When performing the numerical simulation, this 1/4 section of the tunnel boundary was uniformly covered by 24 two-node linear elements. The amplitudes for the radial and axial displacements at the tunnel surface obtained by the analytical methodology outlined in Lu and Jeng (2006) and the 2.5-D BEM approach are illustrated in Fig. 2. The amplitudes for the radial displacement, the axial displacement and the pore pressure at the interior points with the radial coordinate r = 3.5 m for the two cases are shown in Fig. 3. Note that when implementing the analytical approach, the space domain solution is also recovered by the FFT method from the closed-form wavenumber domain solution. Also, the number of the wavenumber domain sample points and the sample spacings in the space domain and the wavenumber domain are the same as those for the 2.5-D BEM methodology. It follows from Figs. 2 and 3 that the numerical results due to the 2.5-D BEM methodology agree very well with those from the analytical approach, which verifies the proposed 2.5-D BEM. 6.2. A half-circle tunnel embedded in a half-space porous medium and subjected to a moving load In this example, a semi-circular tunnel embedded in a half-space porous medium and subjected to a moving load (Fig. 4) is discussed. Eqs. (29) and (30) for the moving load problem are used to solve the problem in the wavenumber domain and the space domain solution is again recovered by the FFT methodology. The radius of the tunnel is R = 4 m. The depth of the tunnel is h = 8 m (Fig. 4). The parameters for the porous medium are assumed as the same values as in Section 6.1. Due to the symmetry of the structure about the y-axis, only half of the boundaries of the tunnel and the half-space need to be discretized. The discretisation scheme for the tunnel boundaries and the half-space is as follows: the bottom right side of the tunnel is covered uniformly by 8 two-node linear elements; the right semi-circular part of the tunnel boundary is divided

370

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

b

a

0.08

0.42

ω=10 1/s

0.07

urμ/(σnR)

0.36

ω=10 1/s

/( nR) uzμ/(σ

0.06 0.05

0.30

Analytical method The proposed 2.5-D BEM

0.04

0.24 Analytical method The proposed 2.5-D BEM

0.03

0.18

0.02

0.12

0.01

0.06

0.00

0.00

-0.01 0

1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500

250

500

750

1000 1250 1500 1750 2000 2250 2500

z (m)

z (m)

c

d

0.14

0.06

urμ/(σnR)

0.12

ω=100 1/s

ω=100 1/s uzμ/(σnR)

0.05

0.10 0.04

0.08 Analytical method The proposed 2.5-D BEM

0.06

Analytical method The proposed 2.5-D BEM

0.03

0.04

0.02

0.02

0.01

0.00 0.00

-0.02 0

50

100

150

200

250

300

350

400

450

500

0

z (m)

50

100

150

200

250

300

350

400

450

500

z (m)

Fig. 2. Comparison of the amplitudes for the radial displacement and the axial displacement at the tunnel surface determined by the analytical methodology of Lu and Jeng (2006) with those obtained by the proposed 2.5-D BEM: (a) the scaled radial displacement uˆr(R, z, x)l/(Rrn) for the case of x = 10 1/s; (b) the scaled axial displacement uˆz(R, z, x)l/(Rrn) for the case of x = 10 1/s; (c) the scaled radial displacement uˆr(R, z, x)l/(Rrn) for the case of x = 100 1/s; (d) the scaled axial displacement uˆz(R, z, x)l/(Rrn) for the case of x = 100 1/s.

uniformly by 18 linear two-node elements; 25 two-node elements are used to discretize the right boundary of the half space. The moving load is distributed uniformly over four bottom elements along the x-direction of right side of the structure (Fig. 4). Along the z-axis, the distribution of the moving load is considered to be a Bell-shaped function. Thus, the boundary condition for the part of the tunnel bottom which is subjected to the moving load is ryy ðx; y; z  vtÞ ¼ rn e½ðzvtÞzs 

2

=a2l

;

rzy ðx; y; z  vtÞ ¼ 0;

0 6 x 6 2:0 m; y ¼ 8 m

ryx ðx; y; z  vtÞ ¼ 0; ð47Þ

in which v is the velocity of the moving load and t is the time; al = 6 m and zs = Dz(N  1)/2. The number of sample points in the wavenumber domain is set as 801. The sample spacing in the space domain is equal to Dz = 2.0 m, and the sample spacing in the wavenumber domain is given by (44). Moving load velocities of 50, 100 and 150 m/s are considered. Note that the rest of the tunnel surface and the half space boundary are assumed to be stress free. Also, the boundaries of the half space and the tunnel are assumed to be permeable, so the pore pressure vanishes along the boundaries.

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

a

371

b

0.42

0.07

urμ/(σnR)

0.35

uzμ/(σnR)

0.06

ω=10 1/s

ω=10 1/s

0.05

0.28

0.04

0.21

Analytical method The proposed 2.5-D BEM

0.03

Analytical method The proposed 2.5-D BEM

0.14 0.02 0.07 0.01 0.00

0.00

1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500

0

250

500

750

1000 1250 1500 1750 2000 2250 2500

z (m)

z (m)

c

d

0.16

0.12 0.14

p/σn

0.12

urμ/(σnR)

ω=10 1/s

0.10

ω=100 1/s

0.10

0.08

0.08 0.06

Analytical method The proposed 2.5-D BEM

0.06

Analytical method The proposed 2.5-D BEM

0.04

0.04 0.02

0.02

0.00

0.00

-0.02 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500

0

50

100

150

200

250

300

350

400

450

500

z (m)

z (m)

f

e 0.07

1.05

uzμ/(σnR)

0.06 0.05

p/σn

0.90

ω=100 1/s

0.75

0.04

Analytical method The proposed 2.5-D BEM

ω=100 1/s

0.60

0.03

0.45

Analytical method The proposed 2.5-D BEM

0.02 0.30 0.01 0.15 0.00 0.00 -0.01 0

50

100

150

200

250

z (m)

300

350

400

450

500

0

50

100

150

200

250

300

350

400

450

500

z (m)

Fig. 3. Comparison of the amplitudes for the radial displacement, the axial displacement and the pore pressure at the interior points with the radial coordinate r = 3.5 m obtained by the analytical methodology of Lu and Jeng (2006) with those obtained by the proposed 2.5-D BEM: (a) the scaled radial displacement uˆr(R, z, x)l/(Rrn) for the case of x = 10 1/s; (b) the scaled axial displacement uˆz(R, z, x)l/(Rrn) for the case of x = 10 1/s; (c) the scaled pore pressure ^pðR; z; xÞ=rn for the case of x = 10 1/s; (d) the scaled radial displacement uˆr(R, z, x)l/(Rrn) for the case of x = 100 1/s; (e) the scaled axial displacement uˆz(R, z, x)l/(Rrn) for the case of x = 100 1/s; (f) the scaled pore pressure ^ pðR; z; xÞ=rn for the case of x = 100 1/s.

372

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

y free surface

x

P5 o

P6

P8

P7 half-space porous medium

h

P4

P3 moving load

semi-circle tunnel

R

P2

P1

Fig. 4. A semi-circular tunnel embedded in a half-space porous medium and subjected to a moving load.

a

b

0.10 0.05

0.02

uyμ/(σnR)

0.01

uzμ/(σnR)

0.00 0.00 -0.05

Point 1

-0.01

Point 1

-0.10 -0.02 -0.15

v=50 m/s v=100 m/s v=150 m/s

-0.20

-0.03 -0.04

-0.25 -0.30 -300 -240

v=50 m/s v=100 m/s v=150 m/s

-0.05 -180 -120

-60

0

z (m)

60

120

180

240

300

-300 -240 -180 -120

-60

0

60

120

180

240

300

z' (m)

Fig. 5. Numerical results at Point 1 for the semi-circular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR).

Numerical results at eight points are given (Fig. 5). Point 1 and Point 2 are located at the bottom of the tunnel; Point 3 and Point 4 are located at the semi-circular tunnel boundary; two points (Point 5 and Point 6) are on the surface of the half space; the other two points (Point 7 and Point 8) are in the calculation domain (Fig. 5). The coordinates for the Points 1–8 are as follows: (x1 = 0 m, y1 = 8 m), (x2 = 2.5 m, y2 = 8 m), (x3 = 0 m, y3 = 4 m), (x4 = 2.83 m, y4 = 5.17 m), (x5 = 0 m, y5 = 0 m), (x6 = 10.36 m, y6 = 0 m), (x7 = 0 m, y7 = 2 m), (x8 = 6 m, y8 = 2 m). Figs. 5–10 illustrate the results for the displacements uy, uz for Points 1–6, respectively. Figs. 11 and 12 plot the numerical results for the displacements uy, uz and the pore pressure p at Point 7 and Point 8, respectively. To have a clear image for the domain with the most significant response, only the variables within the range 300 6 z 0 6 300 m are plotted in the figures, where z 0 = z  vt  zs. It follows from Figs. 5–12 that for the lower velocity case (v = 50 m/s), the response at all the points is almost symmetrical with respect to the axis of the moving load. However, for the cases of v = 100 m/s and v = 150 m/s, the symmetry is lost due to the occurrence of the Mach cone corresponding to the shear wave of the porous medium. Note that the shear wave velocity of the porous medium with the above parameters for an angular frequency of 100 1/s is approximately 98 m/s; since the moving load velocity of v = 100 m/s is very close to the shear wave velocity of the porous medium, the response due to the load with v = 100 m/s is the largest of the three cases. Furthermore, the response of the saturated porous medium for this case appears to be oscillatory in nature, which increases the risk of structural damage to the tunnel or any overlying structure. Figs. 5–12 also show that with increasing velocity, the increase in the axial displacement uz is much greater than that of the vertical displacement uy.

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

0.05

0.02

uyμ/(σnR)

373

uzμ/(σnR)

0.01

0.00

0.00 -0.05

Point 2 -0.01

-0.10

v=50 m/s v=100 m/s v=150 m/s

Point 2

-0.02

-0.15

-0.03

-0.20

-0.04

v=50 m/s v=100 m/s v=150 m/s

-0.05 -0.25 -300 -240

-180 -120

-60

0

60

120

180

240

300

-300 -240 -180 -120

-60

0

60

120

180

240

300

z' (m)

z' (m)

Fig. 6. Numerical results at Point 2 for the semi-circlular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR).

0.05

uyμ/(σnR)

0.030

0.00

0.015

-0.05 -0.10

uzμ/(σnR)

0.000 v=50 m/s v=100 m/s v=150 m/s

Point 3

-0.015

Point 3

-0.15

-0.030

-0.20

-0.045

-0.25

-0.060

-0.30

-0.075

-300 -240

-180

-120

-60

0

60

120

180

240

300

-300 -240

v=50 m/s v=100 m/s v=150 m/s

-180 -120

-60

0

60

120

180

240

300

z' (m)

z ' (m)

Fig. 7. Numerical results at Point 3 for the semi-circlular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR).

0.08 0.04

uyμ/(σnR)

0.030

uzμ/(σnR)

0.015

0.00 -0.04

0.000

Point 4 -0.08

-0.015 v=50 m/s v=100 m/s v=150 m/s

-0.12

v=50 m/s v=100 m/s v=150 m/s

Point 4 -0.030

-0.16 -0.045 -0.20 -0.060

-0.24

-0.075

-0.28 -300 -240

-180

-120

-60

0 '

z (m)

60

120

180

240

300

-300 -240 -180

-120

-60

0

60

120

180

240

300

'

z (m)

Fig. 8. Numerical results at Point 4 for the semi-circular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR).

374

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

0.08

uyμ/(σnR)

0.08

0.00

uzμ/(σnR)

0.04

-0.08 -0.16

0.00 -0.04

Point 5

v=50 m/s v=100 m/s v=150 m/s

-0.24

-0.12

-0.40

-0.16

-180

-120

-60

0

60

120

180

v=50 m/s v=100 m/s v=150 m/s

-0.08

-0.32

-300 -240

Point 5

240

300

-0.20 -300 -240

-180

-120

-60

0

60

120

180

240

300

z' (m)

z' (m)

Fig. 9. Numerical results at Point 5 for the semi-circular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR).

0.10 0.05

0.08

uyμ/(σnR)

uzμ/(σnR)

0.04

0.00 0.00

-0.05 -0.10

Point 6

-0.04

v=50 m/s v=100 m/s v=150 m/s

-0.15

-0.08

-0.20

-0.12

Point 6

v=50 m/s v=100 m/s v=150 m/s

-0.25 -0.16

-0.30 -0.20

-300 -240

-180

-120 -60

0

z' (m)

60

120

180

240

300

-300

-240

-180

-120

-60

0

60

120

180

240

300

z ' (m)

Fig. 10. Numerical results at Point 6 for the semi-circular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR).

Figs. 5 and 6 show that from Point 1 to Point 2, the variation of the axial displacement uz is not significant, and that the vertical displacement due to moving loads of v = 50 m/s and v = 150 m/s decreases more quickly than that when v = 100 m/s. Comparison of Figs. 7 and 8 with Figs. 5 and 6 suggests that for the cases when v = 50 m/s and v = 150 m/s, the vertical displacements at Point 3 and Point 4 are smaller than those at Point 1 and Point 2; however, the vertical displacement and the axial displacement uz for the case v = 100 m/s increases slightly when the observation point moves from Points 1–2 to Points 3–4. It follows from Fig. 9 that compared with Points 1–4, both the vertical displacement and the axial displacement due to v = 100 m/s at Point 5 increase dramatically. Fig. 10 shows that vertical displacement at Point 6 is less significant than at Point 5, conversely with increasing horizontal distance along the free surface, the variation in axial displacement is minimal. Fig. 11 indicates that the vertical displacement for the case when v = 100 m/s at Point 7 is the largest among of all the points studied. However, the axial displacement at Point 7 due to v = 100 m/s is smaller than that at Point 5. Fig. 11 also shows that the pore pressure for the case when v = 50 m/s is always positive, yet for the cases v = 100 m/s and v = 150 m/s, negative pore pressure occurs. Moreover, for the case of v = 100 m/s, the pore pressure is oscillatory, which will increase the potential for liquefaction of the surrounding saturated porous medium. Fig. 12 demonstrates that the magnitudes of the negative pore pressure and the axial displacement for the case when v = 100 m/s at Point 8 are larger than those at Point 7.

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

0.06

uyμ/(σnR)

0.030

0.00

375

uzμ/(σnR)

0.015

-0.06 0.000 -0.12

Point 7

-0.015

-0.18

v=50 m/s v=100 m/s v=150 m/s

-0.24

-0.030

-0.30

-0.045

-0.36

-0.060

-0.42

-0.075

-0.48 -300 -240 -180 -120

-60

0 ,

60

120

180

240

Point 7

v=50 m/s v=100 m/s v=150 m/s

-0.090 -300 -240 -180 -120

300

z

-60

0

,

60

120

180

240

300

z

0.12 0.10

p/μ

0.08 0.06

Point 7

v=50 m/s v=100 m/s v=150 m/s

0.04 0.02 0.00 -0.02 -0.04 -0.06 -300

-240 -180 -120

-60

0

,

60

120

180

240

300

z

Fig. 11. Numerical results at Point 7 for the semi-circular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR); (c) the nondimensional pore pressure p/rn.

7. Conclusions The 2.5-D boundary integral equation for a porous medium in the frequency-wavenumber domain has been derived systematically by means of a 3-D boundary integral equation for a saturated porous medium. Specific 2.5-D boundary integral equations for the wave scattering problem and the moving load problem, which are important in earthquake and civil engineering design, have also been presented. The desingularisation of the obtained wavenumber domain 2.5-D boundary integral equation is realised by introducing an auxiliary problem composed of a plane strain problem, an antiplane elastic problem and a 2-D harmonic potential problem. Thus, the resulting desingularised 2.5-D wavenumber domain boundary integral equation only involves integrable singularity, which facilitates the numerical solution of the 2.5-D boundary integral equation. Discretisation of the wavenumber domain 2.5-D boundary integral equation is achieved by introducing isoparametric boundary elements. Using the discrete 2.5-D wavenumber domain BEM formulation, yields the wavenumber domain solution and the space domain solution can be further obtained by synthesizing the wavenumber solution via FFT. The new 2.5-D BEM proposed in this study is validated by comparing numerical results for a problem with those obtained via an analytical method. To demonstrate the new method, a numerical example for the moving load problem is given. In summary, the 2.5-D BEM for porous media has been systematically developed to facilitate the analysis of various dynamic problems encountered in earthquake engineering, civil engineering and geophysics. Consideration of the analogy between poroelasticity and thermoelasticity will allow this new methodology to be further applied in the solution of similar dynamic thermoelastic problems.

376

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

0.16

0.030 u y μ/( σnR)

uz μ/( σ R) n

0.015

0.08

0.000

0.00

-0.015 -0.08

Point 8 -0.16

v=50m/s v=100 m/s v=150 m/s

-0.24

v=50m/s v=100 m/s v=150 m/s

-0.045 -0.060

-0.32

-0.075

-0.40

-0.090

-300 -240 -180 -120

Point 8

-0.030

-60

0 , z

60

120

180

240

-300 -240

300

-180 -120

-60

0 , z

60

120

180

240

300

0.125 0.100

p/μ

0.075

Point 8

v=50m/s v=100 m/s v=150 m/s

0.050 0.025 0.000 -0.025 -0.050 -0.075 -300 -240

-180 -120

-60

0 , z

60

120

180

240

300

Fig. 12. Numerical results at the interior Point 8 for the semi-circular tunnel subjected to a moving load with a velocity of 50 m/s, 100 m/s and 150 m/s: (a) the nondimensional vertical displacement luy/(rnR); (b) the nondimensional axial displacement luz/(rnR); (c) the nondimensional pore pressure p/rn.

Acknowledgments This project is supported by the National Natural Science Foundation of China with grant number No. 50578071 and ARC Discovery Grant DP#0450906 (2004-2007). The first author would also like to acknowledge the financial support received from the Chinese Education Department. References Banerjee, P.K., 1981. Boundary Element Methods in Engineering Science. McGraw-Hill-Book Co.(UK), New York. Beskos, D.E., 1987. Boundary element methods in dynamic analysis. Appl. Mech. Rev. 40, 1–23. Beskos, D.E., 1997. Boundary element methods in dynamic analysis: Part II 1986–1996. Appl. Mech. Rev. 50, 149–197. Biot, M.A., 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid, I, low frenquency range. J. Acoust. Soc. Am. 28, 168–178. Biot, M.A., 1962. Mechanics of deformation and acoustic propagation in porous media. J. Appl. Phys. 33, 1482–1498. Bonnet, G., 1987. Basic singular solutions for a poroelastic medium in the dynamic range. J. Acoust. Soc. Am. 82, 1758–1762. Bouchon, M., 2003. A review of the discrete wavenumber method. Pure Appl. Geophys. 160, 445–465. Burke, M., Kingsbury, H.B., 1984. Response of poroelastic layers to moving loads. Int. J. Solids Struct. 20, 499–511. Cheng, A.H.D., Badmus, D.E., Beskos, D.E., 1991. Integral equations for dynamic poroelasticity in frequency domain with BEM solution. J. Eng. Mech., ASCE 117, 1136–1157. Dominguez, J., 1991. An integral formulation for dynamic poroelasticity. J. Appl. Mech., ASME 58, 588–591. Dominguez, J., 1992. Boundary element approach for dynamic poroelastic problems. Int. J. Numer. Methods Eng. 35, 307–324.

J.-F. Lu et al. / International Journal of Solids and Structures 45 (2008) 359–377

377

Fryba, L., 1999. Vibration of Solids and Structures under Moving Loads. Thomas Telford, London. Gaffet, S., Bouchon, M., 1989. Effects of two-dimensional topographies using the discrete wavenumber-boundary integral-equation method in P-SV cases. J. Acoust. Soc. Am. 85, 2277–2283. Hall, W.S., 1989. Integration methods for singular boundary element integrands. In: Brebbia, C.A. (Ed.), Boundary Elements, X, vol. 1. Computational Mechanical Publications, Southampton, pp. 219–236. Jin, B., Yue, Z.Q., Tham, L.G., 2004. Stresses and excess pore pressure induced in saturated poroelastic halfspace by moving line load. Soil Dyn. Earthq. Eng. 24, 25–33. Kattis, S.E., Beskos, D.E., Cheng, A.H.D., 2003. 2D dynamic response of unlined and lined tunnels in poroelastic soil to harmonic body waves. Earthq. Eng. Struct. Dyn. 32, 97–110. Kobayashi, S., 1987. Boundary Element Methods in Mechanics. Elsevier Science Publisher B.V.. Lu, J.F., Jeng, D.S., 2006. Dynamic analysis of an infinite cylindrical hole in a saturated poroelastic medium. Arch. Appl. Mech. 76, 263– 276. Lu, J.F., Jeng, D.S., Williams, S., submitted for publication. A 2.5-D dynamic model for a saturated porous medium: Part I. Green’s function. Int. J. Solids Struct, doi:10.1016/j.ijsolstr.2007.07.025. Luco, J.E., Wong, H.L., De Barros, F.C.P., 1990. Three-dimensional response of a cylindrical canyon in a layered half-space. Earthq. Eng. Struct. Dyn. 19, 799–817. Maeso, O., Aznarez, J.J., Garcia, F., 2005. Dynamic impedances of piles and groups of piles in saturated soils. Comput. Struct. 83, 769– 782. Papageorgiou, A.S., Pei, D., 1998. A discrete wavenumber boundary element method for study of the 3-D response of 2-d scatterers. Earthq. Eng. Struct. Dyn. 27, 619–638. Pedersen, H., LeBrun, B., Hatzfeld, D., Campillo, M., Bard, P.Y., 1994. Ground-motion amplitude across ridges. Bull. Seismol. Soc. Am. 84, 1786–1800. Schanz, M., 2001. Application of 3D time domain boundary element formulation to wave propagation in poroelastic solids. Eng. Anal. Bound. Elem. 25, 363–376. Sheng, X., Jones, C.J.C., Thompson, D.J., 2005. Modelling ground vibration from railways using wavenumber finite- and boundaryelement methods. Proc. R. Soc. A 461, 2043–2070. Sneddon, I.N., 1951. Fourier Transforms. McGraw-Hill, New York, NY, USA. Stamos, A.A., Beskos, D.E., 1996. 3-D seismic response analysis of long lined tunnels in half-space. Soil Dyn. Earthq. Eng. 15, 111–118. Stroud, A.H., Secrest, D., 1966. Gaussian Quadrature Formulas. Prentice-Hall, Englewood Cliffs, NJ. Tadeu, A.J.B., Kausel, E., 2000. Green’s functions for two-and-a-half-dimensional elastodynamic problems. J. Eng. Mech., ASCE 126, 1093–1097. Zhang, L.P., Chopra, A.K., 1991. 3-Dimensional analysis of spatially varying ground motions around a uniform canyon in a homogeneous half-space. Earthq. Eng. Struct. Dyn. 20, 911–926. Zimmerman, C., Stern, M., 1993. Boundary element solution of 3-D wave scatter problems in a poroelastic medium. Eng. Anal. Bound. Elem. 12, 223–240.