See

discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/282133303

Numerical study for soil-spring stiffness of pile

group

Conference Paper · September 2015

CITATIONS READS

0 156

3 authors, including:

Lassaad Hazzar Mourad Karray

Université de Sherbrooke Université de Sherbrooke
26 PUBLICATIONS 46 CITATIONS 109 PUBLICATIONS 145 CITATIONS

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Development of new procedure for liquefaction study in Eastern Canada View project

All content following this page was uploaded by Mourad Karray on 24 September 2015.

The user has requested enhancement of the downloaded file. All in-text references underlined in blue are added to the original document
and are linked to publications on ResearchGate, letting you access and read them immediately.
Numerical study for soil-spring stiffness of pile group
(a) (a) (b)
L. Hazzar , M. Karray & A. Pasic
(a)
Department of Civil Engineering, University of Sherbrook, Sherbrooke,
Québec, Canada
(b)
Dessau Inc. - Bridge and Structures, Québec, Québec, Canada

ABSTRACT
Many approaches have been developed to determine the soil-spring stiffness and typical values have been proposed for
different types and densities of soil. However, these typical values ignore the effect of the depth and the degradation of
elastic parameter (G or E) as a function of strain. In this paper, a series of 3D numerical analyses are conducted to
compute the equivalent soil-spring stiffness’s according to pile distortion for a pile group located under the central pier of
bridge structure above the river Sault-au-Mouton (Longue-Rive, Quebec) and subjected to multi-loading conditions. In
these simulations the degradation of the shear modulus is incorporated in order to account for soil nonlinearity. The idea
is to develop equivalent springs that can be adapted to the lateral deformation of the pile. The stiffness of these springs
thus varies depending on the distortion and may be adapted in an iterative process according to the pile deformation at
each depth. This method can be compared to the linear equivalent method used in dynamic analysis where an
equivalent shear modulus is adapted according to the shear distortion. In fact, This type of linear equivalent spring is
very useful for structural engineers who want to incorporate the effect lateral capacity of soil in their models.

RÉSUMÉ
Des nombreuses approches ont été développées pour déterminer la rigidité du sol-ressort et des valeurs typiques ont
été proposées pour différents types et densités du sol. Toutefois, ces valeurs typiques ne tiennent pas compte de l'effet
de la profondeur et de la dégradation de paramètre élastique (G ou E) en fonction de la déformation. Dans cet article,
une série d’analyse numériques 3D ont été réalisées pour déterminer les rigidités équivalentes des ressorts dans le sol
en fonction de la distortion du pieu pour un groupe de pieux situé sous le pilier central de la structure d’un pont au-
dessus de la rivière Sault-au-Mouton (Longue-Rive) et soumis à plusieurs chargements. La dégradation du module de
cisaillement est prise en compte afin de tenir compte de la non-linéarité du sol. L'idée est de développer des ressorts
équivalents qui peuvent être adaptés à la déformation latérale du pieu. Les rigidités de ces ressorts varient en fonction
de la distorsion et peuvent être adaptés dans un processus itératif selon la déformation du pieu à chaque profondeur.
Cette méthode peut être comparée à la méthode linéaire équivalente utilisée dans l'analyse dynamique où un module de
cisaillement équivalent est calculé en fonction de la déformation de cisaillement. En réalité, Ce type de ressort équivalent
linéaire est très utile pour les ingénieurs en structure qui veulent intégrer l’effet de la capacité latérale du sol dans leur
modèles.

