SciELO - Scientific Electronic Library Online

vol.13 issue8A New Approach for Severity Estimation of Transversal Cracks in Multi-layered BeamsEffect of an Opening on Reinforced Concrete Hollow Beam Web Under Torsional, Flexural, and Cyclic Loadings author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Latin American Journal of Solids and Structures

Print version ISSN 1679-7817On-line version ISSN 1679-7825

Lat. Am. j. solids struct. vol.13 no.8 Rio de Janeiro Aug. 2016 


On the Validation of a Numerical Model for the Analysis of Soil-Structure Interaction Problems

Jorge Luis Palomino Tamayoa 

Armando Miguel Awruchb 

aCenter of Applied Mechanics and Computational (CEMACOM), Engineering School of Federal University of Rio Grande do Sul, Av. Osvaldo Aranha 99-3o Floor, 90035-190, Porto Alegre, RS, Brazil,

bDepartment of Civil Engineering, Engineering School Federal University of Rio Grande do Sul, Av. Osvaldo Aranha 99-3o Floor, 90035-190, Porto Alegre, RS, Brazil, +55(51)-33083450,


Modeling and simulation of mechanical response of structures, relies on the use of computational models. Therefore, verification and validation procedures are the primary means of assessing accuracy, confidence and credibility in modeling. This paper is concerned with the validation of a three dimensional numerical model based on the finite element method suitable for the dynamic analysis of soil-structure interaction problems. The soil mass, structure, structure's foundation and the appropriate boundary conditions can be represented altogether in a single model by using a direct approach. The theory of porous media of Biot is used to represent the soil mass as a two-phase material which is considered to be fully saturated with water; meanwhile other parts of the system are treated as one-phase materials. Plasticity of the soil mass is the main source of non-linearity in the problem and therefore an iterative-incremental algorithm based on the Newton-Raphson procedure is used to solve the nonlinear equilibrium equations. For discretization in time, the Generalized Newmark-β method is used. The soil is represented by a plasticity-based, effective-stress constitutive model suitable for liquefaction. Validation of the present numerical model is done by comparing analytical and centrifuge test results of soil and soil-pile systems with those results obtained with the present numerical model. A soil-pile-structure interaction problem is also presented in order to shown the potentiality of the numerical tool.

Keywords: Soil-structure interaction; fully coupled analysis; porous media; finite elements


The soil-structure interaction in general has been a concern and there is a need to further understand and better model this interaction. Civil structures are commonly supported on reinforced concrete shallow or deep foundations. When the rock basement of the soil deposit is far from the soil surface and the shear resistance of the soil is adequate, deep foundations made of concrete piles are typically used. In this way, proper numerical analyses of these structures are of importance for the civil and geotechnical engineering community in order to improve designs in terms of safety and economy. The main function of piles is to transferred axial an lateral loads from the superstructure to a reliable soil. Thus, several studies have been advocated to the study of soil-pile interaction problems under lateral and axial loads by using simplifying or complex approaches (Taiebat and Carter, 2001, Gu et al., 2016, Khoadir and Mohti, 2014, Chatterjee et al., 2015). Other research groups have focused on the development of simplified expressions for evaluating pile displacements and bending moments along the pile axis (Khodair and Mohti, 2014) by using the subgrade reaction approach (Valsamis, 2008, Valsamis et al., 2012, Chaloulos et al., 2013b). In this approach, the pile is treated as an elastic laterally loaded beam and the soil is idealized as a series of springs. Nevertheless, the nature of pile-soil interaction is three-dimensional and to complicate further the soil is a nonlinear and anisotropic medium (Khodair and Mohti, 2014). It is common practice to use deep foundations in areas where liquefaction may occur in surface layers. Thus, soil liquefaction is a close related topic because the soil mass surrounding the pile provides lateral support for the pile. This is the reason why is very important to study the behavior of deep foundation under lateral spreading since a small inclination of the ground surface can cause large deformations and thus, have a devastating effect on civil structures.

Liquefaction and associated shear deformation is a major cause of earthquake-related damage to piles and pile-supported structures. Pile foundation damage due to lateral spreading induced by liquefaction is documented in numerous reports and papers (Tasiopoulou et al., 2015a,b, Chaloulos et al., 2013a, Ou and Chan, 2006, Maheshwari and Sarkar, 2011, McGann et al., 2012, Valsamis et al., 2012, Cuellar 2011, Tamayo 2015). Liquefaction can take place not only during seismic excitation, but also some minutes later, thus consolidation phenomenon is also of interest in this work. The recognition of the importance of lateral ground displacement on pile performance has led to the development of analytical models capable of evaluating the associated potential problems. Modeling lateral ground displacement and pile response involves complex aspects of soil-structure interaction and soil behavior under large strains (Lu et al., 2004).

The soil mass should be modeled explicitly in the finite element mesh in order to simulate all the involved phenomena. In this way, soil-structure interaction effects are automatically accounted for. The numerical modeling of soil-structure interaction problems is rather complex and there are two main methods, namely, substructure and direct methods. Although the substructure method is less expensive than the direct method in terms of computational memory and time, the latter is preferred since its robustness for treating all nonlinearities together (Jeremic, 2004). In the direct method, the soil-structure model, which can be composed of shell, solid and beam-column elements, is analyzed in a single time step or load increment. Thus, the direct method is preferred in this work.

