MAT EC Web of Conferences 16, 020 05 (2014)

DOI: 10.1051/matecconf/ 201 4 16 020 05


C Owned by the authors, published by EDP Sciences, 2014

Axial Non-linear Dynamic Soil-Pile Interaction - Keynote

A. Holeyman1 and V. Whenham2
1
Université catholique de Louvain, iMMC, GeoMEM, 1348 Louvain-la-Neuve, Belgium
2
Fugro Geoconsulting Belgium, 1150 Brussels, Belgium

Abstract. This keynote lecture describes recent analytical and numerical advances in the modeling of the
axial nonlinear dynamic interaction between a single pile and its embedding soil. On one hand, analytical
solutions are developed for assessing the nonlinear axial dynamic response of the shaft of a pile subjected to
dynamic loads, and in particular to vibratory loads. Radial inhomogeneity arising from shear modulus
degradation is evaluated over a range of parameters and compared with those obtained by other authors and by
a numerical radial discrete model simulating the pile and soil movements from integration of the laws of
motion. New approximate non linear solutions for axial pile shaft behaviour developed from general
elastodynamic equations are presented and compared to existing linear solutions. The soil non linear behaviour
and its ability to dissipate mechanical energy upon cyclic loading are shown to have a significant influence on
the mechanical impedance provided by the surrounding soil against pile shaft movement. The limitations of
over-simplified modelling of pile response are highlighted.

1 Introduction

1.1 Typical situations

Piles are used to support civil engineering structures

whenever loads are concentrated relative to the soil
bearing capacity and thus strength. They consist of
elongated structural elements that are embedded in the
soil to a depth that will allow bearing layers to develop a
safe resistance. Typical situations where piles are used
are illustrated on figure 1. In the most common case
referred to as “active” pile, the pile head receives the load
from the superstructure and transmits it to the resisting
soil.
Although all 6 load components (forces components
and moments about 3 orthogonal axes) need to be Fig.1. Typical use of piles to support wind turbines: Tripod,
considered, the pile axial component generally governs Jacket, and Monopile [1]
the design of the pile and in particular its embedment
depth . In addition to static loads resulting from gravity, harmonic axial load. Similar contributions can be
operational transient or periodic loads may warrant developed on the pile lateral response, torsional response
special design requirements. In response to the quickly and would justify a discussion on the coupling. Finally,
expanding market of offshore renewable energy [1], piles most of the discussion will revolve about the friction
have been recently used in several configurations to component of the pile resistance, leaving out the
support wind turbines, as illustrated in figure 1. singularity brought about by the pile bottom end.
Extreme loading cases have also to be considered due Whatever the mode considered, the engineering
to the fact that piles have to be installed at depth: such picture needs to involve three main ingredients: the
cases involve pile impact driving as well as vibratory complete definition of all actions on the pile, the
driving. structural characterization of the pile, and the
Due to space and time constraints, this lecture will characterization of the embedding medium, i. e. the soil.
only focus on the axial response of a pile subjected to a

This is an Open Access article distributed under the terms of the Creative Commons Attribution License 3.0, which permits unrestricted use,
distribution, and reproduction in any medium, provided the original work is properly cited.
Article available at http://www.matec-conferences.org or http://dx.doi.org/10.1051/matecconf/20141602005
 MATEC Web of Conferences

1.2 Pile characterization resistance versus the pile base displacement. figure 2b
schematizes such simplified modelling of the soil-pile
Piles used for large projects are commonly made of interaction.
concrete or steel. Generally concrete is used for onshore By extending Winkler approach initially developed in
applications while steel is used for offshore applications. the static domain, the dynamic response of pile shafts
In the later case, steel pipes are used to construct “pipe embedded in an elastic medium can be studied by
piles”. The advantage of that geometry is that the volume replacing the soil surrounding the pile with a series of
of soil to be displaced to accommodate pile insertion into independent springs and dashpots. Notably, Smith [2]
the soil is a very small fraction of the gross pile volume used that approach to model soil resistance to pile
that will govern the geotechnical capacity of the pile. For driving, leading to the emergence of soil-pile dynamic
large diameter piles, the inertia of the soil volume inside interaction parameters known as “quake” and “damping”,
the pipe will prevent the soil core from moving down in as illustrated as Q and J on figure 3.
unison with the pile during driving, leading to what is
termed a “coring” mode of driving.
Most tubular piles are driven by impact using special
pieces of equipment (hydraulic pile hammers) set on top
of the pile during installation. More rarely, piles can be
vibrated into the ground using vibratory hammers. Once
installed to an appropriate depth, a pile will develop its
bearing capacity over time, moving from its end-of-
installation capacity towards its long-term capacity. The
pile geotechnical capacity will come from contact
stresses generated by the soil along two interfaces: a
shear stress along the pile lateral surface (called shaft)
and a normal stress against the pile end bearing area
(called toe).
Fig. 2. (a) Global load-settlement curve at top of pile
The axial bearing performance of a pile can be (b) Embedded pile as continuously supported column
characterized by its response under axial static loading.
As illustrated in figure 2a, the pile load-settlement curve
provides the overall relationship between the applied load
F and the pile settlement s. One can notice that under
limited loads, the pile responds linearly, but endures non
recoverable displacements under larger loads, with the
ultimate limit state being defined by unlimited
displacements under an asymptotic load. This reflects the
non linear behaviour of the soil while the pile material
remains well within its elastic realm. This forces us to
address an essential feature of the system, namely the soil
behaviour that will be characterized in Section 2.

Fig. 3. Smith visco-elasto-plastic load-deformation curve for

1.3 Simplified Soil-Pile interaction local slice of pile shaft (adapted from [2])

Engineering methods that are commonly applied to assess

For a pile undergoing harmonic axial motion,
pile response under a static load applied at the head of the
coefficients of the Winkler springs and dashpots depend
pile treat the pile as a column collecting reactions along
on the frequency. Such coefficients can be obtained by
its shaft and at the pile toe. Assuming that the embedding
considering the elastodynamic problem of an infinitely
medium is elastic, it would appear possible to model each
long pile subjected to harmonic vertical displacements.
local reaction as proportional to the local displacement.
An alternative approach consists in modelling fully
Further assuming that each equivalent spring is
coupled 3D pile-soil interactions, for example by means
independent from its neighbours allows one to treat each
of the finite element method. However because of the
soil layer with its own properties.
complexity of the problem, especially when the pile is
This simplification known as a “Winkler” model has been
subjected to high strain loading conditions, this approach
extended to cope with non linear behaviour, leading to
makes it difficult to properly simulate the essential
what is known in the geotechnical jargon as “t-z” and “Q-
phenomena at play and is faced with the challenge to
z” curves. A “t-z” curve models the non linear
harness adequate model parameters. Practical use of the
development of the local shear stress along the pile shaft
full 3D finite element approach is further hampered by its
versus local vertical displacement while a “Q-z” curve
high demand for computer resources.
models the non linear development of the pile toe
While structural engineers are quite fond of
simplifications of the pile behaviour boiling it down to a