1 INTRODUCTION approaches (Terzaghi, 1955; Vesic, 1961; Broms 1964a,

1964b; Ménard et al., 1969; Davisson 1970; Matlock,
The behavior of laterally loaded piles is a complex soil 1970; Poulos, 1971; Gilbert 1995) have been developed
structure interaction problem that has received a to determine the soil-spring stiffness and typical values
considerable amount of attention over the last four has been proposed for different type of soils and for
decades mainly in the field of offshore engineering or different densities. However, these typical values
earthquake geotechnical engineering due to large stakes generally ignore the effect of depth and the degradation of
involved. The primary function of a pile or a group of piles elastic parameter of the (G or E) as a function of strain.
is to transfer the external loads from the superstructure to In Québec, Broms’s method is still used in the current
the surrounding soil medium without causing excessive design practice of piles under lateral loads (CGS, 2013).
deflections at the top of the pile or piles. The analysis of This empirical method is not entirely appropriate for the
such a problem is complicated because of the complex design of complex structure such as bridge, which needs
stress-strain behavior of the soil surrounding the piles. In sophisticated three dimensional (3D) numerical analysis.
practice, analysis of laterally loaded piles is carried out In this 3D modeling, the soil behavior should be simulated
using beams on non-linear Winkler springs model (often by a series of springs that surround the pile foundation.
known as p-y method) due to its simplicity, low Thus, several questions were asked by the engineers
computational cost and the ability to model layered soils. (essentially structural engineers):
In this approach, soil-pile interaction along the depth is - What are the appropriate values of the horizontal
characterized by a set of discrete non-linear springs stiffness of the springs in the soil around the piles?
represented by p-y curves where p is the pressure on the - What are the values of the stiffness of the soil during
soil that causes a relative pile deformation of y. A number seismic loading?
of different methods are in use to construct the p-y curves - What are the group reduction factors that must apply to
(Matlock, 1970; Cox et al., 1974; Reese et al., 1974; reflect the pile group effects?
Reese et Van Impe, 2001) usually constructed based on In view of the above stated issues and to answer
semi-empirical correlations. Subsequently, many these previous questions, this paper describes and
discusses the results of a series of 3D finite difference FD be variable not only with the soil type and mechanical
3D
analyses using the program FLAC (Itasca 2009). The properties, but also with stress level and the geometry of
3D numerical analyses are carried out to investigate the the pile. In the absence of better information, the
soil-spring stiffness of pile group located under the central coefficient of horizontal subgrade reaction may be
pier of bridge structure above the river Sault-au-Mouton estimated by several methods (Hazzar, 2014).
(Longue-Rive, Quebec) and subjected to multi-loading It is evident that the Winkler model (Fig. 2) cannot fully
conditions (vertical load, lateral loads and moments). capture the 3D aspect of soils. The fundamental
assumption on which the technique of constructing a p-y
2 COEFFICIENT OF SUBGRADE REACTION curve is based is the similarity between the load
deformation pattern of pile head and the stress-strain
The Winkler method, or sometimes known as the behaviour of the interacting soil from carefully chosen
subgrade reaction method, currently appears to be the element testing (e.g., triaxial tests).
most widely used in a design of laterally loaded piles. The
method was first introduced by Winkler (1867) to analyze
the response of beams on an elastic subgrade by
characterizing the soil as a series of independent linearly-
elastic soil springs. Since then, this concept has been
extensively employed for the laterally loaded pile problem.
The concept of this method is graphically illustrated in
Figure 1.

Figure 2. Single pile under lateral loading: (a) real vertical

pile; (b) Winkler idealization; (c) p-y curves for lateral
Winkler springs.

3 GROUP PILES UNDER MULTI-LOADS

3.1 Description of the bridge project

The project consists of a bridge structure above the river

Sault-au-Mouton (Longue-Rive, Quebec). Figure 3 shows
the location of the study site.
Figure 1. Implementation of Winkler spring concept for
laterally loaded pile problem.

The term of subgrade reaction indicates the lateral

pressure, p, per unit area of the surface of the contact
between a loaded beam or slab and the subgrade on
which it rests and on to which it transfers the loads. The
coefficient of subgrade reaction, k, is the ratio between
the soil lateral pressure, p, at any given point of the
surface of contact and the lateral deflection, y, produced
by the load application at that point:
Figure 3. Site location.
p
k [1]
y The bridge deck contains two spans, each 75 m in
length and 14.0 m in width. According to the architectural
The soil-spring stiffness or the modulus of subgrade plan of the bridge, the footings will be positioned at
2
reaction, ES (kN/m ) is related to the pile diameter, B by elevations following:
the following equation: - West Abutment: top of the footing directly below the
ground surface at the elevation of 37.78 m;
- Central pier: above the sole directly under the surface of
ES  k  B [2]
the riverbed elevation at 27.265 m;
- Abutment: above the sole to 38.495 m elevation in the
Though simple in its definition, the modulus of approach embankment.
subgrade reaction has proved to be a very difficult The central pier is supported by a footing that is
parameter to evaluate. This is because it cannot be supported by a group of 6 drilled shaft piles (figure 4).
measured in laboratory tests, but must be back-calculated
from full-scale field tests. Investigations have shown it to
 Figure 4. Bridge P-15705 in Longue-Rive, 3D modeling.