Another important aspect is the consideration of the dissipation of the excess pore pressure generated during the earthquake. Therefore, in order to count on with a reliable numerical tool for the complete analysis of this problem, it is necessary to validate extensively the numerical accuracy of the model (Tasiopoulou et al., 2015a, b). In this study, emphasis is given to the numerical modeling of soil-structure interaction problems based on deep foundations inserted in saturated soils under seismic actions. Hence, a general three-dimensional numerical model is coded by using the Fortran programming language. The inelastic behavior of the soil mass, soil-pile interface and concrete piles (Tamayo et al., 2013) can be included in the present numerical model. Nevertheless, for the validation examples presented in this work, the influence of slipping at the soil-pile interface in the global response of the soil-pile system was found to be negligible and this effect can be omitted with safety. This fact is also supported in (Chaloulos et al., 2013a) and (Cuellar, 2011) where it is stated that cohesionless soils can move together with piles, therefore none opening occur at the soil-pile interface. It is assumed that concrete piles have high strengths and they can be modeled as linear elastic materials. This last fact is acceptable since concrete piles are usually designed to remain elastic because subsurface damage is difficult to assess or repair. The adjacent saturated soil mass is considered to be fully saturated with water and its modeling is done by using the theory of porous media of Biot (Zienkiewicz et al., 1980, Lewis and Schrefler, 1998). The stress-strain behavior of the soil is represented by a plasticity-based effective stress constitutive model namely PZ-Mark III model (Pastor et al., 1990), which is suitable for simulating the behavior of cohesionless soils, including shear-induced pore-pressure generation (dilatancy) and cyclic mobility (Kumari and Sawant, 2013). The complete soil-structure interaction problem is solved within the framework of an incremental-iterative procedure of the Newton-Raphson type, while the Generalized Newmark-β scheme is used for solving the equilibrium equations in time.


The pile foundations and the structure above the ground level are modeled as one-phase materials (solids) with linear elastic laws and their formulations can be found in any book related to the finite element method. Otherwise, the saturated soil system is modeled as a two-phase material based on the theory of porous media of Biot. A numerical formulation of this theory, known as u-p formulation (in which displacement of the soil skeleton u, and pore pressure p, are the primary unknowns) was implemented. This implementation is based on the assumptions of small deformation and rotation and negligible fluid acceleration relative to the solid. The formulation presented here has been validated extensively with several experimental results in various works (Ou and Chan, 2006, Lewis and Schrefler, 1998, Tasiopoulou et al. 2015a,b, Kumari and Sawant, 2013). For more details about this section the reader is referred to the work of (Tamayo, 2015).

2.1 Governing Equations

The coupled set of equations of Biot that governs the behavior of a saturated porous media is given in the following way:

Equilibrium of mixture is defined in the following manner:


where σ ij is the total stress tensor of the skeleton (tensile positive), ü i is the acceleration of the solid skeleton, bi is the body acceleration per unit mass, ρs , ρf and ρ are the densities of the solid grain, fluid and mixture, respectively, with ρ = (1 - n) ρs+ n ρf and n being the porosity of the porous media and ai represents the component of the fluid acceleration relative to the solid in the xi direction. When the deformed and undeformed configuration of the porous media are almost identical (the case of small deformations and rotations), the concepts of Lagrangian and Eulerian porosity are the same. The porosity and the specific mass of the porous media can vary significantly along time during a consolidation analysis and therefore these variables must be updated continuously. On the other hand, in earthquake related problems, the duration of interest is very short (5-20 sec.) and during the shaking, it is possible to have flow of water within the soil mass, however it is not expected that the porosity change in this short period can be of any significant order (Leung, 1984).

Equilibrium of fluid:


where R represent the viscous drag forces which, assuming the Darcy seepage law can be written as k'ijRj = wi , k' = k / ρ'fg', where ρ'f and g' are the fluid density and gravitational acceleration at which the permeability is measured and k is the permeability of soil with dimensions of [m]/[sec.].

Conservation of mass for fluid phase


where wi,i is the flow divergence in the unit volume, is the increased volume due to a change in strain, n/Kf is the additional volume stored by compression of void fluid due to the fluid pressure increase, (1 - n) ṗ / Ks is the additional volume stored by the compression of grains by the fluid pressure increase and KT ( + ṗ / Ks ) / Ks is the change in volume of the solid phase due to a change in the inter-granular effective contact stress. The mass conservation equation can be further expressed by using the definition of and Q in the following way:


where KT is the average bulk modulus of the solid skeleton, Ks is the average material bulk modulus of the solid components of the skeleton and Kf is the bulk modulus of the fluid, with 1/Qn/Kf + ( - n)/Ks n/Kf + (1 - n)/Ks and = 1 - KT / Ks .Combining (eqs. 2)-(4), together with (eq. 1), neglecting the underlined terms which are generally small, the governing equations of the u-p formulation can be expressed in the following way:



For most soils, ≈ 1, that is, the incompressibility of the soil grains is considered (Ks >>KT ).

2.2 Discretization of the Governing Equation in Space

For the spatial discretization of the governing equations, the finite element method is used. The variables u and p are interpolated by suitable shape functions Nu and Np , respectively, in the following manner:



where û and are the nodal displacement vector and the nodal pressure vector, respectively. The definition of Biot effective stress (in vectorial form) is given in the following form:


where m is the vectorial form of the delta of Kronecker. The governing (eqs. 5)-(6) may be expressed in the following finite element matrix form (Lewis and Schrefler, 1998):












