Dynamic Reissner–Sagoci problem for a transversely isotropic half-space containing finite length cylindrical cavity

Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

Contents lists available at ScienceDirect

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

Dynamic Reissner–Sagoci problem for a transversely isotropic half-space containing finite length cylindrical cavity Azizollah Ardeshir-Behrestaghi a, Morteza Eskandari-Ghadi b,n, Bahram Navayi neya a, Javad Vaseghi-Amiri a a b

Faculty of Civil Engineering, Babol Noshirvani University of Technology, P.O. Box 47148‐71168, Babol, Iran School of Civil Engineering, College of Engineering, University of Tehran, P.O.Box 11165-4563, Tehran, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 5 May 2014 Received in revised form 28 June 2014 Accepted 15 July 2014

A transversely isotropic linear elastic half-space containing a circular cylindrical cavity of finite length with a depth-wise axis of material symmetry is considered to be under the effect of a mono-harmonic torsional motion applied on a rigid circular disc with the same radius of the cavity and welded at the bottom of the cavity. With the aid of Fourier sine and cosine integral transforms, the mixed boundary value problem is reduced to a generalized Cauchy singular integral equation for the unknown shear stress. The Cauchy integral equation involved in this paper is analytically investigated such that the solution is written in the form of a known singular function multiplied by an unknown regular function. The regular part of the shear stress is numerically determined with the use of Gauss–Jacobi integration formula. Series representation of the stress and displacement are obtained, and it is shown that their degenerated form to the static problem of isotropic material is coincide with the existing solutions in the literature. To investigate the effects of material anisotropy and the length of cavity, the tangential displacement and the shear stress in between the rigid disc and the bottom of cavity are numerically evaluated and illustrated, where some differences are distinguished. With the differences illustrated in this paper for different length of cavity, it is recognized that the effect of length of cavity cannot be neglected in analysis and design. Different results for different degrees of anisotropy shows that the anisotropy of the material is a normal behavior, which should be considered in this kind of medium. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Reissner–Sagosi Torsion Rigid disc Cylindrical cavity Transversely isotropic Half space Cauchy integral equation

1. Introduction The problem of determination of stresses and displacements in an elastic medium due to torsional static and dynamic rotation of a rigid disc in welded contact has been an interesting subject in both theoretical and practical mechanics. In this topic, evaluation of the stresses and displacements of an elastic half-space with a finite open cylindrical cavity is a considerable interest in both mathematical and engineering points of view. In the context of engineering, a study of this kind is relevant to foundation drilling, structural and mechanical designs, and borehole technology [10]. The first investigation in this class of problems is due to Westergaard [17], who studied an infinite isotropic solid containing an infinite cylindrical cavity with hydrostatic pressure acting on the wall of the cavity and presented some approximate results. A rigorous investigation of this problem was presented by Tranter [16] and Jordan [5], then treated the dynamic problem of a suddenly applied pressure over a finite interval of the infinite

n

Corresponding author. Fax: þ98 21 88632423. E-mail address: [email protected] (M. Eskandari-Ghadi).

http://dx.doi.org/10.1016/j.soildyn.2014.07.009 0267-7261/& 2014 Elsevier Ltd. All rights reserved.

cavity. Because of the complexities encountered in the problem, numerical results were presented only at large distance away from the location of pressure. As a fundamental example of mixed boundary value problems in the theory of elastic wave propagation, the problem of forced torsional rotation of an elastic halfspace has been first attracted by Reissner and Sagoci [14] and Sogoci [15]. That is why this class of problems is often referred to as Reissner–Sagoci problem. Some researchers, after Reissner and Sogoci, have investigated either the static or dynamic interaction of rigid disc and elastic isotropic half-spaces, among which Collins [2]; Williams [19]; and Pak and Sophores [9] are mentioned. The attempt for the simulation of the mechanical features of different realistic materials in the analyses, however, requires expressions for anisotropic media. The form of anisotropy with the most common application is perhaps the case of transverse isotropy. Rajapakse and Gross [13], for example, have been interested in the influence of transverse isotropy on the axisymmetric response of a borehole. Rahimian et al. [12] have investigated the Reissner–Sagoci problem for transversely isotropic half-space. The solutions of the generalized problem associated with a finite cylindrical cavity in a half-space would be of even greater both engineering and mathematical interest and challenge. It has

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

Nomenclature The following symbols are used in this paper: a Cs H ðpÞ n Im i Jm Km l ^l P nðγ ;βÞ r

radius of cylindrical cavity and rigid disc shear wave velocity in the isotropic plane Hankel function of first and second kinds (p ¼ 1; 2) of order n modified Bessel function of the first kind and mth pffiffiffiffiffiffiffiffiorder 1 Bessel function of the first kind and mth order modified Bessel function of the second kind and mth order depth of cylindrical cavity normalized depth of cylindrical cavity Jacobi polynomials radial coordinate

been found that the additional stiffness of the medium below the bottom of the hole can lead to a noticeable change of the response in the upper region. Pak and Abedzadeh [10] investigated rigorously the problem of torsional shear static traction acting on an open finite cylindrical cavity in an isotropic half-space in detail, and found the corresponding fundamental solution. They also extracted mathematically the resulting load-induced as well as shape-induced singularities in the response. Eskandari-Ghadi et al. [4] have extended this problem for both dynamic case and transversely isotropic material, where they have shown that the anisotropy cannot be neglected. Pak and Abedzadeh [11] has been interested in the static torsion of a rigid disc bonded to the bottom of a finite length cylindrical cavity exist in an elastic half-space, where they have used a combination of Fourier sine and cosine integral transforms to reduce the related mixed boundary value problem to a pair of integral equations, one of which possesses a generalized Cauchy singular kernel. They have applied the theory of analytic functions and Gauss–Jacobi integration formula to evaluate the solution of the mixed boundary value problem. This paper is concerned with the elastostatic and elastodynamic forced torsion of a transversely isotropic linear elastic halfspace containing a finite open circular cylindrical cavity. To attack the mixed boundary value problem, a transversely isotropic halfspace containing a circular cylindrical cavity of finite length is considered as the domain of the problem in such a way that the material symmetry of the half-space is assumed to be both depthwise and parallel to the axis of cylindrical cavity. A rigid disc welded on the bottom surface of the cavity is considered to be oscillatory moved with a mono-harmonic torsional motion. As Pak and Abedzadeh [11] did, both Fourier sine and cosine integral transforms are used to reduce the in hand mixed boundary value problem to a generalized Cauchy singular integral equation for the unknown shear stress. The Cauchy integral equation involved in this paper is analytically investigated and transformed to an equation, which is numerically solved, after which the shear stress, circumferential displacement, and the impedance function are determined in a straight manner. Excellent satisfaction of the boundary conditions, and the excellent agreement between the impedance function determined from this paper for the simpler case of surface rigid disc attached on an isotropic half-space and the existing results prove both the validity and accuracy of the solution reported in this paper. It is shown that neglecting either the degree of anisotropy or the length of cavity results in some wrong results, which means that none of these two parameters can be neglected.

253

r^ t uðr; z; tÞ Ym z z^

