On soil-structurae interaction in large non-slender partially buried structures

On soil-structure interaction in large non-slender partially buried structures * J. Vega1 ,J. J. Azn´arez2 ,A. Santana2 ,E. Alarc´on3 ,L. A. Padr´on2 ,J. J. P´erez2 ,O. Maeso2 1Centro de Modelado en Ingenier´ıa Mec´anica (CEMIM-F2I2), 28006 Madrid, Spain 2Instituto Universitario de Sistemas Inteligentes y Aplicaciones Num´ericas en Ingenier´ıa (SIANI) Universidad de Las Palmas de Gran Canaria, 35017 Las Palmas de Gran Canaria, Spain 3Escuela T´ecnica Superior de Ingenieros Industriales, Universidad Polit´ecnica de Madrid, 28006 Madrid, Spain Abstract This paper addresses the seismic analysis of a deeply embedded non-slender structure hosting the pumping unit of a reservoir. The dynamic response in this type of problems is usually studied under the assumption of a perfectly rigid structure using a sub-structuring procedure (three-step solution) proposed specifically for this hypothesis. Such an approach enables a relatively simple assessment of the importance of some key factors influencing the structural response. In this work, the problem is also solved in a single step using a direct approach in which the structure and surrounding soil are modelled as a coupled system with its actual geometry and flexibility. Results indicate that, quite surprisingly, there are significant differences among prediction using both methods. Furthermore, neglecting the flexibility of the structure leads to a significant underestimation of the spectral accelerations at certain points of the structure. Keywords: Soil-structure interaction, Direct method, Three-step solution, Buried structures, Boundary Element Method 1 Introduction It is a well known fact that soil-structure interaction (SSI) should be considered in the dynamics of deeply embedded non-slender structures under ground motion. When analysing the dynamic response of such an structure (for instance large diameter caissons or shaft foundations, and pumping stations) the hypothesis of infinite stiffness is usually assumed. More precisely, embedded structures are considered to be rigid when the relation between depth of embedment and width is smaller than 2, 3 or 4 depending on the authors and application (Mylonakis,2001;Varun et al,2009), therefore many of these structures are seismically designed following that assumption of rigidness. The basic methods for the analysis of dynamic soil-structure interaction problems may be classified in direct and sub-structuring methods. Direct methods allow soil and structure to be modelled and analysed simultaneously in a single step. On the contrary, substructuring methods study soil and structure separately and, if the assumption of a perfectly rigid structure is adopted, the method leads to a three-step solution (Kausel and Roesset,1974;Kausel et al,1978;Villaverde, 2009) that allows parametric studies on the factors known to influence the response. However, when the hypothesis of rigidity of the structure foundation does not apply, the use of sub-structuring or direct methods requires soil and structure to be modelled with their actual geometry using more powerful numerical methods, which would imply higher computational effort in both cases. The literature is rich in examples of models addressing soil-structure interaction problems (Kausel , * This is the peer reviewed version of the following article: On soil-structure interaction in large non-slender partially buried structures. Bull Earthq Eng (2013) 11(5):1403–1421. The final publication is available at Springer via http://dx.doi.org/10.1007/s10518-013-9433-8. 1 2010) that use either the sub-structuring or the direct approaches in combination with Finite Element Method (FEM) or Boundary Element Method (BEM). This paper addresses the seismic analysis of a real massive deeply embedded structure, a structural typology traditionally analysed with the three-step approach, the objective being the assessment of the motion at discrete points of the structure consistent with a given free field ground motion. The structure is a part of a big industrial installation and hosts equipment that could be seriously affected during an earthquake. The relation between the embedded length and width of the structure is smaller than 2, approximately 1.7, thus it may be considered as a non-slender structure. This study has been performed using sub-structuring and direct methodologies in order to assess the validity of the rigid body assumption, which is the main goal of this paper. In the first part of the paper the problem at hand is presented. Then, the three-step and direct methods are briefly reviewed. Finally, results obtained with both methods are presented, compared and discussed. Attention will be paid not only to final results, but also to partial ones. The boundary element method is used when applying both methodologies. In each case complex valued relations between response at the points of interest and at free field are calculated. Frequency domain results are later transformed to the time domain using Fourier transform. 2 Problem description Figure 1shows a geometrical description of the structure under study, which is an 80 m high almost cylindrical concrete structure, with its lower 50 m embedded in the soil. The exterior diameter of the non-embedded part of the structure is 28 m. In the embedded length, the thickness of the excavation walls must be added, resulting in a exterior diameter of 30 m. The bottom 20 m of the structure are composed by a great number of slabs and stiffeners that provide a large rigidity. Therefore, this bottom part, defined as domain 2 (domain 1 corresponds to the rest of the structure) has been considered solid, with the stiffness of the concrete but an equivalent density. Table 1shows the values of the different properties of both domains of the structure, with µbeing the shear modulus, νthe Poisson’s ratio, ρthe density, ξthe internal damping coefficient and cs the S wave propagation velocity. The total mass Mof the structure is 61.3·106kg, the inertia IGis 3.51 ·1010 kg·m2and the distance hgfrom the base to the center of gravity is 28.36 m. The structure hosts the pumping unit of a large reservoir and other critical devices, such as a crane, that could be seriously affected during an earthquake. For this reason, part of the seismic analyses will focus on the natural frequencies of such devices, estimated to be around 3.89, 7.69 and 12.05 Hz. Table 1: Elastic properties of structure and soil layers domain 1 domain 2 soil layer 1 soil layer 2 soil layer 3 ν0.2 0.2 0.3 0.3 0.3 ρ2685.85 kg/m32253.42 kg/m32000 kg/m32100 kg/m32200 kg/m3 ξ0.05 0.05 0.05 0.05 0.05 µ1.15·1010 N/m21.15·1010 N/m25·108N/m21.029 ·109N/m22.2·109N/m2 cs2069.23 m/s 2259.06 m/s 500 m/s 700 m/s 1000 m/s The soil is composed of three layers. The top layer is a stiff clay formation and reaches 37 m depth. The second layer is formed by conglomerates and is 9 m depth. Finally, the bottom layer is a half space of strongly consolidated detrital sedimentary rock. Only static tests were carried out in the characterization of the different soil layers, so no measurements of S and P wave velocities were available. For this reason, some correlations and empirical rules, especially Ishihara (1996) and Seed and Idriss (1970), were used to estimate S wave velocities and soil dynamic shear stiffness. A sensitivity analysis was carried out in order to find out the most conservative soil profile among several ones, all compatible with the field tests results. Having defined the design accelerograms (described below) at the free field surface, the corresponding acceleration response 2 G O ∅25 m ∅28 m ∅30 m 30 m 30 m 20 m hg +0.0 m -37.0 m -46.0 m domain 1 domain 2 soil surface 2.5 m 1.5 m layer 1 layer 2 layer 3 x z M= 61.3·106kg ; IG= 3.51 ·1010 kg ·m2;hg= 28.36 m Figure 1: Geometrical description of the structure. spectra were computed at -50.0 m (using SHAKE (Schabel et al,1972) to obtain the accelerogram at such depth). The profile in figure 2was finally selected because it leads to higher (conservative) values for the range of frequencies of interest in regard to the hosted equipment. All results shown in this paper will be based upon it. Table 1summarizes the properties finally used to characterize the different soil layers. x z +0.0 m -37.0 m -46.0 m 500 m/s 700 m/s 1000 m/s Figure 2: Shear wave velocity profile for the soil. As previously stated, the aim of the dynamic study is the assessment of the motion at discrete 3 points of the structure consistent with a given free field ground motion. The latter has been defined in terms of the response spectrum shown in figure 3a, with a horizontal peak ground acceleration equal to 0.17 g, and 10 s conventional duration, as defined by the client’s specifications. In order to obtain some statistical significance, three different signals were generated using SIMQKE (Gasparini,1976) and then modulated using the scheme represented in figure 3b (Jennings et al, 1968). The resultant horizontal accelerograms are shown in figure 4, while the corresponding acceleration response spectra are presented in figure 5together with the design acceleration response spectrum. 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Sa/PGA Period (s) 0 1 0 2s 10s 18s ∼t 22∼e−0.268(t−10) ✐ a✐ b Figure 3: (a) Design acceleration response espectrum, (b) Modulation scheme. -0.3 -0.2 -0.1 0 0.1 0.2 0.3 acceleration (g) accelerogram 1 0.173 g -0.3 -0.2 -0.1 0 0.1 0.2 0.3 acceleration (g) accelerogram 2 0.172 g -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 2 4 6 8 10 12 14 16 18 acceleration (g) t (s) accelerogram 30.172 g ❜ ❜ ❜ Figure 4: Horizontal accelerograms 3 Methodology Once geometry, stratigraphy and seismic action have been defined, the seismic response of the structure is going to be computed by using two methodologies: a three-step approach and a direct approach (c.f. figure 6), as briefly described next. The objective is to simultaneously obtain the motion at the surface and at discrete points of the structure, under vertically plane harmonic incident S waves of fixed frequency. The complex valued relationships may be used to transform a given free field motion into the motion a the reference point, in the time domain. 4 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 Sa/PGA Period (s) Accelerogram 1 Accelerogram 2 Accelerogram 3 Design response spectrum Figure 5: Free-field elastic acceleration response spectra üg üg üki ki ¨ F M üki ki ¨ ü  ¨ ü  ¨ Kxx Kx Kx K K = Total Solution Kinematic Interaction 1 Impedance Matrix of Soil 2 3 Direct Methods 3 Steps Solution seismic excitation seismic excitation rigid interface K vs. rigid structure real structure Figure 6: Direct method versus three-step solution Substructuring methods in the hypothesis of a perfectly rigid structure provide the response following a three-step procedure (Kausel and Roesset,1974;Kausel et al,1978): i) computation of the kinematic interaction factors; ii) computation of stiffness and damping coefficients; and finally iii) estimation of inertial interaction (see figure 6). The second alternative is the use of direct methods. They model all domains defining soil and structure simultaneously. Therefore they allow taking into account the interactions among the different elements in a more rigorous way. The more popular numerical techniques used in this approach are FEM and BEM. In this work, the Boundary Element Method (Dom´ınguez,1993) has been used for the direct approach because it is specially suitable to deal with the dynamic analysis of unbounded domains such as the soil layers in this problem. When applying the three-step method, simplified expressions available in the bibliography are usually used to estimate kinematic interaction factors and impedance functions (Mylonakis et al, 2006;Kausel et al,1978). However, accurate enough expressions are not available for the layered soil profile of the problem at hand. For this reason, again BEM is also used to assess the kinematic 5 interaction coefficients and the soil impedances corresponding to the first and second steps of the sub-structuring methodology. When using the BEM to solve the problem following a direct approach, or in order to compute kinematic interaction factors and impedances for substructuring approach, all domains defining the geometry of the problem (soil layers and concrete walls) are modelled as linear homogeneous isotropic viscoelastic regions, and welded conditions are assumed between the different domains. Figures 7and 16 show two boundary element meshes used for the sub-structuring problems (step 1 and 2) and direct approach respectively. It can be appreciated that the boundary element code allows the use of quadratic triangular (6 nodes) and quadrilateral (9 nodes) elements. This particular geometry also includes corner problems where an interface cuts the embedded structure. Here, corner problems are solved by means of a non-nodal collocation strategy, which also allows using non-conforming meshes, as can be seen in the figures. The element size must be smaller than the half-wave length at the corresponding region for the highest frequency of analysis, in this case 25 Hz, although in general, all mesh parameters, such as free-surface extension and number of elements, are defined by performing convergence analyses of the variables of interest for different meshes. It is also worth noting that, even though the figure shows three quarters of the problem geometry, only a quarter of the total geometry is really meshed, as the code is able to take into account the symmetry properties of the problem. More details of the used boundary element code can be found in Maeso et al (2002,2004,2005). Next section provides the kinematic interaction coefficients, while the soil impedances are presented in section 5. These results are later used in section 6to compute the dynamic response of the whole system by solving the inertial interaction step of the sub-structuring method. They are also interesting because they enable understanding the sensitivity of the final response to different characteristics of the system. Finally, and following a direct approach, the real problem is solved modelling all elements simultaneously in section 7. Results are presented compared to those obtained through the three-step method. 4 Kinematic Interaction In this section, the displacements and rotations of the massless and infinitely rigid structure are obtained considering vertically incident S waves. The problem is well known (Roesset,1977), as well as its approach through the boundary element method (Dom´ınguez,1993;´ Alvarez Rubio et al, 2005). Figure 7shows the boundary element mesh used for this problem. The free surface and the interfaces among strata are displayed with different colours. The rigid interface among soil and structure is depicted in red. Figure 8presents absolute values of the translational Iu=u/uffand rotational IΦ=θR/uff kinematic interaction factors of the problem at hand, being uand θthe horizontal displacements and rotations at the base of the structure, Rthe radius of the structure and uffthe free-field horizontal motion at the ground surface. The horizontal axis represents the dimensionless frequency ao=ωR/cs, being ωthe angular frequency of the excitation and csthe velocity of the S waves in the top layer (500 m/s). The figure presents, together with the results obtained using boundary elements, the transfer functions obtained through the simple rules proposed by Elsabee et al (1977) for homogeneous soils (see also Kausel et al (1978); Mylonakis et al (2006)). The agreement between them is good at low and high non-dimensional frequencies, while discrepancies appear between 0.4 and 2.0. These differences would have a significant impact on the seismic assessment since the key components hosted by the structure have their fundamental frequencies within this range. Thus, only kinematic interaction factors obtained using BEM will be used in section 6. 5 Soil impedances. Equivalent springs and dashpots The second step in the three-step method is the determination of the soil impedances (c.f. figure 6). Impedance functions corresponding to a wide variety of cases can be found in the literature (Wolf, 6 Figure 7: Boundary elements mesh used for the kinematic interaction and soil impedance problems. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 1 2 3 4 5 Iu=u/uff Dimensionless frequency (ao) −50.0 m BEM Kinematic interaction Simple rules 0.0 0.1 0.2 0.3 0.4 0.5 0 1 2 3 4 5 IΦ=θR/uff Dimensionless frequency (ao) BEM Kinematic interaction Simple rules Figure 8: Kinematic interaction factors at the base of the structure for vertically incident S waves. 1985;Mylonakis et al,2006). However, no such functions are available for this particular configuration. Therefore, in this work, horizontal and rocking stiffness and damping coefficients have been obtained using the above mentioned boundary element code by prescribing unitary harmonic rigid-body horizontal displacements or rotations at the soil-structure interface. For each harmonic frequency the impedances are obtained by integration of the relevant stresses in the interface. It is worth noting that rotations induce horizontal stresses in the cylinder walls and, reciprocally, horizontal motions produce a momentum at the base. Therefore, with reference to the central point at the base of the structure, horizontal and rocking impedances are coupled. The coupling term becomes more important the deeper the structure is embedded. Impedance functions can then be arranged in a stiffness matrix Krelating forces and moments Fapplied at the center of the base to the resulting displacements and rotations uin a relationship of the type F=Ku that condenses the dynamic properties of the soil profile and that can be written as F M=Kxx Kxθ Kθx Kθθ  u θ(1) where each term Klm depends on the dimensionless frequency aoand represents the force (or moment) at the base of the structure needed to obtain a harmonic unitary displacement (or rotation). 7 −80 −40 0 40 80 120 160 200 0 1 2 3 4 5 kxx/(GR), kθθ/(GR3), kxθ/(GR2) Dimensionless frequency (ao) kxx kθθ kxθ=kθx −80 −40 0 40 80 120 160 200 0 1 2 3 4 5 cxx/(GR), cθθ/(GR3), cxθ/(GR2) Dimensionless frequency (ao) cxx cθθ cxθ=cθx Figure 9: Horizontal and rocking stiffness and damping functions As forces and displacements are out of phase, terms Klm are complex valued, in the form Klm =klm +i aoclm (2) where klm and clm are the stiffness and damping coefficients, respectively. Figure 9presents the stiffness and damping functions of the problem at hand, computed with the mesh shown in figure 7, and normalized using the shear stiffness of the top soil layer and the radius of the embedded structure. It is worth noticing that resonant frequencies associated to the shear deformation of the top two soil strata do not become apparent in these impedance functions. This is due to the fact that the related shear modes do not get excited because the structure goes all the way through the soil layers and does not produce significant vertical S waves. 6 Inertial interaction. Dynamic response of the structure. In previous steps the impedance functions, that condense the dynamic properties of the soil profile taking into account the embedded structure, and the seismic excitation including the effects of kinematic interaction have been assessed. In the third step, the dynamic response of the system is assessed under the assumption of the structure being infinitely rigid. Thus, the rigid body motion equations that govern the dynamic behavior of the structure can be written in the frequency domain as Kxx −ω2M Kxθ +ω2Mhg Kxxhg+Kθx Kxθhg+Kθθ −ω2IG u θ= Kxx Kxθ Kxxhg+Kθx Kxθhg+Kθθ  uki θki (3) where uand θare the horizontal displacement and rotation at the base of the structure; and uki and θki are the input horizontal displacement and rotation due to the seismic excitation, which are obtained from the kinematic interaction factors. Equation 3can be solved numerically for any given frequency to obtain the structural response of the rigid body in terms of its displacements uand rotations θ. This means that, once impedance functions and kinematic interactions factors have been computed, the cost of computing the system response is very low. This allows to perform parametric analyses and study the influence of different aspects on the final response at low cost, which is one key benefit of the three-step approach. Figure 10 presents the dynamic response of the structure, computed making use of this substructuring methodology, as a function of the dimensionless frequency ao. The translational kinematic interaction factors presented in section 4are also plotted in this figure in order to allow 8 0 0.2 0.4 0.6 0.8 1 1.2 u/uff +0.0 m Kinematic interact. Inertial interact. −7.5 m Kinematic interact. Inertial interact. 0 0.2 0.4 0.6 0.8 1 1.2 u/uff −15.0 m Kinematic interact. Inertial interact. −22.5 m Kinematic interact. Inertial interact. 0 0.2 0.4 0.6 0.8 1 1.2 0 1 2 3 4 5 u/uff Dimensionless frequency (ao) −33.5 m Kinematic interact. Inertial interact. 0 1 2 3 4 5 Dimensionless frequency (ao) −50.0 m Kinematic interact. Inertial interact. Figure 10: Frequency response functions for the normalized horizontal displacements from inertial interaction analysis and translational kinematic interaction factors for vertically incident S waves assessing the relative importance of kinematic interaction and inertial interaction on the final response of the system. In both cases, the transfer functions relate the horizontal displacement uat the base of the structure to the free-field motion uffat the ground surface. Under the rigid body motion assumption, horizontal displacements at any point of the structure can be easily obtained, and results for 6 different depths are presented. The results presented in figure 10 indicate that the presence of the structure filters out an important part of the seismic signal, mainly for frequencies ao≥2. In order to illustrate the effects of this filtering, figure 11 shows acceleration records evaluated on the cylinder at depths 0 and -50 m when the system is subjected to the first design accelerogram (c.f. figure 4). These responses have been computed as the inverse Fourier transform of the product between the frequency response function obtained for step 3 (inertial interaction) and the discrete Fourier transform of accelerogram 1, using the Fast Fourier Transform algorithm. The free-field accelerogram 1 is also presented in figure 11 in order to provide a baseline. The input seismic signal is filtered by the presence of the structure, as the peak acceleration is reduced from 0.173 to 0.087 g (∼ −50%) at -50.0 m depth, while the rotation induced by the incident wave-field increases the acceleration to 0.114 g at 0.0 m. The latter acceleration value is still 35% smaller than the free field peak acceleration of the input seismic signal. Elastic response spectra can also be obtained using these filtered signals at different depths as the excitation of one degree-of-freedom systems. Figure 12 presents the pseudo-acceleration 9 0 0.1 0.2 0.3 0.4 0.5 Sa (g) +0.0 m Free−field (0.0 m) Rigid model Flexible model Design response spectrum −7.5 m 0 0.1 0.2 0.3 0.4 0.5 Sa (g) −15.0 m −22.5 m 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Sa (g) Period (s) −33.5 m 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Period (s) −50.0 m Figure 18: Pseudo-acceleration response spectra at different levels obtained from the direct methodology (flexible model) and three-step solution (rigid model) methodologies. cannot be considered to be a rigid one. Also, the sub-structuring analysis have been useful to understand the key role of damping factors and rocking stiffness on the seismic response of the structure. It has been shown that, when damping factors are high, kinematic interaction is an extremely important element in the estimation of the structural response, and could even be used as a very good approximation to the inertial interaction step. This approximation will be even better if the natural frequencies of the system do not lie within frequency content of the seismic excitation. Acknowledgments This work was supported by the Subdirecci´on General de Proyectos de Investigaci´on of the Ministerio de Econom´ıa y Competitividad (MINECO) of Spain and FEDER through research project BIA2010-21399-C02-01 and also by the Agencia Canaria de Investigaci´on, Innovaci´on y Sociedad 16 −50 −40 −30 −20 −10 0 0 0.2 0.4 0.6 0.8 1 Deep (m) PA/PGA Rigid model Flexible model Figure 19: Validation of the hypothesis of rigidity of the structure. Maximum accelerations with depth under accelerogram 1. Vertically incident S waves. −50 −40 −30 −20 −10 0 0 10000 20000 30000 40000 50000 60000 Deep (m) Peak Normal Stress (N/m2) Rigid model Flexible model Figure 20: Validation of the hypothesis of rigidity of the structure. Maximum normal stresses with depth under accelerogram 1. Vertically incident S waves. de la Informaci´on (ACIISI) of the Government of the Canary Islands and FEDER through research project ProID20100224. A. Santana is recipient of the FPI research fellowship BES-2009-029161 from the MINECO. The authors would like to thank the engineer Mr. El´ıas Fern´andez who, by describing the problem as well as his concerns on the applicability of the current methods to solve it, motivated the research presented in the paper. References ´ Alvarez Rubio S, Benito JJ, S´anchez-Sesma FJ, Alarc´on E (2005) The use of direct boundary element method for gaining insight into complex seismic response. Computers & Structures 83:821– 835 Dom´ınguez J (1993) Boundary Elements in Dynamics. Computational Mechanics Publication: Southampton and Elsevier Applied Science: New York Elsabee F, Morray JP, Roesset JM (1977) Dynamic behavior of embedded foundations. Rep. No. R77-33, Massachusetts Institute of Technology, Cambridge MA Gasparini DA (1976) SIMQKE: A program for artificial generation. Rep. No. R76-4. Massachusetts Institute of Technology, Cambridge MA 17 Gerolymos N, Gazetas G (2006) Winkler model for lateral response of rigid caisson foundations in linear soil. Soil Dynamics and Earthquake Engineering 26(5):347–361 Ishihara K (1996) Soil Behaviour in Earthquake Geotechnics. Oxford Science Publications Jennings PC, Housner GW, Tsai NC (1968) Simulated Earthquake Motions. Caeltech Kausel E, Roesset JM (1974) Soil-structure interaction problems for nuclear containment structures. In: Electronic Power and Civil Engineer Conf Paper Power Div Specially Conf, Boulder, Colorado, pp 469–498 Kausel E, Whitman RV, Morray JP, Elsabee F (1978) The spring method for embedded foundations. Nuclear Engineering and Design 48:377–392 Kausel E (2010) Early history of soil-structure interaction. Soil Dyn Earthquake Eng vol. 30, no. 9, pp. 822–832 Maeso O, Azn´arez JJ, Dom´ınguez J (2002) Effects of space distribution of excitation on seismic response of arch dams. Journal of Engineering Mechanics (ASCE) 128 (7):759–768 Maeso O, Azn´arez JJ, Dom´ınguez J (2004) Three-dimensional models of reservoir sediment and effects on the seismec response of arch dams. Earthquake Engineering and Structural Dynamics 33:1103–1123 Maeso O, Azn´arez JJ, Garc´ıa F (2005) Dynamic impedances of piles and groups of piles in saturated soils. Computers & Structures 83:769–782 Mylonakis G (2001) Elastodynamic model for large-diameter end-bearing shafts. Soils and foundations 41(3):31–44 Mylonakis G, Nikolaou S, Gazetas G (2006) Footings under seismic loading: Analysis and design issues with emphasis on bridge foundations. Soil Dynamic and Earthquake Engineering 26:824– 853 Roesset JM (1977) Chapter 19: Soil Amplification of Earthquakes. In CS Desai & JT Christian Eds. Numerical Methods in Geotechnical Engineering. McGraw-Hill: NY Saitoh M, Watanabe H (2004) Effects of flexibility on rocking impedance of deeply embedded foundation. Journal of Geotechnical and Geoenvironmental Engineering 130 (4):435–445 Schabel B, Lysmer J, Seed HB (1972) SHAKE: A computer program for analysis of horizontally layered sites. Rep. No. EERC/72-12, University of California, Berkeley Seed HB, Idriss IM (1970) Soil Moduli and damping factors for dynamic response analysis. Report No 70-1. EERC Varun D, Assimaki D, Gazetas G (2009) A simplified model for lateral response of large diameter caisson foundations: Linear elastic formulation. Soil Dynamics and Earthquake Engineering 29(2):268–291 Villaverde R (2009) Fundamental Concepts of Earthquake Engineering. CRC Press, Taylor & Francis Wolf JP (1985) Dynamic Soil-Structure Interaction. Prentice Hall: Englewood Cliffs, NJ 18