1 INTRODUCTION
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.
2 FINITE ELEMENT FORMULATION AND IMPLEMETATION
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, b_{i} 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 a_{i} represents the component of the fluid acceleration relative to the solid in the x_{i} 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'_{ij}R_{j} = w_{i} , k' = k / ρ'_{f}g', 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 w_{i,i} is the flow divergence in the unit volume, is the increased volume due to a change in strain, nṗ/K_{f} is the additional volume stored by compression of void fluid due to the fluid pressure increase, (1 - n) ṗ / K_{s} is the additional volume stored by the compression of grains by the fluid pressure increase and K_{T} ( + ṗ / K_{s} ) / K_{s} 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 K_{T} is the average bulk modulus of the solid skeleton, K_{s} is the average material bulk modulus of the solid components of the skeleton and K_{f} is the bulk modulus of the fluid, with 1/Q ≡ n/K_{f} + ( - n)/K_{s} ≅ n/K_{f} + (1 - n)/K_{s} and = 1 - K_{T} / K_{s} .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 (K_{s} >>K_{T} ).
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 N^{u} and N^{p} , 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}):
with,
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, f_{u} and f_{p} 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 Cû, 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 dσ and D^{ep} are the incremental stress and elastoplastic constitutive material tensors, respectively. The elastoplastic constitutive material is defined in the following way:
where D^{e} , n, n_{gL/U} and H_{L/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:
and
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; M_{g} is slope of the critical state line; α _{f} and α _{g} are constants; P_{c} and P_{g} 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 dε^{p} _{v} and dε^{p} _{s} are incremental plastic volumetric and deviatoric strains, respectively. M_{g} 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 n_{gL} is given in the following way:
The non-associated flow rule is used, and then the loading direction is expressed as:
with
M_{f} maintains a constant ratio with M_{g} . (^{Pastor et al., 1990}) assumed this ratio to be dependent on relative density (D_{r} ) suggesting a relation for M_{f} in the following manner:
In the P-Z Mark III model, the plastic modulus for loading (H_{L} ) is obtained as:
with
where H_{o} , β _{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 n_{gU} and the unloading plastic modulus H_{u} .
η _{u} is called the unloading stress ratio given by
H_{uo} and γ _{u} are specified material constants.
3 NUMERICAL EXAMPLES
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 N^{o} 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 D_{r} = 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.
Property | Loose sand (D_{r} = 40%) | Dense sand (D_{r} = 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} | K_{evo} = | 770 | 2000 | kPa |
Shear modulus at p' _{o} | K_{eso} = | 1155 | 2600 | kPa |
Reference pressure p' _{o} | p' _{o} | 4 | 4 | kPa |
Critical state line | M_{g} = | 1.15 | 1.32 | |
State line for loading | M_{f} = | 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 | H_{o} = | 600 | 750 | kPa |
Plastic modulus for unloading | H_{uo} = | 4000 | 40000 | kPa |
Parameter for plastic unloading | γ _{u} = | 2 | 2 | |
Other properties | ||||
Specific mass of the soil | ρ = | 2.089 | kN.s/m^{4} | |
Specific mass of the fluid | ρ _{f} = | 0.98 | kN.s/m^{4} | |
Volumetric modulus of the solid particle | K_{s} = | 1017 | kPa | |
Volumetric modulus of the fluid | K_{f} = | 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/s^{2} |
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)V_{s} /f in order to permit good wave transmission within the model, where V_{s} 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.
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.
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 D_{p} = 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.
Material | Property | Units | |
---|---|---|---|
Soil | |||
Specific weight (γ _{s} ) | 17 | kN/m^{3} | |
Elasticity modulus (E_{s} ) | 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/m^{3} | |
Concrete | |||
Specific weight (γ _{c} ) | 23 | kN/m^{3} | |
Elasticity modulus (E_{c} ) | 30x10^{6} | 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 T_{v} = k(1 - v_{s} )E_{s}t/y_{f} (1 - 2v_{s} )(1 + v_{s} )D ^{2} _{p} , load rate ω = d(H/y_{f}D ^{3} _{p} )/dTv and pile diameter D_{p} .
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,T_{v} 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.
Secondly, in a series of elasto-plastic analysis, the total lateral load was varied from H = 5y_{f}D ^{3} _{p} to H = 60y_{f}D ^{3} _{p} , where y_{f} is the specific weight of the water . In each case the total load was applied during a time interval of T_{v} = 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 = 15y_{f}D ^{3} _{p} . The response of the pile in elastic soil is also presented.
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.
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 / (y_{f}D_{p} ).The distributions of the dimensionless excess pore pressures at the end of rapid loading for H = 15y_{f}D ^{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.
3.3 Soil-Pile System Under Harmonic Load (^{Abdoun, 1997})
In the centrifuge test reported by (^{Abdoun, 1997}), a single pile model (called model N^{o} 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, D_{r} = 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 D_{r} = 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.m^{2}, and is free at the top. The model has an inclination angle of 2^{o} (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.
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 2^{o}, 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 r_{u} approaching to 1.0 (r_{u} = (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.
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.
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 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}).
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 r_{u} 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 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 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.
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/m^{3}. In order to show liquefaction, the soil mass is considered to be composed of a uniform liquefiable sand layer with a relative density D_{r} = 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.
In Figure 25 is shown the evolution of the liquefaction potential factor r_{u} 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.
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.
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.
4 CONCLUSIONS
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 (r_{u} 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.