02005-p.2
 CSNDD 2014

single spring and sometimes a dashpot, they tend to where e is the soil void ratio and σ ' is the effective
overlook that a pile is a complex infrastructural system
that interacts with the soil surrounding and with the confining mean stress usually calculated as:
superstructure. Furthermore soil is a medium that is far σ v' + 2σ r' , 0
from behaving linearly, as can be summarized below. σ = '
(2)
3
Where σ v' is the easily calculated effective initial
2 Soil characterization vertical stress, σ r' , 0 = kσ v' is the effective horizontal

Soil is a multiphase medium made of solid particles stress, with k being the coefficient of horizontal stress in
whose composite behaviour depends on many factors: the soil. For a “wished-in-place” pile assumption, an at
attributes of particles, fluid filling the voids left between rest coefficient is estimated herein using Jacky’s formula
the solids, stress history, just to name a few. While the k = 1 − sin ϕ ' wherein ϕ ' is the soil internal friction
particles can be characterized by their nature, size, and angle.
shape, their overall behaviour with respect to the water
content can be characterized by their plasticity index (PI).
The PI of sand is zero while clay minerals can exhibit 2.1.2 Strain hardening
values in excess of 50, with silts having intermediate
Thanks to numerous forms of soil testing, the relationship
values.
between shear stress and shear strain has been shown to
Volume variations and distortion of soil depend solely
deviate from the initial tangent value Gmax defined above
on the soil “effective” stress, i.e. the stress between
as shear strain increases, as shown in figure 4. This can
particles, the pore pressure having no intrinsic role other
be viewed as “strain hardening” since the shear stress
than taking a part of the total stress. One rather unique
increases beyond an “elastic” limit to be identified. In soil
feature of soil behaviour is its tendency to change volume
mechanics, this feature is preferably described in terms of
when sheared. Loosely packed soils tend to contract
shear modulus degradation [4] inasmuch the secant
while densely packed soils tend to dilate. Soil shear
modulus Gs degrades with strain. Two examples of
strength essentially comes from friction, which is
models commonly accepted to characterize the shear
controlled by effective stress while the effective internal
modulus degradation are discussed in Sections 2.2 and
friction angle generally assumes a value close to 30°.
2.3. Densely packed soils, such as stiff clays can also
When saturated with water, a contractant assemblage
exhibit some form of strain softening, which will not be
of particles can only modify its volume inasmuch water
covered in this lecture.
has the time to drain away from it. This means that low
At large strains, the soil reaches a “critical state”
permeability soils or soils undergoing fast loading have to
shear strength characterized by a constant volume and
deform without changing their volume, which implies a
mean stress. That ultimate limit state will be
substantial change in their effective stress. Such a
volumetric constraint explains why the strength of a soil characterized by the value of τ max shown on figure 4.
depends so drastically upon its loading rate. At one
extreme, loading is so slow that volume changes can be
accommodated without any interference from the pore
fluid, and the soil will behave as “drained”. At the other
extreme, loading is so fast that volume remains constant
and the material will behave as undrained.
Whether drained or undrained, the latter case being
more common under dynamic loading, soil beahvior
exhibits several features that are characterized in the
following section, namely, stiffness, strain hardening and
yield criteria, implying material damping upon cyclic
loading.

2.1. Key attributes of soil behavior Fig. 4. Hyperbolic shear stress – shear strain model for soils [4]

2.1.1 Small strain stiffness 2.1.3 Material damping and viscous equivalent

An initial (or maximum) shear modulus can be calculated According to Masing rules [5] which will be deemed
for rounded grained sands using the relation of Hardin applicable, if the loading curve is characterized by the
and Black [3]: relationship f ( F , δ ) = 0 with F the load and δ the
(2.17 − e)² ' 0.5 (δ c , Fc )
Gmax = 6908 σ [kPa] (1) displacement, then for a cycle between points
1+ e and ( −δ c ,− Fc ) , the loading-unloading curves are:

02005-p.3
 MATEC Web of Conferences

δ − δ c* F − Fc* 2.1.4 Cyclic degradation and liquefaction

f( , )=0 (3)
2 2 Beyond strain related shear modulus “degradation”
discussed above, soil is subject to fatigue whenever strain
where (δ c , Fc ) is the point of loading inversion. figure
* *
cycles of sufficient amplitude accumulate. In loose
5 provides an illustration of Masing’s rules, showing how granular soils that are saturated, the cyclic degradation
the backbone loading curve can be expanded and rotated can be compounded by the onset of “liquefaction”. This
to generate the unloading and reloading curves. phenomenon occurring in contractive materials involves
The loop developed within this stress-strain plane the increase of the pore pressure to the point that effective
highlights the dissipation of mechanical energy during a stresses vanish. Since soil strength is intrinsically related
complete loading cycle. Such a material dissipation is to friction, the removal of any effective stress actually
characterized by the soil “damping coefficient”, a relative transforms the soil into a medium unable to resist to shear
measure of the dissipated energy ΔW within one cycle to stress, i.e. a fluid.
the maximum accumulated elastic energy. If τc is the
amplitude of the shear stress and γ c the amplitude of the
2.2. Hardin and Drnevich [4] model
shear strain, the soil damping coefficient is defined as
ξ = ΔW /(4πW ) , with W = γ c .τ c /2 . The relationship between shear stress and shear strain
under undrained conditions can be assumed to follow the
It should be noted that for a given τ max the material soil model suggested by Hardin & Drnevich [4] and
damping depends on γc but not on the frequency of a based on Kondner [6] formulation.
potentially considered harmonic movement. Harmonic
γ γ
displacements prescribed by the pile generate cyclic τ= = τ max . (4)
deformations and stresses within the influenced soil zone 1 / Gmax + γ / τ max γ +γr
that can be conveniently expressed by:
τ = τ c .eiωt = G.(1 + 2iξ ).γ c .eiωt = G*.γ c .eiωt where τ γ is the shear strain, τ max is
is the shear stress,
*
wherein G is the complex shear modulus. the maximal shear stress (shear strength), Gmax is the
initial shear modulus and γ r = τ max /Gmax is the
reference strain. Degradation of the secant shear modulus
with the shear strain can therefore be defined by:

τ τ max γ § τ ·
G= = = Gmax. r = Gmax.¨¨1− ¸¸ (5)
γ γ +γr γ +γr © τ max ¹

Applying Masing’s rules to the Hardin & Drnevich [4]

loading curve leads to the following expression of the
Fig. 5. Masing’s rules [5] applied to stress-strain curves damping coefficient as a function of γ c :

Such an expression postulates that energy losses can be

attributed to an out of phase term treating the soil as 2§ γ § γ · γ ·
ξ = .¨¨ 2 r2 .(γ r + γ c ).ln¨¨ r ¸¸ + 2 r + 1¸¸ (6)
visco-elastic. In that case, a so called “linear equivalent” π © γc ©γr +γc ¹ γc ¹
soil model is invoked, although strictly speaking, it is
linear only under static conditions.
Assuming hysteretic energy losses can be handled at Illustrations of the Hardin & Drnevich [4] formulation
a given frequency by an equivalent viscous-elastic shear and Masing’s rules [5] are provided in figures 4 to 6.
modulus G = G.(1 + 2iξ ) , the equivalent soil viscosity
*

can be expressed under harmonic conditions at frequency

2.3. Ishibashi and Zhang [7] model
ω as η = ∂τ ∂γ = 2Gξ / ω , with γ being the shear
strain rate. Although it does not respect the shape of the Alternatively, the stress-strain relation suggested by
stress-strain loop, the assumption of the “linear Ishibashi and Zhang [7] can be considered to characterize
equivalent” medium is generally accepted because of its soil nonlinearity, especially in the weakened zone close
mathematical convenience, as will be shown in Section to the pile shaft. The degradation of the secant shear
3.3. modulus is expressed as a function of the shear strain γ ,
of the effective confining mean stress σ' , and of the
plasticity index IP:

02005-p.4
 CSNDD 2014

G ­ 0.000102 + n ½ 'm a representative number of cycles and degradability into

= 0.5®1 + tanh(0.492 ln )¾σ (7a) account. That refinement can however be explicitly
Gmax ¯ γ ¿ accounted for in the numerical models developed by
where Holeyman [9], where several degradation laws of
­ 0.000556 ½ −0.0145 IP 1.3
(7b) τ max are implemented according to the local shear strain
m = 0.272®1 − tanh(0.4 ln )¾e
¯ γ ¿ history.
­0.0 for IP = 0
° −6
°3.37 × 10 IP for 0 < IP ≤ 15
1.404
(7c)
n=® −7 2.5. More advanced models
°7 × 10 IP for 15 < IP ≤ 70
1.976

°2.7 × 10 IP
− 5 1 .115
for IP > 70 Many more models have been developed by researchers
¯
attempting to capture various features of the complex soil
The hysteretic damping coefficient resulting behaviour. These can be approached by separating
from equations 7a, b, c can be expressed as [7]: recoverable and non recoverable deformations, the latter
1.3
1 + e −0.0145IP G G being handled through plastic theory. Three-dimensional
ξ= [0.586( )² − 1.547 + 1] (8)
representations of the yield function in the stress space
6 Gmax Gmax
and the choice of flow rules in the strain space are then
The asymptotic value of ξ depends on the plasticity necessary, requiring the knowledge of up to tens of
1.3 parameters that are difficult to determine experimentally
1 + e −0.0145 IP for many engineering projects.
index ( ξ IP =
max
), reaching a maximum of The separation between elastic and plastic domains
6
33% for IP=0 and 18.3% for IP=50. This feature makes can be circumvented by the use of so-called
the Ishibashi and Zhang’s model more flexible than the “Hypoplasticity”, which appears to gain popularity
basic hyperbolic model [4] and applicable to several thanks to a more reasonable number of parameters. While
types of soils according to their plasticity attributes. the initial hypoplastic model [10] required 8 parameters,
those incorporating the intergranular concept and cyclic
features [11] can go up to 13. That number contrasts with
the more manageable 2 or 3 parameters necessary to
2.4. Experimental evidence and complicating understand what are believed to be the essentials of pile
factors response used in the remainder of this invited lecture.
Figure 6 presents a comparison between experimental Moreover the basic parameters used (e and PI or Gmax
curves established by Vucetic and Dobry [8] for soils of and τ max ) are standard geotechnical parameters widely
varying PI and the above theoretical formulations of the
available from customary site characterization.
damping coefficient (6 and 8), for various reference
strains. It can be noted that experimental values of ξ can
typically range between 0 and 0.4.
In soils that are subject to cyclic degradation, such as 3 Problem statement and linear solution
loose sands or sensitive clays, the maximum shear stress
(shear strength) τ max evolves as cycles accumulate. Such 3.1 Idealized conditions
an evolution is not explicitly accounted for in the
analytical models presented hereafter; rather it is The problem considered in this keynote lecture involves a
accounted in the choice of an equivalent τ max that takes vertical cylindrical floating pile shaft of infinite length
and rigidity, embedded within an infinite homogeneous
soil medium. The pile shaft is subjected to a purely
harmonic axial displacement prescribed by
w0 = w0c . cosωt where w0c is the amplitude of
displacement of the pile shaft, and ω is the circular
frequency, and t the time. The examined system is a unit
slice of the problem as shown on figure 7, isolating a
single pile shaft segment and associated unit thickness
soil layer of infinite radial extent. Plane strain conditions
prevail across any horizontal slice because of the infinite
extent of the considered problem in the axial direction
and uniformity of the prescribed movement along the
vertical direction. The layer outside the pile can be
viewed as an infinite shear plate with a circular hole
Fig. 6. Soil damping coefficient ξ as a function of cyclic shear about which the harmonic vertical motion is prescribed.
strain γc: experimental curves vs. models [4] and [7] The prescribed dynamic displacement generates
cyclic deformations and stresses within the analyzed soil

02005-p.5
 MATEC Web of Conferences

∂ ² w § G * ∂G * · ∂w ∂²w
G *
+ ¨¨ + ¸¸. = ρ. (10)
∂r ² © r ∂r ¹ ∂r ∂t ²

Extending the Winkler concept introduced in Section 1.3,