normalized radial coordinate time variable displacement component in θ-direction Bessel function of the second kind and mth order vertical coordinate normalized vertical coordinate μ shear modulus in the plane normal to the axis of symmetry μ0 shear modulus in planes normal to the plane of transverse isotropy Θ rotation angle of the rigid disc about the z-axis θ angular coordinate ρ material density τij ði; j ¼ r; θ; zÞ shear stress tensor ω angular frequency ω0 non-dimensional frequency ξ Fourier parameter

2. Formulation of the problem A transversely isotropic homogeneous linear elastic half-space containing circular cylindrical cavity with radius a 40 and depth lZ 0 is considered in a cylindrical coordinate system ðr; θ; zÞ, with a depth-wise z-axis, in such a way that the material axis of symmetry of the medium is parallel to both the z-axis and the axis of cylindrical cavity. A rigid disc of radius a bounded on the medium at the bottom of the flat-ended cavity, as depicted in Fig. 1, is considered to be affected by a time-harmonic torsional rotation, Θ e iωt , with Θ and ω, respectively, being the angle and circular frequency of the motion. Because of axial symmetry of the boundary value problem, the displacement vector has only one non-vanishing component, i.e. uθ ¼ uðr; z; tÞ. Following Pak and Abedzadeh [10,11] (see also Eskandari-Ghadi et al. [4]), it is convenient to define two different regions as indicated in Fig. 1 and find the response of each region with satisfying the boundary and continuity conditions. These two regions are defined as Region 1 ¼ fðr; θ; zÞjr 4 a; 0 o θ r 2π ; z 4 0g;

ð1Þ

Region 2 ¼ fðr; θ; zÞj0 r r o a; 0 o θ r 2π ; z 4 lg:

ð2Þ

The second region defines an open cylindrical region under the cavity and the first region is the remaining part of the whole halfspace, which is also an open region. In the absence of the body

Fig. 1. A rigid disc on an indented transversely isotropic half-space.

254

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

force, the time-harmonic equations of motion in terms of the nonzero displacement, in each region, are written in the form of ∂2 u1 1 ∂u1 u1 ∂ 2 u1 ω2  2 þ α 2 2 ¼  2 u1 ; þ 2 r ∂r r ∂r ∂z Cs

r 4 a;

z40

ð3Þ

∂2 u2 1 ∂u2 u2 ∂ 2 u2 ω2  þ α 2 2 ¼  2 u2 ; þ ∂r 2 r ∂r r 2 ∂z Cs

r o a;

z4l

ð4Þ

where the indices 1 and 2 show the Region 1 and 2, respectively. In ffi pffiffiffiffiffiffiffiffi addition, C s ¼ μ=ρ is the shear wave velocity in any plane perpendicular to the z-axis, commonly called as isotropic plane, α2 ¼ μ0 =μ, ρ the material density, μ the shear modulus of the material in the isotropic plane, and μ0 is the same function in any plane parallel to the material axis of symmetry. α2 is used to show the degree of anisotropy. It equals unity for isotropic material. The relevant stress-displacement relations in a transversely isotropic material for non-zero stresses are [6]   ∂u u ∂ u  ¼μr ; ð5Þ τrθ ¼ μ ∂r r ∂r r

τzθ ¼ μ0

  ∂u ∂z

ð6Þ

The stresses and displacement boundary conditions for the problem may be written as

τrθ1 ða; z; ωÞ ¼ 0;

0 oz o l

ð7Þ

τzθ1 ðr; 0; ωÞ ¼ 0;

r 4a

ð8Þ

u2 ðr; l; ωÞ ¼ Θ r;

r oa

ð9Þ

∂u1 ðr; 0; ωÞ ¼ 0; ∂z

r 4a

ð10Þ

r

∂ u1  ¼0 ∂r r

r ¼ a;

0 o z ol

ð11Þ

where Θ r stands for the angular displacement field of the solid. Moreover, the radiation condition for both Region 1 and 2 that guarantees no incoming wave is written as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uq ðr; z; ωÞ-0; ðr 2 þz2 Þ-1; q ¼ 1; 2 ð12Þ Using T l ðωÞ to denote the applied torque required to sustain the rotation of the rigid disc about the z-axis, one may write the loading condition at the base of the cylindrical cavity as Z a T l ð ωÞ ¼  2 π r 2 τ zθ2 ðr; l; ωÞdr ð13Þ 0

The continuity conditions at the common boundary of Region 1 and Region 2 may be written as [10]

τrθ1 ða þ ; z; ωÞ ¼ τrθ2 ða  ; z; ωÞ; ∂u1 þ ∂u2  ða ; z; ωÞ ¼ ða ; z; ωÞ: ∂z ∂z

z4l

ð14Þ

z 4l

ð15Þ

where f ða 7 Þ ¼ lim7 f ðxÞ. Because of the boundary conditions, it is x-a convenient to use integral transforms for the solution of the boundary value problems involved here. In view of the boundary conditions (9) and (10), it is appropriate to define the depth-wise Fourier cosine and sine transforms respectively for Region 1 and 2 as [11] Z 2 1 f~ 1C ðr; ξ; ωÞ ¼ f 1 ðr; z; ωÞ cos ðξzÞ dz; for Region 1 ð16Þ

π

2 f~ 2S ðr; ξ; ωÞ ¼

π

0

Z

1 l

f 2 ðr; z  l; ωÞ sin ξðz lÞ dz:

for Region 2 ð17Þ

With these integral transforms, the partial differential Eqs. (3) and (4) lead to   2 d u~ 1C 1 du~ 1C 1 2 2  u~ 1C ¼ 0; þ α λ þ ð18Þ r dr r2 dr 2   2 d u~ 2S 1 du~ 2S 1 2 2  α2 λ þ 2 u~ 2S ¼  α2 ξΘr; þ 2 r dr π r dr

ð19Þ

2 2 2 2 where λ ðξ; ωÞ ¼ ξ  ξp , ξp ¼ ω2 =ðα2 C s 2 Þ, and u~ 1C ðr; ξ; ωÞ and u~ 2S ðr; ξ; ωÞ are respectively the Fourier cosine and sine transforms of u1 ðr; z; ωÞ and u2 ðr; z; ωÞ as defined in (16) and (17). The general solutions of Eqs. (18) and (19) are

u~ 1C ðr; ξ; ωÞ ¼ A1 ðξ; ωÞK 1 ðαλrÞ þ B1 ðξ; ωÞI 1 ðαλrÞ; u~ 2S ðr; ξ; ωÞ ¼ A2 ðξ; ωÞK 1 ðαλrÞ þB2 ðξ; ωÞI 1 ðαλrÞ þ

ð20Þ 2 α 2 Θr ξ ; 2

ð21Þ

πλ

