A 3D time domain numerical model based on half-space Green's function for soil-structure interaction analysis

Computational Mechanics manuscript No. (will be inserted by the editor) A 3D time domain numerical model based on half-space Green’s function for soil-structure interaction analysis. P. Galv´ın ·J. Maestre ·A. Romero Received: date / Accepted: date Abstract This paper presents a numerical method based on a three dimensional boundary element-finite element (BEM-FEM) coupled formulation in the time domain. The proposed model allows studying soil-structure interaction problems. The soil is modelled with the BEM, where the radiation condition is implicitly satisfied in the fundamental solution. Halfspace Green’s function including internal soil damping is considered as the fundamental solution. An effective treatment based on the integration into a complex Jordan path is proposed to avoid the singularities at the arrival time of the Rayleigh waves. The efficiency of the BEM is improved taking into account the spatial symmetry and the invariance of the fundamental solution when it is expressed in a dimensionless form. The FEM is used to represent the structure. The proposed method is validated by comparison with analytical solutions and numerical results presented in the literature. Finally, a soil-structure interaction problem concerning with a building subjected to different incident wave fields is studied. Keywords Half-space fundamental solution ·Internal soil damping ·BEM-FEM coupling · Soil-structure interaction ·Time domain 1 Introduction Soil-Structure Interaction (SSI) is a field of interest that involves structural analysis considering flexibility and damping due to the soil. Induced vibrations by rotatory machines, dynamic effects due to railway and road traffic, seismic problems, and dynamic behaviour of foundation systems are examples where SSI is an important issue. The effects of SSI are important and cannot be neglected in this type of problems [28]. Kausel [15] presented an exhaustive review of the main developments in this topic. Recently, Clouteau et al. [6] also reviewed different proposed numerical models to study the dynamic behaviour of structures on elastic media. Escuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla, Camino de los Descubrimientos, 41092 Sevilla, Spain Tel.: +34 955487293 (94) Fax: +34 954487296 E-mail: [email protected] 2 P. Galv´ın et al. Numerical models based on the Boundary Element Method (BEM) and the Finite Element Method (FEM) allow to study SSI problems rigorously. The BEM [8] is especially suited for the analysis of wave propagation in soils. Sommerfeld radiation condition [9] is satisfied implicitly and its semi-infinite character is well considered. The FEM is very useful to analyse the dynamic behaviour of structures accounting for nonlinear effects [29]. BEM models in time domain are usually based on full-space fundamental solutions due to their simplicity. Among other authors, Rizos and Karabalis [23] presented a BEM model to study general elastodynamic problems using the solution to Stokes’ problem within B- Spline polynomials. Schanz [26] proposed a BEM formulation based on a fundamental solution in the Laplace domain which is transformed at each time step using the convolution quadrature method. Ara´ujo et al. [1] developed a BEM model based on the time-marching Wilson θ -method. Marrero and Dom´ınguez [18] presented a BEM formulation for the solution of transient problems in bounded domains using an explicit full-space fundamental solution [9]. Later, Galv´ın and Dom´ınguez [10] extended this model to study wave propagation in unbounded regions. SSI problems can be solved by BEM models based on a full-space fundamental solution. However, it requires an additional discretization of the free field around the studied soil area to avoid cut-off errors related with the mesh boundaries. Therefore, the computational effort to solve SSI problems using a full-space solution could dramatically increase. This problem does not appear using a half-space fundamental solution since the free field condition is implicitly satisfied. Soil discretization is limited to those regions differing either from the half-space geometry or the free field condition. In SSI problems, soil discretization only concerns to soil-structure interfaces. However, a fully general half-space fundamental solution is not available in explicit form. Some solutions were proposed assuming some simplifications. Firstly, Pekeris [22] presented the half-space response for a vertical load and Poisson’s ratio 1/4. Afterwards, Chao [4] extended Pekeris’ solution for horizontal forces. Also, Mooney [19] proposed expressions for computing the half-space response excited by a vertical load and arbitrary Poisson’s ratio. A general half-space fundamental solution was proposed by Johnson [14]. In that work, the complete solution in implicit form to the Lamb’s problem was derived using the Cagniard-de Hoop method by applying Laplace transform with respect to the time and the spatial coordinate. Recently, Kausel [16] has provided a complete set of expressions for the half-space solution. The number of works that use the time formulation of the BEM and a half-space fundamental solution is quite reduced. Two works were found in the literature review. Firstly, Triantafyllidis [27] presented an approach using Johnson’s solution [14]. This work provides a physical meaning of the Green’s function terms. Afterwards, Bode et al. [3] proposed a BEM-FEM model to study SSI problems. The present paper proposes a time domain BEM formulation using the half-space fundamental solution presented by Johnson [14]. The innovative of this work is the effective treatment of the fundamental solution and its implementation in the BEM. The Green’s function has two kind of singularities that should be adequately solved to avoid unstable behaviour. A detailed integration procedure is presented. Moreover, internal soil damping is included in the fundamental solution. The paper is organized as follows. Firstly, the half-space fundamental solution is discussed, singularities are properly treated, and the internal soil damping is added. Afterwards, the BEM formulation based on the previous solution is described emphasising in the time integration of the Green’s function. Finally, a comprehensive SSI problem is solved using the proposed methodology. A 3D time domain numerical model based on half-space Green’s function for SSI analysis 3 2 Half-space fundamental solution In this work, the half-space fundamental solution proposed by Johnson [14] is analysed. Figure 1 shows schematically the problem geometry, where θ is the direct emission angle, ris the distance between the observation point and the source, and φ is the azimuth angle between the reference direction xand the projection of ron the plane z=0. The solution for displacements is a second order tensor, g(x,t;xi), that relates displacements at a point xdue to a Heaviside load applied at xi. The tensor component glk represents the kdisplacement component due to a load applied according to direction l. xi x θ r x y z φ Fig. 1 Problem geometry. The Green’s function for the half-space can be expressed as follows [14]: g(x,t;xi) = P(x,t;xi) +S(x,t;xi)+PP(x,t;xi) +SS(x,t;xi)+PS(x,t;xi) +SP(x,t;xi)(1) where Pand Sare tensors containing direct P and S wave fields, and PP,SS,PS and SP concern to P and S reflected waves in the half-space surface. The solution of Equation (1) involves the evaluation of expressions as: I=Re  Zpw 0 R σ  1 qt r2−1/c2 γ −p2  Ht−r/c γ dp  (2) where c γ is either the P-wave propagation velocity, cp, or the S-wave propagation velocity, cs, σ is the Rayleigh function, and the tensor Raccounts for the contribution of different waves. The integration is done from p=0 up to p=pw. The upper integration limit depends on the wave propagation velocities and the distance between the receiver point and the source. The physical meaning of these terms can be found in Reference [14]. Equation (2) presents two kinds of singularities. Firstly, a weak singularity is found at the poles of q(t/r)2−1/c2 γ −p2. A simple change of variable allows an effective treatment of this singularity [14]: p=(t/r)2−1/c2 γ 1/2+w2 γ =p,s(3) 4 P. Galv´ın et al. where wis the new integration variable. The second singularity occurs when the direct emission angle is θ = π /2, the source and the observation points are located at the half-space surface and, moreover, the time is longer than the arrival time of the Rayleigh waves. This singularity depends on the Rayleigh function poles. The Rayleigh function is expressed as: σ (p) = 1/c2 s−2(q2−p2)2+4(q2−p2)q1/c2 s−(q2−p2)q1/c2 p−(q2−p2)(4) being q(p) = −tsin( θ )/r+icos( θ )q(t/r)2−1/c2 γ −p2(5) where i=√−1. The homogeneous solution of the Equation (4) yields to the Rayleigh poles, pR: p2 R=cot2( θ ) c2 γ +(−1/c2 R+t2/r2)csc2( θ )± 2tcos( θ )csc2( θ )iqc2 γ /c2 R−1 c γ r(6) where cRis the Rayleigh wave propagation velocity. A strong singularity appears at Rayleigh poles defined by the Cauchy Principal Value. Johnson [14] proposed a refinement of the integration around pRusing the Romberg integration scheme to evaluate the integral. This procedure increases the computational effort and it is unstable for a direct emission angle close to π /2. This work proposes an alternative procedure to evaluate integrals as defined by Equation (2) based on the Residue Theorem. The integrals can be solved using a closed complex path to avoid the strong singularity (Figure 2). The integration path is divided into a real path defined between p=0 and p=pw, a deformed branch that excludes the non-analytical points due to P and S wave propagation velocities, pp=q(t/r)2−1/c2 pand ps=p(t/r)2−1/c2 s, and a complex path that closes the integration contour. The Equation (2) is rewritten as follows: 0 Im(p) Re(p) arc Deformed branch pw pp ps pR Fig. 2 Complex integration path. I=ReIb Rd p −Zde f.branch b Rd p −Zarc b Rdp= Re2 π iResb R,pR−Zde f.branch b Rd p −Zarc b Rdp(7) A 3D time domain numerical model based on half-space Green’s function for SSI analysis 5 where b Rencloses the integrand in Equation (2), and Resb R,pRis the residue of b Rat p=pR. Equation (7) is simplified taking into account that the kernels over the complex path are conjugated: ReZde f .branch b Rdp=Re2Zpw ps Im{b R}dp=0 (8) The residue also vanishes when p=pR. This fact is shown by separating b Rin two parts, the first one containing the Rayleigh function, σ , and the other one, N, accounts for the rest of parameters. If the integration variable tends to the Rayleigh pole pRthe following identity is obtained: Ren2 π iResb R,pRo=Re2 π ilim p→pR (p−pR)N/ σ =Re2 π iN/ σ ′(9) where σ ′at the Rayleigh pole, pR, is: σ ′=d σ (p) dp pR =4pR−2c4 Rc2 s+3c2 Rc4 s+c2 p−4c4 s+2c4 Rq(1/c2 p−1/c2 R)c2 sq1−c2 s/c2 R+ c2 Rc2 s3−4q(1/c2 p−1/c2 R)c2 sq1−c2 s/c2 R/c2 pc4 Rc2 sq(1/c2 p−1/c2 R)c2 sq1−c2 s/c2 R (10) The expression between the brackets in Equation (9) is purely imaginary since Nand σ ′are real values. Therefore, the real part of the residue vanishes. Finally, Equation (2) can be written as follows: I=Re      −Zarc R σ     2 r−w2−2qt r2−1/c2 γ    Ht−r/c γ dw       (11) Equation (11) does not show any singularity and it is integrated by a simple quadrature rule. In this paper, the internal soil damping in addition to the radiation soil attenuation is considered. An hysteretic damping model depending on the angular frequency, ω , the distance between the source and the observation point, r, and the S-wave propagation velocity, cs, is considered: gd(x, ω ;xi) = g(x, ω ;xi)e−i ξω r cs(12) where gd(x, ω ;xi)denotes the soil response in the frequency domain considering the internal damping. ξ is the hysteretic damping ratio. Equation (12) is expressed in the time domain by applying the inverse Fourier transform: gd(x,t;xi) = g(x,t;xi)∗csr ξ π (r2 ξ 2+c2 st2)(13) where ∗means time convolution process. Equation (13) represents the half-space fundamental solution considering internal soil damping in the time domain. The proposed procedure to evaluate the half-space fundamental solution is low time consuming due to the analytical treatment of the singularities. The traction solution, fd, is evaluated using the constitutive law for a linear elastic, homogeneous and isotropic solid by computing the spatial derivatives of the Green’s function. In the next sections, subscript dwill be omitted for a clearer notation. 6 P. Galv´ın et al. 2.1 Numerical verification The proposed solution is verified by comparison with known analytical solutions. Pekeris [22] presented an explicit displacement solution for a Heaviside vertical load applied at the half-space surface. Later, Chao [4] extended Pekeris’ solution for a horizontal load. Both solution are only valid for a half-space with Poisson’s ratio ν =1/4. Half-space properties were: P-wave propagation velocity cp=519.6m/s, S-wave propagation velocity cs=300m/s, and density ρ =1900kg/m3. Three different internal damping ratios were considered: ξ 1=0, ξ 2=0.02, and ξ 3=0.06. The response at a point located at r=10m from the source was analysed. Numerical results were computed by the non-stabilised integrals defined in Equation (2) and by the methodology proposed in Equation (13). Figure 3 shows the dimensionless displacement, e g(x, τ ;xi) = g(x,t;xi) µ r π /P, where µ =c2 s/ ρ is the transversal shear modulus and Pis the load amplitude. The response is represented versus dimensionless time τ =tcs/r. The response shows an initial perturbation at the arrival time of the P- wave ( τ =0.57). The maximum displacement is reached between the arrival time of S- wave ( τ =1) and the arrival time of Rayleigh wave ( τ =1.08). The agreement between the computed results and those obtained from the analytical solution [4,22] is quite good. The internal soil damping induces a smoother behaviour around the different arrival time of the waves. Numerical solution computed by the non-stabilised procedure shows an unstable behaviour after dimensionless time τ =1.08. These oscillations are eliminated using the proposed integration path. 3 Boundary element formulation This section presents the three-dimensional (3D) boundary element formulation in time domain based on the half-space fundamental solution previously developed. The integral representation of the displacement uifor a point ion the soil surface at time t, with zero body forces and zero initial conditions can be written as [8]: ci lkui k(xi,t) = Zt+ 0Z Γ u∗ lk(x,t− τ ;xi)pk(x, τ )d Γ d τ −Zt+ 0Z Γ p∗ lk(x,t− τ ;xi)uk(x, τ )d Γ d τ (14) where ukand pkstand for the kcomponent of the displacement and traction, respectively; u∗ lk(x,t− τ ;xi)and p∗ lk(x,t− τ ;xi)are the half-space fundamental solution displacement and traction tensors, respectively, at point xdue to an impulsive load applied at xiin ldirection. The integral-free term ci lk depends only on the boundary geometry at point xi, and its value is unitary at the soil surface for a smooth boundary. For non-smooth boundaries ci lk is computed as presented in Reference [17]. Half-space fundamental solution for an impulsive load is derived from the solution for a Heaviside load as: u∗ lk(x,t;xi) = ∂ glk(x,t;xi) ∂ t(15) p∗ lk(x,t;xi) = ∂ flk(x,t;xi) ∂ t(16) A 3D time domain numerical model based on half-space Green’s function for SSI analysis 7 0 0.5 1 1.5 2 −0.5 0 0.5 Dimensionless time [−] Dimensionless displacement [−] (a) 0 0.5 1 1.5 2 −0.5 0 0.5 Dimensionless time [−] Dimensionless displacement [−] (b) 0 0.5 1 1.5 2 −0.5 0 0.5 Dimensionless time [−] Dimensionless displacement [−] (c) 0 0.5 1 1.5 2 −0.5 0 0.5 Dimensionless time [−] Dimensionless displacement [−] (d) Fig. 3 Soil surface dimensionless displacement, e g(x, τ ;xi), for a point located at r=10m from the source due to a Heaviside load: (a) eg11, (b) eg31, (c) eg22, and (d) eg33. Analytical solution [4,22] (circles), non-stabilized solution (dashed black line), and proposed method for internal damping ξ 1=0 (black solid line), ξ 2=0.02 (dark grey solid line), and ξ 3=0.06 (light grey solid line). Displacements and tractions over the boundary are approximated at the time step t=n ∆ t from their nodal values um j kand pm j k, using the space interpolation functions φ j(r)and ψ j(r), and the time interpolation functions η m( τ )and µ m( τ )as: ci lkui k(xi,t) = n ∑ m=1 Q ∑ j=1Z Γ jZ ∆ tm u∗ lk(x,t− τ ;xi) µ md τψ jd Γ pm j k −Z Γ jZ ∆ tm p∗ lk(x,t− τ ;xi) η md τφ jd Γ um j k(17) where Qis the total number of boundary nodes and Γ jrepresents the elements to which node jbelongs. Time kernels in Equation (17) are analytically integrated using constant and linear piecewise time interpolation functions µ m( τ )and η m( τ )for tractions and displacements, respec- 8 P. Galv´ın et al. tively: Unm lk =Z ∆ tm u∗ lk(x,t− τ ;xi) µ md τ =glk(x,[n−m+1] ∆ t;xi)−glk(x,[n−m] ∆ t;xi) (18) Pnm lk =Z ∆ tm p∗ lk(x,t− τ ;xi) η md τ =1/ ∆ tFlk(x,[n−m+1] ∆ t;xi) −2Flk(x,[n−m] ∆ t;xi)+Flk(x,[n−m−1] ∆ t;xi)(19) where F(x,t;xi)is the time integral of the traction tensor f(x,t;xi). This integral is easily computed considering that the tensor f(x,t;xi)presents implicitly a time derivative [14]. After interpolating the boundary variables, the integral equation (17) becomes: ci lkui k(xi,t) = n ∑ m=1 Q ∑ j=1Z Γ j Unm lk ψ jd Γ pm j k−Z Γ j Pnm lk d τφ jd Γ um j k(20) Equation (20) is written in a more compact form as: ci lkuni k= n ∑ m=1 Q ∑ j=1hGnmi j lk pm j k−b Hnmi j lk um j ki(21) Once the integral-free term ci lk is included in the system matrix, the integral representation for point iat time t=n ∆ tbecomes: Hnnun=Gnnpn+ n−1 ∑ m=1 [Gnmpm−Hnmum](22) where Hnmi j lk collects for ci lk when i=jand n=m. In this work nine-node rectangular quadratic boundary elements are used. Identical space interpolation shape functions φ j(r)and ψ j(r)are assumed. The spatial integration is done only in those parts of an element under the effects of the fundamental solution waves, according to the causality condition of each term of the fundamental solution [18]. Each side of the element is divided into equal parts in the natural coordinate domain yielding an element subdivision. Equation (22) yields a system of equations that is solved step by step to obtain the time variation of displacements and tractions at the boundary. After boundary unknowns are solved, the radiated wave field un rat any internal point or at the free field is computed by means of the integral representation of the Somigliana identity: un r= n ∑ m=1Gnm rpm−b Hnm rum(23) where Gnm rand b Hnm rare the boundary element matrices computed considering only domain point as source points. The BEM performance is enhanced taking into account dimensionless expressions for the fundamental solution, e g(x, τ ,xi) = g(x,t,xi) µ r π /P, and the time vector, τ =tcs/r. Then, the fundamental solution is computed considering cylindrical coordinate system as: e gcyl(x, τ ;xi) = µ r π P egrr 0egrz 0eg φφ 0 egzr 0egzz  (24) A 3D time domain numerical model based on half-space Green’s function for SSI analysis 9 The fundamental solution in the cartesian coordinate system is obtained transforming the Equation (24) for each integration element: e g= Ω Ω Ω e gcyl Ω Ω Ω T(25) where Ω Ω Ω is the transformation matrix from cylindrical to cartesian coordinate system: Ω Ω Ω =  cos φ sin φ 0 −sin φ cos φ 0 0 0 1  (26) Equation (24) represents an antisymmetrical second order tensor. This dimensionless fundamental solution does not depend on the distance from the source to the observation point when the direct emission angle is θ = π /2. In this case, Equation (24) is computed just once and it is transformed for each collocation point and integration element. Several test have shown that the evaluation of the fundamental solution should be extended, at least, three times after the arrival time of the Rayleigh wave, i.e., τ =4.32. 3.1 Numerical validation The soil response due to a vertical Heaviside point load is studied to validate the proposed BEM formulation. The computed results are compared with analytical solution presented by Pekeris [22], and with numerical results presented by Triantafyllidis [27] and Bode et al. [3]. This benchmark problem has been previously used by Rizos and Karabalis [24] and by Romero et al. [25] to validate different boundary element formulations based on full space fundamental solutions. A half-space with the following properties was studied: cp=519.6m/s, cs=300m/s, and ρ =1900kg/m3. Three different internal soil damping ratios were considered: ξ 1=0, ξ 2=0.02, and ξ 3=0.06. The soil surface was discretized by one boundary element with a characteristic length ∆ l=0.5√2/2m (Figure 4). A load, p(t) = 1H(t)N, was applied at the central node of the element. Soil response was computed at a semicircular grid of 5151 receivers located in the free field and inside the soil. x y z p(t) l=0.5 p(t) t 1N Fig. 4 Half-space representation excited by a Heaviside load. The time step was set at ∆ t=4×10−4s in order to avoid stability problems [18]. The selected time step and element size give a value of the stability parameter β =cs ∆ t/ ∆ l= 0.34. 16 P. Galv´ın et al. peaks moved to lower frequencies when the softest soil was considered. In this case, the resonance frequencies were f3=0.6Hz and f4=7.0Hz. Figure 13.(b) shows the amplitude at the third resonance frequency at each floor for the different studied soil. The building response at the ground and at the first floor presented a significant difference because of the stiffness of the foundation. The soil flexibility caused a global translational displacement on the structure. The floor responses increased with the level. Figure 14 shows the building response at two resonance frequencies for the medium soil. The deflection shape presented a large translational displacement at the foundation, specially at the resonance frequency f3=4.0Hz. 0 5 10 15 20 25 30 10−6 10−4 10−2 100 Frequency [Hz] Displacement [m/Hz] (a) 01234 x 10−3 0 2 4 6 8 10 12 Amplitude [m/Hz] Storey level (b) Fig. 13 (a) Frequency content of the vertical displacement at the observation point Alocated at the top floor and (b) amplitude at the resonance frequency f3, for a building founded on stiff soil (black line), medium soil (dark grey line), and soft soil (light grey line). Axial and bending floor stiffness are considered. The previous results show that the proposed methodology allows to describe the dynamic SSI due to incident wave fields. In this case, the building response due to a transversal wave field was determined by the translational mode shapes. The structural behaviour due to vertical waves is highly influenced by the soil flexibility. 5.2 Dynamic response due to a moving load This example concerns with the building response due to a vertical load travelling at constant speed. the ground-borne vibrations were computed from the analytical solution proposed by Barber [2]. The structure was subjected to this incident wave field. This procedure is only valid if the structure is far enough from the load and its influence is negligible. The load moved according to a straight line located at y=26m from the building (Figure 15) at speed v=97.2m/s. Load travelled from xi=−696m to xf=1274m, and its value was P=170×103N. The medium soil with internal damping ξ =0.06 (Table 1) was considered. Figure 16 shows soil and building responses at three time steps: t1=6.88s, t2=6.98s, and t3=7.08s. These times correspond with load passage at x1=−16m, x2=0m, and A 3D time domain numerical model based on half-space Green’s function for SSI analysis 17 x y z (a) (b) 1.00 0.71 0.52 0.32 0.23 0.13 Fig. 14 Building response due to an incident wave field applied in zdirection: (a) f3=4.0Hz and (b) f4= 10.0Hz. Pv x y A Fig. 15 Geometry of the building excited by a moving load. x3=16m. Soil and building displacements are represented in different scale for showing soil and building deformed shapes. The load travelled in sub-Rayleigh regime (v<cR). Therefore, the induced wave field was symmetric around the load. The structure response was delayed with the load passage. The maximum displacement occurred at the top floor and at the third time step considered. Figure 17.(a) shows the time history for the vertical displacements at the observation point Aat the ground, and at the sixth, and at the twelfth floors. Maximum vertical displacement occurred at time t=7.1s when the load was in front of the observation point. The structural response was damped after the load pass, and exhibited oscillations according to 18 P. Galv´ın et al. (a) (b) (c) 4.0×10−3 3.4×10−3 2.8×10−3 2.2×10−3 1.6×10−3 1.3×10−3 Fig. 16 Soil and building displacements due to a moving load travelling at v=97.2m/s at time steps: (a) t1=6.88s, (b) t2=6.98s and (c) t3=7.08s. Colorbar is valid for building response. its first characteristic period, T1=1.19s. Maximum displacements amplitude was reached at the top floor. Figure 17.(b) compares the mobility at different floor levels. The structural behaviour due to the load passage was dominated by the first mode shape at frequency f1=0.86Hz, and it also appeared a peak at fR=21.0Hz. The characteristic wavelength of the Rayleigh propagation velocity (cR=141m/s) at frequency fRmatches almost exactly with the half of the length of the foundation ( λ R=v/fR=6.71m). Moreover, a low level peak was detected at 12.75Hz related with the Rayleigh wavelength. These wavelengths induced a periodic excitation in the building. The wavelengths did not match exactly with the building size due to the Doppler effect: the response at fixed receivers due to a harmonic moving load travelling at speed vis determined in a frequency range defined by a characteristic phase velocity, c, as [ ω (1−v/c), ω (1+v/c)] [7]. 0 5 10 15 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0x 10−5 Time [s] Displacement [m] (a) 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 1.2 x 10−11 Frequency [Hz] Mobility [m/sN] (b) Fig. 17 (a) Time history and (b) mobility of the vertical displacement due to a moving load travelling at v=92.7m/s at observation point Alocated at ground (solid black line), and 6th (solid dark grey line), and 12th (solid light grey line) floors. The mobility of the incident wave field (dashed black line) is also represented. A 3D time domain numerical model based on half-space Green’s function for SSI analysis 19 5.3 Dynamic response due to El Centro (1940) earthquake Finally, the dynamic behaviour of the building due to El Centro (1940) seismic accelerogram is evaluated. The y-direction was considered coincident with the North-South (NS) earthquake component. Earthquake accelerogram was described as an incident wave field. The accelerations were twice integrated to compute the incident displacement. The computed signal was filtered with a third order Chebyshev highpass digital filter, with a cut off frequency fc=1Hz. The medium soil with internal damping ratio, ξ =0.06, was also studied in this example. The earthquake accelerogram induced the maximum amplification at the top floor as can be seen in Figures 18.(a,b). The time histories of the floor responses showed a maximum peak acceleration around 6m/s2in the transversal direction, while low level acceleration was found at the vertical direction. Figures 18.(c,d) show the response spectra of the acceleration for El Centro earthquake and the building. The response spectra was computed for a period response from Ti=0.01s to Tf=10s, with a time step ∆ T=5×10−3s, and damping oscillator ζ =0.05 [5]. El Centro accelerogram was characterized by energy distribution over a frequency range below 20Hz. The structural response amplification was large at low frequencies, and it became much lower at higher frequencies. Amplifications of the horizontal response occurred at frequencies around f1=0.86Hz and f2=3.03Hz corresponding with the first and second mode shapes, respectively. A resonance was also found at the vertical building response at the fourth mode shape, f4=10Hz. Both acceleration components were highly amplified with the floor level. 6 Conclusions This paper has developed a time domain BEM formulation using Green’s function for the half-space as the fundamental solution. This methodology was used together with the FEM to study soil-structure interaction problems. The half-space fundamental solution is advantageous because of soil discretization is limited to soil-structure interface. However, the solution has two kind of singularities related to the arrival time of the Rayleigh waves. An effective treatment of this singularities that allows a low time consuming evaluation of the solution has been proposed. Visco-elastic soil behaviour is accounted for using an approximation based on an hysteretic damping model. The proposed fundamental solution has been verified with a benchmark problem. A BEM formulation using the previously analysed half-space solution has been presented. The fundamental solution can be written in dimensionless form to reduce the computational effort of the BEM. It is only evaluated once during the analysis. The proposed methodology was used to study the dynamic response of a building subjected to different sources: a horizontal and vertical uniform wave fields, a moving load, and an earthquake defined by El Centro (1940) accelerogram. Next conclusions were obtained: –Horizontal incident wave field excited low frequencies mode shapes related to the deformation of the global structure. Vertical incident wave field induced local deformation at the floors at higher frequencies. In both cases, maximum amplitudes were found at the top floor. –SSI modified the structural behaviour. The resonance frequencies moved to lower values, and the vertical mode shapes changed for softer soils due to the influence of the soilfoundation system. 20 P. Galv´ın et al. 0 5 10 15 20 25 −6 −4 −2 0 2 4 6 Time [s] Acceleration [m/s2] (a) 0 5 10 15 20 25 −6 −4 −2 0 2 4 6 Time [s] Acceleration [m/s2] (b) 0 5 10 15 20 25 30 0 5 10 15 20 25 Frequency [Hz] Acceleration [m/s2/Hz] (c) 0 5 10 15 20 25 30 0 5 10 15 20 25 Frequency [Hz] Acceleration [m/s2/Hz] (d) Fig. 18 (a,b) Time history and (c,d) response spectra of the (a,c) ycomponent and (b,d) vertical component at the observation point Alocated at ground (black line), and 6th (dark grey line), and 12th (light grey line) floors due to El Centro (1940) earthquake (dotted black line). –Building response induced by a moving load was determined by its resonance frequencies. In addition, Rayleigh phase wavelengths contributed to the building response when they matched with the structure foundation characteristic length. Acknowledgments This research was funded by the Spanish Ministry of Economy and Competitiveness (Ministerio de Econom´ıa y Competitividad) through research project BIA2010-14843. Financial support is gratefully acknowledged. The support given by the Andalusian Scientific Computing Centre (CICA) is also gratefully. References 1. Ara´ujo, F.C., Mansur, W.J., Nishikava, L.K.: Linear θ time-marching algorithm in 3D BEM formulation for elastodynamics. Engineering Analysis with Boundary Elements 23(10), 825–833 (1999) 2. Barber, J.R.: Surface displacements due to a steadily moving point force. Journal of Applied Mechanics, Transactions ASME 63(2), 245–250 (1996) 3. Bode, C., Hirschauer, R., Savidis, S.A.: Soil-structure interaction in the time domain using halfspace Green’s functions. Soil Dynamics and Earthquake Engineering 22(4), 283–295 (2002) A 3D time domain numerical model based on half-space Green’s function for SSI analysis 21 4. Chao, C.C.: Dynamical response of an elastic half-space to tangential surface loadings. Journal of Applied Mechanics, Transactions of the ASME 27(3), 559–567 (1960) 5. Clough, R.W., Penzien, J.: Dynamic of Structures. McGraw-Hill, New York (1975) 6. Clouteau, D., Cottereau, R., Lombaert, G.: Dynamics of structures coupled with elastic media - A review of numerical models and methods. Journal of Sound and Vibration 332(10), 2415–2436 (2013) 7. Ditzel, A., Herman, G.C., Drijkoningen, G.G.: Seismograms of moving trains: a comparison of theory and measurements. Journal of Sound and Vibration 248(4), 635–652 (2001) 8. Dom´ınguez, J.: Boundary elements in dynamics. Computational Mechanics Publications and Elsevier Aplied Science, Southampton (1993) 9. Eringen, A.C., Suhubi, E.S.: Elastodynamics, Volume 2, Linear theory. Academic Press, New York, USA (1975) 10. Galv´ın, P., Dom´ınguez, J.: Analysis of ground motion due to moving surface loads induced by high-speed trains. Engineering Analysis with Boundary Elements 31(11), 931–941 (2007) 11. Galv´ın, P., Romero, A., Dom´ınguez, J.: Fully three-dimensional analysis of high-speed train-track-soil- structure dynamic interaction. Journal of Sound and Vibration 329(24), 5147–5163 (2010) 12. Galv´ın, P., Romero, A.: A MATLAB toolbox for soil-structure interaction analysis with finite and boundary elements. Soil Dynamics and Earthquake Engineering 57, 10–14 (2014) 13. Heidebrcht, C., Stafford-Smith, B.: Approximate analysis of tall wall-frame structures. Structural Division 99, 199–221 (1973) 14. Johnson, L.R.: Green’s Function for Lamb’s Problem. Geophysical Journal of the Royal Astronomical Society 37(1), 99–131 (1974) 15. Kausel, E.: Early history of soil-structure interaction. Soil Dynamics and Earthquake Engineering 30(9), 822–832 (2010) 16. Kausel, E.: Lamb’s problem at its simplest. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 469(2149) (2013) 17. Mantic, V.: A new formula for the C-matrix in the Somigliana identity. Journal of Elasticity 33(3), 191–201 (1993) 18. Marrero, M., Dom´ınguez, J.: Numerical behavior of time domain BEM for three-dimensional transient elastodynamic problems. Engineering Analysis with Boundary Elements 27(1), 39–48 (2003) 19. Mooney, H.M.: Some numerical solutions for Lamb’s problem. Bulletin of the Seismological Society of America 64(2), 473–491 (1974) 20. Newmark, N.M.: A method of computation for structural dynamics. ASCE Journal of the Engineering Mechanics Division 85(1), 67–94 (1959) 21. Oliveto, G., Santini, A.: A simplified model for the dynamic soil-structure interaction of planar framewall systems. Engineering Structures 15, 431–437 (1992) 22. Pekeris, C.L.: The seismic surface pulse. Proceedings of the National Academy of Sciences of the United States of America 41, 469–480 (1955) 23. Rizos, D.C., Karabalis, D.L.: An advanced direct time domain BEM formulation for general 3-D elastodynamic problems. Computational Mechanics 15(3), 249–269 (1994) 24. Rizos, D.C., Karabalis, D.L.: A time domain BEM for 3-D elastodynamic analysis using the B-spline fundamental solutions. Computational Mechanics 22(1), 108–115 (1998) 25. Romero, A., Galv´ın, P., Dom´ınguez, J.: 3D non-linear time domain FEM–BEM approach to soil– structure interaction problems. Engineering Analysis with Boundary Elements 37(3), 501–512 (2013) 26. Schanz, M.: A boundary element formulation in time domain for viscoelastic solids. Communications in Numerical Methods in Engineering 15(11), 799–809 (1999) 27. Triantafyllidis, T.: 3-D time domain BEM using half-space Green’s functions. Engineering Analysis with Boundary Elements 8(3), 115–124 (1991) 28. Wolf, J.P: Dynamic soil-structure interaction. Prentice Hall, Englewood Cliffs (1985) 29. Zienkiewicz, O.C, Taylor, R.L.: The Finite Element Method, Volume I: The Basis. McGraw-Hill, New York (1967)