the dynamic reaction of the soil surrounding the pile shaft
may be expressed with reference to an equivalent spring-
dashpot system anchored to a stationary point, as
illustrated by figure 8. The dynamic soil reaction
opposing the prescribed pile shaft movement can then be
expressed as

Pz ( r, t ) = Cz .w ( r, t ) + K z .w( r, t ) (11)

Fig. 7. Unit layer considered within infinite pile shaft and withPz the soil reaction per unit length of shaft, Cz the
embedding soil
damping coefficient and K z the stiffness coefficient.
layer shown on figure 7 that can be represented Since the problem has been stated within a unit thickness
by τ c . cos(ωt − θ ) where τ c is the stress at radial soil layer, it should be noted that Pz , C z , and K z are
distance r and θ is the phase difference with respect to expressed per unit length of pile shaft, and thus typically
the displacement applied at the soil-shaft interface. in the following respective units [ kN / m ] , [kPa / s ] ,
In practice, stress anisotropy induced due to the and [kPa ] .
weight of the soil will result in a specific distribution of
the shear modulus with depth. Furthermore, soil layering Assuming Pz ( r, t ) is harmonic, we can define the
is not homogeneous as the pile can endure axial unit (lineal) shear impedance of the soil against the pile
compression, making the infinite extent of the pile and shaft movement in the z direction as:
surrounding soil assumption less legitimate. Averaging of
soil properties along the depth of the pile shaft should be Pzc
considered prior to using a single layer model. Iz = = (C z .i.ω + K z )
Relationship between shear stress and shear strain wc
under undrained conditions will be assumed to follow the (12)
2π .r0
soil model suggested by Hardin & Drnevich [4] and = G s 0 (C za .i + K za ) = .τ c ( r0 )
based on Kondner [6] formulation. wc (r0 )
with Gs 0 the shear modulus at the pile shaft-soil
3.2 Soil impedance to pile shaft movement interface ( r = r0 ) , and

K za = ℜ{I z }/ Gs 0
Under the assumption of small deformations and absence
(13)
of slippage at the pile shaft-soil interface, and provided
radial deformations as well as the pile mass effect can be Cza = ℑ{I z }/(ω.Gs 0 )
neglected in the analysis, the differential equation
describing the vertical motion w(r,t) of a floating rigid the dimensionless stiffness and damping parameters.
pile shaft embedded in a homogeneous isotropic elastic
soil medium of shear modulus G and volumetric mass
ρ is given by:

∂ ² w § G ∂G · ∂w ∂²w
G +¨ + ¸. = ρ. (9)
∂r ² © r ∂r ¹ ∂r ∂t ²

Let us consider further that the vertical movement is

harmonic and stationary; it can be characterized by the
iωt
following relationship: w( r , t ) = wc .e , where wc is
the amplitude of the soil displacement that solely depends
on the radial distance r. Assuming hysteretic energy
losses can be handled at a given frequency by an
equivalent viscous-elastic shear modulus Fig. 8. Equivalent Winkler spring-dashpot soil model
G = G.(1 + 2iξ ) , the equation of movement becomes: supporting the pile shaft
*

02005-p.6
 CSNDD 2014

3.3 Analytical solution for equivalent medium

I z = −2π .a .G s (1 + 2iξ ).
* ( )
H 1( 2 ) a *
Assuming that the shear modulus is independent of the
radial distance r to the pile shaft, a unique shear wave
H (2)
0 (a )
*
(19)

velocity can be defined as: Vs = Gs / ρ for a purely Dimensionless impedance parameters defined by (12) and
(13) are depicted in figure 9, emphasizing the influence
elastic medium (or Vs* = Gs* / ρ for the equivalent of the hysteretic damping coefficient. In the absence of
visco-elastic medium) (10) can be thus expressed as: viscous damping ( ξ = 0 ) only radiation or geometrical
damping prevails. In that case, the stiffness parameter
∂ ² wc 1 ∂wc (real part K z , a of impedance) tends toward π for
+ . + ( k ) .wc = 0
* 2
(14)
∂r ² r ∂ r increasing frequencies ( a → ∞ ) , per (9). Except for

where k = ω / Vs is the shear wave number, the asterisk

low a values, the dimensionless damping term C z , a
indicating that the viscous soil behavior (characterized by linearly increases with frequencies at a rate of
the Kelvin-Voigt formulation) is taken into consideration, ΔC z ,a = 2πΔa . An increased hysteretic damping
i.e. k = ω / Vs .
* *
coefficient enhances the quasi-linear increase of the total
Defining the dimensionless frequency for the pile (hysteretic and radiation) damping term, but decreases the
in-phase stiffness component. Figure 9 also shows that
shaft-soil interface as a = ω.r0 / Vs or a * = ω.r0 /Vs* , K z , a can become equal to 0 at particular values of a
the general solution of (14) is given by
( a ( K z , a = 0 ) ), implying a potential for some form of
§ r·
wc = w0c . H ¨¨ a*. ¸¸ H 0( 2) a *
( 2)
0 ( ) (15) ‘resonance’ effect. Influence of the hysteretic damping on
© r0 ¹ the infinite annular shear plate apparent ‘resonant’
frequency can be appreciated in figure 10. The physical
with H υ = J υ − i.Yυ the Hankel function and J ν ,
( 2)
reason for that effect remains unclear, despite the fact that
Yν , Bessel functions of order υ of the first and second similar phenomena have been experimentally observed
type, respectively. Based on the on model tests [13].

τ (r , t ) = G . ∂w(r , t ) ∂r
*
relationship, the solution can
also be expressed in terms of stress amplitudes:
a* .w0c § r·
τ c = Gs .(1 + 2iξ ). ( )
. H1( 2) ¨¨ a* . ¸¸ H 0( 2 ) a*
r
(16)
r0 © 0 ¹

Particular solutions can be obtained by applying

adequate boundary conditions. A first boundary condition
corresponds to the imposed displacement at the pile shaft-
soil interface, i.e. wc = w0 c for r = r0 . The second
boundary condition is a radiation condition imposing that
the wave should only propagate away from the vibration
source at the outer boundary ( w( r → ∞, t ) , also known
as “Sommerfeld” condition [12]). Fig. 9. Impedance parameters for homogeneous shear modulus

→ wc (r ) = ℜ 2 {wc (r )} + ℑ2 {wc (r )}
§ r· § r· (17)
J 02 .¨¨ a . ¸¸ + Y02 .¨¨ a . ¸¸
= w0c © r0 ¹ © r0 ¹
J 0 .(a ) + Y0 .(a )
2 2