where M is the consistent mass matrix, Q is the coupled matrix, H is the permeability matrix, S is the compressibility matrix, G is the seepage matrix, fu and fp are the volume forces that act on the surface Г for the solid and fluid phase, respectively, L is a matrix operator of derivate, B is the usual strain-displacement matrix and is the prescribed traction on boundary and is the prescribed influx. Viscous damping is also incorporated into the dynamic equation of the solid phase (eq. 10) in the form of , where


is called the Rayleigh damping matrix (Kumari and Sawant, 2013). The coefficients α1 and α 2 can be obtained by selecting a damping ratio ξ n and a certain frequency ωn such that


When consolidation analysis is of interest, it is only necessary to eliminate all inertial terms in the above formulation. Otherwise, nonlinear elasticity and theory of generalized plasticity are used to determine the relationships between incremental stresses and strains. The incremental stress-strain relationship is expressed in the following form:


where and Dep are the incremental stress and elastoplastic constitutive material tensors, respectively. The elastoplastic constitutive material is defined in the following way:


where De , n, ngL/U and HL/U are the elastic constitutive tensor, loading direction vector, flow direction vector under loading or unloading conditions, and loading or unloading plastic modulus, respectively.

2.3 Discretization in Time

In this work, (eqs. 10)-(11) must be integrated in time using the single-step Generalized Newmark-β (GNpj) method. Using GN22 for the displacements u and GN11 for the pore pressure p, the following expressions are obtained:







where δ = 0,50, β=0.25 and θ = 1.0 are used for unconditional stability of the integration scheme and t refers to current time. For a time interval ∆t, the second term on the left hand side of (eq. 10) can be expressed as with .

2.4 Constitutive Model for Sands

Soil behavior under cyclic loading is complex. Hence, the constitutive model used in a numerical code should be able to capture important features of soil behaviour under cyclic loading such as permanent deformation, dilatancy, and hysteresis loops to obtain reliable solutions of displacements and pore water pressure. For this study the constitutive model described by (Pastor et al., 1990) was used for the sand. The P-Z Mark III model is a generalized plasticity-bounding surface-non associative type model (Kumari and Sawant, 2013). The model is described by means of yield surfaces and potential surfaces which are described by the following equations:



in which, p is mean confining stress; q is deviatoric shear stress; Mg is slope of the critical state line; α f and α g are constants; Pc and Pg are size parameters. The dilatancy of the sand in the P-Z Mark III model is approximated using the linear function of the stress ratio η = q / p as (Kumari and Sawant, 2013):


and p v and p s are incremental plastic volumetric and deviatoric strains, respectively. Mg is related to the angle of friction Φ' by the Mohr-Coulomb relations in the following way:


value of S is ± 1 based on compression or extension. The plastic flow direction under loading ngL is given in the following way:


The non-associated flow rule is used, and then the loading direction is expressed as:




Mf maintains a constant ratio with Mg . (Pastor et al., 1990) assumed this ratio to be dependent on relative density (Dr ) suggesting a relation for Mf in the following manner:


In the P-Z Mark III model, the plastic modulus for loading (HL ) is obtained as:




where Ho , β o and β1 are model parameters; and dξ is plastic deviatoric strain increment. The undrained triaxial test predict rapid pore pressure build up on unloading. This highlights the necessity to predict plastic strains on unloading in a constitutive model. The P-Z Mark III model uses the following expression for the plastic flow direction ngU and the unloading plastic modulus Hu .




η u is called the unloading stress ratio given by


Huo and γ u are specified material constants.


3.1 Sand Deposit Submitted to Harmonic Loading at Base (Taboada and Dobry, 1993)

For verification of the developed code towards liquefaction analysis, the class A prediction of the experiment No 1 of VELACS (Verification of Numerical Procedures for the Analysis of Soil Liquefaction Problems) project is considered. The experiment carried out by (Taboada and Dobry, 1993) consists of a 0.20 m high, horizontal, uniform Nevada sand layer, which is placed in a laminar box at a relative density of about 40% (loose sand). The purpose of the laminar box is to simulate the response of a semi-infinite loose sand layer during shaking. A sketch of the laminar box and the instrumentation used for this experiment is presented in Figure 1 (a). Material properties are listed in Table 1 (see loose sand column with relative density Dr = 40%). The experiment was carried out at 50-g acceleration, where g is the gravity acceleration. The soil deposit is submitted to a lateral movement at its base according to harmonic function shown in Figure 2 with a maximum peak acceleration of 0.2g.

Figure 1 Geometry and mesh of the finite element model 

Figure 2 Horizontal input motion at bottom. 

Table 1 Soil properties. 