where the last term in (21) is a particular solution for u~ 2S ðr; ξ; ωÞ, and I 1 ðxÞ and K 1 ðxÞ are the first-order modified Bessel functions of the first and second kinds, respectively. To be consistent with the radiation condition, one must also define a Riemann surface for λ so that it will appropriately be single-valued everywhere. This can be achieved by specifying a branch cut for λ on the complex ξ-plane with a branch point emanating from ξp ¼ ω=α C s such that the real part of λ is always non-negative. For imaginary values of λ, 2 which happens when ξ o ω2 =ðα2 C s 2 Þ, the modified Bessel functions I 1 ðαλrÞ and K 1 ðαλrÞ are transformed to the ordinary Bessel function J 1 ð  αλrÞ and the Hankel function  ð1=2Þπ H ð2Þ 1 ðαλrÞ, respectively. With the conditions specified here, B1 ðξ; ωÞ should be identically zero as I 1 ðαλrÞ is unbounded when r approaches infinity. This is also compatible with the radiation condition. Likewise, A2 ðξ; ωÞ should be zero for the displacement u2 ðr; z; ωÞ in the Region 2 to be bounded at r ¼ 0. The displacement u1 ðr; z; ωÞ can be obtained by means of the inversion theorem for the Fourier cosine transform, while u1 ðr; z; ωÞ is determined by implementing inverse of Fourier sine transform as follows Z 1 A1 ðξ; ωÞK 1 ðαλrÞ cos ðξzÞ dξ; r 4 a; z 4 0 ð22Þ u1 ðr; z; ωÞ ¼ 0

u2 ðr; z; ωÞ ¼

Z

1 0

B2 ðξ; ωÞI 1 ðαλrÞ sin ξðz  lÞdξ þ Θre  iξp ðz  lÞ :

z 4l; r o a

ð23Þ Based on the continuity conditions (14) and (15) on the cylindrical boundary r ¼ a, one may write the shear stress τrθ1 in Region 1 as ∂ u1  ¼ μχ ðz; ωÞ τrθ1 ðr ¼ a þ ; z; ωÞ ¼ limþ μr z40 ð24Þ r-a ∂r r where χ ðz; ωÞ ¼ 0 for 0 o z o l and χ ðz; ωÞ ¼ τðz; ωÞ for z Z l. In this expression, τðz; ωÞ is an unknown function to be determined from the solution. In terms of the function τðz; ωÞ, an equivalent statement of the loading condition in (13) can also be given by Z 1 T ðωÞ τðζ ; ωÞdζ ¼  l 2 ð25Þ 2πμa l On the other hand, τrθ2 in Region 2 at r ¼ a is written as ∂ u2  ¼ μ χ ðz; ωÞ; τrθ2 ðr ¼ a  ; z; ωÞ ¼ r-a lim μ r zZl ∂r r

ð26Þ

Substituting the displacements from Eqs. (22) and (23) into Eqs. (24) and (26) and implementing the inverse theorem for Fourier cosine and sine transforms, one may respectively find Z 1 2 A1 ðξ; ωÞ ¼ χ ðz; ωÞ cos ðξzÞ dz; ð27Þ π αλ K 2 ðαλaÞ 0

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

B 2 ð ξ ; ωÞ ¼

2 π αλ I2 ðαλaÞ

Z

1 l

χ ðz; ωÞ sin ξðz  lÞ dz;

ð28Þ

Z 0

where the following recurrence formulas have been used [1]  ∂  αλ ½K 1 ðαλrÞ ¼ K 0 ðαλrÞ þ K 2 ðαλrÞ ; ∂r 2 2 K ðαλrÞ; K 2 ðαλrÞ ¼ K 0 ðαλrÞ þ αλr 1  ∂ αλ ½I 1 ðαλrÞ ¼ I 0 ðαλrÞ þ I 2 ðαλrÞ ; ∂r 2 2 I 2 ðαλrÞ ¼ I 0 ðαλrÞ  I ðαλrÞ: αλr 1

απ Z

0

½φk3 ðr; ζ þ z; ωÞ þ φk3 ðr; ζ  z; ωÞχ ðζ ; ωÞ dζ ;

ð29Þ

ð30Þ

ð31Þ

for Region 1, and Z 1 1 ½φI1 ðr; ζ  z; ωÞ  φI1 ðr; ζ þ z  2l; ωÞχ ðζ ; ωÞ dζ u2 ðr; z; ωÞ ¼

απ l þ Θr e  iξp ðz  lÞ Z μ0 1 τzθ2 ðr; z; ωÞ ¼ ½φI2 ðr; ζ  z; ωÞ þ φI2 ðr; ζ þ z  2l; ωÞχ ðζ ; ωÞ dζ απ l  iμ Θ r ξp e  iξp ðz  lÞ ;

τrθ2 ðr; z; ωÞ ¼

μ απ

Z

1 l

½φI3 ðr; ζ  z; ωÞ  φI3 ðr; ζ þ z  2l; ωÞχ ðζ ; ωÞ dζ ;

for Region 2, where Z 1 K 1 ðαλrÞ cos ðξdÞ dξ; φk1 ðr; d; ωÞ ¼  λK 2 ðαλaÞ 0 Z 1 I 1 ðαλrÞ cos ðξdÞ dξ; φI1 ðr; d; ωÞ ¼ λI2 ðαλaÞ 0 Z 1 ξK 1 ðαλ rÞ sin ðξdÞ dξ; φk2 ðr; d; ωÞ ¼ λK 2 ðαλ aÞ 0 Z 1 ξI1 ðαλ rÞ φI2 ðr; d; ωÞ ¼ sin ðξdÞ dξ; λI2 ðαλ aÞ 0 Z 1 K 2 ðαλrÞ cos ðξdÞ dξ; φk3 ðr; d; ωÞ ¼ K 2 ðαλaÞ 0 Z 1 I 2 ðαλrÞ φI3 ðr; d; ωÞ ¼ cos ðξdÞ dξ: I 2 ðαλaÞ 0

ð32Þ

1



0

sin ðξdÞ dξ:



ξI1 ðαλ rÞ  ξða  rÞ e sin ðξdÞ dξ λI2 ðαλ aÞ

roa

ð36Þ

With the aid of an asymptotic analysis for the quotients K 1 ðαλrÞ=K 2 ðαλaÞ and I 1 ðαλrÞ=I 2 ðαλaÞ as r approaches a, it can be shown that [4]  Z 1 Z 1 ξK 1 ðαλ a þ Þ ξK 1 ðαλ a þ Þ 1 sin ðξdÞ dξ ¼  1 sin ðξdÞ dξ þ ; d λK 2 ðαλ aÞ λK 2 ðαλ aÞ 0 0 ð37Þ Z 0

1

ξI1 ðαλ a  Þ sin ðξdÞ dξ ¼ λI2 ðαλ aÞ

Z 0

1





ξI1 ðαλa  Þ 1 1 sin ðξdÞ dξ þ d λI2 ðαλaÞ ð38Þ

Using Eq. (38) in φI2 ða ; z  ζ ; ωÞ and φI2 ða ; z þ ζ 2l; ωÞ, and Eq. (37) in φk2 ða þ ; z  ζ ; ωÞ and φk2 ða þ ; z þ ζ ; ωÞ, one may express (34) in the form of

Z 1 2 1 1 þ  τðζ ; ωÞdζ ζ  z ζ þ z  2l ζ þz l



Z 1 Z 1 ξI1 ðαλa  Þ þ  1 ½ sin ξðζ zÞ þ sin ξðz þ ζ  2lÞdξ τðζ ; ωÞdζ λ I ð αλ aÞ 2 l 0 Z 1 Z 1 ξK 1 ðαλa þ Þ þ 1½ sin ξðζ  zÞ λK 2 ðαλaÞ l 0 





 sin ξðζ þzÞ dξ τ ðζ ; ωÞdζ ¼ f ðz; ωÞ;

zZl

ð39Þ