→ τ c ( r ) = ℜ 2 {τ c (r )} + ℑ2 {τ c (r )}
§ r· § r· (18)
J 12 .¨¨ a . ¸¸ + Y12 .¨¨ a . ¸¸
= Gs .
a.w0 c
. © r0 ¹ © r0 ¹
r0 J 0 .(a ) + Y0 .(a )
2 2

Fig. 10. Dimensionless “resonant” frequency of the soil annular

shear plate vs. hysteretic damping coefficient
Using [Equ.8] we deduce the equivalent impedance:

02005-p.7
 MATEC Web of Conferences

Typical results expressed in terms of dimensionless 4 Non linear aspects

displacement and stress amplitudes profiles are presented
in figure 11 for two dimensionless frequency values and
two damping coefficients. The dimensionless 4.1. Literature review
displacement (stress) amplitude is the ratio between the Nonlinear models of axial pile-soil vibration started with
displacement (stress) amplitude wc (τ c ) at a given the works of Novak and Sheta [15], Mitwally and Novak
radial distance to that w0c (τ 0 c ) at the pile shaft-soil [16], Han and Sabin [17], and El Naggar and Novak
([18],[19]) who suggested distinguishing two separate
interface. radial soil zones around the pile shaft: an inner zone with
For ξ =0 and high a values, displacement reduced shear stiffness and an outer zone where the
amplitudes attenuate according to the inverse of the elastic solution is considered.
To eliminate undulations in the impedance
square root of the radial distance, while for ξ = 0 and
functions due to wave reflections from the interface
low a values, shear stress amplitudes attenuate between the two media, some researchers proposed a
according to the inverse of the radial distance. These continuously increasing modulus with radial distance to
results can be demonstrated by considering the the pile shaft. Gazetas and Dobry [20] and Veletos and
asymptotic behavior of the Bessel functions for a → ∞ Dotson [21] suggested schemes in which the modulus
(17) and (18). It can also be noted that material damping increased unboundedly. Han and Sabin [17] formulated
enhances the radial attenuation of both displacement and impedances based on a parabolic variation of the medium
stress amplitudes. properties so that the inner zone has properties smoothly
approaching those of the outer zone.
These contributions however address the problem
of lateral soil heterogeneity with only qualitative
reference to the non-linear soil response, since the
variations of soil properties invoked are merely
hypothetical. To aid practical applications, Michaelides et
al. ([22],[23]) utilized experimental data (e.g. Vucetic &
Dobry [8]) characterizing the dependence of the secant
shear modulus and hysteretic damping of soil on the
shear strain amplitude and the nature of the soil (the latter
represented by the plasticity index PI). The variation of
modulus and damping is then related to the magnitude of
the applied load through the amplitude of the shear
strains induced within a succession of co-axial cylinders.
Such an approach involves assumptions related to the
shear stress distribution and implies the use of an iterative
procedure to calculate the variation of modulus as a
function of the distance to the pile shaft.
Some modifications to the Michaelides et al.’s
model have been proposed by Holeyman et al. [24] to
simplify definitions of the model parameters as well as
calculation procedures. Using the modified method, a
refined soil discretization can be achieved based on more
rigorous soil behavior description and without a priori
assumptions about the shear modulus or shear stress
radial distributions.
The following sections describe analytical solutions
Fig. 11. Radial distribution of displacement and shear stress assuming various theoretical radial variations of shear
amplitudes modulus. The results are evaluated over a range of
parameters and compared with those obtained from the
Impedance parameters for use in the Winkler semi-analytical model derived from [22] and [23], and
approach have been studied by many researchers in the from a radial discrete model simulating the pile shaft and
past. First models (e.g. Novak [14]) were based on the soil movements by integrating of the laws of motion
assumptions that the soil behavior is governed by the ([25], [26], [27] (pile driving model), [9], [28], and [29]
laws of (viscous-) elasticity and the soil is perfectly (Vibratory pile driving)).
bonded to the pile. In practice however the soil region
immediately adjacent to the pile can undergo a large
degree of straining which causes the soil-structure system
to behave non-linearly and even degrade under cyclic 4.2 Radius-dependent shear modulus
loading. Slippage can also occur about the contact area.
Assuming that the stress developed into the soil
attenuates according to the inverse of the radial distance,

02005-p.8
 CSNDD 2014

which is exact in the static case, the following can be m

§r·
established: τ / τ 0 = r0 / r . If we further consider that G = G .¨¨ ¸¸ *
s0 (22)
the stress τ0 is a fraction f of the shear strength τ max , i.e. © r0 ¹
with m values decreasing with the distance to the pile
τ / τ max = f .r0 / r , (5) becomes:
shaft. This discretization also allows taking into account
damping coefficients ξ varying with r, albeit through a
§r ·
G s = G max .(1 − f .¨ 0 ¸) (20) piecewise approximation.
©r¹ Michaelides et al. assume that the shear stress
where f can be viewed as a “loading factor” or soil distribution is independent from the shear modulus
distribution, in order to alleviate interdependence
strength mobilization ratio at r = r0 : it is actually the
between Gs , τ c and γ c . Based on this assumption and
inverse of the factor of safety of the pile shaft capacity.
Because of the variation of the shear modulus with the use of empirical rules for the shear modulus
the radial distance, two extreme shear wave velocities can distribution, they proposed following equation:
be distinguished: Vs 0 = Gs 0 / ρ at the pile shaft-soil § ­ r ½
0.72
·
Gs = Gmax.¨1 − ®Λ. 0 .h(ar )¾ ¸ (23)
interface and V ff = Gmax / ρ in the free field at the ¨ ¯ r ¿ ¸
© ¹
furthest distance away from the pile. Dimensionless
where Λ = τ τ max is a loading intensity factor and
frequencies can thus be defined respectively for the pile
shaft-soil interface: a0 = ω.r0 / Vs 0 (or h(ar ) a shape function. Since the method imposes the
use of iterations to determine h(ar ) values, Michaelides
a 0 = a 0 / 1 + 2iξ = ω.r0 / Vs 0 ) and the free-field:
* *
limited the number of radial increments to four. That
a ff = ω.r0 / V ff (or a = a ff / 1 + 2iξ ). Defining
* approximation has been further enhanced by Bertin [30],
ff
who developped a special routine able to iterate on many
ς = r / r0 > 1 the dimensionless distance to the pile more radial increments.
shaft, the general equation of movement becomes: Because of the limitations of the Michaelides et
al.’s method, Bertin [30] suggested another approach
∂ ² wc 1 ζ ∂w ζ based on analytical elements to discretize the radial
+ . . c + (a *ff ) 2 . .wc = 0 (21)
∂ζ ² ζ ζ − f ∂ζ ζ −f coordinate. The shear modulus distribution is still
described by a set of parabolas, but based on the Hardin
Analytical solutions of (21) are presented in [24] and & Drnevich [4] functions and Masing rules [5].
discussed below. Influence of the ‘f’ parameter on the Furthermore, a higher number of radial steps, up to 270,
shear modulus distribution is depicted in figure 12. The is considered. Because values of m and ξ are different
reader is referred to Bertin [30] for solutions based on
other assumptions of shear modulus distributions. for each element, displacements and stresses are obtained
by assuring continuity of displacement and stress
equilibrium. For each discretization step, the following
equations (24) are used:

­ ª A .J (
. χ .λ .ζ ( i −1) º
°ζ m0 / 2 .« (i −1) χ ( i −1) −1 (i −1) 0( i −1) (i −1)
1/ χ
) » if r ≤ ri
° (i −1) «+ B .Y
° ¬ ( i −1) χ ( i − 1 ) −1 . χ ( (
i −1) .λ 0 ( i − 1)
.ζ ( i −1) )
1 / χ ( i −1 )
»
¼
wc (r ) = ®
ª χ (
° m1 / 2 (i ) χ ( i ) −1 (i ) 0( i ) (i )
A . J . .λ .ζ
1 / χ( i )
) º
° ζ (i ) .« » if ri ≤ r ≤ ri +1
°
¯
« + B .Y . χ ( .
¬ (i ) χ ( i ) −1 (i ) (i ) (i ) ¼λ .ζ
1 / χ (i )
)»

The four integration constants A( i −1) , A(i ) , B( i −1)

and B (i ) are deduced from continuity and equilibrium
Fig.12. Influence of ‘f’ parameter on the shear modulus conditions, adopting previously described boundary
distribution, assuming stress attenuation according to r-1 conditions (imposed displacement and outer radiation).
As illustrated by figure 13, Bertin analytical
4.3 Semi-analytical solutions elements approximation [30] emulating Michaelides et
al.’s concept is correct for low values of the
Michaelides et al. [22],[23] suggested the use of a radial dimensionless frequency, but is more questionable for
discretization and approximation of the shear modulus Gs higher values.
distribution within each zone using the following
expression

02005-p.9
 MATEC Web of Conferences

Fig. 14. Numerical model geometry (adapted from [28])

The program was further developed using Matlab®

routines [30] to produce the results presented herein with
a view to compare them with the above described semi-
analytical methods. Strain rate effects as well as cyclic
degradation effects which are accounted for in the
Hipervib-II program were however disabled to produce
results compatible with the assumptions adopted in the
other methods discussed in this paper.
When comparing Hipervib-II with the (semi-)
analytical methods, following considerations are made:
(a) only results corresponding to the steady state (after a
few second of simulations) are presented herein, (b)
impedance parameters are calculated as follows:
2π /ω

Fig. 13. Shear stress radial distributions for different values of I z = 2 π . r0

³
. π0
τ ( r 0 , t ) dt
(25)
2 /ω
(22) m exponent compared to [30] enhanced approximation
³ 0
w ( r 0 , t ) dt
with w(t ) imposed at the pile shaft boundary and τ (t )
calculated by Hipervib-II at the pile shaft-soil interface.
In order to compare the numerical results with the
4.4 Numerical solution
analytical ones presented herein for the shaft only, the
base resistance that can be modeled in Hipervib-II has
4.4.1 Model been set to zero.
Holeyman [26],[27],[28],[29] has suggested the use of a
radial discrete model (see figure 14) to calculate the
vertical shear waves propagating away from the pile 4.4.2 Radial distribution of shear modulus
shaft. The pile is considered as a rigid body and the soil is
represented by discretizing the medium into concentric Relevant results are related to the evolution of the shear
rings that have their own individual masses and that modulus as a function of the radial distance, by reference
transmit forces to their neighboring ones. The movement to its value at the pile shaft-soil interface. Figure 15
of the pile and the rings is calculated from time shows results obtained with three methods: the analytical
integration of the law of motion: the equations of solution to (21), Michaelides et al. original approximation
movement are integrated for each cylinder based on their (4 increments, labelled as “M-4”), and Hipervib-II
dynamic shear equilibrium in the vertical direction. An program.
energy absorbing boundary condition in accordance with
Table 1. Reference parameters for comparative study
plane-strain elasticity theory [14] limits the lateral extent
of the model. Soil density ( ρ ) 1.8 T/m³
The model makes use of constitutive relationships Maximum shear modulus (Gmax) 60 MPa
representing the large-strain, dynamic and cyclic shear Free field shear wave celerity (Vff) 183 m/s
stress-strain strength behavior of the medium surrounding
the pile shaft. Initially implemented for vibratory driving Shear strength ( τ max ) 0.15 MPa
modeling in a Basic computer code “Hipervib-II” [9], it Reference shear strain ( γ r ) 1.25 10-3
applied the hyperbolic Kondner law (4) and Masing rules
[5] to model the shear force-displacement relationships Hysteretic damping value ( ξ ) 0 -
between successive rings. Pile shaft radius (r0) 0.5 m

02005-p.10
 CSNDD 2014

Fig. 15. Shear modulus distributions according to different

approaches for same acceleration amplitude of 9.9 m/s2
Fig. 16. Stiffness and damping parameters deduced from the
Reference parameters used to that end are summarized in
semi-analytical method (see Table 1 for reference parameters)
Table 1, for two different combinations of imposed
displacement amplitude and frequency leading to a given
acceleration amplitude of 9.9 m/s2.
For the same imposed acceleration amplitude
( w0 c .ω ² ), all three methods indicate a more heavily
degraded shear modulus at the pile shaft-soil interface for
a larger displacement amplitude. When however
compared to the Michaelides original method, the
analytical results [30] better match those of the numerical
model.

4.5 Comparison of calculated impedance

