BEM-FEM coupling model for the dynamic analysis of piles and pile groups

BEM-FEM coupling model for the dynamic analysis of piles and pile groups ∗ L.A. Padr´on, J.J. Azn´arez, O.Maeso Instituto Universitario de Sistemas Inteligentes y Aplicaciones Num´ericas en Ingenier´ıa (SIANI) Universidad de Las Palmas de Gran Canaria Edificio Central del Parque Cient´ıfico y Tecnol´ogico Campus Universitario de Tafira, 35017, Las Palmas de Gran Canaria, Spain {lpadron,jjaznarez,omaeso}@siani.es 1 september 2006 Abstract This paper shows a BEM-FEM coupling model for the time harmonic dynamic analysis of piles and pile groups embedded in an elastic half-space. Piles are modelled using Finite Elements (FEM) as a beam according to the Bernoulli hypothesis, while the soil is modelled using Boundary Elements (BEM) as a continuum, semi-infinite, isotropic, homogeneous or zoned homogeneous, linear, viscoelastic medium. It is assumed that the soil continuity is not altered by the presence of the piles, and the tractions at the pile-soil interface are considered as a load applied within the half-space. The formulation is exposed in detail. In order to validate the model, selected numerical results of time harmonic impedances of different pile groups configurations are evaluated and contrasted with other reference values taken from the literature. 1 Introduction The frequency domain dynamic analysis of piles and pile groups embedded in a half-space has been treated numerically by several authors using different boundary integral formulations for the soil, and FEM (Finite Elements Method) for the piles, considered as monodimensional beam elements [1, 2, 3, 4, 5, 6, 7, 8, 9]. Also, the coupling of framed structures approached by FEM with threedimensional bodies represented by BEM (Boundary Elements Method) in time domain has been presented in [10, 11], where piles would be approximated using a special cylindrical boundary element. Using BEM for both soil and piles, more versatile and rigorous numerical models have been developed: vibration isolation by a row of piles has been analyzed in [12, 13] and dynamic impedances of pile groups has been studied by two of the authors in [14, 15]. High computational cost is a disadvantage of these models. The present paper, in order to reduce the number of degrees of freedom in the problem, presents a different BEM-FEM model for the time harmonic dynamic analysis of piles and pile groups embedded in an elastic half-space, taking advantage of the particular characteristics of each one of the methods. This model is based on the idea of a previous static model developed in [16, 17, 18] ∗This is the peer reviewed version of the following article: L.A. Padr´on, J.J. Azn´arez, O.Maeso, BEM-FEM coupling model for the dynamic analysis of piles and pile groups, Engineering Analysis with Boundary Elements 31 (2007) 473484, which has been published in final form at doi:10.1016/j.enganabound.2006.11.001. This work is released with a Creative Commons Attribution Non-Commercial No Derivatives License. 1 where it is assumed that the continuity of the soil is not altered by the presence of the piles and where the tractions in the pile-soil interface are considered as a load applied within the half-space in the boundary integral representation of the soil. Piles are modelled using FEM as beams according to the Bernoulli hypothesis, while soil is modelled using BEM as a continuum, semi-infinite, isotropic, homogeneous or zoned homogeneous, linear, viscoelastic medium. Welded boundary contact condition at the pile-soil interfaces are assumed. The formulation allows the analysis of problems including soil strata, rigid rocky beds and any topography for the soil surface. Furthermore, as the pile boundary does not need to be discretized, low computing times and memory requirements are needed. First, the formulation will be exposed in detail. After, in order to validate the model, numerical results for vertical, horizontal and rocking impedances are evaluated and compared to other reference values taken from the literature. 2 Pile FE equations The behaviour of a pile submitted to dynamic loads can be described by the following differential equation M ¨u(t) + C ˙u(t) + K u(t) = f(t) (1) where M,Cand Kare the mass, damping and stiffness matrices of the pile, u(t) is the vector of nodal displacements and f(t) the vector of nodal forces over the pile. It will be assumed now that the pile is subjected to a harmonically varying load. In this case, the vectors of nodal displacements and forces can be expressed as u(t) = upeiωt ;f(t) = Feiωt (2) where upis the vector of nodal translations and rotations amplitudes, Fis the vector of nodal forces amplitudes, ωthe circular frequency of the excitation and i = √−1. Then, and considering a pile with zero internal damping, Eq.(1) becomes (K−ω2M)up=F(3) Piles are modelled by FEM as vertical beams according to the Bernoulli hypothesis, and are discretized using a three-nodes element, shown in Fig. 1, that has been defined in order to be able to adjust the deformed shape accurately with a scarce number of elements. There are 13 degrees of freedom defined on it: two lateral displacements and a vertical displacement on each node, and two rotations θon each one of the extreme nodes, one about x1axis and another one about x2. The lateral displacements u1and u2along the element are approximated by a set of fourth degree shape functions, while vertical displacements u3are approximated by one of second degree. Thus ui=ϕ1uki+ϕ2θki+ϕ3uli+ϕ4umi+ϕ5θmi;i= 1,2 (4) u3=φ1uk3+φ2ul3+φ3um3(5) where 2 Figure 1: Finite element definition ϕ1=ξ(−3 4+ξ+1 4ξ2−1 2ξ3) ϕ2=1 4ξ(−1 + ξ+ξ2−ξ3) ϕ3= 1 −2ξ2+ξ4(6) ϕ4=ξ(3 4+ξ−1 4ξ2−1 2ξ3) ϕ5=1 4ξ(−1−ξ+ξ2+ξ3) and φ1=1 2ξ(ξ−1) φ2= 1 −ξ2(7) φ3=1 2ξ(ξ+ 1) being ξthe elemental dimensionless coordinate varying from −1 to +1. Using the principle of virtual displacements and the shape functions defined above, the stiffness sub-matrix for the lateral behaviour of this element (denoted by super-index l) can be obtained as (see e.g. reference [19]) kl ij =ZL ϕ′′ iEIϕ′′ jdx3;i, j = 1, ..., 5 (8) and the one for the axial behaviour (denoted by super-index a) as ka ij =ZL φ′ iEAφ′ jdx3;i, j = 1,2,3 (9) 3 where Eis the Young’s Modulus for the pile, Aand Iare the area and the moment of inertia of the section of the pile and Lis the element length. Here, super-indexes ′denote derivative referring to x3coordinate. Finally, the matrices obtained are           fki mki fli fmi mmi           =EI 5L           316 L2 94 L −512 L2 196 L2 −34 L 94 L36 −128 L 34 L−6 −512 L2 −128 L 1024 L2 −512 L2 128 L 196 L2 34 L −512 L2 316 L2 −94 L −34 L−6128 L −94 L36                     uki θki uli umi θmi           ;i= 1,2 (10) and     fk3 fl3 fm3     =EA 3L    7−8 1 −8 16 −8 1−8 7        uk3 ul3 um3     (11) Similarly, the mass influence coefficients for an element, that represents the inertia force opposing the acceleration experimented by a certain degree of freedom, can be evaluated by a similar procedure as ml ij =ZL ϕi¯mϕjdx3;ma ij =ZL φi¯mφjdx3(12) Using the same functions that were used for calculating the stiffness matrix, the result obtained is the consistent-mass matrix. Thus, considering a beam with uniformly distributed mass ¯m, the matrices obtained for the lateral and axial behaviours are, respectively Ml=L¯m           13 63 L 63 4 63 −23 630 L 180 L 63 L2 630 2L 315 −L 180 L2 1260 4 63 2L 315 128 315 4 63 −2L 315 −23 630 −L 180 4 63 13 63 −L 63 L 180 L2 1260 −2L 315 −L 63 L2 630           ;Ma=L¯m 15     2 1 −1 2 181 −1 21 2     (13) The considered forces acting over the pile are punctual forces and moments applied at the top of the pile and distributed forces arising from pile-soil interaction. An axial force at the tip of the pile has been taken into account also. So, the vector of nodal forces Fcan be decomposed as F=Fext +Feq =Ftop +Fp+Feq (14) where Fext include the forces at the top Ftop and the axial force at the tip of the pile Fp; and Feq is the vector of the equivalent nodal forces from the pile-soil interaction, that can be calculated as Feq =Q·qp(15) where Qis the matrix that transforms nodal force components to equivalent nodal forces. The external forces defined over the generic element are schematized in Fig. 2. The tractions qpalong the pile-soil interface are approximated by the set of shape functions defined by Eq.(7) as 4 Figure 2: External punctual forces (left) and tractions along the pile-soil interface defined on the generic element qi=φ1qki+φ2qli+φ3qmi;i= 1,2,3 (16) Again, using the principle of virtual displacements, the coefficients of matrix Qfor lateral forces can be obtained as ql ij =ZL ϕiφjdx ;i= 1, ..., 5 ; j= 1,2,3 (17) and the ones for axial forces as qa ij =ZL φiφjdx ;i, j = 1,2,3 (18) This way, one can obtain the following matrices for lateral and axial equivalent nodal forces respectively           feq ki meq ki feq li feq mi meq mi           =           23L 140 11L 105 −L 28 L2 84 L2 105 −L2 210 4L 105 16L 35 4L 105 −L 28 11L 105 23L 140 L2 210 −L2 105 −L2 84               qp ki qp li qp mi     ;i= 1,2 (19) and     feq k3 feq l3 feq m3     =L 30     4 2 −1 2 16 2 −1 2 4        qp k3 qp l3 qp m3     (20) 5 Once all elemental matrices have been obtained for the whole pile, one can write, for each one, the following expression ¯ K up=Fext +Q qp(21) where ¯ K=K−ω2M. As each pile will be discretized using as many elements as necessary to follow its deformed shape accurately, matrices ¯ Kand Qare global matrices, obtained as usual from the elemental ones. 3 Soil BE equations The soil is modelled by BEM as a linear homogeneous isotropic viscoelastic unbounded region with complex valued shear modulus µof the type µ=Re[µ](1 + 2iξ), where ξis the damping coefficient. The boundary integral equation for a time-harmonic elastodynamic state defined in the domain Ω with boundary Γ can be written in a condensed and general form as ckuk+ZΓ p∗udΓ = ZΓ u∗pdΓ + ZΩ u∗XdΩ (22) where ckis the local free term matrix at collocation point ‘k’, which is a diagonal matrix with a 1 on rows corresponding to internal points and a 0.5 on rows corresponding to boundary points where the boundary is smooth; Xare the body forces in the domain Ω; uand pare the displacement and traction vectors; and u∗and p∗are the elastodynamic fundamental solution tensors on the boundary Γ due to a time-harmonic concentrated load at point ‘k’: u=   u1 u2 u3   ;p=   t1 t2 t3   ;u∗=  u∗ 11 u∗ 12 u∗ 13 u∗ 21 u∗ 22 u∗ 23 u∗ 31 u∗ 32 u∗ 33  ;p∗=  t∗ 11 t∗ 12 t∗ 13 t∗ 21 t∗ 22 t∗ 23 t∗ 31 t∗ 32 t∗ 33  (23) The fundamental solution used in this work corresponds to the complete space because there is no closed form expression for the half-space fundamental solution, which depends on boundless integrals and would require approximate procedures for its evaluation. Therefore, the soil free surface should be discretized. In practice, however, only a small region around the foundation has to be included in the model to reach accurate results. Generally, body forces Xare considered to be zero in most of elastodynamic problems. Nevertheless, in this work, pile-soil interaction is replaced by internal punctual forces at the piles tip and load-lines placed along the piles axis, while it is assumed that the soil continuity is not altered by the presence of the piles, which is the main hypothesis that supports this model. The load-lines qsjwithin the soil and its relationship with tractions qpare represented in Fig. 3 together with the internal punctual forces Fpjat the tip of the piles. Then, Eq.(22) can be written as ckuk+ZΓ p∗udΓ = ZΓ u∗pdΓ + np X j=1 "ZΓpj u∗qsjdΓpj−Υj kFpj#(24) where Γpjis the pile-soil interface of pile j,npis the total number of piles and Υj kis a 3 components vector that represents the contribution of the axial force Fpjat the tip of the jth pile, when the concentrated load is applied on point ‘k’. The boundary surface Γ is discretized into quadratic elements of triangular and quadrilateral shapes with six and nine nodes, respectively (see reference [20]). Once the boundary has been 6 Figure 3: Load-lines representation. discretized and Eq. (24) has been written for all nodes in Γ, the equation can be expressed in matrix form as Hssus=Gssp+ np X j=1 Gspjqsj− np X j=1 ΥsjFpj(25) where usis the vector of nodal displacements on the surface, Hss and Gss are matrices obtained by integration over Γ of the 3-D elastodynamic fundamental solution times the shape functions of the boundary elements, and Gspjis the matrix obtained by integration over Γpjof the 3-D elastodynamic fundamental solution times the interpolation functions defined in (7), when the unit load is applied over Γ. Assuming free traction surface (p= 0), Eq. (25) becomes Hssus− np X j=1 Gspjqsj+ np X j=1 ΥsjFpj= 0 (26) Furthermore, Eq. (24) will be also applied on the piles internal nodes, so, for a certain pile i, one can write upi k+Hpisus− np X j=1 Gpipjqsj+ np X j=1 ΥpijFpj= 0 (27) where upi kis the vector of nodal displacements at the node kof the pile iwhere the unit load is applied, Hpisis the matrix obtained by integration over Γ of the 3-D elastodynamic fundamental solution times the shape functions of the boundary elements, and Gpipjis the matrix obtained by integration over Γpjof the 3-D elastodynamic fundamental solution times the interpolation functions defined in (7), when the unit load is applied over a pile i. Besides, as an axial force at the pile tip is considered, an extra equation needs to be written. To do so, the unit load must be applied in the x3direction at any no-nodal point. The point with elemental dimensionless coordinate ξ=−0.5 in the bottom element of the pile has been chosen because of its nearness to the tip. This way, the extra equation is 1 83ubk 3+ 6ubl 3−ubm 3+ZΓ ˆ p∗udΓs= np X j=1 "ZΓpj ˆ u∗qsjdΓpj−Υj b3Fpj#(28) where ubk 3,ubl 3and ubm 3are the vertical displacements of nodes k,land mof the bottom element, ˆ p∗={p∗ 31, p∗ 32, p∗ 33}and ˆ u∗={u∗ 31, u∗ 32, u∗ 33}. In a matrix form, Eq. (28) can be written as 7 D upi b+Hpis eus− np X j=1 Gpipj eqsj+ np X j=1 Υpij b3Fpj= 0 (29) where upi bis the vector of nodal displacements at the bottom element nodes of the pile iwhere the unit load is applied, Hpis eis a vector obtained by integration over Γ of the 3-D elastodynamic fundamental solution times the shape functions of the boundary elements, and Gpipj eis a vector obtained by integration over Γpjof the 3-D elastodynamic fundamental solution times the interpolation functions defined in (7), when the unit load is applied on the extra point of the pile i.Dis the vector 1/8{0,0,3,0,0,6,0,0,−1}. The integrals of u∗qsjin (24) and ˆ u∗qsjin (28) over Γpjare calculated as a monodimensional integral extended to a load-line, defined by the pile axis, when the collocation point is outside pile j. However, this integrals have a singularity at the collocation point when applied on the integrated pile. In this case, in order to avoid this singularity, the integrals are evaluated over a cylinder which radius Rpis pπ/A. This way, let us consider the pile-soil interface Γp(whatever the kind of section) as a cylinder of radius Rpwhere tractions σps are applied. The last term of Eq. (22) includes integrals of the type ZΓp σpsu∗dΓp=ZΓp qs 2πRp u∗dΓp=1 2πRp Ne X e=1 X i=k,l,m qs iZΓpe φiu∗dΓp(30) where Neis the number of elements in which the load-line has been discretized and Eq. (16) has been used to express qsalong each element. The elastodynamic fundamental solution, that gives the displacement in a kdirection when the load is applied in the ldirection, can be written as u∗ lk =1 4πµ [ψδlk −χr,kr,l] ψ=−c2 c121 z2 1r2−1 z1rez1r r+1 z2 2r2−1 z2r+ 1ez2r r(31) χ=−c2 c123 z2 1r2−3 z1r+ 1ez1r r+3 z2 2r2−3 z2r+ 1ez2r r where µis the shear modulus, δlk is the Kronecker delta function, c1and c2are the wave velocities for P and S waves and zj=−ikj, being k1and k2the wave numbers for P and S waves, respectively. Then, the integrals in the last term of Eq. (30) can be solved, in cylindrical coordinates (see Fig. 4), as ZΓpe φiu∗dΓp=Zxr 3Zθ 1 4πµ [ψδlk −χr,kr,l]φiRpdθ dx3= Rp 4πµ Zxr 3 φih2πψδlk −π r2χRlkidx3(32) where xr 3=x3−xk 3and Rlk =  R2 p0 0 0R2 p0 0 0 2(xr 3)2 (33) 8 Figure 4: Integration over pile-soil interface when the collocation point belongs to the pile Now, the integral of Eq. (30) can be written as ZΓp σpsu∗dΓp= Ne X e=1 RpLe 8µX i=k,l,m qs jZ1 −1 φih2ψδlk −χ r2Rlkidξ (34) Figure 5: No-nodal collocation strategy Nevertheless, to compute integrals over Γpjfrom the same pile, a no nodal collocation strategy could also be carried out. This will lead to a procedure that allows the reinterpretation of the previous equation. This way, in order to avoid breaking the problem symmetries, at least four collocation points, symmetrically placed around the pile, should be chosen (see Fig. 5). A single equation can be obtained by adding these four equations divided by four, so that the arising coefficients are of 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 ao kG xx/(N·ks xx) Miura et al. BEM−FEM 4 x 4 2 x 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 3.5 4 ao cG xx/(N·ks xx) Miura et al. BEM−FEM 4 x 4 2 x 2 Figure 15: Horizontal impedances of 2x2 and 4x4 pile groups. Comparison with Miura’s solution. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 ao kG zz/(N·ks zz) Miura et al. BEM−FEM 4 x 4 2 x 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 ao cG zz/(N·ks zz) Miura et al. BEM−FEM 4 x 4 2 x 2 Figure 16: Vertical impedances of 2x2 and 4x4 pile groups. Comparison with Miura’s solution. (a) Shear forces due to horizontal excitation (b) Axial forces due to vertical excitation Figure 17: Forces on pile caps a 4x4 pile group Piles are modelled using Finite Elements (FEM) as a beam according to the Bernoulli hypothesis, while the soil is modelled using Boundary Elements (BEM) as a continuum, semi-infinite, isotropic, homogeneous, linear, viscoelastic medium. In the soil equations, the presence of piles is replaced by internal punctual forces at the tip and load-lines placed along the pile axis, while it is assumed that 16 the soil continuity is not altered by the presence of the piles. Then, the associated unknowns at internal points (displacements and tractions) are linked to piles variables through equilibrium and compatibility equations. The main advantage of this model is the capacity of computing dynamic behaviour of piles with low computing times and low memory requirements in comparison to other methods that need to discretize the pile surface or volume. This way, pile groups with a big number of members can be analyzed without difficulty. Besides, once the surface (not necessarily flat) has been discretized, it does not need to be changed to analyze different sets of piles. Several results of stiffness have been presented and compared to well known values taken from the literature, obtaining an excellent agreement. More cases than these presented in this work have been tested, all of them with the favourable conclusions. Although all results have been given for homogeneous soils, different configurations including soil strata and rigid rocky beds can be easily taken into account. Furthermore, other variables such as displacements or internal axial and shear forces and bending moments along the pile can be obtained without difficulty. Acknowledgements This work was supported by the Ministry of Education and Science of Spain through research project BIA2004-03955-C02-02. L.A. Padr´on is recipient of the FPU research fellowship AP-2004-4858 from the Ministry of Education and Science of Spain. Authors would like to thank for this support. References [1] A.M. Kaynia, “Dynamic stiffness and seismic response of pile groups”, Research Report R8303. Massachusetts Institute of Technology. Cambridge, Mass. 1982. [2] R.S.Y. Pak, P.C. Jennings, “Elastodynamic response of the pile under transverse excitation”, J Eng Mech, ASCE 1987;113(7):1101-16. [3] R.K.N.D. Rajapakse, A.H. Shah, “On the longitudinal harmonic motion of an elastic bar embedded in an elastic half-space”, Int J Solids Struct 1987;23(2):267-85. [4] R.K.N.D. Rajapakse, A.H. Shah, “On the lateral harmonic motion of an elastic bar embedded in an elastic half-space”, Int J Solids Struct 1987;23(2):287-303. [5] A.M. Kaynia, E. Kausel, “Dynamics of piles and pile groups in layered soil media”, Soil Dyn Earthq Eng 1991;10:386-401 [6] G. Gazetas, K. Fan, A.M. Kaynia, E. Kausel “Dynamic interaction factors for floating pile groups”, J Geotech Eng, ASCE 1991;117:1531-48. [7] A.M. Kaynia, M. Novak, “Response of pile foundations to Rayleigh waves and obliquely incident body waves”, Earthq Eng Struct Dyn 1992:21;303-18. [8] K. Miura, A.M. Kaynia, K. Masuda, E. Kitamura, Y. Seto, “Dynamic behaviour of pile foundations in homogeneous and non-homogeneous media”, Earthq Eng Struct Dyn 1994;23:183-92. [9] A.M. Kaynia, S. Mahzooni, “Forces in pile foundations under seismic loading”, J Engng Mech, ASCE 1996:122,46-53. 17 [10] H.B. Coda, W.S. Venturini, M.H. Aliabadi “A general 3D BEM/FEM coupling applied to elastodynamic continua/frame structures interaction analysis”, Int J Numer Meth Eng 1999;46:695712. [11] H.B. Coda, W.S. Venturini, “On the coupling of 3D BEM and FEM frame model applied to elastodynamic analysis”, Int J Solids Struct 1999;36:4789-804. [12] S.E. Kattis, D. Polyzos, D.E. Beskos, “Vibration isolation by a row of piles using a 3-D Frequency domain BEM”, Int J Numer Meth Eng 1999;46:713-28. [13] S.E. Kattis, D. Polyzos, D.E. Beskos, “Modelling of pile wave barriers by effective trenches and their sreening effectiveness”, Soil Dyn Earthq Eng 1999:18:1-10. [14] F. Vinciprova, J.J. Azn´arez, O. Maeso, G. Oliveto, “Interaction of BEM analysis and experimental testing on pile-soil systems.”, in “Problems in structural identification and diagnostic: General aspects and applications.”, C. Davini, E. Viola (Editors), Springer-Verlag, 195-227, 2003. [15] O. Maeso, J.J. Azn´arez, F. Garc´ıa, “Dynamic impedances of piles and groups of piles in saturated soils”, Comput Struct 2005;83:769-82. [16] A.V. Mendon¸ca, J.B. de Paiva, “A boundary element method for the static analysis of raft foundations on piles”, Eng Anal Boundary Elem 2000:24;237-47. [17] A.V. Mendon¸ca, J.B. Paiva, “An elastostatic FEM/BEM analysis of vertically loaded raft and piled raft foundation”, Eng Anal Boundary Elem 2003:27;919-33. [18] R. Matos Filho, A.V. Mendon¸ca, J.B. Paiva, “Static boundary element analysis of piles submitted to horizontal and vertical loads”, Eng Anal Boundary Elem 2005;29:195-203. [19] R.W. Clough, J.Penzien, “Dynamics of structures”, McGraw-Hill, 1982. [20] J. Dom´ınguez, “Boundary elements in dynamics”, Southampton, New York, Computational Mechanics Publications & Elsevier Applied Science; 1993. [21] H.B. Li, G.M. Han, H.A. Mang, “A new method for evaluation singular integrals in stress analysis of solids by the Direct Boundary Element Method”, Int J Numer Meth Eng 1985;21:2071-98. [22] F. Chirino, O. Maeso, J.J. Azn´arez, “Una t´ecnica simple para el c´alculo de las integrales en el sentido del valor principal en el MEC 3D”, Revista Internacional de M´etodos Num´ericos para C´alculo y Dise˜no en Ingenier´ıa 2000;16(1):77-95. 18