Boundary element approach to the dynamic stiffness functions of circular foundations

BOUNDARY ELEMENT APPROACH TO THE DYNAMIC STIFFNESS FUNCTIONS OF CIRCULAR FOUNDATIONS E. ALARCON AND J. J. CANO AND J. DOMINGUEZ INTRODUCTION In the past twenty years numerous investigations which deal with the computation of dynamic stiffness (also called impedance) matrices of massless rigid foundations with different shapes have been reported. Taking advantage of the simplifications derived from axial symmetry, circular and cylindrical foundations have been studied in many cases. Analytical solutions for stiffness functions of circular foundations on a uniform half-space were presented for horizontal and rocking motions by Veletsos and Wei1 and for vertical and torsional motions by Luco and Westman.2 The development of absorbing boundaries permitted the use of the finite element method to study the dynamic response of axisymmetric foundations on soils that consist of one or several layers based on a rigid rock.3"5 An alternative approach to the dynamic analysis of rigid foundations is the use of boundary integral equations which allow for the modelling of many different geometries and soil profiles. Both the direct and the indirect integral equation formulations have been used in dynamic soil-structure interaction problems.6"8 Dynamic stiffness functions of rigid circular foundations on uniform or layered soils have been computed by Apsel9 and Apsel and Luco7 using the indirect formulation in combination with Green's function of a ring load in a uniform or layered halfspace, derived by the same authors.9"11 This Green's function is written in terms of integral forms that include products of two Bessel functions. Owing to that, it has to be evaluated segmentally and its computation becomes rather involved. An auxiliar boundary where the unknown sources are located has to be defined.7'9,12 The boundary element method, based on the direct integral equation formulation, has been used to compute dynamic stiffness functions of rectangular6 and strip foundations in the frequency domain13,14 as well as in the time domain.15 In the present paper the boundary element method is formulated for time harmonic axisymmetric problems using the full-space point load fundamental solution. The approach is applied to the computation of the dynamic stiffness functions of rigid circular foundations on layered viscoelastic soils. The formulation of the BEM for axisymmetric elastostatic problems was presented by Kermanidis16 and Cruse et al11 They used the ring load fundamental solution. However, a closed form expression for the ring load does not exist in elastodynamics18 and when the same formulation is followed the computations become involved because of the difficulties in the evaluation of the fundamental solution at each integration point. On the contrary, the point load fundamental solution is simple and may be integrated numerically along the azimuthal coordinate without any particular difficulty. As an example of the range of sites to which the proposed approach can be applied, three different sites are considered; a uniform half-space, a soil layer on a half-space, and a soil consisting of four horizontal layers and a compliant half-space. The numerical results computed by the proposed approach are compared with results obtained by different procedures. BOUNDARY ELEMENT FORMULATION FOR ELASTODYNAMIC AXISYMMETRIC PROBLEMS Following the idea of Chapel,19 the full-space point load fundamental solution is used to solve time harmonic axisymmetric problemsJ The axisymmetric representation of the geometry and field variables is maintained. The boundary of the generating surface of the body is discretized into line elements and the point load collocated at each node. The fundamental solution, in terms of cylindrical co-ordinates, is integrated along the boundary elements of the generating surface and along the azimuthal co-ordinate. The basic BEM equation for zero body forces can be written in Cartesian co-ordinates as usual: CV+ Tudr= UtdT (1) where u and t stand for the displacements and traction vectors, respectively; U and T are, respectively, the matrices representing displacements and tractions produced by the point load applied at point / and C is the independent coefficient matrix such that Ckl = (1/2)<5W when the surface is smooth at the point i. The relation between Cartesian and cylindrical co-ordinates may be written for vectors u and t at a certain pointy as u = Qu (2) t=Qt cylindrical (Figure 1) and the transformation Q COS& sind 0 sind 0 cos 0 0 0 1 (3) By substitution of equation (2) into equation (1) and premultiplying by Q,vr the following equation is obtained: QiT C Q' u^ + J (Q*'-T TQ)ucdr (Q*'-TUQ)tcdr (4) r where QJ is the transformation matrix for the collocation point i. Equation (4) can now be written as C^+ jrTc ucdr uc tc dr (5) r which is the same as equation (1) written in cylindrical co-ordinates. [ The 9 co-ordinate of the collocation point j is fixed during the integration process and can be set equal to any desired value; for instance $' = 0. In such a case T i.T Q UTQ rncos0+ 7\2sin$ !T21cos9 + r22sin9 r31cosd+r32sin0 rnsind+r12cos3 Tx r21sin»+r22cos» T r31sin9+r32cos9 T 3 23 33 (6) where the terms Tkl are the well-known Cartesian co-ordinate components of the tractions in the / direction at the integration point, due to a unit point load in the k direction at the point i. The same equation can be written for Uc. x i x 2 Figure 1. Description of co-ordinate systems When a circular foundation under vertical or torsional motion is considered, the field variables are axisymmetric. Assuming that the boundary is smooth at i, equation (5) becomes 1 2 u i p 3 t4 + 2 P r r^cS+r^sd o r31cd + r32s0 o T2ls$+T22c& 0 7-13 0 ^33 d& up "3 "z dr 2 P l/nc&+t/12s3 0 U3! cd + U32 sd 0 U21sa+L/22c9 0 l>13 0 t/33 dO ', '3 t2 dr (7) where T* is the boundary of the generating surface (Figure 2), s and c stand for sin and cos, respectively and a zero has been placed for the skewsymmetric terms, where the integrals around the azimuthal axis are null. It should be mentioned that the torsional and vertical radial motions are uncoupled in equation (7) and that the kernels to be integrated along d are known functions. When a circular foundation is under horizontal motion along Xx or rocking around X2, the field variables are of the form (Figure 3) u u 3 ulpcos9 ; "id sin » ; t rlpcosd t 3 tld sin $ (8) u2 = ulz cos 3 ; tz = tXz cos $ where ulp, ulv, ulr, rlp, tx</> and fj2 are independent of 9. Instead of equation (7) the following P (a) (b) Figure 2. Generating surface of an axisymmetric domain: (a) genera] problem; (b) surface foundation on a half-space X 2 X 1 X 3 X 2 Horizontal X 1 X 1 Rocking Figure 3. Horizontal and rocking displacement equation is obtained: 1^ 2 u IP 0 u it + 2 P 2 rnc2d + r12s&cd rns2o-r12sdcd o 0 2 r31c29 + r32s9c9 r31s2d- r32sdcd r, 3 c» 0 r33cd dd "IP «1S "l* dr 2 P 2 l/n^S+l/jaSacd (/nS29-l/12sdc9 l/i 3 0 0 0 2 l/31c2d+l/32s9c0 l/31s29-U32sdcd U 33 c& c9 dd h. h* *i. dr (9) Equation (9) can be used to represent the radial and vertical displacements but not the azimuthal displacement for which an identity would be obtained. The amplitude of a function such as ud = ux a sin 9 cannot be computed with a collocation point where the function is zero. Because of this, a different collocation point where the shape function —sin 3' has a unit amplitude (»' 7t/2) is selected for the azimuthal equation: T o ,«.T Q'-'TQ T°2lc$-T°22sB r^ca+r^sa r?,cs r?,s9+r?2c» T°3is$+T°32c$ n 23 7^3 (10) where the superindex '0' indicates that the collocation point is at 8' represents the three components is 7t/2. The equation that 1 2 id if «u «i. + 2 P r* J (7ilC2s+r12 s9c9)dS 72 (rj^s+rjjsScSjdd -«/2 2 P r* I (7-3, c29+T32s9c9)d9 c29+l/I2s9c9)d9 72 (C/?1c2»+{/o2s8c8)d8 -«/2 ["(U" c29+l/32s9c9)d9 J> s23-r12s&cS)d» «/2 (r?jsa»-r?288o»)d» -«/2 Jo s29-7"32s9c9)d9 s29-l/,2s9c9)d9 72 0 „2 (l/?,s29-l/?2s9c9)d9 -«/2 T If*"" s2»-l/32$3c9)d3 c9)d3 :/2 (r?3c»)ds -«/2 (733cS)d» u ip «i» u Xx dr r ^ c9)d9 n (l/?3c»)dS -*/2 (l/33c»)d» r IP 'ia 'i. dr (11) In general, non-axisymmetric boundary conditions can be analyzed using a plane model by representing the field variables by a Fourier series along the azimuthal co-ordinate. The series is of the form 00 U P (u*np cos n§ + ulp sin nS) « = o 00 uB = £ (-149 sin n9 + u^cos nO) (12) n = 0 ao U (u*zcos nO + wJJr sin nO) Fl = 0 terms and a boundary equation may be written for the amplitude of the symmetric and the antisymmetric parts of each mode. The equations for the amplitude of the terms of the form sin n3 are obtained using a point at 0' fin as the collocation point except element that contains the collocation point. In the next section, an integration procedure for constant elements contained in planes x3 = constant is presented. This kind of element is the only one needed for the analysis of circular foundations on uniform or layered soils. INTEGRATION OVER THE BOUNDARY ELEMENTS Constant boundary elements with one node per element are considered (Figure 4). The boundaries of the circular foundation problems analyzed are perpendicular to the X3-axis and the line boundary co-ordinate T* coincides with the radial cylindrical co-ordinate p. Equations (7) and (11) can be written after discretization as 4 N «£ + j-i p Pi [JT.d»]d,>. USAM*™ (13) t Massless foundation 1 In tegrat ion point Collocat ion • • Collocation point Figure 4. Boundary element discretization where N is the number of boundary elements, pi indicates the segment that forms the element j9 and Tc and Uc stand for the matrices in equations (7) or (11). The double integrals extend over one half of the circular crowns represented by the line elements (Figure 4). The collocation point for equation (13) can be either at 9l =0 or at & n12 and for each term the integration domain is either the element to which the collocation point belongs (i=y) or a different one (i #7). The latter case is analysed first. The integrals are done numerically using a Gaussian quadrature formula after a special co-ordinate transformation. The integration domain is shown in Figures 5(a) and 5(b), the latter being for the case when the collocation point is at S1 = — it/2. In order to apply the Gaussian quadrature the domain is transformed into a square in the dimensionless co-ordinates (— 1 ^ £ <; 1, — 1 ^ */ < 1). The radius p is transformed into r\ by 9 1 2 l(Rx+R2) + ri(R2 1 dp = ^(R2-Rx)dti *i)] (14) + A linear transformation of the same kind was tested for the co-ordinate 0; however, the accuracy of the computed values of the integrals was poor even when a large number of integration points was used. The accuracy was improved, and satisfactory results were obtained, by using a quadratic transformation formula that increases the number of integration points in the vicinity of the X j-axis or the X2-axis in the cases of Figures 5(a) or 5(b), respectively.20 The following (a) X 1 (c) * 5 (b) Figure 5. Integration domains and coordinates when the collocation point is outside the integration domain: (a) 9' = 0, (b) $' = - 7t/2 parabolic transformation has been applied; n a =^<i + 02 n (15) dd=-(l + £)d£ when the collocation point is at $' = 0 (Figure 5(a)) and 4 n (16) d9=-(l+£)d£ when the collocation point is at 01 7c/2 (Figure 5(b)). In this way one half of the total number of integration points are located in the quarter of the domain closer to the collocation point. \ The number of quadrature points is two for the radial co-ordinate and twenty for the azimuthal co-ordinate. Thus, the integrals in equation (13) become H*u P pi 2 20 dSJdp- X Z(7T/4)[(/?1+/?2) + nm(R 2 = 1 n=l Rl)1Tc(tn.'lm)lR2-RllV+U">*Wm 2 20 G'j = | p Ucd» dp- £ ZinMti^ + Rt) Pj L JS J m~lm*l + ^/?2-^i)]Uc(^,0[«2-^i][l+^]>vll>vm (17) where Tc(^, fjj and UJ&,, r\m) represent the values of the matrices in equations (7) and (11) computed at the (n, m) integration point and wn, wm are the weights of the Gaussian quadrature formula. When the integration is carried out over the element which contains the collocation point, a singularity exists in the fundamental solution and the integration scheme is different. First of all, the part containing the singularity is separated from the rest of the fundamental solution to be integrated. This is done by subtracting the static fundamental solution from the dynamic one. The difference is non-singular at any point of the integration domain and can be integrated by means of the same kind of numerical quadrature formula, combined with the quadratic transformation of the co-ordinate $, as before. The coefficients of equation (13) for which i==/ can be written as where Hft" and G|J are the static counterparts of H*" and G", respectively, and H&'= f p{ f rTc«Tc(static)jddjdp= f pj [Tc,dinda|dp (18) 2 20 X (rt/4)l(Ri + R2) + rim(R2 ~ *i)]TCldtaK„ rjm) = 1 n=l [«2-«l3[l+^]wnW G&= P< [Uc-Uc(static)]d3Vdp = pi Uc,djnd9Vdp Pi IJB 2 20 l.(n/4fmRl+-R2) + fim(R2-Rln = 1 11=1 Uc,din(^,^)[/?2-K1][l+^]wnw (19) where TCidin(<^„, r\m) and UK<6in(£H, rjm) stand for the values at £„, rjm of matrices of the same form of those in equations (7) and (11) but with the difference between the dynamic and the static fundamental solution instead of the dynamic fundamental solution (Tkldin = Tkl — !rt, (static) and *Aj.din = Uki — (A/(static) instead of Tkl and Ukl, respectively). Figure 6 shows the integration domains and the collocation points for this case. The process is completed with the integration of the static fundamental solution. The coefficients Tu (static) and Ukl (static) are known in terms of the distance r and its derivatives (Figure 6): 2 . D2 r = J{p2 + R2 - 2 Rpcos 3) R Rt + R2 2 /-, = (/? — pcos$)/r r2 = (psin9)/r (20) Figure 13. Circular foundation on a multilayered soil on a half-space. Boundary element discretization CONCLUSIONS Dynamic stiffness functions of rigid circular foundations on a uniform or layered viscoelastic soil have been computed by a boundary element approach that makes use of the complete space point load fundamental solution. This solution is simple and easy to integrate over the boundary elements. The price paid for using such a simple fundamental solution is that the soil free surface and the layer interfaces have to be discretized. However, the discretization of the soil surface and the layer interfaces can be reduced to a limited zone sorrounding the foundation and is automatically done by the code. Approaches based on Green's functions for a layered half-space7,9-11 have the advantage of requiring only the discretization of the soil foundation interface but have the disadvantage of dealing with a fundamental solution that is not known in closed form but in terms of integral forms that include products of two Bessel functions and have to be evaluated segmentally by elaborated numerical integration procedures. The proposed approach is easier to implement than those other 1.0 k hh 0.5 0.0 1.0 hh 0.5 0.0 V r 1 ^—J p^4 i.d krr 0.5 0.0 | • j >s>^ | "«».• ^ .2 .4 .6 .8 1 1.2 1.4 .2 .4 .6 .81.1.21.4 1 .0 k vv 0.5 • 0.0 1 .5 c vv 1.0 0.5 .2 .4 .6 .8 1. 1.2 1.4 Luco B.E.M Figure 14. Stiffness coefficients for circular foundation on a multilayered soil on a half space approaches. The amount of computer time required to solve a particular problem is smaller for one case or the other depending on the number of layers in the soil profile. To illustrate the capabilities of the proposed approach, dynamic stiffness functions for a uniform half-space and two layered soils have been obtained. The results obtained for the halfspace have been validated by comparison with analytical solutions; those for the layered soils by comparison with the results of an approach based on numerical solution of sets of integral equations and, in the case of a single soil layer on a half-space, also by comparison with other boundary element results. The comparisons indicate a good degree of accuracy. ACKNO WLEDGEM ENTS The work described herein is partially supported by a grant provided by the local government of Andalusia, Spain. The authors express their sincere appreciation for this support. APPENDIX Integrals of the static part of the fundamental solution over the element where the point load is located, are given in this Appendix. Terms corresponding to vertical or torsional motion of the foundation (equation (7)): J Pi Jo 2/4 +^3-cos0+—cosa-^3-(l+cos29) pdOdp 2/4 7,cos9d8+ /.cos&d& + /?2 J,cos&d9 VI- R j /3(l+cos2d)d3 2 0 (24) where A = (l/16);rp(l - v) and B = 3 - 4v. G,2,2c,s,= 2| I l/22c>$lpd3dp = 2,4j J !-_+£_ JpdSdp 2/4 B /jCOsSdd + K /3sin2Odd (25) G^C.«=2J f,t/33c.,,pdddp=2/l [ [*p<&dp = 2AB\ 1,6$ (26) P i J O J Pi G il _ fi%i _ fiu _ fin __ fiu _ fin _ (\ 12c,st— °13c,st — u23c,st — u21c,«t — u31c,st — u32c,*t — u (27) «?!*..«-2 J j r13c>4tpd3dp = 2^J fr^(pcosd-/?)Jpdadp 2AB\ | /3cos9d9-.R |/2d9 (28) #5iU=2J fr31c,$tpdddp = 2,4 J rr-^(p-Rcosd)"|pd»dp 2^B| I Add-/? fYcosOdS ,• ] £7*ii _ H*" — *7*" — 77*" — 77*" —£7*" _ 77*" 13 llcst — /722c,st — "33c,st — A 12c,st ~ "21c,st — "23c,$t — "32c,st Terms corresponding to horizontal or rocking motion (equation (11)): H- G?M = 2A\B\ /1cos2»d9+ /4cos2»dd 0 JO + R2\ J2cos23-R J3(l+cos29)cos3d9 Gf2Citl = 2i4 B \Ixsm2&d$ + R2 f/2sin2dd»-/? | /3sin29cos9dd L Jo Jo Jo G 22c, st — 2-4 j^fl j"\ cos2 a da + 7* f 73sin2acosada G*2!ic,st = 2i4 73 71sin2ada+ 74sin2ada —7? 73 sin2 a cosa da G3''3c,st = 2/IB 7lCos»dd u13c,sl — u23c,st — u31c,st — °32c.st — u tf?3c..t = 2i4£ /3cos2&d9-K /2cos»d9 0 JO /!&,„ = 2,4 J? /3 sin2 Odd J o Hfi^=-2AB\ ^cosOdO-K J2cos20d3 /*« #*"c.st = 2/4B/? 72sin29d9 77*ii 77*" 77*" 77*" 77*" A 11 licst — /322c,si — /733c,st "** " 12c.st — n 21c,st — u