where f ðz; ωÞ ¼ i π α Θ aξp e  iξp ðz  lÞ :

ð40Þ

The Eq. (39) is a high-singular integral equation that needs some especial attention to be solved for τðζ ; ωÞ. For further analysis, it is convenient to introduce the dimensionless parameters ^l ¼ l ; a

ζ ζ^ ¼ :

ð41Þ

a

^l

ð33Þ

ð34Þ

Eq. (34), which is a generalized Cauchy singular integral equation (see [3, 8]) is an integral equation for determining τðζ ; ωÞ. For further analysis of this integral equation, it is useful to write the functions φI2 ðr; d; ωÞ and φK2 ðr; d; ωÞin (33) in the following forms [4]:  Z 1 Z 1 ξK 1 ðαλ rÞ ξK 1 ðαλ rÞ  ξðr  aÞ sin ðξdÞ dξ ¼ e sin ðξdÞ dξ λK 2 ðαλ aÞ λK 2 ðαλ aÞ 0 0 Z 1 þ e  ξðr  aÞ sin ðξdÞ dξ; r 4 a ð35Þ 0

e

 ξða  rÞ

Z

With the use of these parameters, the functions τðz; ωÞ and f ðz; ωÞ are denoted as τ^ ðz^ ; ωÞ and f^ ðz^ ; ωÞ, respectively, and the governing integral Eq. (39) can be written as # Z 1" 2 1 1 þ  τ^ ðζ^ ; ωÞdζ^ ^l ζ^  z^ ζ^ þ z^ 2^l ζ^ þ z^ Z 1 þ ½k1 ðζ^  z^ ; ωÞ þ k1 ðζ^ þ z^  2^l; ωÞ

l

zZl

1

z z^ ¼ ; a

By virtue of (31) and (32) the continuity condition (15) can be stated in the form of Z 1 ðφI2 ða  ; z  ζ ; ωÞ  φI2 ða  ; z þ ζ  2l; ωÞÞτðζ ; ωÞ dζ l Z 1 þ ðφk2 ða þ ; z  ζ ; ωÞ þ φk2 ða þ ; z þ ζ ; ωÞÞτðζ ; ωÞ dζ ¼  i π α Θ r ξp e  iξp ðz  lÞ :

Z

0

0

1

ξI 1 ðαλ rÞ sin ðξdÞ dξ ¼ λI2 ðαλ aÞ

þ

Substituting A1 ðξ; ωÞ and B2 ðξ; ωÞ from (27) and (28) into (22) and (23), and also in the stress–displacement relationships, the displacements and stresses are written as Z 1 1 u1 ðr; z; ωÞ ¼ ½φk1 ðr; ζ þ z; ωÞ þ φk1 ðr; ζ  z; ωÞχ ðζ ; ωÞ dζ ; απ 0 0 Z 1 μ τzθ1 ðr; z; ωÞ ¼ ½φk2 ðr; z  ζ ; ωÞ þ φk2 ðr; ζ þ z; ωÞχ ðζ ; ωÞ dζ ;

μ τrθ1 ðr; z; ωÞ ¼ απ

1

255

þ k2 ðζ^  z^ ; ωÞ  k2 ðζ^ þ z^ ; ωÞτ^ ðζ^ ; ωÞdζ^ ¼ f^ ðz^ ; ωÞ; where k1 ðd; ωÞ ¼ k2 ðd; ωÞ ¼

Z

1



0

Z 0

1



z^ Z ^l

ð42Þ



ξI1 ðαaλÞ 1 sin ðξdÞ dξ; λI2 ðαaλÞ

ð43Þ



ξK 1 ðαaλÞ  1 sin ðξdÞ dξ: λK 2 ðαaλÞ

ð44Þ

The unknown in the singular integral Eq. (42) is the shear stress on the surface of r ¼ a for z 4l. To find the solution for (42), it is important to note that τ^ ðz^ ; ωÞ is apt to be singular at z ¼ l owing to both the geometry of the domain and the rigidity of the disc attached on the domain at the bottom of the cavity.

3. Evaluation of Bessel's integrals Before proceeding to the analysis of (42), it is appropriate to seek first an effective method to evaluate the two kernel functions k1 ðd; ωÞ and k2 ðd; ωÞ. As they are odd functions of d, it suffices to consider these for d Z0. For an accurate evaluation of k1 ðd; ωÞ, it is

256

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

useful to define the integrand function as  8 ξI ðαiλ^ Þ > < iλ^ I1 ðαiλ^ Þ  1 sin ðξdÞ ξ o ξp  2  f ðd; ξ; ωÞ ¼ Þ > : ξλII1 ððαλ αλÞ  1 sin ðξdÞ ξ 4 ξp

ð45Þ

2

where λ^ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R ξp 2  ξ2 , α ¼ αa;and k1 ðd; ωÞ ¼

1 0

f ðd; ξ; ωÞ dξ. On

the other hand, with the use of the relation between Bessel function and modified Bessel function, I υ ðyÞ ¼ e  ðπ =2Þυi J υ ðyeðπ =2Þi Þ, one may find ! ! ξI 1 ðαiλ^ Þ  ξJ 1 ðαλ^ Þ  1 sin ðξdÞ ¼  1 sin ðξdÞ ð46Þ λ^ J ðαλ^ Þ iλ^ I ðαiλ^ Þ 2

2

This shows that the singular J 2 ðyÞ ¼ 0. If we show these then the singular points qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ξRpn ¼ ξ2p  ðy2n =α2 Þ; ξRpn 40.

points are the roots of λ^ ¼ 0 and roots by yn for n ¼ 1; 2; : : : ; 1, of the integrand are ξp and With the definition of ξ ¼ ω2 = 2 p

ðα C s Þ, these singular points may be written as ξpn ¼ α1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ða2 ω2 =C 2s Þ yn 2 , which shows that they exist only in dynamic 2

R

2

case and also they varies in terms of degree of anisotropy. To investigate the issue in detail, it is useful to define   ξI1 ðαλÞ 1 eiξd f 1 ðd; ξ; ωÞ ¼ ð47Þ λI2 ðαλÞ and f 2 ðd; ξ; ωÞ ¼





ξI1 ðαλÞ þ1 eiξd λI2 ðαλÞ

To evaluate the infinite series in (49), one may note that since the Bessel function J 2 ðyÞ is equivalent to cosine function for large y, the zeros of Bessel function J 2 ðyÞ behave asymptotically as   3 yn - n þ π ; n-1: ð50Þ 4 By virtue of (50), it can be shown that (49) can be reduced to R 1 π 1 π 1 e  iξp d þ ∑ e  idξpn þ ∑ e  dξpn α αn¼1 αn¼1 I 1 π 1 1 dZ0 þ ∑ ½e  dξpn e  dξpn   ; d αn¼1

k1 ðd; ωÞ ¼

2π

where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðn þ ð3=4ÞÞ2 π 2 1 2 ξpn ¼  ξp ; 2

ð51Þ

ξ1 pn 4 0

α

ð52Þ 1

For the static case, (52) can be reduced to ξpn ¼ ðn þð3=4ÞÞðπ =αÞ and (51) is identical to that given by Eskandari-Ghadi et al. [4]  # 1 1 ðn þ 3=4Þπ d e  7π d=4α  yn d  α α k1 ðdÞ ¼ 2π  þ π þπ ∑ e e ; α d 1 e  π d=α n¼1 1