parameters Fig. 17. Comparison of analytical and semi-analytical solutions
for woc/ro=0.4 10-3 (see Table 1 for reference parameters)
Impedance parameters deduced from application of the
semi-analytical method using Table 1 parameters are observation could be extended to other values of woc/ro
presented in figure 16 with a view to emphasize the and ξ , it would mean that the simplest form of stress
influence of the imposed displacement amplitude. Figure
attenuation (corresponding to the static case) could be
16 represents the real and imaginary parts of the
adopted for the dynamic cases.
impedance versus the dimensionless frequency aff for
Combined influences of frequency and imposed
three displacement amplitudes.
displacement amplitudes are represented in figure 18 for
Comments made for the homogeneous model in
the analytical method. Both real and imaginary parts of
Section 3.3 are also applicable to the non-homogeneous
the impedance decrease as the imposed displacement
models, as far as the evolution of the impedance curves
increases, whatever the frequency. Inspection of (12) and
as functions of the dimensionless frequency is concerned,
figure 19 shows how the stress value at the pile shaft-soil
when trading the hysteretic damping coefficient for the
interface naturally tends towards the shear strength with
amplitude of displacement. A comparison between semi-
increasing displacements, evidencing another
analytical and analytical results is presented in figure 17,
manifestation of stiffness degradation.
indicating some quantitative variations in the impedance
The evolutions of the pile shaft friction mobilization
results that mainly depend on the free field dimensionless
frequency. In the above example, analytical results are ratio (τ oc / τ max ) versus dimensionless displacement
quite close to the Bertin [30] approximation. If that (woc/ro ) shown in figure 19 can be compared to so-called

02005-p.11
 MATEC Web of Conferences

“t-z curves” published in the literature to evaluate the pile convergence problems. By contrast, the semi-analytical
shaft load-displacement behavior under axial static method suffers from numerical limitations arising from
( a ff = 0 ) conditions. As an example, static friction the radial soil discretization. The semi-analytical method
correctly models a radial variation of the soil hysteretic
mobilization curves adopted by Holeyman [25] based on damping, contrary to the analytical method which
an extension of the influence radius Rm approach assumes a homogeneous hysteretic damping.
suggested by Randolph and Wroth [31] to incorporate a Similar approaches can be used to address the
hyperbolic stress-strain law are also plotted in figure 19 dynamic non linear response of piles under a lateral mode
for two values of the boundary radius Rm. Confirmation of deformation, as well as coupling effects between the
of the effective reduction of the apparent “quake” value axial and lateral modes of deformation [32].
as velocity increases, observed by Holeyman [7] can be
found in figure 19 for increasing values of a ff .

Fig. 18. Impedance parameters as functions of the imposed

displacement (analytical model for Table 1 parameters). Fig. 19. Shear stress developed at the pile shaft-soil interface as
a function of the imposed cyclic displacement amplitude
(analytical model with Table 1 parameters vs. typical static t-z
5 Conclusions curve)

The dynamic axial response of pile shafts can be

approached using the concept of a continuously
distributed mechanical impedance replacing the 6 Notations
embedding medium. Equivalent impedance parameters
can be defined to characterize the equivalent in-phase The following symbols are used in this paper:
‘spring’ and the equivalent out-of-phase ‘dashpot’. When
considering a pile shaft undergoing axial oscillations, a = dimensionless frequency ω.r0 / Vs [-]
shear strains (and shear stresses) are induced in the
surrounding soil, the amplitude of which attenuates a * = complex value of dimensionless frequency
radially away from the pile shaft. Analytical equivalent
linear methods and a numerical method are shown to a / 1 + 2iξ = ω.r0 / Vs* [-]
adequately derive those amplitudes of shear stresses and
shear strains as functions of the radial distance r to the a0 = dimensionless frequency for the pile shaft-soil
pile shaft. The corresponding dynamic impedance interface ω.r0 / Vs 0 [-]
components are then readily determined.
The analytical solution proposed herein a ff = dimensionless frequency in the free field
accommodates a continuous variation of soil properties
alleviating wave reflections and avoiding numerical ω.r0 / V ff [-]

02005-p.12
 CSNDD 2014

C z = damping coefficient [kPa.s-1] woc = amplitude of imposed displacement at the pile

Cza = ℑ{I z }/(ω.Gs 0 ) = dimensionless damping shaft [m]
parameter [-] Yυ = Bessel function of order υ of second type
f = “loading factor” or soil strength mobilization ratio at γ = shear strain [-]
r = r0 [-]
γ = shear strain rate [s-1]
G = shear modulus [kPa]
γc = amplitude of shear strain [-]
G = complex shear modulus G.(1 + 2iξ ) [kPa]
*
γr = reference strain τ max / Gmax [-]
Gmax = initial (maximal) shear modulus [kPa]
ξ = damping coefficient [-]
Gs = secant shear modulus [kPa] ρ = soil density [T/m³]
Gs 0 = shear modulus at the pile shaft-soil interface τ = shear stress [kPa]
( r = r0 ) [kPa] τc = amplitude of shear stress [kPa]

H n = Hankel function J ν − iYν

2
τ0 = shear stress at the pile shaft-soil interface [kPa]

I 0 = modified Bessel function of order 0 of first type τ 0c = amplitude of shear stress at the pile shaft-soil
interface [kPa]
I z = unit (lineal) shear impedance of the soil against the
pile shaft movement in the z direction [kPa] τ max = maximal shear stress (shear strength) [kPa]

J υ = Bessel function of order υ of first type ω = circular frequency [rad.s-1]

k = shear wave number ω / Vs [m-1] ς = dimensionless distance to the pile shaft r / r0 > 1 [-]

k = complex shear modulus ω / Vs [m-1]