3.2 Soil properties

Table 2. Values of applied forces.
A total of eight boreholes, identified F25 to F32, were Case n° 1 2 3
made in 2004 to investigate the soil conditions and to Lateral load, Px (kN) 123 2.280 7.010
determine the roc position. The soil profile is constituted
by two sandy layers with groundwater at the surface. Axial load, Py (kN) 39.152 28.100 27.800
Table 1 summarizes the properties of these soil layers. Lateral load (kN) 1.950 3.150 1.340
Bending moment, Mx 30.000 53.010 22.200
Table 1. Soil properties after geotechnical tests (kN.m)
Depth (m) 0.0-6.0 6.0-20.0 Torque, My (MN.m) 1.800 120 350
Unit weight, w (kN/m )
3 Bending moment, Mz 1.900 37.100 113.500
21.00 21.00
(kN.m)
Angle of friction,  (°) 35.00 40.00
Shear velocity, Vs (m/s) 200.00 230.00 4 FINITE DIFFERENCES MODELLING
Poisson’ ratio,  0.33 0.33
4.1 Finite differences mesh and boundary conditions
3.3 Piles properties 3D
The 3D FD program FLAC (Itasca, 2009) was employed
to study the behavior of a single pile and pile groups
A total of 6 drilled shaft piles are installed. The piles have
(lateral displacement, lateral resistance and stiffness of
a diameter B of 2.0 m and a length L of 16 m.
each pile) under several loads. Figures 5 and 6 show 3D
Configuration of 2  3 piles has been proposed with
finite difference grids used in the numerical analyses.
spacing of 5.0 m  10.0 m. the mechanical characteristics A mesh generator subroutine was implemented using
of the pile are: a Poisson’s ratio p of 0.20 and a modulus the FISH built-in programming language providing the
of elasticity Ep of 25 GPa. possibility of mesh refinement and geometry variation.
The bottom elevation and the lateral sides of the
3.4 Types of loads applied to the bridge computational domain were taken far enough from the
group to avoid any significant boundary effect. Based on
The bridge is located in seismic performance zone 3. the experience gained through this and previous
Referring to the specification of the CAN/CSA-S6-06 numerical works, a mesh refinement around the piles
(2006) code, this structure is classified in the category of leads to a more accurate distribution of stresses and soil
lifeline bridge. Subsequently, a seismic analysis with yielding. In order to make sure that the zones size has no
multi-mode spectral method is considered and 3 loading effect on the response of the characteristic piles, trial
cases have been analyzed (Table 2): analyses have been carried out to optimise mesh
discretization. More specifically, the response of each
- Case n°1: applied forces without taking into account characteristic pile has been established carrying out an
the seismic excitation; analysis in which only this pile was activated. The
- Case n°2: applied forces taking into account a 100% responses of these piles have then been compared and
transversal seismic excitation and a 30% longitudinal re-adjustment of the mesh was done if necessary.
seismic excitation ;
- Case n°3: applied forces taking into account the 100%
longitudinal seismic excitation and the 30%
transversal seismic excitation;
 Table 3. Soil parameters according to M.C.
Depth  Gmax Kmax cu(kPa) 
(m) (kg/m³) (MPa) (MPa)
0.0-6.0 2100 5.40 11.70 0 35
6.0-20.0 2100 12.60 27.30 0 40

Standard elastic/plastic constitutive laws such as

Mohr-Coulomb (M.C.) can also produce the reduction of
shear-modulus that can be evaluated from degradation
curves as follows: consider the M.C. model with a
Figure 5. Three dimensional single pile model. constant shear modulus, Gmax, and a constant yield
stress, m, corresponding to cyclic shear strain of
amplitude, . In pre failure phase, the secant shear
modulus, G, is assumed to be equal to Gmax. For a cyclic
excitation that involves failure, the secant modulus is
given by:

m
G [4]


The shear modulus-reduction curve relates the ratio

G/Gmax to the amplitude of shear strain,  ; it is simply
obtained from Eq. 5 by setting m = m /Gmax , one obtains,
for >m,

G 
Figure 6. Finite differences mesh for a 53 pile group: (a)  m [5]
Gmax 
Global model and (b) pile group with footing.

4.2 Soil model For numerical simulation, the recognised shear-

modulus degradation curve proposed by Seed and Idriss
The Mohr-Coulomb’s law (M.C.) constitutive relation was (1970) for sands soils is adopted (Figure 7).
used to model soil layers. In fact, it allows us to define
shear failure of a soil mass based on its state of stress.
This elastic perfectly-plastic model requires five basic
input parameters, namely a maximum elastic bulk
modulus “Kmax”, maximum elastic shear modulus “Gmax”,
mass density “”, undrained shear strength “cu” and the
friction angle “”.
For sandy soils, Gmax, is estimated from the shear
wave velocity, Vs, and the mass density, , by the
correlation:

Gmax  Vs2 [2]

The maximum elastic bulk modulus, Kmax is related to

Gmax and the Poisson’s ratio,  by:

2 1   
Figure 7. Adopted curve of shear modulus degradation for
K max  Gmax [3] sandy soils (Seed and Idriss, 1970).
 3  6 
4.3 Pile model
Therefore, the material properties adopted in the
analyses for the two layers (referred to the properties The pile is modelled as linear-elastic material. Three
deduced from geotechnical tests) are presented in Table parameters are required to define the pile material
3. behavior. These parameters are the elastic bulk modulus,
Kp, the elastic shear modulus, Gp, and the mass density,
p.
4.4 Soil-pile interface model deflection of the pile, y, first the coordinates of all the
nodes of the vertical pile axis were stored in the memory
The interface elements are modeled by the linear of the computer. Once the model reached equilibrium, the
Coulomb shear-strength criterion that limits the shear same nodes were identified and their displacements
force acting at an interface node. The shear-strength calculated by subtracting their initial and final horizontal
criterion that limits the shear force acting at an interface coordinate.
node for sandy soils is given by Eq. (6).
4.6 Soil lateral pressure
Fs max   Fn  p1 A tgi [6]
The soil lateral pressure, p, can also be calculated by
summing the forces in the relevant direction acting on the
Where Fsmax is the limiting shear force at pile–soil
soil-pile interface nodes at the same depth. A schematic
interface, Fn is the normal force, i is the angle of friction of the pile-soil system is presented in Fig. 8. Each
of the interface surface, p1 is pore pressure (interpolated interface node is associated with a normal force and a
from the target face) and A is the contact area between shear force (Hazzar, 2014).
pile and soil.
The value of friction angle of the interface surface
corresponds to critical state and is reduced compared to
the friction angle of the surrounding soil (Ortigao, 1995).
Separation is able to cause a significant increase in
displacements (Poulos and Davis, 1980) and therefore the
interface elements are allowed to separate if tension
develops across the interface and exceeds the tension
limit of the interface. Once gap is formed between the
pile-soil interfaces, the shear and normal forces are set to
zero.
The normal and shear forces at the interface nodes
are determined by the following equations:

Fn(t t )  knun A   n A [7]

F ( t t )
si F (t )
si
( t  0.5 t )
 ks u
si   si A [8] Figure 8. Schematic plot of pile-soil system with interface
between them (Hazzar, 2014).
Where Fn and Fsi are the normal and shear force,
respectively, kn and ks the normal and shear stiffness, The x-component of the lateral pile-soil pressure is
respectively,usi the incremental relative shear summed over all the interface nodes to calculate the
displacement vector, un the absolute normal penetration lateral soil pressure, p per unit length along the pile at a
of the interface node into the target face, n the additional particular pile section and is expressed as:
normal stress added due to interface stress initialization,