"

α

ð53Þ For an accurate evaluation of k2 ðd; ωÞ, it is useful to define the integrand function as  8 ξK ðαiλ^ Þ > < iλ^ K1 ðαiλ^ Þ  1 sin ðξdÞ ξ o ξp  2  ð54Þ gðd; ξ; ωÞ ¼ Þ > : ξλKK 1 ððαλ αλÞ  1 sin ðξdÞ ξ 4 ξp 2

ð48Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 in the complex ξ-plane. With the definitions of λðξ; ωÞ ¼ ξ  ξp qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 and λ^ ¼ iλ ¼ ξp  ξ , it is clear that f 1 ðd; ξ; ωÞ based on (47) is changed to f 1 ðd; ξ; ωÞ ¼ ðð  ξJ 1 ðαλ^ Þ=λ^ J 2 ðαλ^ ÞÞ 1Þeiξd for ξ o ξp , however it is unchanged for ξ Z ξp . In the same way, f 2 ðd; ξ; ωÞ based on (48) is changed to f ðd; ξ; ωÞ ¼ ðð  ξJ ðαλ^ Þ=λ^ J ðαλ^ ÞÞ þ 1Þ 2

1

2

eiξd for ξ o ξp , and it is unchanged for ξ Z ξp . By virtue of the

relation between I v and J v it is evident that for ξ Z ξp both f 1 and f 2 have an infinite number of simple poles corresponding to the roots of J 2 ðyÞ. If we show these roots by yn for n ¼ 1; 2; : : : ; 1, then the singular points of the integrand function are collapsed on the imaginary axis in the complex ξ-plane in such a way that qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 I 2 I iξpn ¼ i yn2  ξp with ξpn 40. One can show with the aid of α

R 1 where k2 ðd; ωÞ ¼ 0 gðd; ξ; ωÞ dξ. λ, based on Eq. (54), is asymptotically equivalent to ξ, when ξ approaches infinity. In addition, lim K j ðyÞ-0 for j ¼ 1 and 2, however K j ðyÞ may be replaced by y-1  1=2 h i 2 2  12 12 Þð4j2  32 Þ π K j ðyÞ  2y e  y 1 þ 4j1!ð8yÞ þ ð4j 2!ð8yÞ þ :::: when y approa2 Þ ches infinity [18]. So that lim ξλKK 1 ððαλ αλÞ-1. On the other hand, by

ξ-1

2

^ the definition for λ^ , one readily determine that lim ξ^K 1 ðαiλ^Þ -0. For

ξ-0iλK 2 ðαiλÞ

the other values of ξ, the integrand function gðd; ξ; ωÞ is regular. With the use of the relation between Bessel function and ð2Þ Hankel function, K υ ðiyÞ ¼  ðπ =2Þi1  υ H ð2Þ υ ðyÞ and H n ðxÞ ¼ J n ðxÞ  iY n ðxÞ, one may write gðd; ξ; ωÞ as  8 > ξ½J 1 ðαλ^ ÞJ 2 ðαλ^ Þ þ Y 1 ðαλ^ ÞY 2 ðαλ^ Þ  ið2=παλ^ Þ >  1 sin ðξdÞ ξ o ξp < ^λ½ðJ ðαλ^ ÞÞ2 þ ðY ðαλ^ ÞÞ2  2 2 gðd; ξ; ωÞ ¼   > ξK 1 ðαλÞ > : λK ðαλÞ 1 sin ðξdÞ ξ 4 ξp 2

contour integration of f 1 and f 2 over C 1 and C2 on the upper half-plan (see Fig. 2) that k1 ðd; ωÞ ¼

2π

α

e  i ξp d þ

dZ0

π 1  idξRpn π 1  dξIpn 1 ∑ e þ ∑ e  ; d αn¼1 αn¼1

ð55Þ therefore

d Z0

k2 ðd; ωÞ ¼ k21 ðd; ωÞ þ k22 ðd; ωÞ

ð56Þ

ð49Þ where k21 ðd; ωÞ ¼

k22 ðd; ωÞ ¼

Z ξp

0 h

ξ J 1 ðαλ^ ÞJ 2 ðαλ^ Þ þY 1 ðαλ^ ÞY 2 ðαλ^ Þ  ið2=παλ^ Þ

@

λ^ ½ðJ 2 ðαλ^ ÞÞ2 þ ðY 2 ðαλ^ ÞÞ2 

0

Z

1

ξp





ξK 1 ðαλÞ  1 sin ðξdÞ dξ ¼ λK 2 ðαλÞ

Z 0

1

i

1  1A sin ðξdÞdξ; !

ξ~ K 1 ðαλ~ Þ  1 sin ðξ~ dÞdξ λ~ K 2 ðαλ~ Þ

ð57Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ~ ~ ~ and λ ¼ ξ  ξp ; ξ ¼ ðξ þ ξp Þ. For the treatment of k22 ðd; ωÞ, one may consider the integral of !

hðξÞ ¼ Fig. 2. Contours of integration for f1 and f2.

~ ξ~ K 1 ðαλ~ Þ  1 eiξ d ; λ~ K 2 ðαλ~ Þ

d Z 0;

ð58Þ

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

over a contour similar to C 1 on the upper half-plane (see Fig. 2). By taking the radius of C 1 to be infinity, the contour integral yields k22 ðd; ωÞ ¼ Im Z

Z

0 1

¼ Im 0

B @

1 0

h i ~ ~ ~ ~ ~  ðξ  iξp Þ J 1 ðαλ~ ÞJ 2 ðαλ~ Þ þ Y 1 ðαλ~ ÞY 2 ðαλ~ Þ  ið2=παλ~ Þ ~

~

~

λ~ ½ðJ 2 ðαλ~ ÞÞ2 þ ðY 2 ðαλ~ ÞÞ2 

1 C  iAe  ðξ  iξp Þd dξ

ð59Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ξðξ 2iξp Þ. By setting ω ¼ 0 it is easy to show that for the static case, (56) is identical to those given by Eskandari-Ghadi et al. [4] ! Z 1 2 k2 ðdÞ ¼ dZ0  1 e  ξd dξ; παξ½J 22 ðαξÞ þ Y 22 ðαξÞ 0 ~ where λ~ ¼

ð60Þ The kernel 2=ðζ^  z^ Þ þ 1=ðζ^ þ z^  2^lÞ, in the singular integral Eq. (42), is a generalized Cauchy kernel [3]. In this kernel, the terms 2=ðζ^  z^ Þ and 1=ðζ^ þ z^  2^lÞ become unbounded when either ζ^ approaches z^ or both ζ^ and z^ approach the end point of the cavity, ^l. To investigate the singularity of the shear stress, it is useful to write 0 o ReðβÞ o 1

ð61Þ

F 2 ðz^ ; ωÞ ¼

1

π 1

π

Z

1

^l

Z ^l

1

" # x 2 1 1 þ  v x v v þ x  2vx x þ v þ 2vxð^l  1Þ       1 1 1 1 1  ; ω þ k1 þ  2^l; ω þ 2 k1 v x v x v    