* *
The asterisk indicates that the viscous soil behavior
(characterized by the Kelvin-Voigt formulation) is taken
k *ff = complex shear modulus in the free field a *ff / r0 into consideration.
[m-1]
K 0 = modified Bessel function of order 0 of second type
K z = the stiffness coefficient [kPa] References
K za = dimensionless stiffness ℜ{I z }/ Gs 0 [-] 1. Upwind – Integrated Wind Turbine Design. Project
funded by the European Commission under the 6th
r = distance to the pile shaft [m] (EC) RTD Project No. 019945 (SE6), (2010)
r0 = pile shaft radius [m] 2. Smith E. A. L., Pile driving analysis by the wave
equation. Journal of Soil Mechanics and Foundation.
Rm = influence radius of shear strain field 86(4): 35-61, (1960)
t = time [s] 3. Hardin, B0, Black, WL. Vibration modulus of
normally consolidated clay. Journal of the Soil
Vs = shear wave velocity Gs / ρ [m.s-1] Mechanics and Foundations Division ASCE 1968;
92(2): 353-369
4. Hardin, B. O. and Drnevich, V. P., Shear modulus
Gs / ρ
*
Vs* = complex value of shear wave velocity and damping in soils: design equations and curves
[m.s-1] Journal of the Soil Mechanics and Foundations
Vs 0 = shear wave velocity at the pile shaft-soil interface Division, American Society of Civil Engineers 98 7,
pp. 667–692, (1972)
Gs 0 / ρ [m.s-1] 5. Masing, G., Eigenspannungen und Verfeistigung
beim Messing Proc. 2nd Int. Congress of Appl.Mech.,
V ff = shear wave velocity in the free field Gmax / ρ pp. 332-335 (1926)
-1 6. Kondner, R. L., Hyperbolic Stress-Strain Response:
[m.s ] Cohesive Soils. Journal of the Soil Mechanics and
w = displacement [m] Foundations Division, ASCE, 89, No. SM1, pp.115-
143 (1963)
wc = time-independent value of displacement [m] 7. Ishibashi, I., & Zhang, X., Unified dynamic shear
wo = displacement at the pile shaft [m] moduli and damping ratios of sand and clay. Soils
and Foundations, 33(1) , 182-191 (1993)

02005-p.13
 MATEC Web of Conferences

8. Vucetic, M. & Dobry, R., Effect of soil plasticity on http://dx.doi.org/10.1016/j.soildyn.2012.09.006, 44

cyclic response J. Geotech. Engng, ASCE 117, No. 115-126 (2013)
1, 89-107 (1991) 25. Holeyman, A., Contribution à l'étude du
9. Holeyman, A., HIPERVIB-II, A detailed numerical comportement transitoire non-linéaire des pieux
model proposed for Future Computer pendant leur battage. Doctoral thesis, Université
Implementation to evaluate the penetration speed of Libre de Bruxelles, April, 1984, 584 p. (1984)
vibratory driven sheet Piles Research report for 26. Holeyman, A., Dynamic non-linear skin friction of
BBRI, EarthSpectives, Irvine, Ca, USA, 54p (1993) piles, Proceedings of the International Symposium
10. Bauer, E., Calibration of a Comprehensive on Penetrability and Drivability of Piles, San
Hypoplastic Model for Granular Materials. Soils and Francisco, 10 August, Vol. 1, pp. 173-176 (1985)
Found. 36(1), 13-26 (1996). 27. Holeyman, A., Technology of Pile Dynamic Testing,
11. Von Wolffersdorff, P. A., A hypoplastic relation for in: Application of Stress-Wave Theory to Piles,
granular materials with a predefined limit state edited by F. Barends, Balkema, Rotterdam, pp. 195-
surface. Mechanics of Cohesive-Frictional Materials 215 (1992)
(1),251–271 (1996) 28. Holeyman, A., and Legrand, C., Soil Modeling for
12. Sommerfeld, A., Mechanics, Lectures on Theoretical Pile Vibratory Driving, U.S. FHWA International
Physics, volume I. Academic press Inc., 4 edition Conference on Design and Construction of Deep
(1952) Foundations, Orlando, Florida, December 1994, Vol.
13. Storz, M., Chaotic Motion in Pile-Driving. First II, pp. 1165-1178 (1994)
International Conference Soil Dynamics and 29. Holeyman, A., Vibratory Pile Driving, in: Quality
Earthquake Engineering (pp. 503-512). Assurance on Land and Offshore Piling, Edited by S.
Southhampton, England: Computational Mechanics Nyyama and J. Beim, Balkema Publishers,
(1991) Rotterdam, 2000, pp.479-494 (2000)
14. Novak M., Dynamic stiffness and damping of piles 30. Bertin R., Modélisation de l’interaction axiale sol-
Can. Geotech. J. 11, No. 4, 574-598 (1974) pieu - Détermination des paramètres d’impédance
15. Novak, M. and Sheta, M., Approximate approach to d’un milieu non-homogène Master’s thesis,
contact effects of piles In Special technical Université Catholique de Louvain (2009)
publication on dynamic response of pile foundations: 31. Randolph, M.F. and Wroth, C.P., Analysis of
analytical aspects (eds M. W. O'Neill & R. Dobry). deformation of vertically loaded piles, J. Geotech.
New York: ASCE (1980) Eng. Div. ASCE 104(GT12): 1465-1488 (1978)
16. Mitwally H and Novak M , Pile driving analysis 32. Malek, A. and Holeyman, Flexural Analysis in
using shaft and FEM. Proceedings of the third Dynamic Pinned Head Pile Testing, Geotech Geol
international conference on the application of stress Eng 32:59–70 (2014)
wave theory to piles; Bitech Publishers, Vancouver 33. Malek, A. and Holeyman, A., Numerical evaluation
(1988) of nonlinear lateral pile vibrations on nonlinear axial
17. Han Y.C. and Sabin G.C., Impedances for radially response of pile shaft Soils and foundations,
inhomogeneousviscoelastic soil media J. of Engrg. http://dx.doi.org/10.1016/j.sandf.2013.04.002, (2013)
Mechanics, Vol.121, No.9, pp.939-947 (1995)
18. El Naggar, M. & Novak, M., Nonlinear axial
interaction in pile dynamics. Journal of Geotechnical
Engineering, 120(4): 678-696 (1994)
19. El Naggar, M. & Novak, M., Nonlinear model for
dynamic axial pile response. Journal of Geotechnical
Engineering, 120(2): 308-329 (1994)
20. Gazetas, G. and Dobry, R., Simple radiation damping
model for piles and footings. J. Engrg. Mech.,
ASCE, 10(6), 937-956 (1984)
21. Veletsos, A. S. & Dotson, K. W., Vertical and
torsional vibration of foundations in inhomogeneous
media J. Geotechnical Engineering Division, ASCE
114, No. 9, 1002-1021(1988)
22. Michaelides O., Gazetas G., Bouckovalas G. and
Chrysikou E, Approximate nonlinear dynamic axial
response of piles Geotechnique 48, No.1, 33-53
(1997)
23. Michaelides O., Bouckovalas G. and Gazetas G.,
Non-linear soil properties and impedances for axially
vibrating pile elements Soils and foundations,
Vol.38, No.3, 129-142 (1988)
24. Holeyman, A., Bertin, R., and Whenham, V.,
Impdedance of pile shafts under axial vibratory
loads, Soil Dynamics and Earthquake Engineering,

02005-p.14