Property Loose sand (Dr = 40%) Dense sand (Dr = 60%) Units
Elastic linear analysis
Elasticity modulus 30000 30000 kPa
Coefficient of Poisson 0.3 0.3
Non-linear analysis with PZ-Mark III
Compressibility modulus at p' o Kevo = 770 2000 kPa
Shear modulus at p' o Keso = 1155 2600 kPa
Reference pressure p' o p' o 4 4 kPa
Critical state line Mg = 1.15 1.32
State line for loading Mf = 1.035 1.3
Dilatancy parameter α g = 0.45 0.45
Dilatancy parameter α f = 0.45 0.45
Shear hardening parameter β o = 4.2 4.2
Shear hardening parameter β1 = 0.2 0.2
Plastic modulus for loading Ho = 600 750 kPa
Plastic modulus for unloading Huo = 4000 40000 kPa
Parameter for plastic unloading γ u = 2 2
Other properties
Specific mass of the soil ρ = 2.089 kN.s/m4
Specific mass of the fluid ρ f = 0.98 kN.s/m4
Volumetric modulus of the solid particle Ks = 1017 kPa
Volumetric modulus of the fluid Kf = 1.092x106 kPa
Porosity n = 0.363
Permeability (prototype scale) k = 3.3x10-3 m/s
Permeability (model scale) k = 6.6x10-5 m/s
Gravity acceleration g = 9.81 m/s2