1 1 1 1  ; ω  k2 þ þ 2ð^l  1Þ; ω ; þ k2 v x v x   1 ^ þ l  1; ω : gðx; ωÞ ¼ iπαΘaξp e  iaξp ðð1=xÞ  1Þ ; τ^ 0 ðv; ωÞ ¼ τ^ v

τ^ ðζ^ ; ωÞ ^ 1 dζ ¼ π ζ^  z^ τ^ ðζ^ ; ωÞ

dζ^ ¼

ζ^ þ z^  2^l

Z

1

^l

1

π

v¼

1

ζ^  ^l þ 1

x¼

1 z^  ^l þ 1

:

ð71Þ

Utilizing an even extension of τ^ 0 ðv; ωÞ with respect to the origin changes the limits of integration in (69) to be a symmetric interval as ½  1; 1. Thus, (69) can be written as Z 1 1 Gðjxj; jvjÞτ^ 0 ðjvj; ωÞdv ¼ gðjxj; ωÞ: 1rxr1 ð72Þ 2 1 Recognizing the degree of singular behavior of τ^ ðz^ ; ωÞ as mentioned in Eqs. (61) and (68), the unknown and fundamental solution of (72) may be expressed as

τ^ 0 ðx; ωÞ ¼ Tðx; ωÞð1  x2 Þ  ð2=3Þ

ð73Þ

Z ^l

2

ηðζ^ ; ωÞ dζ^ ; ^ ðζ  z^ Þðζ^  ^lÞβ 1

ηðζ^ ; ωÞ

ðζ^ þ z^  2^lÞðζ^  ^lÞβ

ð62Þ

dζ^ ;

ð63Þ

as in [8], where ηðζ^ ; ωÞ represents the bounded part of the function τ^ ðζ^ ; ωÞ. It is readily observed that both functions F 1 ðz^ ; ωÞ and F 2 ðz^ ; ωÞ are singular at z^ ¼ ^l. One may decompose the function F n ðz^ ; ωÞ, n ¼1 and 2, into a singular and a regular part as Pak and Abedzadeh [11] and write F 1 ðz^ ; ωÞ ¼ ηð^l; ωÞcotðπβ Þeπ iβ =ðz^  ^lÞβ þ F 1 n ðz^ ; ωÞ;

ð64Þ

F 2 ðz^ ; ωÞ ¼ ηð^l; ωÞeπ iβ = sin ðπβÞðz^  ^lÞβ þ F 2 n ðz^ ; ωÞ;

ð65Þ

where F n1 and F n2 are regular and thus bounded with the following conditions  α  n  F ðz^ ; ωÞ oc1 = z^  ^l 1 ; ð66Þ 1  α  n  F ðz^ ; ωÞ oc2 = z^  ^l 2 : 2

ð67Þ

Here c1 ,c2 , α1 and α2 are real positive constants such that α1 and α2 are smaller than ReðβÞ. Substituting (64) and (65) into (42), one may find the characteristic equation for the exponent β as 2 cotðπβ Þ þ 1= sin ðπβÞ ¼ 0:

interval  1 rx r1, and ð1  x2 Þ  3 is analogous to a weighting function. Since the latter is the same as the weighting function for Jacobi polynomials P nðγ ;βÞ ðxÞ for γ ¼ β ¼  2=3 [3], one may use a numerical procedure based on Gauss-Jacobi integration formula to solve the integral Eq. (72). With the use of (73), (72) is written as Z 1 1 Tðjvj; ωÞ Gðjxj; jvjÞ dv ¼ gðjxj; ωÞ; 1rxr1 ð74Þ 2 1 ð1  v2 Þ2=3 With the aid of the Gauss–Jacobi numerical integration rule and the collocation method [3], one may change the integral Eq. (74) to a system of linear algebraic equations as         1 2N ∑ W Gðxj ; vk ÞTðvk ; ωÞ ¼ gðxj ; ωÞ; 2k¼1 k

4. Numerical solution of singular integral equation

ð1=3;1=3Þ

P 2N  1 ðxj Þ ¼ 0;

ð75Þ

j ¼ 1; :::; 2N  1

ð76Þ

and ð  2=3;  2=3Þ

P 2N

ðvk Þ ¼ 0;

k ¼ 1; :::; 2N

ð77Þ

respectively. In addition, the weighting function is given by [3] Wk ¼

 ð4N þ γ þ β þ 2ÞΓ ð2N þ γ þ 1ÞΓ ð2N þ β þ 1Þ2γ þ β ; ðγ ;βÞ ðγ ;βÞ ð2N þ 1Þ!ð2N þ γ þ β þ 1ÞΓ ð2N þ γ þ β þ 1ÞP 2N þ 1 ðvk ÞdP 2N ðvk Þ=dv

ð78Þ for γ ¼ β ¼  2=3. Since the roots xj and vk are symmetric with respect to the origin, (75) can be reduced to a simpler system of equations as N

Because of its complex kernel function, the governing integral Eq. (42) cannot be solved analytically. For its evaluation, it is convenient to first write the integral Eq. (42), as in Pak and Abedzadeh [10], in the form of Z 1 Gðx; vÞτ^ 0 ðv; ωÞ dv ¼ gðx; ωÞ; 0rvr1 ð69Þ

j ¼ 1; :::; 2N  1

where 2N is the number of collocation points, and xj and vk are the roots of Jacobi's functions

ð68Þ

whose root is β ¼ 2=3 as in the isotropic case [11].

0

;

ð70Þ

where Tðx; ωÞ is regular and analytic with respect to x in the

and define F 1 ðz^ ; ωÞ ¼

where

Gðx; vÞ ¼

hðξÞdξ

ηðζ^ ; ωÞ τ^ ðζ^ ; ωÞ ¼ ^ β ; ðζ  ^lÞ

257

∑ W k Gðxj ; vk ÞTðvk ; ωÞ ¼ gðxj ; ωÞ:

k¼1

j ¼ 1; :::; N

ð79Þ

Eq. (79) provide N equations to determine N unknowns Tðvk ; ωÞ at N collocation points vk ; k ¼ 1 to N. The last equation in (79) involves the root xN ¼ 0, which corresponds to infinity for z, and it is automatically satisfied. On the other hand, the condition u1 ða; l; ωÞ ¼ Θ a may be satisfied by the following equation

258

Z

1 l

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

½φk1 ða; ζ þ l; ωÞ þ φk1 ða; ζ  l; ωÞτðζ ; ωÞ dζ ¼ α π Θ a

ð80Þ

With the use of Eqs. (41), (70) and (71), Eq. (80) is changed to N

∑ W k GN ðvk Þ Tðvk ; ωÞ ¼ α Θ π ;

ð81Þ

k¼1

where GN ðvÞ ¼

   

1 1 1 φk1 a; þ 2^l  1; ω þ φk1 a;  1; ω : 2 v v v

ð82Þ

Eqs. (79) and (81) are a system of linear algebraic equations, whose solution can be computed numerically. Then, with the use of Tðv; ωÞobtained from (79) and (81), the displacement and stresses are written as   a N 1 ∑ W k φk1 r^ ; þ ^l  1 þ z^ ; ω u1 ðr; z; ωÞ ¼ v απ k ¼ 1  

1 ^ Tðvk ; ωÞ þ φk1 r^ ; þ l  1  z^ ; ω ; v v2k   aμ0 N 1 τzθ1 ðr; z; ωÞ ¼ ∑ W k φk2 r^ ; z^   ^l þ 1; ω απ k ¼ 1 v^  