p    i nx ,i   xy ,i ny ,i   Ai 
n
and si the additional shear stress vector due to interface [9]
stress initialization. i 1
In many cases, particularly when linear elastic analysis
is performed, values for interface stiffness are assigned to Where,
simulate the nonlinear behavior of a physical system. In xi
the present study, where nonlinear analysis is carried out, nx ,i  cos  
the value for the interface stiffness should be high xi2  yi2
enough, in comparison with the surrounding soil, in order yi
to minimize the contribution of those elements to the n y ,i  sin  
accumulated displacements. According to the results of x  yi2
2
i
9
numerical analyses a value of 10 MPa/m for both kn and σi : normal stress at the interface node at point Q;
ks was sufficient to ensure that no additional deflections τxy,i : shear stress at the interface node at point Q;
were attributed to the pile due to the deformation of the xi : x-coordinate of the interface node at point Q;
springs representing the interface. The use of yi : y-coordinate of the interface node at point Q;
considerably higher values is tempting as it could be Ai : representative area of interface node.
considered as more appropriate, but in that case the
solution convergence would be very slow. 5 NUMERICAL RESULTS AND DISCUSSIONS

4.5 Lateral deflection of pile Based on the constitutive model parameters described
previously, the response of the laterally loaded single pile
When a pile is laterally loaded, the nodes move along the is presented in terms of p-y curves stiffness of equivalent
direction of applied load in large strain mode. Therefore, it soil-spring versus pile distortion. For group piles, the
is not possible to identify a particular node after loading response is presented in terms of pile group effects.
using its original coordinate. To calculate the lateral
5.1 Stiffness of equivalent soil-spring
According to the definition of modulus of subgrade
The most widely used nonlinear analysis for laterally reaction or spring-soil stiffness in section 2, the variation
loaded piles is the p-y curves. In this paper p-y curves are of the stiffness of equivalent springs, with the distortion of
obtained based the methods of prediction of lateral pile are plotted for at several depths in Fig. 10. Figure 10
deflection, y and lateral pressure, p, described in sections indicates that the stiffness’s of spring-soils are not
4.5 and 4.6 respectively. constants as they have already mentioned the most
Figure 9 shows numerical p-y at seven depths of the methods considered in practice (Hazzar, 2014). Table 4
pile. It is clear that the lateral pressure of soil increases as shows the spring-soil stiffness predicted by the current
the depth does. This can be explained by a decrease in numerical analysis and those given by several methods
lateral deflection with depth. adopted in the practice. According to Table 4, these
methods give reasonable values but ignore the effect of
depth and the pile distortion.

5.2 Piles group effects

The lateral response (profiles of lateral deflection and

lateral pressure) of the pile group 2×3 are studied with
3D
FLAC for the several loads conditions (Table 3). The
results are compared with the lateral response of single
pile. Consequently group reduction factors are determined
for the equivalent spring-soil stiffness along the depth
(Table 4). It can be seen that these reduction factors
increase with depth and its values depend of the location
of each pile in the group.

Figure 9. p-y curves between 4 and 10 m of depth.

Figure 10. Stiffness of equivalent spring-soil versus distortion of pile between 4 and 11 m of depth.

Table 4. Modulus of subgrade reaction, ES (kN/mm): Comparisons between the current analysis and several methods.
Depth (m) Current Broms Ménard et Poulos Gilbert
analysis (1964b) al. (1969) (1971) (1995)
0.0-6.0 4 to 33 22.0 25.0 11.8 26.96
6.0-20.0 2 to 40 22.0 32.0 27.5 30.9

6 CONCLUSIONS Seed and Idriss (1970) to adjust the value of shear

modulus of soil is considered.
For prediction the spring-soil stiffness of piles group under The results obtained are interpreted rationally to
multi-loads (lateral load, axial load and torque), a series of conclude, initially, that the stiffness’s of equivalent springs
rigorous 3D numerical analysis based on finite difference vary depending on the pile distortion and may be adapted
technique has been performed. The Mohr Coulomb in an iterative process according to the pile deformation at
criterion is used to model the soil parameters and each depth. This process can be compared to the linear
recognised shear-modulus degradation curve given by equivalent method used in dynamic analysis where an
equivalent shear modulus is adapted according to the
 shear distortion. This type of linear equivalent spring is Ménard, L., Bourdon, G., Gambin, M. 1969. Méthode