Numerical modeling is done in model scale using a three dimensional formulation with a plane-strain condition. For this purpose, lateral faces at the xz plane are allowed only to move in its plane. Due to symmetry, only half of the model is considered. The finite element mesh is composed of 5120 coupled hexahedral finite elements with 8-node for pore pressure and 8-node for solid displacements (called here 8-8 node elements). The number of degrees of freedom to solve is 24276. The mesh is regular and uniform as shown in Figure 1 (b). The laminar box is modeled with the constraint of the lateral tied nodes. The displacements of nodes located at the two ends of the soil at the same level are restrained to have the same value. The base nodes are fixed in both horizontal and vertical directions. Dissipation of pore pressure is allowed only through the top surface of the layer; the lateral boundaries and the base are kept impermeable. The maximum size of a finite element in the mesh is limited to (L = (1/5 - 1/8)Vs /f in order to permit good wave transmission within the model, where Vs is the shear wave velocity of the soil and f is the predominant frequency of loading. The limitation in the xy plane is more flexible and the subdivision in this plane can be chosen between 3(L - 5(L (Tamayo, 2015). The Rayleigh damping matrix was used with a damping coefficient of ξ = 5% and with a circular frequency of loading ω = 2 πf. First a static analysis due to application of gravity (model's own weight) is performed before seismic excitation. The resulting fluid hydrostatic pressures and stress-states along the soil mass are used as initial conditions for the subsequent dynamic analysis (Ou and Chan, 2006).The magnified deformed mesh and excess pore pressure at the end of the analysis are shown in Figures 3 (a)-(b), respectively. In (Figures 4)-(5) are compared the development of the excess pore pressure at points P1, P2 and P3, P4, (see locations in Figure 1 (a)), respectively, as predicted by the present numerical model and those recorded in the experiment. In Figure 6, the lateral displacement at locations LVDT3 and LVDT4 are also depicted. As it can be seen, a reasonable good agreement between numerical and experimental results is shown.

Figure 3 Deformation and excess of pore pressure after 16.38 s. 

Figure 4 Excess pore pressure histories for P1 and P2. 

Figure 5 Excess pore pressure histories for P3 and P4. 

Figure 6 Lateral displacement for LVTD3 and LVDT4. 

Cycles of shear stress-strain histories for different depths in the soil mass are shown in the left part of Figure 7. The shear stress τ h refers to the shear stress component τ xz acting in the plane perpendicular to the longitudinal movement due to shaking, while γ h is the associated shear deformation component γ xz . As it may be observed, the shear deformation reaches a maximum value around 1.5%. Similarly, the right part of Figure 7 shows the associated shear stress-effective vertical stress paths in the soil mass. As it may be observed, soils at depths 2.81 m and 5.31m are in a liquefied state because the effective vertical stresses has almost reduced to zero at the end of shaking.

Figure 7 Cycles of shear stress-strain histories (left) and shear stress-effective vertical stress paths (right) in the soil at 2.81 m, 5.31 m, and 7.81 m depths 

3.2 Soil-Pile System Under Consolidation Load (Taibet and Carter, 2001)

In (Taibet and Carter, 2001) is studied the time-dependent behavior of a vertical pile inserted in a saturated soil mass and which is submitted to a lateral loading H at its head. In the mentioned reference, a semi-analytical finite element method based on discrete and continuous Fourier transformations was used for the analysis of the saturated porous media. The pile diameter is Dp = 2.0m and this is inserted in a non-cohesive saturated soil deposit that follows the Mohr-Coulomb law. The influence of the dilatancy angle ψ in the soil and pile response was studied by the authors considering a constant porosity value along time. The initial effective stresses in the soil mass are determined with a lateral earth coefficient at rest of 0.5. The geometry of the problem, boundary conditions and the finite mesh used in this work are shown in Figure 8. It is important to comment that in (Taibet and Carter, 2001), a semi-analytical plane finite element mesh was used and this mesh served as a basis for the discretization of three-dimensional mesh shown in Figure 8 (b). The pile is considered to be elastic and material properties are listed in Table 2.

Figure 8 Geometry, boundary conditions and finite element mesh 

Table 2 Material properties. 

Material Property Units
Specific weight (γ s ) 17 kN/m3
Elasticity modulus (Es ) 30000 kPa
Coefficient of Poisson (νs) 0.3
Cohesion c 0.0 kPa
Friction angle (Ф') 30
Dilatancy angle (Ψ) 0 (non-associative) or 30 (associative)
Permeability (k) 1x10-4 m/s
Fluid specific weight (γ f ) 10 kN/m3
Specific weight (γ c ) 23 kN/m3
Elasticity modulus (Ec ) 30x106 kPa
Coefficient of Poisson (νc) 0.2

The finite element mesh is formed by 996 20-8 node hexahedral finite elements (20 nodes for the solid phase and 8 nodes for the fluid one) for simulating the saturated soil and 84 20-node hexahedral elements for modeling the concrete pile. The number of equations to solve is 16032. In order to provide a direct comparison with the results provided by (Taiebat and Carter, 2001), all obtained results are expressed in terms of a dimensional time Tv = k(1 - vs )Est/yf (1 - 2vs )(1 + vs )D 2 p , load rate ω = d(H/yfD 3 p )/dTv and pile diameter Dp .

An elastic analysis of the problem was first conducted to evaluate the accuracy of the present model. A horizontal load H was applied rapidly to the pile head, thereafter the load was held constant with time (ramp load). The predicted lateral displacement of the pile head in the direction of the applied load is plotted against the dimensionless time,Tv in Figure 9. Results of the analysis using the discrete and continuous Fourier series method suggested by (Taiebat and Carter, 2001) are also shown for comparison. As it may be observed, the results of the analysis using the discrete Fourier approach are in close agreement with the results obtained in the present analysis.

Figure 9 Lateral displacement at pile head 

Secondly, in a series of elasto-plastic analysis, the total lateral load was varied from H = 5yfD 3 p to H = 60yfD 3 p , where yf is the specific weight of the water . In each case the total load was applied during a time interval of Tv = 0.0001, with a loading rate of ω= 100000. This loading rate was sufficiently high to approximate an initial undrained loading. Thereafter, the load was held constant in time and the analysis was continued, allowing excess pore pressures to dissipate, and thus for the soil to consolidate. The time-dependent lateral displacements at the pile head predicted by the elasto-plastic analyses with both associated (dilatancy and friction angles are equal) and non-associated (null dilatancy angle) flow rules are plotted in Figure 10 for the particular case of H = 15yfD 3 p . The response of the pile in elastic soil is also presented.

Figure 10 Lateral displacement at pile head 

The largest pile head displacement is predicted by the elasto-plastic soil model with a non-associated flow rule. This displacement value is almost twice the value predicted by using an associated flow rule. The stiffer behavior of a pile in soil with associated flow rule can be attributed to the dilative characteristics of the soil after failure. Expansion of the soil after failure increases confining pressures which in turn increase soil resistance, causing stiffer behavior in comparison to the behavior of soil with non-associated flow rule (Taiebat and Carter, 2001). The predicted load-horizontal displacement curves at the pile head under fully drained (load is applied slowly) and undrained conditions are presented in Fig. 11.

Figure 11 Lateral displacement at pile head versus applied lateral force H 

In order to permit a direct comparison with the results provided by (Taiebat and Carter, 2001), the excess pore pressures p are expressed in a non-dimensional form p / (yfDp ).The distributions of the dimensionless excess pore pressures at the end of rapid loading for H = 15yfD 3 p are shown in Figures 12 (a)-(b) for the cases of soil with an associated and non-associated flow rule, respectively. The interested reader can compare these distributions with those provided in the work of (Taiebat and Carter, 2001) and good agreement can be inferred.

Figure 12 Dimensionless excess pore pressure close to the pile 

3.3 Soil-Pile System Under Harmonic Load (Abdoun, 1997)

In the centrifuge test reported by (Abdoun, 1997), a single pile model (called model No 3 in the experimental work) was embedded in a saturated soil deposit and submitted to a lateral movement at its base. The experiment was conducted using the rectangular, flexible-wall laminar box container shown in Figure 13 (a). The soil profile consists of two layer of fine Nevada sand saturated with water: a top liquefiable layer of relative density, Dr = 40% and 6 m prototype thickness, and a bottom slightly cemented non-liquefiable sand layer with a thickness of 2 m. According to the experimental work, the cemented non-liquefiable sand layer has similar properties of a non-liquefiable sand layer with Dr = 60%. Material properties are listed in Table 1. The prototype single pile is 0.6 m in diameter, 8 m in length, has a bending stiffness, EI = 8000 kN.m2, and is free at the top. The model has an inclination angle of 2o (model scale) and is subjected to the predominantly 2Hz harmonic base excitation shown in Figure 14 with a peak acceleration of 0.3g, where g is the gravity acceleration.

Figure 13 Geometry and mesh of the finite element model. 

Figure 14 Movement at base (Abdoun, 1997). 

The centrifuge test was simulated using the above-described three-dimensional finite element model. As it may be observed in Figure 13 (b), the soil domain and the single pile were discretized using 3D 8-8 node and 8-node brick finite elements, respectively. The finite element size in the mesh is limited to a maximum value in order to permit a good wave transmission (see example of section 3.1). A half mesh configuration is used due to symmetry considerations. The number of degrees of freedom in the mesh is 6056. The boundary conditions were (i) dynamic excitation was defined as the recorded base acceleration, (ii) at any depth, displacement degrees of freedom of the downslope and upslope boundaries were tied together (both horizontally and vertically) to reproduced a 1D shear beam effect, (iii) the soil surface was traction free, with zero prescribed pore pressure, and (iv) the base and lateral boundaries were impervious. A static application of gravity (model own weight) was performed before seismic excitation. The resulting fluid hydrostatic pressures and soil stress-states served as initial conditions for the subsequent dynamic analysis. The Rayleigh damping matrix was used with a damping coefficient of ξ = 5%.

With a mild inclination of 2o, model 3 attempts to simulate an infinite slope subjected to shaking parallel to the slope. However, it was noted that, in the centrifuge test, the soil surface gradually lost the slope and became level during the shaking phase. To simulate such behavior of losing the surface slope, a horizontal component of gravity varying with time was applied to the finite element simulation. The load-time history of the applied horizontal gravity component was calculated based on the recorded lateral displacement at ground surface (Lu et al., 2004). Figures 15 (a)-(b) displays the computed lateral displacements and pore pressure ratios, respectively, at the end of the analysis. As it can be seen, liquefaction was reached down to a depth of 5.0 m (see Figure 15 (b)), as indicated by the pore-pressure ratio ru approaching to 1.0 (ru = (p/σ'v where (p is the excess pore pressure and σ'v is the initial effective vertical stress). The Nevada sand layer remained liquefied until the end of shaking and beyond. Thereafter, excess pore pressure started to dissipate.

Figure 15 Finite element results at the end of shaking. 

The mild inclination of model 3 also imposed a static shear stress component (due to gravity), causing accumulated cycle-by-cycle lateral deformation. The recorded and computed short-term and long-term excess pore pressure histories for two control points (PP1 and PP2) at depths of 1 m and 5 m are compared in Figure 16 and Figure 17, respectively. Both computed and recorded results displayed a number of instantaneous sharp pore pressure drops after initial liquefaction. The numerical results obtained in the works of (Lu et al., 2004) and (Valsamis, 2008) are also depicted for comparison.

Figure 16 Short-term excess pore pressure time histories of PP1 and PP2 

Figure 17 Long-term excess pore pressure time histories of PP1 and PP2 

The permanent lateral displacement of the ground surface after shaking is approximately 100 cm. All lateral displacement occurred in the top 6.0 m within the liquefiable sand layer. The top graph in Figure 18 shows the recorded and computed pile lateral displacement at the soil surface during and after shaking. The computed pile lateral displacement increased to 40 cm and decreased to approximately 8 cm at the end of shaking. The bottom slightly cemented sand layer, as indicated in the bottom graph, did not slide with respect to the base of the laminar box. Lateral displacements at some other control points in the soil mass are also shown in Figure 18.

Figure 18 Lateral displacement at various monitoring points 

Figure 19 shows the profile of pile lateral displacements obtained with the present numerical model for the time in which the maximum lateral displacement occurs at the pile head. In the same graph are also depicted the pile displacement predictions according to some simplifying approaches based on the works of (Valsamis et al., 2012), (Brandenberg, 2005), (Cubrinovski et al., 2006), (Tokimatsu and Asaka, 1998), American Petroleum Institute (API, 2005), High Pressure Gas Safety Institute of Japan (HPGS, 2000) and Railway Technical Research Institute (RTR, 1999). All these references have used a methodology based on the subgrade reaction approach (p-y method). As it may be observed, there is a great difference among all methodologies. The closer predictions to the experimental values are due to the works of (RTR, 1999), (API, 1995) and (Brandenberg, 2005).

Figure 19 Profile of pile lateral displacements 

Because it is true that major liquefaction does not necessarily takes place at the end of the shaking, in Figure 20 is depicted the liquefaction evolution measured by the ru factor at different time instants. As it may be observed, a dilatation zone appears close to the pile head (blue color zones). At the end of the shaking (after 25 seconds of analysis), liquefaction almost took place in the top liquefiable sand layer, thus reflecting a similar condition as established in the experimental work (Lu et al., 2004).

Figure 20 Evolution of liquefaction in the soil mass measured by the liquefaction potential factor ru

Figure 21 shows the cycles of shear stress-strain histories in the soil mass for depths of 1.5 m, 3.5 m and 5.5 m, respectively, in the soil zones close to the pile (near field) and at the edge of the laminar box (free-field). As it may be observed, the maximum shear deformation component expressed in percentage occurs for a depth of 1.5 m (around 22%) in the near field, while it is around 23% in the free-field for a depth of 3.5 m. As it was stated in the introduction section, liquefaction and associated shear deformation in the present example are considerable.

Figure 21 Cycles of shear stress-strain histories in the soil at 1.5 m, 3.5 m, and 5.5 m depths for model 

Figure 22 shows the shear stress- effective vertical stress paths during shaking for the soil zones located in the near field and free-field. As it may be inferred from the graphs, all soil samples have an initial effective vertical stress and a static shear stress value as a result of the initial conditions of the soil (gravity load) and due to the surface inclination. During shaking the effective vertical stress is almost reduced to zero for various soil depths due to soil liquefaction.

Figure 22 Shear stress-vertical effective stress paths in the soil at 1.5 m, 3.5 m, and 5.5 m depths 

3.4 Soil-Structure Interaction Examples Under Harmonic Load

In order to show the potentiality of the present numerical model, a soil-structure interaction problem is studied. A steel frame building is supported by a group of concrete piles which are square and have a length and side length of 10.5 and 0.5 m, respectively. The number of piles in the group is 35 (configuration 7x5, seven piles along the x direction and five piles along the y direction with a space ratio s/d = 3, where d is the side length of the pile and s is the distance among piles) and they are joined together by a concrete cap. The soil domain is defined by a region of 41.5 m x 30.5 m x 14m along the x, y and z directions, respectively. The typical story height of the steel frame is 3 m and this story is composed of three and two spans of 3.25 m each one along the x and y directions, respectively. The finite element mesh used in the analysis is shown in Figure 23 and it is formed by 1380 8-node hexahedral finite elements for modeling the pile and cap, 15816 8-8 node hexahedral finite elements for modeling the soil domain, 230 8-node quadrilateral zero-thickness contact elements for modeling the soil-pile interface and 228 2-node truss elements for modeling the steel frame. A half mesh configuration is used due to symmetry considerations and the number of degrees of freedom to solve is 77980. The structural steel properties used in the computations correspond to a W13x426 section with a specific weight of 7.8 kN/m3. In order to show liquefaction, the soil mass is considered to be composed of a uniform liquefiable sand layer with a relative density Dr = 40% (see Table 1). The load is applied at the base of the model by using the E-W component of the Centro earthquake (1940) as shown in Figure 24. The predominant frequency of the earthquake is 2Hz and the damping ratio for the Rayleigh damping matrix is 5%. The boundary conditions used in the example of section 3.3 are also used here.

Figure 23 Finite element mesh for soil-steel frame system 

Figure 24 The Centro Earthquake (1940, E-W component) 

In Figure 25 is shown the evolution of the liquefaction potential factor ru for different time steps. As it may be observed, liquefaction starts after 4 seconds of analysis in the soil zones among piles. Thereafter, the liquefaction mechanism spreads out to the lateral edges of the soil domain. Figure 26 shows the lateral displacement histories at the pile cap and stories of the steel frame. As it may be observed, the maximum lateral displacement at the cap is around 0.10 m. Because the soil zones among piles are almost liquefied and there is not exist a non-liquefiable surface layer, the soil zones surrounding the piles do not provide a rigid lateral support, thus pile lateral displacements increase considerably.

Figure 25 Evolution of liquefaction potential according to ru factor 

Figure 26 Horizontal displacements at different levels of the structure 

In order to study the soil-structure effect in the problem, Figure 27 shows the axial force histories in truss elements for the case of the steel frame clamped to the ground directly (fixed base) and when the soil domain is included in the analysis (flexible case). It may be observed that both cases yielded considerably different axial forces in the elements with the absolute higher values associated to the flexible base condition. Precisely, Figure 28 shows the axial force histories in two representative elements namely element 1 and element 2 (see Figure 23) of the steel frame for both supporting conditions.

Figure 27 Axial force histories in truss elements (kN) 

Figure 28 Axial force histories in elements 1 and 2 

As it may be observed, the axial forces in these elements are significantly different for the cases of rigid and flexible base support. Thus, for the present example, the soil deposit must be included in the analysis. This result was expected since the soil mass is represented by a loose sand layer, which is not rigid at all.


In this study, a computer program for the three dimensional finite element analysis of soil-structure systems was developed. The numerical model uses the coupled dynamic field equation with the u-p formulation of Biot's theory for modeling saturated soils, while concrete caps and piles are modeled as monophasic materials. The superstructure can be modeled by using shell, solid and beam-column finite elements. The numerical model is firstly validated with some benchmarks found in the technical literature. Detailed comparison between numerical and experimental results for soil and soil-pile systems showed acceptable matching. Also, the numerical model was able to predict dilatation zones close to the pile heads, which characterize soil-pile systems involving liquefaction. Only after validation processes have been successfully completed, the following parametric studies can be done with the aim of the present numerical tool. Numerical modeling and simulation can be used to predict the behavior of piles or pile groups embedded in fully saturated soils and thus used to improve design in terms of safety and economy. When liquefaction is involved in the analysis, a complex constitutive model capable of simulating dilatancy and cyclic mobility in the soil mass was proved to perform well. Potentiality of the numerical tool is shown by modeling a steel frame building supported by a group of concrete piles under seismic load. The obtained results showed that the soil zones among piles have a high potential of liquefaction (ru close to one). Also, axial forces in the steel frame elements are underestimated when the frame is considered to be clamped to the ground directly. Studies are in progress in order to improve the modeling of hysteretic damping in the soil.


The authors gratefully acknowledge the financial support provided by CAPES and CNPq.


Abdoun T. (1997). Modeling of seismically induced lateral spreading of multi-layered soil and its effect on pile foundations. Ph.D. Thesis, Rensselaer Polytechnic Institute, USA. [ Links ]

API (1995). Recommended practice for planning, designing and constructing fixed offshore platform, Washington, DC: American Petroleum Institute [ Links ]

Brandenberg SJ. (2005). Behavior of pile foundations in liquefied and laterally spreading ground. Ph.D. Thesis, University of California, Davis, USA. [ Links ]

Chaloulos Y. K., Bouckovalas G.D, Karamitros D. K. (2013a). Pile response in submerged lateral spreads: Common pitfalls of numerical modeling techniques, Soil Dynamics and Earthquake Engineering 55: 275-287 [ Links ]

Chaloulos Y. K., Bouckovalas G.D, Karamitros D. K. (2013b). Analysis of liquefaction effects on ultimate pile reaction to lateral spreading, Journal of Geotechnical and Geoenviromental Engineering 140: 1-11 [ Links ]

Chatterjee K.,Choudhury D., Poulos H.G. (2015). Seismic analysis of laterally loaded pile under influence of vertical loading using finite element method, Computers and Geotechnics 67: 172-186 [ Links ]

Cubrinovski M., Kokusho T., Ishihara K. (2006). Interpretation from large-scale shake table tests on piles undergoing lateral spreading in liquefied soils, Soil Dynamics and Earthquake Engineering 26: 275-286 [ Links ]

Cuellar, P. (2011). Pile foundations for offshore wind turbines: numerical and experimental investigations on the behavior under short-term and long-term cyclic loading. Ph.D. Thesis. Berlin:Technischen Universitat, Germany [ Links ]

Gu Q., Yan Z., Peng Y. (2016). Parameters affecting laterally loaded piles in frozen soils by an efficient sensitivity analysis method, Cold Regions and Technology 121:42-51 [ Links ]

High Pressure Gas Safety Institute of Japan (2000). Design method of foundation for level 2 earthquake motion (In Japanese) [ Links ]

Jeremic B. (2004). Domain reduction method for soil-foundation-structure interaction analysis. Technical Report UCD-CGM 01-2004. [ Links ]

Khodair Y., Mohti A. (2014). Numerical analysis of pile-soil interaction under axial and lateral loads, International Journal of Concrete Structures and Materials 8: 239-249 [ Links ]

Kumari S., Sawant V. A. (2013). Use of infinite elements in simulating liquefaction phenomenon using coupled approach, Coupled Systems Mechanics 2: 375-387 [ Links ]

Leung K. H. (1984). Earthquake response of saturated soil and liquefaction, Ph.D. Thesis, University College of Swansea, UK [ Links ]

Lewis R. W., Schrefler B. A. (1998). The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media (2th Edition). John Wiley and Sons Ltd (New York). [ Links ]

Lu J., He L., Yang Z., Abdoun T., Elgamal A. (2004). Three dimensional finite element analysis of dynamic pile behavior in liquefied ground. In: Doolin D, Kammerer A, Nogami T, Seed R, eds. Proceedings of the 11th International Conference on Soil Dynamics and Earthquake Engineering(11ICSDEE), Berkeley, United States [ Links ]

Maheshwari B. K., Sarkar, R. (2011). Seismic behavior of soil-pile structure interaction in liquefiable soils: parametric study, International Journal of Geomechanics 11: 335-347 [ Links ]

McGann C. R., Arduino P., Mackenzie-helnwein P. (2012). Stabilized single-point 4-node quadrilateral element for dynamic analysis of fluid saturated porous media, Acta Geotechnica 7: 297-311 [ Links ]

Ou J. H., Chan A. H. C. Three dimensional numerical modeling of dynamic saturated soil and pore fluid interaction. In: Topping B H V, Montero G, Montenegro R, eds. Proceedings of the Fifth International Conference on Engineering Computational Technology, Stirlingshire, United Kingdom, 2006 [ Links ]

Pastor M., Zienkiewicz O.C., Chan A. H. C. (1990). Generalized plasticity and the modeling of soil behavior, International Journal for Numerical and Analytical Methods in Geomechanics 14: 151-190 [ Links ]

Railway Technical Research Institute (1999), Earthquake resistant design code for railway structures, Maruzen Co. (In Japanese) [ Links ]

Taboada V. M, Dobry R. (1993). Experimental results and numerical predictions of model No 1. In: Arulanandan K, Scott R F, eds. Proceedings of the International Conference on the Verification of Numerical Procedures for the analysis of soil liquefaction problems, California, United States. [ Links ]

Taiebat H.A., Carter J.P. (2001). A semi-analytical finite element method for three-dimensional consolidation analysis, Computers and Geotechnics 28:55-78 [ Links ]

Tamayo J. L. P. (2015). Numerical simulation of soil-pile interaction by using the finite element method. Ph.D. Thesis. Federal University of Rio Grande do Sul, Brazil. [ Links ]

Tamayo J. L. P.,Morsch I.B, Awruch A. M. (2013). Static and dynamic analysis of reinforced concrete shells, Latin American Journal of Solids and Structures 10: 1109-1134 [ Links ]

Tasiopoulou P., Taiebat M., Tafazzoli N. (2015a). On validation of fully coupled behavior of porous media using centrifuge test results, Coupled Systems Mechanics 4:37-65 [ Links ]

Tasiopoulou P., Taiebat M., Tafazzoli N. (2015b). Solution verification procedures for modeling and simulation of fully coupled porous media: static and dynamic behavior, Coupled Systems Mechanics 4:67-98 [ Links ]

Tokimatsu K., Asaka Y. (1998). Effects of liquefaction-induced ground displacements on pile performance in the 1995 Hyogoken-Nambu earthquake, Soil and Foundations, Special Issue: 163-177 [ Links ]

Valsamis A. I (2008). Numerical simulation of single pile response under liquefaction induced lateral spreading. Ph.D. Thesis. Department of Geotechnical Engineering, School of Civil Engineering, National Technical University of Athens, Greece. [ Links ]

Valsamis A. I, Bouckovalas G. D, Chaloulos Y. K. (2012). Parametric analysis of single pile response in laterally spreading ground, Soil Dynamics and Earthquake Engineering 34: 99-110 [ Links ]

Zienkiewicz O. C., Chang C. T., Bettes P. (1980). Drained, undrained, consolidating and dynamic behavior assumptions in soils, Geotechnique 30: 385-395 [ Links ]

Received: September 07, 2015; Revised: March 31, 2016; Accepted: April 03, 2016

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License