1 ^ Tðvk ; ωÞ þ φk2 r^ ; þ l  1 þ z^ ; ω ; v v2k aμ N 1 τrθ1 ðr; z; ωÞ ¼ ∑ W k φk3 ðr^ ; þ ^l  1 þ z^ ; ωÞ v απ k ¼ 1

1 ^ Tðvk ; ωÞ þ φk3 ðr^ ; þ l  1  z^ ; ωÞ ; ð83Þ v v2k

K TT ðl; ωÞ is evaluated for isotropic material with no cavity i.e. l=a ¼ 0, and compared with the results reported by Luco and Mita [7]. Fig. 3 illustrates the real and imaginary parts of dimensionless impedance function in terms of dimensionless frequency for different values of N. As it is observed, there exist some errors if the numerical results are truncated for N ¼20, and the larger the frequency the larger the error results in. The errors in a wide range of the frequencies are decreased when N ¼ 40, but the curve is not collapsed on the solution given by Luco and Mita [7]. Moreover, there are some large errors for larger frequencies. However, if N is selected as 60 or larger than 60, then very accurate numerical evaluations are achieved for a very wide range of frequency, say up to ω0 ¼ 4:0, and it is not very well for ω0 4 4:5. The error at ω0 4 5:0 for N ¼80 is about 3.43 percent for real part and 1.35 percent for imaginary part, which is accepted in this paper, and thus, the numerical evaluations of the rest of this paper is done with N ¼80. To prove that the results are in good accuracy, two more comparisons are made one of which is for static analysis of an isotropic half-space containing a cylindrical cavity and the other is for a static forced displacement of a transversely isotropic half-space by a rigid disc attached on the free surface of the halfspace. Fig. 4 illustrates the function TðxÞ=T l for an isotropic halfspace with a finite length cylindrical cavity, which is going to displaced statically with a rigid disc. The function TðxÞ=T l has been plotted for l=a ¼ 0:001; 0:01; 0:05; 0:1 and 1:0, from this

for Region 1, and

  N 1 ∑ W k φI1 r^ ; þ ^l  1  z^ ; ω v απ k ¼ 1  

1 ^ Tðvk ; ωÞ ^ þ Θ a r^ e  iξp aðz^  lÞ ;  φI1 r^ ;  l  1 þ z^ ; ω v v2k   aμ0 N 1 τzθ2 ðr; z; ωÞ ¼ ∑ W k φI2 r^ ; þ ^l  1  z^ ; ω v απ k ¼ 1  

1 ^ Tðvk ; ωÞ ^  iμ Θar^ ξp e  iξp aðz^  lÞ ; þ φI2 r^ ;  l  1 þ z^ ; ω v v2k   aμ N 1 τrθ2 ðr; z; ωÞ ¼ ∑ W k φI3 r^ ; þ ^l  1  z^ ; ω v απ k ¼ 1  

1 ^ Tðvk ; ωÞ ; ð84Þ  φI3 r^ ;  l  1 þ z^ ; ω v v2k u2 ðr; z; ωÞ ¼

a

for Region 2, where r^ ¼ r=a. In addition, the torsion impedance function related to the rigid disc is defined as: K TT ðl; ωÞ ¼

T l ðωÞ

Fig. 3. Comparison of the results of this study with different number of Jacobi's points with that of Luco and Mita [7].

ð85Þ

Θ

where T l with the use of Eq. (25) is determined as N

T l ðωÞ ¼  2πμ a2 ∑ W k k¼1

Tðvk ; ωÞ ; vk 2

ð86Þ

5. Computational results The procedure for numerical evaluations has been given in previous section. With the use of this procedure, we are going to investigate the effect of both the degree of anisotropy and the length of cavity on the displacement, shear stress and the impedance functions. It needs to be pointed out that all numerical results presented here are dimensionless, with a non-dimensional pffiffiffiffiffiffiffiffiffiffi frequency defined as ω0 ¼ aω ρ=μ0 . First of all, the value of N is determined in such a way that the numerical results to be obtained accurately. For this reason, the impedance function

Fig. 4. Comparison of TðxÞ in isotropic elastic half-space with Pak and Abedzadeh [11] for l=a ¼ 0:001; 0:01; 0:05; 0:1 and 1.

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

259

research and the paper published by Pak and Abedzadeh [11]. As seen, very accurate results are observed. In addition, it is seen that although the function TðxÞ is a regular function, its gradient at the vicinity of the edge of the rigid disc is different from surface, l=a ¼ 0:001, to the depth l=a ¼ 1:0. Fig. 5 compares the

rotation/displacement of a radial line at different depths, say z=a ¼ 0; 0:4; 0:8; 1:2; 1:6 and 2:0, of a transversely isotropic half-space under Reisner-Sagoci condition due to a rigid disc attached on the surface of a the half-space evaluated from this

Fig. 5. Comparison of displacement solutions at z=a ¼ 0; 0:4; 0:8; 1:2; 1:6 and 2:0 versus r=a in transversely isotropic elastic half-space with Rahimian et al. [12].

Fig. 8. Real and Imaginary parts of torsional impedance functions versusω0 for different depth of cavity indented isotropic half space.

Fig. 6. Real parts of solution TðxÞ for torsional rotation of rigid disc with amplitude Θ and ω0 ¼ 0:5.

Fig. 7. Imaginary parts of solution TðxÞ for torsional rotation of rigid disc with amplitudeΘ and ω0 ¼ 0:5.

Fig. 9. Real and Imaginary parts of torsional impedance functions versus ω0 for different depth of cavity indented transversely isotropic half space.

Fig. 10. Real parts of normalized angular displacement at z ¼ l in terms of horizontal distance for ω0 ¼ 0:5.

260

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

study and the research reported by Rahimian et al. [12]. Again, an excellent agreement can be observed from this figure. Figs. 6 and 7 show the real and imaginary parts of the regular part of the shear stress, Tðx; ω; lÞ, in the terms of x for different degrees of anisotropy and different length of cavity when the dimensionless frequency is ω0 ¼ 0:5. We should notice that x ¼ 0

corresponds to infinity and x ¼ 1 is related to the bottom of cavity. As it is observed, the function Tðx; ω; lÞ in dynamic case treats as static counterpart, which means the smaller the depth of cavity the larger the gradient of the function Tðx; ω; lÞ infers. Furthermore, both the degrees of anisotropy and the length of cavity affect

Fig. 11. Imaginary parts of normalized angular displacement at z ¼ l in terms of horizontal distance for ω0 ¼ 0:5.

Fig. 14. Real parts of normalized shear stress τzθ at z ¼ l in terms of horizontal distance for ω0 ¼ 0:5.

Fig. 12. Real parts of normalized angular displacement along r ¼ a in terms of normalized cavity depth for ω0 ¼ 0:5.

Fig. 15. Imaginary parts of normalized shear stress τzθ at z ¼ l in terms of horizontal distance for ω0 ¼ 0:5.

Fig. 13. Imaginary parts of normalized angular displacement along r ¼ a in terms of normalized cavity depth for ω0 ¼ 0:5.

Fig. 16. Frequency variation of real part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼0.5 and l ¼ 2a.

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

