green's function for tangentialy loaded HORIZONTALY LAYERED HALF-SPACE TOMAŽ PLIBERŠEK and ANDREJ UMEK About the authors Tomež Pliberšek University of Maribor, Faculty of Civil Engineering, Smetanova 17, 2000 Maribor, Slovenia E-mail: tomaz.plibersek@uni-mb.si Andrej Umek University of Maribor, Faculty of Civil Engineering, Smetanova 17, 2000 Maribor, Slovenia E-mail: umek@uni-mb.si Abstract The topic of this paper is a novel evaluation of the integral representation of the surface Green's function for a layered half-space, loaded on its surface by a harmonic tangential point force. The equations of motions are reduced to wave equations by the introduction of wave potentials. The Hankel transform is applied to them and they are consecutively solved leading to the integral representation of surface displacements. They are consecutively evaluated by the proposed three step procedure. First the singularity is extracted. It is further noted that so obtained integrals, after suitably chosen branch cuts and analytic continuation of integrands are introduced, can be evaluated by contour integration for an arbitrary number of layers. They are, therefore, expressed by number of residues at the poles of integrand and the integrals along finite portions of the branch cuts. The latter ones can easily be evaluated to any desired accuracy leading to a closed form solution. Some numerical results corroborating the presented approach are given. Keywords elastodynamics, elastic wave propagation, Green's function, horizontaly layered half-space, horizontal point load 1 INTRODUCTION Modeling the elastodynamic characteristics of soil is required in a number of engineering problems, e.g. dynamic soil-structure interaction, dynamically loaded foundations etc. The soil is geometrically in prevalent number of cases modeled as a half space, which is, to be more realistic, endowed with some structure. As a starting point to determine the elastodynamic characteristics of soil one can use the fundamental solution or the Green's function. The use of the fundamental solution, which is known from the literature, results in the integrals over the whole surface of the half-space, interface between the soil and the superstructure and over the interfaces between the materials with different elastodynamic properties. Some of them are of infinite or semi-infinite extend. Their evaluations, which can be performed only numerically, represent the major difficulty of this approach. For practical evaluation the integration area has to be made finite, what results into introduction of the fictitious boundary, where the radiation conditions should be satisfied. The latter represent a demanding and not completely resolved problem e.g. Premrov [1]. The Green's function approach leads us to the integrals extending across the interface surface between the soil and the superstructure only. Their evaluation is straightforward and relatively easy to perform, ones the Green's function is given. The difficulty of the problem is transferred to the determining and evaluating the Green's function itself. The problematic of the homogeneous as well as the layered half-space has drawn the attention of considerable number of authors, not all of them can be mentioned here. The first elasticity solutions for whole-and half-space problems, static as well as dynamic ones have been obtained by Kelvin [2], Boussinesq [3], Cerruti [4], Lamb [5], Mindlin [6] and the others. The basic results on layered media were presented by Ewing et al. [7]. In more recent times the authors sought the solutions by two basically different approaches. On one side there are the approximate methods e.g. Luco's ray method [8, 9], Kausel's [10] thin layer method as well as BEM and FEM methods, their combinations e.g. ACTA GeOTeCHNICA SLOVENICA, 2008/i 37. T. PLIBERŠEK & A. UMEK: GREEN'S FUNCTION FOR TANGENTIALY LOADED HORIZONTALY LAYERED HALF-SPACE Gaitanaros et al. [11], Triantafyllidis [12], Auersch [13] and refinements e.g. Aubry et al. [14]. On the other hand we have methods leading to exact solutions in the form of integrals of semi- or infinite extend e.g. Vostroukhov [15] and Jin [16]. Their evaluation in the concept of FFT [17] concludes the latter approach. It, however, appears that the use of FFT like integration is a successful method to evaluate the Hankel transforms inversion integrals in the cases, where there are no singularities in the resulting displacement field. It must be however noted that in the case of the Green's function the vital information about its singularity comes from the portion of integration path, where the integration variable is very large or infinite. This fact makes in the case of the Green's function the FFT like evaluation of integrals less efficient. Kobayashi et al. [18] considered the homogeneous, elastic half-space and succeeded to reduce the semi-infinite integrals representing the Green's function to a part containing the singularity, the residue at the Rayleigh pole and the integrals of finite extend along the properly chosen branch cut, which can be easily evaluated with any desired accuracy. Štrukelj et al. [19] succeeded to modify the Kobayashi's approach, so that it could be applied to the problem of vertically loaded layered half-space. The authors have first derived the Green's function for a single layer [20] and have demonstrated that under the assumption of infinite thickness of the layer their solutions lead to the Kobayshi's ones. Our decision to focus on the layer has also been motivated by the investigations of the mechanical properties of soils e.g. Žlender et al. [21] and [22], which are given at a point and can be easily generalized to a layer. This paper continues and completes the method to determine the surface Green's function for a layered half-space loaded with the harmonic force acting in any direction. It is motivated by the previous works by Kobayashi et al. [18, 23], Štrukelj et al. [19] and Pliberšek et al. [20]. The load acting in a general direction is decomposed into normal and tangential components with respect to the surface of the half-space, so that the problem of the latter one has to be dealt with only. First the general equations of motion for a single layer in Hankel transform domain are derived by adapting the Vostroukhov's [15] approach and then transformed back into the geometrical domain. These single layer solutions are then combined into the solution for a layered half-space making use of the continuity conditions on the interfaces between the layers and boundary conditions on the surface of the half-space. It is proven that the basic mathematical properties of these solutions do not depend on the number of layers. They are the same for a homogeneous half-space as well as for the half-space with any number of layers. The inverse Hankel transform integrals appearing in the Green's function can be therefore for any number of layers expressed with part proportional to r— , integrals of finite extend along the properly chosen branch cut and some residues at the poles of the integrand. The closed form solutions, obtained by our approach, for the Green's function of the tangentially loaded layered halfspace are consecutively presented graphically for some selected number of cases. 2 METHOD OF ANALYSIS 2.1 GOVERNING EQUATIONS FOR A LAYER We consider a horizontally layered half-space, which consist of n layers on a half-space, as shown on the Fig. 1. The material properties of each layer and of the underlying half-space are assumed to be isotropic and homogeneous. Shear modulus E, Poisson's ratio vi, mass density p and the dumping coefficient E are the material constants of i-th layer. The homogeneous halfspace is labeled as the layer number H. The global, cylindrical co-ordinate system and the local cylindrical co-ordinate systems having their origins at the top of each layer are introduced. The model is subjected on its free surface to a concentrated tangential point load, which varies harmonically in time. Without loss of generality, it is assumed that it acts in the direction of positive x-axis. The governing equation for each homogeneous, elastic layer is the well known [24] reduced Navier's equation of motion in frequency domain: Elv2^ + VV• u, = -w2ui ; i e[ 1,2,...,n,H]. (1) Pi Pi The boundary conditions on the free surface of the halfspace can be in the cylindrical coordinates written as: 4 (r,w) = a„ (r,0,w) • cos(^) = -9MJW • cos(^) (2) zj (r, ê, z, u) = aêz (r, z, u) ■ sin (ê) = 2 ■ n ■ r Q(")■ S (r) . 2 ■ n ■ r azzJ (r,ê,0,u) = 0 . (4) sin(ê) (3) where S (r ) is the Dirac delta function. The continuity equation along the interfaces of consecutive layers, where perfect bonding is assumed, and radiation conditions for r complete the definition of boundary value problem under consideration. 52. ACTA GEOTECHNICA SLOVENICfl, 2008/l T. PLIBeRSeK & A. UMeK: GREEN's FUNCTION FOR TANSeNTIALY LOADeD HORIZONTALY LAYeReD HALF-SPACe r Figure 1. A horizontally layered half-space subjected to a surface horizontal harmonic point load. The geometry of the above defined boundary-value problem is axi-symmetric. Therefore the 9-dependence of the problem is governed by the 9-dependence of the loading only. Due to very simple 9-dependence as a first step in the solution procedure reduced stresses: arzi (r,',z,w) = arzi (r,z,w)-cos(') ; i g[i,2,...,n,H] (5) J (r,z,w) = t (r,z,w)- sin(') ; i e[1,2,...,n,H] (6) azzj (r, z, w) = azzj (r, z, w)-cos(') ; i s[l,2,...,n,H] (7) and reduced displacements: ur t (r,z,w) = ur t (r,z,w)- cos(') ; i e[ 1,2,...,n,H] (8) u,ffJ (r,z,w) = u^ (r,z,w)- sin(') ; i e[ 1,2,...,n,H] (9) uz. (r,z,w) = uz. (r,z,w) - cos(') ; i e[ 1,2,...,n,H] (10) are introduced. In the second step we make use of the Helmholtz wave potentials [25] to decouple the equations of motion (1): ui = V-i. (r,z,w) = }zj (r,z,w)-sin(') ; i e[1,2,...,n,H]. (16) The substitution of Eqs. (13) to (16) into equations of motion (1) yields: ACTA SeOTeCHNICA SLOVENICA, 2008/l 53. T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LAYERED HALF-SPACE d2& dr2 1 - + -r dr 1 _ 7 ■ & d2& dz2 ) + kl = 0 ; i G [ 1,2,. .,«, H (17) d2 dr2 1 r d dr 2 - + r d2 dz2 ■ Vr,i 2 - +— ■ Vv,, + r -kT,, ■ Vr,, = 0 ; i G [ 1,2,... «, H ] (18) d2 dr2 1 + -r d dr 2 -- + r d2 dz2 ■  2 - + — ■ Vr r ,, + -kl. ■ = 0 ; i G [ 1,2,... «, H ] (19) d2 dr2 1 r d dr 1 - + r d2 dz2 ■ V,, + k 2 T ,i = 0 ; i G [ 1,2,.. ,«, H (20) Vr,." -Âi + r ■ dVr dr 1 dVz,i dz = 0 ; i G [ 1,2,...,«,H] (21) We note that the above equation system contains two decoupled equations (17) and (20) and two still coupled equations (18) and (19). To decouple the latter ones we add them and subtract them and introduce two new reduced wave potentials Xi and Ki as: X, = V,i + i K = V,i _ V (22) where i G [ 1,2,...,«, H] and the equations (18) and (19) are replaced by: d2 +1 Br1 r dr Dz2 d2 +1 ±+ Br2 r dr r2 dz X, + K,, ■ X, = 0 ; iG[ 1,2,...,«,H] (23) K + k2Tii ■ K = 0 ; i G [1,2,...,«,H] (24) and the compatibility condition becomes: 1 d(X, + K ) dV,, k, + r dr dz = 0 ; i G[ 1,2,...,«,H]. (25) The most efficient way to solve the boundary value problem under consideration is by the introduction of Hankel integral transform. 2.2 HANKEL TRANSFORM AND SOLUTIONS IN HANKEL TRANSFORM DOMAIN To solve the equations of motion (17), (20), (23) and (24) for each layer with appropriate boundary conditions (2) to (4) and continuity conditions Hankel integral transform r ^ Ç : m fH" (Ç ) = Hn [f (r )] = f f (r )■ J« (Ç r) r ■ dr (26) and its inverse £ — r : to f (r ) = H- [fHn (£)] = f fH«(£ ) ■ J«(£ r )• £ • d£ , (27) 0 are introduced [26]. n is the order ofthe transform and Jn(£ r j is Bessel function of the first kind and order n. To transform the equations of motion to their canonical form the integral transforms of different orders are employed. Equations (17) and (20) are transformed through Hankel transform of order 1, equation (23) of order 0 and equation (24) of order 2. This yields the following system of equations: d (ç ) dz2 d 2X,-H° (Ç) dz J2£ H _(2 _ kl)H1 (£) = 0 ; i G[ 1,2,...,«,H] (28) _((2 _ )-XH° (Ç) = 0 ; i e[ 1,2,...,«, H ] (29) d2kih2 (Ç) dz2 _((2 _ Kj )■K1 (Ç) = 0 ; i G[ 1,2,...,«,h] (30) dV H (ç) dz2 _(( _ k2,i H/' (Ç) = 0 ; i g[ 1,2,...,«, H ] . (31) The solutions of the above equations can be written as: &Hl = G,,, ■ ea z> + c2,t ■ e_a''z' : ; ÎG[ 1,2,...,«, H ] (32) X,H0 = G3,i ■ efi ' z + C4J ; ÎG[ 1,2,...,«, H ] (33) h2 = C5, ■ e * z + c6i ■ e_* z : ; ÎG[ 1,2,...,«, H ] (34) VH = C7, ■ e ■z'+ CSj ■e_*<■« ; , iG[ 1,2,...,«,H] , (35) where: a =yl Ç2 _ kl, ; = J Ç2 _ k Ij ; iG[ 1,2,...,«, H ] . (36) 2 0 54. ACTA GEOTECHNICA SLOVENICA, 2008/l T. PLIBËRSËK & A. UMEK: GREEN'S FUNCTION FOR TflNGËNTIflLY LOflDËD HORIZONTflLY LflYËRËD HALF-SPACE The constraint condition the transformed wave potentials must satisfy is obtained from equation (25). This yields: dfe (£, w £■ KH2 (£,z,w)-£■ xH" (£,z,w) + 2• i e[ 1,2,...,«, H ] dz = " ; (37) The yet unknown integration constants Cj ; i e[ 1,2,...,8] and j e [ 1,2,...,«,H] will be determined from boundary, continuity and radiation conditions. For this purpose the reduced displacement components are expressed through the transformed wave potentials as: 1 to ri(^zw)=- f£ ■ " ;>2„T H £■0H1 (£,z,w) + £■ AH- (£,z,w)- dKH (£,z,w) dz 2 d2^z,iHl (£,z,w) + £ dz2 Ok*2 (£,z,w) dz ■J"(£r)- -£■vH1 (£,z,w)+ £■ AH (£,z,w)+ (38) ■ J 2 (£ r )[■ d£ ; i e[ 1,2,...,«, H ] w( r, z, w)=2 f£ " U2 T H, £ ■ tH1 (£,z,w)- OkH2 (£,z,w) dz + 2 d2tH (£,z,w) £ ' dz2 Ok"2 (£,z,w) dz +£ ■ A J1 (£, z, w) J " (£ r ) + l-£ ■tH1 (£, z, w) + +£ ■A H (£, z, w] J2 (£r)}■ d£ ; i e[ 1,2,...,«,H] (39) 'z,i(r,z,w) = f £ dTH1 (£, z, w) dA H (£, z ,0 dz dz (4") -£■ k"2 (£,z,w)|- J, (£r)■ d£ ; i e[ 1,2,...,«,H]. From the above expressions strains are derived and introduced into the constitutive equation for linear, homogeneous and isotropic solid. The pertinent reduced stresses are then given as: I TO ëz.i(^z w) = y ■] f £' 2 . dTH- (£ ) 2 d3AH (£,z,w) 2 ■ £--(£, z, w)------ dz ' £ dz3 d2kH2 (£,z,w ^2 Th dz2 - £2 ■ KH2 (£, z, w] J "(£ r )■ (41) +1-2 ■ £ ■ dTH1 (£, z, w) + 2 ■ £ ■ dA H (£, ^ w) dz dz + + £2 ■ ¿H (£, z, w)- J2 (£ r ) [ ■ d£ L ; i e[ 1,2,...,«, H] " ACTA SEOTECHNICA SLOVENICA, 2008/l 57. T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LAYERED HALF-SPACE w(r,z,w) = E•!/£•! -2• £ dy,Hl (£,z,w) 2 l Jo * dz 2 d3i>zH (£,z,w) d'ËH (£,z,w) dz3 +[-2 • £ •ÎL(£,z,w) + , dz + £2 • (£,z, w)l J2 ( / dz2 +£2 • (£, z, w) + J o ( r ) + 0%H2 (X,z,w) + 2 £ dtH (£,z,w) -3-+ 2 •£-- dz dz ; i e[ 1,2,...,«,H] (42) to ëzz,i (r, z, w) = E f £ • o f2 H1 \ + 2 • e d%Hl (X,z,w) dz2 - £2 •r1 (£, z, w; + 2 • £2 • fp*1 (£,z,w) - 2 • £ •d(£,z,w) + dz - 2 d2t Z1 (£, z, w) dz2 • J1 (r)• d£ ; i e[ 1,2,...,«,H]. (43) The equations (38) to (43) permit us to construct a liner, algebraic equation system in integration constants Cij as unknowns. It is worth-while to note that in these equations due to constraint condition (37) only three wave potentials are appearing, what results in six integration constants per layer or underlying half-space. If a branch cut in the complex £-plane is introduced, which makes Re(a; )> 0 and Re(^)> 0 on the real positive £-axis, then the radiation conditions demand: C1,H = C5,H = C7,H = 0 (44) The boundary conditions on the surface of the halfspace are given by equations (2) to (4) and for the continuity condition perfect bonding between the layers and the underlying half-space i.e. the continuity of displacements and stresses across the interfaces, is assumed. Thus the boundary and continuity conditions result in 6 • « + 3 equations for the same number of unknown integration constants Cij. For the reason of in the forthcoming paragraph introduced solution procedure we must study some properties of the matrix of this equation system and its submatrices. To simplify the further mathematical derivation some dimensionless variables and constants are introduced: £ = n • kT Yi =" ; a,I =- ai = kT1 •ai ; P = kT 1 • P t ; - k-T* 1 ' h; .El E (45) i e[ 1,2,...,«,H]. We begin the derivation of the equation system with the boundary conditions defined on the upper surface of the first layer shown on the Fig. 2. Figure 2. Material and geometrical properties of the first layer. 54. ACTA GEOTECHNICA SLOVENICA, 2008/l T. PLIBËRSËK & A. UMEK: GREEN'S FUNCTION FOR TflNGËNTIflLY LOflDËD HORIZONTflLY LflYËRËD HALF-SPACE We introduce the expressions for stresses (41) to (43) into the transformed boundary conditions (2) to (4) and take into account the equation (45). This yields: a • kTi -|-2 • n • a -[ci.i- ci.2 ] + (2 + A2 )-[Ci,3 + Ci,41 + +-• Â3 •Ci, - Ci, l) = 9M n ) n (46) Mi • kTi "I-2 n • ai •fCi.i - Ci.2 ] + (2 + Âi2 WCi.3 + Ci.4 1 + +2 • n • Âi •Ci.5 - Ci.61} = 0 (47) Mi • kTi 1 (Ai + 2 • Mi) 2 2\ , 2 ----(i -n ) + n 2 Mi |Ci.i + Ci.2 I + -n • A • [Ci.3 - Ci.4 ] - Âi2 • Ci.5 + Ci.< (48) = 0 The other equations follow from continuity conditions on the interfaces between layers and the nth layer and the underlying half-space respectively. The interface between the ith and the (i+1)st layer is depicted in Fig. 3. _ Figure 3. Material and geometrical properties of ith and (i+1)st layer. Introduction of equations for the displacements and the stresses (38) to (43) into continuity equations leads to: azzA h - a \Zi =h: zz .i + il Mi.i lz;+, —0 f i (A, + 2 • M, ) •(a,2 - n2 )+n2 12 m, IQ. • ea * + Ci.2 • e~a'* 1 - n • ( • I Ci.3 • e * - Ci.4 • e^h 1 -(,2 I i (( + i + 2 • M, +i ) 2 ^ 2 -Mi+i.i •] 2---(a,+i -n )+n I 2 M'+i C,5 • e•" +Ci6 • e~Â4' n • A +i •[Ci + i.3 Ci + i.4 1 A + Ci + i.5 + Ci + i.6 [Ci + i.i + Ci + i.2] — 0 (49) — -Mi.i • {-2 • n • a, • [C,.i • ea'ti - Ci.2 • e a'h ]+(2 + â [C,.3 • e+ Ci.4 • e-^ ] + + -# •[Ci.5 • e^ - Ci.6 • e^ j[ + m,+u ^-2 n • (50) [Ci + i.i - C, + i.2 ]+(( + Âi+i ) • [Ci + ,3 + C +i.4 ] + - . (i+i • [Ci + i.5 - Ci+i.6 ] [ = 0 ACTA SEOTECHNICA SLOVENICA, 2008/l 57. T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LAYERED HALF-SPACE - = au • { -2 • n • a • \Cy • ea - C,,2 • e-"J + ( + 3, fc,3 • e3'h + C -3, •t, + +2•r•A •[c,,5 • e3-Ct,, 0-A t - a+1,1 •{ -2 • r • c,+i,i - C+1,2J + +(r2 + 3+1) [c,+1,3 + c,+i,4 J + 2 • r • 3,+1 • [c,+1,5 - c,+wj}=0 (51) -r • fc, 3 • e3 4> + c, = { a, •Ic,,1 • ea - Ct,2 • e~a 4> -3 •\Ct,5 • e3-Ct,6 • e-3^ H + -3, t -{ a+1 •[Q+1,1-c+1,2J-r•\c,+1,3 + c,+1,4j-3+1 •[Q+1,5-c,+1,6j } =0 (52) = { r •Ic, J • ea ^ + C ,2 • e~a 4' J- 3 \c,3 • e3<> - C,4 • e-3'<> j + r - - • 32 r fc,5 • e3, ^ + C, -3, t Ir •[C;+1,1 + C,+1,2 J 3+1 '[C,+1,3 C;+1,4 J + ,5 1 ¡,6 2 -T r - - • 3+1 r \c,+1,5 + c,+1,6j =0 (53) M,, - M„ {-r • [c,J • ea^ + C,,2 • e~aj + 3, •[c,,3 • e3* - C,,4 • e-3^ j + r fc,,5 • e3* + + C, ,6 • e -3 t - {-r • \c,+1,1 + c,+1,2 J + 3+1 • [c,+1,3- c,+1,4 j+r • C, +1,5 +C, + 1,6 =0 (54) where , e [ 1,2,...,«,« +1 = H] and the equation (44) has to be accounted for. The above created system of 6 • n + 3 equations can be written in the matrix form: A • C = b (55) The matrix A in the above equation is a band matrix with the bandwidth of maximum 9 terms. C is a properly ordered vector of integration constants C{ j. The right side of the system is a vector, where each term except the first one equals to zero. As we have limited our interest to motion at the surface, only the first six integration constants are needed. They are obtained by Cramer's rule. The value of the determinant of the matrix A is determined by its development along the first row. According to the fact, that only first six terms of the first row are different from zero the expression of the determinant |A has the following form: j=6 A = £(-1)j+1 • «1,j • \Aj\ , (56) where a1,j ; j e [1,2,.. .,6] represent the first six non-zero terms offirst rowofmatrix [ A] and | Aj |; j e [1,2,. ,66 represent corresponding sub matrices. Due to the fact that the vector on the right side of the system equations (55) has the following form: b= Q( n • a • kT ; b, = 0 ; j e[2,3,... ,6n + 2J, (57) the six constants, which determine the surface motion of the half-space, can be written as: cu = (-1) j+1 • A • Q( lAl n • A • kT1 j e[ 1,2,.,66. ; k e [1,2,5,6,7,8] (58) j=1 Introducing the equations (32), (33), (35) and (58) into equations (38) to (40) and evaluating the latter at z = 0 lead to the displacements on the surface of the half-space as: 58. ACTA SEOTECHNICA SLOVENICA, 2008/l T. PLIBERŠEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LflYERED HALF-SPACE u (*» = QfL • f J A -[n -(( - a2| -((+|AI) + n - 2 -(n2-i)J( I-A . A J2 (an)- dri Jo (a n)- dn+J A • -n -(a, |-| aJJ + (59) +Vnr-I-(A3 l+l A41)+n-N1-A6 ^^ (a,0) = f- -'|J(jA|)-[n -(A, 1-1A2I .(A3I+1A41) + 2n - i, cT1 n - - -i)|-(A51-1 A |A| J2 (an)- dri J 0 (a n)- dn+J A --n ( H A2 +Vnr-I-(A3 l+l A41)+n-N1-A6 (60) uz (a,0) = -J n - VT^r-(a, I +| A21)-n-((-|A1) + n - ¡1 CT,1 0 lAl 1 ' 1 ; -V^2-! -((+| A61)1- j, (a n)- dn, (6!) where following change of variable has been introduced: , f - r a = kT,, - r = . (62) CT ,1 For further analysis it is worth vile to note that in the above equations only three distinct integrals appear. They are denoted as I1 , I2 and I3 and are given as: I, = JA r ((1 -A2l)-^-(N+1A41) + n—(-1) n (|-|A61)|- J0 (an)- dri = JB, (n)- J0 (an)- dr, (63) I2 = J |A-[-n-((1H A2 N l+l A41)+n-((A l-l A6 - j2 (a n)-dn = J B2 (n)-J (a n)-dn (64) to - J A -(Ail+1A21)-n-(A3 - A41) + i3 = , 3 ' |A| -4^-1 -((1-|A61) - Ji (an)- dri = JB3 (n)- Ji (an)- dn. (65) The newly introduced functions B{ (n), i = 1,2,3 can be identified from the above equations. They also allow us to write the reduced surface displacements in a more compact form. They are given as: TO TO TO TO 0 TO 0 TO 0 ACTA GEOTECHNICA SLOVENICA, 2008/i 69. T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LAYERED HALF-SPACE ur (a,0)= Q( u) 2n ■ ¡ii u (I, +12 Q( ) to to ■ f b, (n) ■ J0 (a n)■ dv + f b2 (v)^ j2 (a n) dn 2n ■ J J (66) u^ (a,0) = Q( u) 2n ■ ¡ii '-T ,j Q(u) ju_ 2n ■ cT, to to f B, (n)J 0 (a n) dn+f B2 (n)■J 2 (a n )■ dn (67) 0 0 _ , Q(u) u T Q(u) u „ , s . / \ , uz (a,0) = — ----13 ---I B3 (n)■ Jj (a n) ■ dn. (68) n ■ ¡1 CT ,j n ■ ¡1 CT ,j J0 3 eVfiLUfiTION OF THE INVeRSION INTEGRALS As it is known from the literature [27] and it was already demonstrated in the previous paper by the authors [19] the leading term of singularity depends on the local conditions only. Due to this fact the singularity at the surface of the layered half-space can be determined by considering the homogeneous half-space with the material properties of the uppermost layer, what considerably simplifies the analysis. The integrals I, i = 1,2,3 equations (63) to (65) are first reduced to the case of a homogeneous half-space and consecutively their limits are evaluated. This yield: lo lo S,1 = ■ lim f a—0 J n n 2Ct,1 S,j 2Ct,1 a—o\ [vn2-! \pH (n) Jo (an)■ dn oj f11 j n- lim T ,1 o u 3 - 2 7 2 cT1 2 — 2 y \fh (n 2 to 2 V fJ o (a n )■ dn = 1 J (an) ■ dn (69) 2r 0 and in an analogous way: 2c, -IS,2 = T,1 2c, T,1 ■ lim f a—0 J jnF—1 \Fh (n) Jo(ar) dn = t1 2r (70) where: \Fh(n)\ = (2■ n2—i)2 — 4■ n2■Jnr—1 Vn2 — y2 . (7i) The limit of the integral I3 as a — 0 is zero and therefore it remains regular for all values of a. The singulari- ties given by equations (69) and (70) are now subtracted from integrals I1 and I2 respectively and taken under the integral sign. This yields: 58. ACTA SEOTECHNICA SLOVENICA, 2008/l T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TANSENTIALY LOADED HORIZONTflLY LAYERED HALF-SPACE I = A - IS,1 = f A {n '0Ail-I A1 -0^31+1A1) + + n--(n2-1) -(N-|a6|)-(2-vj0(an)-dn to = f Bi (n)- J o (a n)- dn (72) I2 = I2 - is,2 = f A-{-n-(n-i a2 d+vn2-!-((+|a4 |) + +n - (As I-| A61)-vx}-j 2 (a n)-dn to = f B2 (n)- J (a n)- dn . (73) The horizontal displacements components are now given by: - / ^ Q(w) 1 W f-T T • ur (a,0) = ^-~- - + — -( +12 2n - j1 r cT,j _ Q(w) fl w ^ 2n - j1 \r cT1 -+—- f Bi (n)- J0 (an)- dn+f B2 (n)- J2 (an)- dn (74) (a,0)= Q(w) 1-v w 2n - j- Q( w) 2n - j1 +--(-I1 +12 1-v W T,1 0 to to -f B1 (n)- J0 (an)- dn+f B2 (n)- J2 (an)- dn [ • (7s) to 0 to 0 0 0 0 0 We note that the singular terms, which appear in the horizontal displacement components, are given explicitly. The values of the integrals I1 and I2 are bounded for all values of the parameter a or r respectively. And what is also important for a numerical evaluation and of course for our later considerations B1 (n), B2 (n) and B3 (n) tend to zero as n ^to . 3.1 EXTENDING THE RANGE OF INTEGRATION To transform the integrals I{, i = 1,2 and I3 into a form permitting their evaluation by contour integration in a complex n -plane we must make their integrands single valued and extend the range of integration from -to to to . Thus we note that the functions B1 (n), B2 (n) and B3 (n) are not single-valued due to the terms ai and fj. appearing in them. They are made single valued by introducing the branch cuts in the complex n -plane. In the selection of a suitable branch cut we are however limited by the following requirements: imposed radiation conditions, which require that the real parts of ai (n) and ffi (n) are positive on the positive real n -axis; that it does not intersect the big semi-circle in the upper n -half-plane and by the demand that ai (n) and ffi (n) are odd functions of n on the real n -axis. The latter is needed to extend the range of integration from semi-infinite to infinite. The branch-cut fulfilling the above stated requirements is shown in Fig. 4. ACTA GEOTECHNICA SLOVENICA, 2008/i 6l. T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LAYERED HALF-SPACE Figure 4. Branch points of expressions ai and ß with introduced branch cut and the corresponding Rayleigh pole. For greater clarity of the figure some material damping has been assumed and expressed by complex shear module ß. = ß • e"p . We realize that the chosen branch-cut indeed makes ai (n) and fj. (n) odd functions of n on the real n -axis and does not intersect the big semi-circle in the upper n -half-plane. We further note that all terms of matrix A, equation (55), except the exponential functions: s \ ±aiti jn2 —7.2 t, p{ (n) = e . . = e ^ ' ' and ±1 \ ±pih ±Jn2 —«f t, qt (n) = e " = e . . (76) which are neither odd nor even. To make the terms in matrix A determined with respect to evenness or oddness respectively we replace them by their analytic continuations, which do not change their values on positive n -axis and are even functions on the real n -axis. As functions satisfying the stated requirements we have chosen: p± (n) = e n = e±a'h and h - q±(n) = e n = e±-ti . (77) The replacement of the functions p± (n) and q± (n) in equations (49)-(54) through p± (n) and q±(n) leads to the matrix A , which has on the integration path of inverse Hankel transform exactly the same values as the original matrix A. The terms of matrix A have interesting and for further analysis very useful properties. It can be seen from equations (46) to (48) that all the terms in the first three rows of matrix A are even functions of n on the real n -axis. The further rows of matrix A are coming from continuity conditions on interfaces between layers, six rows for each interface. From equations (49) to (54) it is easy to recognize that the first tree of these six rows contain only terms, which are even functions of n , and the next three only terms, which are odd functions of n on the real n -axis. Matrix A therefore has 3 ■(n + 1) rows, where all the terms are even, and 3 ■ n rows, where all the terms are odd and n is the number of layers. This implies that all the determinants |A| and |A,|, i s[ 1,2,...,6] are in the case of odd number of layers n odd functions and on the other hand in the case of even number of layers n are even functions on the real n -axis. We now derive functions Bi (n) exactly the same way as we have derived functions Bi (n) only that we make use of determinants |A| and | A, |, i e[1,2,...,6] , instead ofthe determinants |A| and |A,|, i e[ 1,2,.,6] . This yields: 58. ACTA SEOTECHNICA SLOVENICA, 2008/l t. plibersek & a. umek: green's function for tangentialy loaded horizontaly layered half-space B1 (n) = pA-[n-((-| a2 I)—vn2—1 -((+| A< |) + 2 n—(n2—1) -(( — A 1) n J B (n) = jAj-[—n-((—|a2 |)+Vn2—1 -((+| A41) + +n (I—I A61)]— (78) (79) and 4(n) = -[V n2 — y2-((+|a2 I)—n-((—| A41)—Vn2—1 +|A6|)| . (80) A where all the Bi (n), i' €[ 1,2,3] are even functions of n on the real n -axis and Bt (n) = Bt (n), i € [1,2] and B3 (n) = B{ (n) on the real, positive n -axis. Taking into account these equalities, the equations (72), (73) and (65) can be written as: m m 1 = J B1 (n) -/0 (an) - dn = / B1 (n) -Jo (an) - dn = ¡1 (81) 0 0 I2 = / b2 (n)/ (an)- dn = / B2 (n)/ (an)- dn = ¡2 (82) 00 and m m /3 = J B3 (n)- J1 (an)- dn = / B3 (n)- J1 (an)- dn = ¡3.(83) 00 At this point we make use of integral representations of Bessel functions known from the literature e.g. Gradshteyn et al. [28]: /. (an) = h (an) + h (—an); i = 1,2,3 (84) where: 1 n h0 (an) = ^~ Jelavsin' d' (85) h (a n) = 2 H1 (a n) (86) 1 n h2 (a n) = — J e*(ansm22')d' (87) 2 n and H1 (an) is the Hankel's function of the first order and the first kind. Making use of the relationship (84) equation (81) yields: Where we have made use of the change of variables n ^ — n and the symmetry of the function B1 (n). In a very analogous fashion we obtain for the other two pertinent integrals given by equations (82) and (83) following expressions: m m /2 = ¡2 = J B2 (n) -/2 (an)- dn = J B2 (an)- dn = i2 (89) 0 —m m m 13 = JB5(n)-J1 (an)-dn = JB5(n)-h1 (an)-dn = 4.(90) The equations (88) to (90) clearly show that we have successfully replaced the original inverse Hankel transform integrals with the range of integration from 0 to +m with newly defined integrals having the range of integration from —m to +m . The integrals I., i = 1,2,3 can be evaluated by contour integration in the complex n -plane as it will be shown in the next paragraph. 3.2 EVALUATION OF INTEGRALS BY CONTOUR INTEGRATION Integrals I1, I2 and I3, are finally in the form permitting their evaluation by the contour integration in the complex n -plane. The most suitable contour is shown in Fig. 5. By the residue theorem it can be written: I+ + /r +1— + 4 + 4 + = 2 n i ^ res, ; i = 1,2, 3. (91) It can be easily shown that the value of the integral along the big semi-circle in the upper n -half plane is identi- 1 = ¡1 = J B1 (n)-/ 0 (an)- dn = J B1 (n)- h0 (a n)- dn+J ¡¡1 (n)- h0 (—a n)- dn 0 0 0 m —m m = J B1 (n)- h0 (a n)- dn + J B1 (—n)- h0 (a n)- (—dn) = J B1 (n)- h0 (a n)- dn = I1. (88) m 0 0 ACTA SEOTECHNICA SLOVENICA, 2008/l 63. T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LAYERED HALF-SPACE cally zero. This is due to the behaviour of integrands, semi-circle n can bi given as n = R • e"f , where p takes which are dominated by hi (an) functions. On the big the values from 0 to n . Equations (85) and (87) yield: 1 n lim h0(an)= lim — I exp[iR(cosp + isinp)sin§1 d§ = |n|^m Rim n J lm(n)>0 0 (92) 1 n = — lim I exp(-Rsinp + iRcosp)sin§d§ i 0 n Rim o 1 n lim h2(z) = lim — I exp[iR(cosp + isinp)sin§ — 2i§1d§ = '|im Rim 1 ^ = — lim I exp [(—R sinp + iR cos p)sin § — 2i§1 d§ i 0 . n Rim J ^ Inl Im(n)>0 (93) 0 And the same behavior for h1 (an) = H(1) (an 'j, namely that lim H1(1) (R • e'v) i 0 for 0 < p < n is well documented in the literature [29]. Therefore it can be concluded: itR = 0; i = 1,2,3. (94) Taking in account the above equation the equation (91) can be rewritten as: i = i+ +1- = 2 n i £ res,. — ^ — 4 — ^ ; i = 1,2,3. (95) It is clear from Fig. 5 that all three integrals appearing in the right hand term of the above equation have finite integration path. Therefore by the equation (95) our fundamental goal has been achieved. It can be further noted that if the integrals Iih1 and Im are led along one and the other side of the branch cut and the value of the integral I.r is equal zero. For the numerical calculation it is advantageous to express the integrals ih1 and Ih2 through a sum of integrals of even shorter integration range stretching from one singularity on the branch cut to the other. These singularities are either branch points of functions ai and f defined through the equation (45) or the poles defined by the zeros of the determinant |A , which lie on the branch cut. Introducing the equation (95) and considering the equations (88)-( 90) into equations (74), (75) and (68) yields the surface displacements as: 58. ACTA SEOTECHNICA SLOVENICA, 2008/l T. PLIBERŠEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LflYERED HALF-SPACE ur (a,°) = iM- • | i + — ■ [2ni ■ + ^) - IM - I2bl - Ilb2 - I2b211 2n ■ Mi [r ct ,1 1 JJ Uf) («,0) = ■ i ^ + — ■ [2 n i ■^(-reS1 + reS2 ) + I1bi - I2bi + I1b2 - I2b2 1 [ 2n ■ M1 [ r CT ,1 1 JJ (96) (97) UZ («>0) = ^^ ^ — (2^f^reS3 - I3b1 - Jsi 2 n ■ cT ,1 v (98) The above three equations are the final result of our analysis. In them the surface displacements of a layered half-space due to a tangential point load are given by a singular term, in the components, which become singular as r ^ 0 , a sum of the residues and a sum of integrals with the finite integration path. The latter ones can be evaluated numerically with any desired accuracy. We can further note that the residues represent the surface waves and the integrals are due to body waves. 4 NUMERICAL EXAMPLE As an illustrative example a one-layer half-space with the same geometrical and material characteristics as the one considered by Štrukelj et al. [19] has been chosen. It is shown in Fig. 1 with n = 1. On its surface a horizontal, harmonic point-load is applied. The material properties of this half-space are as follows: ratio of material densities in the underlying half-space and the layer is phIP1 = 1.5. the ratio of shear modules ^h/= 2.0 , the Poisson's ratio for the layer is v1 = 1/3 and the Poisson's ratio for underlying half-space is vH = 1/4. The materials in the layer and in the underlying half-space are assumed to have no material damping. The ratio of the layer thickness h1 and the wave length of shear waves in the layer A1 is taken to be h^A1 = 2.0 . As it can be seen from equations (96) to (98), the displacement components are expressed as sums of several terms. It is worthwhile to note that in the both horizontal displacements components the same terms appear with exception of the singular term. In the vertical displacement, however, completely different terms are forthcoming. The numerical effort can be, therefore, considerably reduced by computing first each of these terms separately and later combine them into displacement components as given by equations (96) to (98). In the Fig. 6 and 7 as an example three such characteristic terms are given. Our choice of the geometrical and material properties of the half-space is based on the fact that it is nearly impossible to obtain the data in the pertinent literature, with which our results could be compared to prove their valid- /MM>(<'■")! 112 0.1 s II L -tMW -u.i + g Figure 6. Real parts of terms 2n i^^ res3 and I, ACTA GEOTECHNICA SLOVENICA, 2008/i 69. T. PLIBERŠEK & A. UMEK: GREEN'S FUNCTION FOR TANGENTIALY LOADED HORIZONTALY LAYERED HALF-SPACE /m («.£>)) (1.05 ■0.05 -0.15 -0.2 > a Figure 7. Imaginary parts of terms 2n i ^^ res3 and 7 ity and accuracy. Therefore we have decided to make use of the principle of reciprocity in elastodynamics, which is well described in the literature e.g. by Achenbach [30]. It can be concluded from this principle, that the vertical displacements along the ray § = 0 in a half-space loaded with a horizontal, harmonic, unit point-load are equal to the radial displacements in an identical half-space loaded with a vertical, harmonic and unit point-load. Therefore we can state that for our above described choice of the layered half-space the displacement component uz (a,0), as given by equation (98), should equal the radial component of the surface displacements presented by Štrukelj et al. [19]. It is however worth vile to note that this equivalence can not be seen from the two expressions before their numerical evaluation. The integrals in equation (98) have integrands, which are based on the determinant and sub-determinants of a 9x9 matrix A. The corresponding integrals are in the case of the vertically loaded half-space based on a 6x6 matrix. The real and imaginary parts of both displacement functions are shown in Fig. 8 and 9. It can be seen from both figures that they are in an excellent agreement with the results presented by Štrukelj et al [19]. The results of the evaluation of the radial displacement component ur (a,0) are shown in the Fig. 10 and 11. In exactly the same way the circumferential displacement component u § (a,0) can be evaluated through a different combination of terms appearing in ur (a,0). The results of this evaluation will not be presented in this paper. Figure 8. The real part of the displacement function uz (a,0) given by the dotted line and the real part of the function ur (a,0) presented through the solid line. 66. ACTA GEOTECHNICA SLOVENICA, 2008/l T. PLIBERŠEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LflYERED HALF-SPACE Figure 9. The imaginary part of the displacement function uz (a,0) given by the dotted line and the real part of the function ur (a,0) presented through the solid line. Re(u,(a,0)) Figure 10. The real part of the displacement function u (a,0) . ^ rUfo.OH Figure 11. The imaginary part of the displacement function u (a,0) . ACTA GEOTECHNICA SLOVENICA, 2008/i 69. T. PLIBERSEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LAYERED HALF-SPACE 5 CONCLUSION REFERENCES The displacement components on the surface of a horizontally layered half-space due to a tangential point-load were expressed through a combination of tree distinct Hankel's inverse integrals and trigonometric functions with the circumferential coordinate as argument. For the evaluation of these Hankel's integrals a novel three step procedure is employed. In the first step the singularity at the generic point from the integrals, where it exists is extracted and the resulting new integrals are made regular. In the second step we replaced the new integrands functions with their suitable analytic continuations, by which we were able to extend the integration range of Hankel's integrals to —< to +<. By this extension of the integration range we were in the last step permitted to evaluate them by contour integration. Through these three steps we were able to transform the Hankel's integrals into sum of three terms. The first one contains the singularity in the form C/r , the second one is given by a sum of the residues of the integrand and finally the third term consists of finite number of integrals along the suitable chosen branch-cut. The latter ones regular and finite in their integration range can be easily evaluated numerically. The results presented in this paper together with the our previous results, Štrukelj et al. [19], constitute, what we believe, a robust and numerically efficient method to evaluate the displacements on the surface of the horizontally layered half-space due to a point force of any direction. The method of evaluation presented in this paper provides us with exact and closed form expressions for the singularities of displacement field, what makes our results very suitable to be used in soil-structure interaction problems. We are convinced that an even more efficient integration line around the branch cut from the one used in this paper can be developed. Before we could come up with a definite recommendation concerning the integration path more numerical research is needed. It is however believed that this problem is beyond the scope of this paper, where we succeeded to demonstrate that the Green's function for a layered half-space can be expressed as a combination of terms, which can be easily, especially in comparison with original Hankel's inversion integrals, and accurately evaluated. [1] Premrov M, Spacapan I. (2004). Solving exterior problems of wave propagation based on an iterative variation of local DtN operators. Appl. Math. Model, 28, 3, 291-305. [2] Kelvin (Lord) (1843). Displacements due to a point load in an indefinitely extended solid. Citation no. 66 in Love's Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover Press, New York,. [3] Boussinesq J. (1878-1885). Citation no. 67 in Love's Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover Press, New York. [4] Cerruti V. (1882). Citation no. 68 in Love's Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover Press, New York. [5] Lamb H. (1904). On the propagation of tremors over the surface of an elastic solid. Phil. Trans. Royal Soc. Lond.; A 203, 1-42. [6] Midlin RD. (1936). Force at a point in the interior of a semi-infinite solid. J. Appl. Phys., 7, 195-202. [7] Ewing MW, Jardetzky WS, Press F. (1957). Elastic waves in layered media. McGraw-Hill Book Company New York. [8] Apsel J, Luco J.E. (1983). On the Green's functions for a layered half-space, Part I. Bull. Seismol. Soc. Am., 73(4), 909-929. [9] Apsel J, Luco J.E. (1983). On the Green's functions for a layered half-space, Part II. Bull. Seismol. Soc. Am., 73(4), 931-951. [10] Kausel E, Peek R. (1984). Dynamic loads in the interior of a layered stratum, An explicit solution. Bull. Seismol. Soc. Am., 74, 1459-1481. [11] Gaitanaros AP, Karabalis DL. (1988). Dynamic analysis of 3D flexible embedded foundations by a frequency domain BEM-FEM. Earthquake Engng. Struct. Dyn., 16, 653-674. [12] Triantafyllidis T. (1991). 3D time domain BEM using half-space Green's functions. Engng. Anal. Bound. Elem., 8 (3), 115-124. [13] Auersch L. (1994). Wave propagation in layered soils, Theoretical solution in wave number domain and experimental results of hammer and railway traffic excitation. J. Sound Vibration, 173(2), 233264. [14] Aubry D., Clouteau D. (1991). A regularized boundary element method for stratified media. In: Cohen G., Halpen L., Joly P., editors: Proceedings of the First International Conference on Mathematical and Numerical Aspects of Wave Propagation Phenomena, 660-668. [15] Vostroukhov A.V., Verichev S.N. (2004). Kok A.W.M., Esveld C., Steady-state response of a stratified half-space subjected to a horizontal arbitrary 58. ACTA SEOTECHNICA SLOVENICA, 2008/l T. PLIBERŠEK & A. UMEK: GREEN'S FUNCTION FOR TflNGENTIflLY LOADED HORIZONTALY LflYERED HALF-SPACE buried uniform load applied at a circular area. Soil Dyna. Earth. Engng., 24, 449-459. [16] Jin B., Lia H. (1999). Exact solution for horizontal displacement at center of surface of an elastic halfspace under horizontal impulsive punch loading. Soil Dyna. Earth. Engng., 18, 495-498. [17] Bringham E.O. (1974). The Fast Fourier Transform. Prentice - Hall, New York. [18] Kobayashi T, Sasaki F. (1991). Evaluation of Green's function on semi-infinite elastic medium. KICT Report No. 86, Kajima Technical Research Institute, Kajima Corporation. [19] Štrukelj A., Pliberšek T., Umek A. (2006). Evaluation of Green's function for vertical point load excitation applied to the surface of a layered semiinfinite elastic medium. Arch. Appl. Mech., 76, 465-479. [20] Pliberšek T., Štrukelj A., Umek A. (2005). Green's function for an elastic layer loaded harmonically on its surface. Acta Geotechnica Slovenica, 2, 4-21. [21] Žlender B., Lenart S. (2005). Cyclic liquefaction potential of lacustrine carbonate silt from Julian alps. Acta Geotechnica Slovenica, 1, 22-31. [22] Žlender B., Trauner L. (2007). The dynamic properties of the snail soil from the Ljubljana marsh. Acta Geotechnica Slovenica, 2, 48-61. [23] Kobayashi T. (1981). Evaluation of response to point load excitation on semi-infinite elastic medium (in Japanese). Tran. Archi. Institute., Japan, 302, 29-35. [24] Achenbach J.D. (1973). Wave propagation in elastic solids. North-Holland, Amsterdam, London. [25] Miklowitz J. (1980). The theory of elastic waves and waveguides. North-Holland, Amsterdam, New York, Oxford. [26] Sneddon I.H. (1972). The use of integral transforms. McGraw-Hill Book Company, New York, Toronto. [27] Pao Y.H., Mow C.C. (1971). Diffraction of elastic waves and dynamic stress concentrations. Crane Russak, New York. [28] Gradshteyn I.S., Ryzhik, I.M. (1994). Table of Integrals, Series, and Products, Academic Press Limited, London. [29] Abramowitz M., Stegun I.A. (1972). Handbook of mathematical functions. New York, Dover publications. [30] Achenbach J.D. (2003). Reciprocity in Elastody-namics, Cambridge University Press, Cambridge. ACTA GEOTECHNICA SLOVENICA, 2008/i 69.