very useful for structural engineers who want to générale de calcul d'un rideau ou d'un pieu sollicité
incorporate the effect of lateral capacity of soil in their latéralement en fonction des résultats
model. Secondly, the reduction factors to reflect the group pressiométriques, Sols–Soils, 22-23: 16-29.
effects depends of the depth and the location of each pile Matlock, H. 1970. Correlations for design of laterally
in the group (consequently the spacing between the loaded piles in soft clay. Journal Soil Mechanics and
piles). Foundations, ASCE. SM5: 63–91.
Ortigao, J.A. 1995. Soil mechanics in the light of critical
Table 4. Factors of reduction to reflect the group effects. state theories. Rotterdam: Balkema.
Depth Factor of Factor of Factor of Poulos, H.G. 1971. Behaviour of laterally loaded piles: I-
(m) reduction / reduction / reduction / single piles. Journal Soil Mechanics and Foundation,
longitudinal transversal transversal ASCE. 5(97): 711–31.
direction direction / direction Poulos, H.G., Davis, E.H. 1980. Pile foundation analysis
middle piles /back piles and design. New York: John Wiley & Sons Ltd.
4 0.49 0.13 0.45 Reese, L. C., Cox, W.R., Kooper, F.D. 1974. Analysis of
th
5 0.48 0.14 0.38 laterally loaded piles in sands. Proceedings of the 6
6 0.51 0.10 0.39 Annual Offshore Technology Conference, Houston,
7 0.51 0.09 0.42 Texas. OTC 2080: 437-483.
8 0.52 0.12 0.49 Reese, L.C., Van Impe, W.F. 2001. Single piles and pile
9 0.53 0.21 0.62 groups under lateral loading. A.A. Balkema,
10 0.53 0.34 0.70 Rotterdam. 508 pages.
11 0.53 0.51 0.90 Seed, H.B., Idriss, I.M. 1970. Soil moduli and damping
12 0.58 0.73 1 factor for dynamic response analyses. Report No.
13 0.72 0.95 1 EERC 70-10, Earthquake Engineering Research
14 0.88 >1 >1 Center, University of California, Brekley, California.
15 >1 >1 >1 Terzaghi, K. 1955. Evaluation of coefficients of subgrade
16 >1 >1 >1 reaction. Géotechnique. 5: 297-326.
17 >1 >1 >1 Vesic, A. 1961. Design of pile foundations. National
18 >1 >1 >1 cooperative highway research program synthesis of
19 >1 >1 >1 highway practice. Report No. 42, Transportation
Research Board, Washington, D.C., 68 pages.
Winkler, E. 1867. Die lehre von der elasticita et und
REFERENCES
festigkeit. Prag, Dominicus
Broms, B.B. 1964a. Lateral resistance of piles in
cohesive soils, Journal of Soil Mechanics. Foundation
Division, ASCE. 90(2): 27-64.
Broms, B.B. 1964b. Lateral resistance of piles in
cohesionless soils, Journal of Soil Mechanics.
Foundation Division, ASCE. 90(SM3): 123-156.
Cox, W. R., Reese, L.C., Grubbs B.R. 1974. Field testing
th
of laterally loaded piles in sand, Proceedings of the 6
Annual Offshore Technology Conference, Houston,
Texas. OTC 2079, 14 pages.
CAN/CSA-S6-06. 2006. Canadian Highway Bridge Design
Code, CSA group, 930 pages.
CGS. 2013. Canadian Foundation Engineering Manual
th
(4 Ed.), Richmond, B.C., Canadian Geotechnical
Society.
Davisson, M.T. 1970. Lateral load capacity of piles,
Highway Research. Record 333: 104-112.
Hazzar, L. 2014. Analyse numérique de la réponse des
pieux sous sollicitations latérales, PhD. Université de
Sherbrooke (Québec-Canada), 214 pages.
Itasca Consulting Group. 2009. FLAC3D: Fast lagrangian
analysis of continua in 3-dimensions 4.1, manual.
Itasca, Minneapolis.
Gilbert, C. 1995. Une nouvelle approche des calculs
d'interaction sol-structure, Révue Française de
Géotechnique. 7: 3-9.

View publication stats