the shear stress significantly, and as a result none of these two parameters can be neglected. Figs. 8 and 9 illustrate, the real and imaginary parts of nondimensional torsion impedance function in terms of dimensionless

Fig. 17. Frequency variation of imaginary part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼ 0:5 and l ¼ 2a.

Fig. 18. Frequency variation of real part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼ 0:5 and l ¼ 5a.

Fig. 19. Frequency variation of imaginary part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼ 0:5 and l ¼ 5a.

261

frequency for isotropic and transversely isotropic material, respectively for α ¼ 1:0 and α ¼ 0:5, however for different values of length of cavity from zero to a very large length (l ¼ 30a). l=a ¼ 0 is

Fig. 20. Frequency variation of real part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼ 1:0 and l ¼ 2a.

Fig. 21. Frequency variation of imaginary part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼ 1:0 and l ¼ 2a.

Fig. 22. Frequency variation of real part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼ 1:0 and l ¼ 5a.

262

A. Ardeshir-Behrestaghi et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 252–262

6. Conclusion

Fig. 23. Frequency variation of imaginary part of displacement/rotation, uða; zÞ=Θa, at different depth either on the wall of cavity or below the cavity for α ¼ 1:0 and l ¼ 5a.

related to surface foundation and l ¼ 30a may be corresponded to the bottom part of a deep pile. As observed there exist significant difference between surface foundation and foundations with cavity, which is due to singular behavior of the cavity. However, since singular behaviors of cavities with different lengths are almost the same, the behavior of cavities with different lengths is almost the same. On the other hand, as it is observed from Fig. 8, the real part of the torsion impedance function treats different compared to other impedance functions like vertical and rocking impedance functions. Figs. 10–13 depict respectively the real and imaginary parts of circumferential displacement in terms of either radial distance from r ¼ 0 to a large distance or depth from z ¼ 0 to z ¼ l, for different lengths of cavities and for different degrees of anisotropy. First, it is seen from Figs. 10 and 11 that the displacement boundary condition has been satisfied very accurately. Second, it is clear from these figures that the effect of anisotropy in vertical direction is more significant than in radial direction, while the effect of length of cavity is significant in both vertical and radial directions. Third, the displacement, for the same dimensionless frequency, is more wavy for larger α than smaller one, which means the less the α  value the larger the stiffness of the system results in. The radial variation of the real and imaginary parts of shear stress for dimensionless frequency of ω0 ¼ 0:5 are depicted in Figs. 14 and 15, where as expected a singular behavior is seen at the vicinity of the edge of the rigid disc. Although the degree of singularity is not affected by the length of cavity, however it treats different for different degrees of anisotropy, and the larger the value of α the more singular the variation of shear stress is. One significant phenomenon that can occur in the excitation of a layered medium is dynamic resonance. This is illustrated in Figs. 16–23, where the primary displacement responses at some points either on the wall of the cavity or below the cavity due to torsional excitation of the rigid plate as a function of the excitation frequency, length of cavity and anisotropy of the half-space are given. First of all, as it is expected, the static part (the zone with short dimensionless frequency) is completely separated from the dynamic part (with frequency content in the range of ω0 4 0:5). As seen both the wave length and the amplitude of the displacement/ rotation is different from l ¼ 2a and l ¼ 5a, from one point to another point, and from one material to another one. As illustrated the amplitude for some points are very large, while the input excitation are the same for them.

The response of a transversely isotropic half-space with a finite length open circular cylindrical cavity to a mono-harmonic motion applied on a rigid circular disc attached on the surface of the bottom of the cavity has been rigorously investigated. By means of Fourier sine and cosine transforms, the dynamic mixed boundary value problem has been transformed to a generalized Cauchy singular integral equation for the traction in between the rigid disc and the bottom of cavity. Investigating the kernel of the generalized Cauchy singular integral equation in detail results in a decomposition for the unknown traction as a multiplication of a known singular function and an unknown regular one. Representing the unknown part of traction in the form of a series with unknown coefficients, enables us to determine the known coefficients with the use of the Gauss–Jacobi procedure and collocation method. Numerical results show an excellent agreement with the existing results for a simpler case of Reisner–Sagoci problem of a transversely isotropic half-space and an isotropic half-space containing finite length cylindrical cavity both in the static case. The effects of material anisotropy and length of cavity have been investigated numerically and depicted graphically. It is shown that both degree of anisotropy and the length of cavity have some significant effects on the displacement and shear stress, which means that none of these two parameters can be neglected.

Acknowledgment The partial support from the University of Tehran through 27840/1/ 05 to M. E.-G. during this work is gratefully acknowledged. References [1] Andrews LC. Special functions of mathematics for engineers. 2nd ed. Oxford University press; 1998. [2] Collins WD. The forced torsional oscillations of an elastic half-space and an elastic stratum. Proc London Math Soc 1962;S3–S12:S226–44. [3] Erdogan F. Mixed boundary value problems in mechanics. Mech Today 1976;4:1–86. [4] Eskandari-Ghadi M, Mahmoodian M, Pak RYS, Ardeshir-Behrestaghi A. Analytical solution of torsion vibration of a finite cylindrical cavity in a transversely isotropic half-space. J ZAMM Z Angew Math Mech 2012;92(7):583–95. [5] Jordan DW. The stress wave from a finite cylindrical explosive source. J Math Mech 1962;11:503–51. [6] Lekhnitskii SG. Theory of elasticity of an anisotropic body. Moscow: Mir Publisher Moscow; 1981. [7] Luco JE, Mita A. Response of a circular foundation on a uniform half-space to elastic waves. Earthq Eng Struct Dynam 1987;15:105–18. [8] Muskhelishvili NI. Singular integral equations Boundary problems of function theory and their application to mathematical physics. 2nd ed. Mineola, New York: Dover Publications; 1992. [9] Pak RYS, Sophers J-DM. Rocking rotation of a rigid disc in a half-space. Int J Solids Struct 1991;28:389–401. [10] Pak RYS, Abedzadeh F. Torsional traction on an open finite cylindrical cavity. Proc R Soc London A 1992;438:133–44. [11] Pak RYS, Abedzadeh F. A torsional contact problem for an indented half-space. J Appl Mech 1996;63:1–6. [12] Rahimian M, Ghorbani-Tanha AK, Eskandari-Ghadi M. The Reissner–Sagoci problem for a transversely isotropic half-space. Int J Numer Anal Methods Geomech 2006;30(Issue 11):1063–74. [13] Rajapakse RKND, Gross D. Traction and contact problems for an anisotropic medium with a cylindrical borehole. Int J Solids Struct 1966;33:2193–211. [14] Reissner E, Sagoci HF. Forced torsional oscillations of an elastic half-space, I. J Appl Phys 1944;15:652–4. [15] Sagoci HF. Forced torsional oscillations of an elastic half-space, II. J Appl Phys 1944;15:655–62. [16] Tranter CJ. On the elastic distortion of a cylindrical hole by a localized hydrostatic pressure. Q Appl Math 1946;4:298–302. [17] Westergaard HM. Theodore von Karman anniversary volume. Caltech 1941:154–61. [18] Watson GN. A treatise on the theory of Bessel Function. Cambridge: University press; 1958. [19] Williams WE. Note on a problem in elastodynamics. Proc London Math Soc 1963;138:119–22.