- Email: [email protected]
Application of a continuum numerical model for pile driving analysis and comparison with a real case
Available online at www.sciencedirect.com
Computers and Geotechnics 35 (2008) 406–418 www.elsevier.com/locate/compgeo
Application of a continuum numerical model for pile driving analysis and comparison with a real case Shahram Feizee Masouleh, Kazem Fakharian
*
Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran 15875, Iran Received 4 February 2007; received in revised form 24 July 2007; accepted 28 August 2007 Available online 29 October 2007
Abstract This paper presents the results of a so-called continuum numerical model for wave propagation analysis and soil-pile dynamic response during pile driving. An axisymmetric finite difference numerical model is developed having solid elements for both pile structure and the soil media surrounding and below the pile. Interface elements are used between the pile shaft and the soil to facilitate the sliding between the two media. The performance of the developed model is verified in two stages. First, a simple rod is subjected to a half sinewave force function at the rod head and the corresponding reflections of force and velocity (multiplied by impedance) are presented for different boundary conditions at the rod tip. The model is then used for signal matching analysis of a real driven pile for which complete information of soil layering, dynamic test signals, and static load test results are available. The signal matching analysis was performed successfully and comparison between several other predicted and measured parameters proved the reasonably good performance of the developed continuum model. 2007 Elsevier Ltd. All rights reserved. Keywords: Pile driving; Stress wave analysis; Continuum numerical model; Finite difference method; Signal matching analysis
1. Introduction Despite the fact that pile design approaches have advanced enormously within the past few decades, but yet, the most fundamental aspect of pile design, that of estimating the axial bearing capacity, relies heavily on empirical correlations as discussed by Randolph [1]. The piles are classified in two main groups in terms of their installation method, including non-displacement piles and displacement piles. In non-displacement piles, in case of utilizing casing or drilling mud as well as proper construction conditions, the initial state of stresses are relatively maintained adjacent to the pile shaft and below tip [2]. In displacement piles such as driven piles, on the other hand, the pile is pushed down by hammer impact or jack. The
*
Corresponding author. Tel.: +98 21 22 90 5633; fax: +98 21 22 27 2000. E-mail address: [email protected] (K. Fakharian). 0266-352X/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2007.08.009
changes in stress state such as degree of compaction around shaft or below tip are difficult to measure, thus making it more difficult to estimate the axial pile capacity. Many attempts have been made for developing an adequate dynamic method to accurately estimate the static axial capacity. The main objective of the dynamic methods has been to establish a relationship between the blow counts and the axial capacity of the pile. Many dynamic formulas have been developed and are still used by some contractors during construction, but the most important innovation was the introduction of one-dimensional wave equation analysis in piles (WEAP) by Smith [3], in which pile is simulated by a number of masses attached to each other by elastic springs, and soil is replaced by a number of springs, sliders and linear viscous dampers to simulate the visco, elasto-plastic response of the soil. The WEAP model, however, requires to estimate some key soil parameters such as quake and damping coefficients as well as hammer efficiency. Such factors significantly influence the analysis results [4–6].
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
407
Nomenclature Notations A pile cross section area (m2) C wave speed in pile (m/s) D pile diameter (m) E elastic modulus (N/m2) F force (N)
Rausche [7,8] proposed the signal matching technique and developed an analysis method, CAPWAP, to overcome the shortcomings of WEAP model. The signal matching technique relies on the pile dynamic test results (PDA). The model used in CAPWAP is very similar to WEAP, except that unloading effects are also considered, though such effects on capacity are insignificant [9]. The signal matching analysis is the most reliable dynamic method in practice. This method, however, uses mass–spring–dashpot system to model the surrounding and below tip soil media, therefore imposes some limitations. Pile is a continuum body having infinite number of degrees of freedom surrounded by a continuum media that is soil. In the mass–spring system, the boundary conditions are different than the real system, causing deviations in the system response. Also the effect of soil inertia (also called ‘‘radiation damping’’ or ‘‘geometrical damping’’) can not be easily considered. To overcome the above shortcomings of mass–spring– dashpot models, different attempts have been made to make use of other numerical solutions. Randolph tried to implement the far field effect of soil using additional dashpot (so-called inertial dashpot representing the surrounding soil inertia or radiation damping), as shown in Fig. 1 [1]. Some other studies have tried to use finite element method so that the soil media around the pile shaft and below the pile tip are directly considered in the model. Coutinho et al. [10] conducted a research on the drivability of piles into the Brazilian coast deposits. They considered the effect of some parameters influencing the pile driving including pile-soil interface and damping. Borja [11] adopted a FEM scheme to model the pile driving system, capable of continuously interpolating displacement, velocity and acceleration profiles along the pile length. As opposed to the previous models that were using the impact velocity of hammer as initial condition, Borja utilized the impact force curve of the hammer with respect to time. This allowed a more realistic simulation of the pile–hammer interaction. Some results were presented on the basis of analyses on an assumed example problem, and no case study was used for comparison with the numerical predictions. Mabsout and Tassoulas [12] and Mabsout et al. [13] also used an axisymmetric FEM model (Fig. 2) for pile driving analysis. The surrounding soil was assumed to be saturated normally-consolidated clay and its stress–strain behavior
Kn Ks L v Z q
interface normal stiffness (N/m3) interface shear stiffness (N/m3) pile length (L) velocity (m/s) impedance (EA/C) mass density (kg/m3)
Fig. 1. Dynamic soil-pile shaft model to account for inertial damping. (After Randlolph, 2003).
was simulated by a non-linear bounding-surface-plasticity model. The slide-line algorithm was used for the interface modeling to enable the soil-pile relative displacements. The at-rest lateral stress coefficient was considered equal to unity. Viscous boundaries were used to absorb the radial stress waves, resulted from hammer impact on the pile. The hammer impact was modeled as a force–time function on the pile head. The transient dynamic analysis was performed using the constant-average-acceleration algorithm. A bored pile with embedment depths of 6, 12 and 18 m was assumed for the driving analysis. The tip and shaft resistances were distinguished and the results showed higher resistances with increase in pile embedment depth. It was also observed that resistance to driving increases with increase in number of blows. The excessive pore pressure at pile tip was shown to dramatically decrease with distance from the pile. No comparison was, however, made between the predictions and a real-life pile driving data. It was suggested, however, as their conclusions to verify the predictions of the FEM model against actual pile driving data. The main objective of this paper is to develop a so-called continuum numerical model to overcome some of the limitations of the mass–spring–dashpot system, and present an improvement in comparison with the previous FEM studies, in that a signal matching procedure to be adopted for comparing the predictions of the model with a well-docu-
408
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
Lagrangian calculation which is adequate for large deformation modeling. One of its advantages is possibility of modeling large systems with a limited memory, as no stiffness matrix is required (as opposed to FEM, for example). In addition, large deformations do not significantly increase the run-time as no stiffness matrix updates are required after each load or time increment. The wave propagation analysis in an elastic rod with different end (tip) boundary conditions is conducted first to verify the model performance and to evaluate the sensitivity of some parameters. Three series of analyses are performed on: (1) a rod with fixed-end (tip) condition (2) a rod with free-end (tip) condition (3) a rod with elastic support at the end (tip). No soil or other supports are considered around the rod shaft yet. Fig. 3 shows the axisymmetric mesh and the boundary conditions of the fixed-end, freeend, and elastic support, in 3a, 3b and 3c, respectively, for a 0.5 m diameter, 20 m long rod. More details of the model and results are presented in the following sections. 2.1. Fixed-end rod
Fig. 2. Continuum soil-pile model for simulation of pile driving. (After Mabsout, 1995).
mented actual driven pile. Among the major advantages of the model is the consideration of soil inertia and/or radiation damping effects. An axisymmetric finite difference model is developed using FLAC-2D with its axisymmetric capability to simplify the 3D model to 2D. For preliminary verification of the model performance, an elastic rod with different boundary conditions at the tip (fixed, free, variable elastic modulus) is considered for wave propagation analysis and response from tip. To further investigate the model capability, a real driven pile that was carefully monitored in Netherlands is considered for modeling, for which dynamic and static test results, the soil classification and CPT profiles are available. A signal matching procedure (similar to CAPWAP) was followed by changing the strength and deformation parameters of soil and interface between pile-soil. The fine-tuned parameters were then used to simulate the static axial load–displacement response to be compared with the real test results. 2. Numerical model and preliminary verification A finite difference scheme is adopted to simulate an axisymmetric model for pile driving analysis, using FLAC-2D version 4. This software is developed on the basis of
In the first series of analyses, a fixed-end round rod is assumed having length, diameter, elastic modulus, Poisson’s ratio, and unit weight of 20 m, 0.5 m, 25,000 MPa, 0.15, and 25 kN/m3, respectively. An axisymmetric mesh of 5 · 5 cm is used for modeling the rod half-section, as shown in Fig. 3. The rod gravity is neglected in the analysis. The hammer impact on the rod head is simulated by a half-sine stress wave with an amplitude of 5 MPa and frequency of 320 Hz. Regarding pffiffiffiffiffiffiffiffiffi that the wave speed in the rod is 3162 m/s ðc ¼ E=qÞ, the wave length of the halfsine signal generated in the rod is 4.94 m. Therefore, the force and velocity records are picked up at 5 m below the rod head, to more clearly represent the end (tip) dynamic responses, and prevent the mixing up of the upward stress wave and downward reflection from the rod head. The analysis results for a rod without damping effects are presented in Fig. 4a. The force (F) and velocity (v) variations with time are plotted in which velocity is multiplied by the rod impedance (Z = EA/C) to be equivalent to force (Zv). This is the same procedure used in PDA test output presentations. The horizontal axis is scaled in L/C, which is the time required for stress wave to travel between the pick up records and the rod end (tip). The force sign conventions are positive and negative for compressive and tensile waves, respectively. Downward displacement of a particle on the rod section shows a positive velocity, and upward displacement shows a negative velocity. Fig. 4a shows that after a time equivalent to 2L/C, that is the time required for stress wave to travel between the record point at 5 m below the rod head to rod tip and return to the same point (L = 15 m), the F wave is suddenly shifted up (positive) and Zv wave is shifted down (negative) and their amplitudes are equal to the generated wave amplitude. This observation is exactly in accordance with the 1D longitudinal wave propagation theory in rods
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
CL
CL
409
CL
Mesh size: 5x5cm
Elements with Variable Elastic Modulus
Viscous boundary
(a)
(b)
(c)
Fig. 3. Axisymmetric mesh for a rod with different boundary conditions at tip: (a) fixed-end; (b) free-end; (c) elastic support at the end.
a
b 1200
1200 F
F 800
Zv
400 0 0
2
4
6
-400
8
10
Force (kN)
Force (kN)
800
Zv
400 0 0
2
4
6
-400 -800
-800
-1200
-1200 L/C
L/C
Fig. 4. Dynamic response of the fixed-end rod: (a) without viscous damping; (b) with viscous damping.
8
10
410
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
without damping. The immediate F and Zv waves shifted down (negative) after tip reflections, are the reflections from the rod free head (top) boundary condition. In fact the downward initial compressive wave is reflected compressive at rod fixed-end (tip) and reflected tensile at the rod free-top (head) boundary conditions. Fig. 4b represents similar results for the same rod and loading, but with a 1% assumption of viscous damping ratio for rod. The ‘‘local damping’’ option of the model is activated that is on the basis of increase or decrease in mass of mesh nodes at specific times during oscillations, such that when velocity sign is changing, mass is increased, and at minimum and maximum velocities, the mass is reduced [14]. The F wave and Zv wave amplitudes are attenuated with time due to the damping effects of rod, as expected. 2.2. Free-end rod In the second series of analyses, a free-end (tip) rod is assumed, having the same geometry and material properties of the previous example. The results without and with damping are presented in Fig. 5a and b, respectively. Similar, but opposite in sign, variations are observed as expected, which is in fact representative of the 1D longitudinal wave propagation in a rod with free-end boundary conditions at rod tip and head. 2.3. Elastic support at rod tip In the third series of analyses, an elastic spring was assumed at the pile tip. The spring was simulated by using solid elements similar to the ones used for the rod itself with the same total cross section area of the rod (0.196 m2), but with a length of 1 m (Fig. 3c). The stiffness and impedance of the 1-m rod (acting as a spring) was changed by varying its elastic modulus from 10 MPa to 10,000 GPa, to see the effect on the dynamic response due to the impact on the rod head. A viscous boundary condition was imposed at the bottom of the 1-m rod to absorb the portion of
b
1200
1200
800
800
400
400
0 0
2
4
6
8
10
-400
Force (kN)
Force (kN)
a
the wave not reflected at the tip. The same half-sine stress wave was applied at the rod head to simulate the dynamic impact. The results are presented in Fig. 6a–h, for E values of the 1-m segment varying (decreasing) between 10,000 GPa to 10 MPa, respectively. For the highest value of E (10,000 GPa) in Fig. 6a, the response is identical to the fixed-end rod of Fig. 4a. For the lowest value of E (10 MPa) in Fig. 6h, on the other hand, the response is the same as the free-end rod of Fig. 5a. A gradual transition is observed between compressive to tensile wave responses, when the elastic modulus varies from high to low magnitudes. It is interesting to note that when the elastic modulus is equal to 30 GPa (Fig. 6d), which is very close to the elastic modulus of the rod, both force and velocity reflections are very close to zero, i.e. no significant reflections are observed. This is due to the fact that at such an elastic modulus, the impedance (EA/C) of both the rod and the segment are equivalent to each other, and therefore, no wave is reflected. It is reminded from the basics of 1D wave propagation that the magnitude of the reflected wave is proportional to the impedance variations along its propagation path. For impedances higher than the rod (Fig. 6a–d), the reflections are compressive, whereas it is the other way around for lower impedances (Fig. 6e–h). The observations very well confirm the performance of the numerical model, as a base for more complicated analyses of a real pile embedded within various soil layers below the pile tip and around the pile shaft. The results of the three series of analyses presented above on a 20 m rod with open surroundings, but different boundary conditions at tip (free, fixed, constant impedance), indicates that the developed model is capable of simulating the pile driving analysis and the wave propagation and reflections from tip and head. Both force (F) and velocity (Zv) signals versus time are identifiable and promising to be used in a signal matching type analysis similar to CAPWAP, for example. To further evaluate the developed model capabilities, a real pile driven into a layered soil is considered for which complete site conditions and dynamic and static test results are available.
0 0
2
4
6
Zv
-1200
-800
Zv
- 1200 L/C
10 F
F -800
8
-400
L/C
Fig. 5. Dynamic response of the free-end rod: (a) without viscous damping, (b) with viscous damping.
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
411
E = 1000 GPa
E = 10000 GPa 1200
1200 F
800
F
800
Zv
400 0 0
1
2
3
4
- 400
Force (kN)
Force (kN)
Zv 400 0 0
- 400
- 800
- 800
-1200
-1200
1
2
3
L/C
L/C
E = 30 GPa
E = 100 GPa 1200
1200 F
800
F
800
Zv
400 0 0
- 400
1
2
3
4
Force (kN)
Force (kN)
Zv 400 0 0
- 400
- 800
- 800
-1200
-1200
1
2
3
E = 1 GPa
E = 10 GPa 1200
1200 F
800
F
800
Zv
400 0 0
1
2
3
4
Force (kN)
Force (kN)
Zv 400 0 - 400
0
1
3
4
-1200
-1200
L/C
L/C
E = 10 MPa
E = 100 MPa 1200
1200 F
800
0 1
2
3
4
- 800
Force (kN)
Zv
400
0
F
800
Zv
- 400
2
- 800
- 800
Force (kN)
4
L/C
L/C
- 400
4
400 0 - 400
0
1
2
3
4
- 800
-1200
-1200 L/C
L/C
Fig. 6. Effect of the support elastic modulus at tip on dynamic response of the rod.
3. Case study pile The case pile is selected from an elaborate pile testing research project along with the ‘‘4th International Conference on the Application of Stress-wave Theory to Piles’’ held in Netherlands, 1992 [15,16]. Several precast concrete piles were driven at the test site behind the auditorium of Delft Technological University. Static and dynamic test results are available for the selected test pile used in this study.
The soil profile of the site together with CPT results are shown in Fig. 7. The pile is precast concrete with 0.25 m dimension and elastic modulus of 36,000 MPa. The 19 m long pile was driven by a hydraulic hammer ICE SC-40, to an embedment depth of 18 m. The blow counts for the last 1 m of driving for 25 cm penetrations are 37, 29, 32 and 41, respectively. The force (F) and velocity (Zv) signals recorded by installing strain gages and accelerometers at 1 m below the pile head and related to the last blow at the end of pile driving are available. The results of a static
412
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
having an average correlation with both the cross section area and the circumferential area of the main pile, as shown in Fig. 8. The pile tip has penetrated about 3 m into a medium dense to dense sand (so called Pleistocene sand). Except a 4 m thick clay and peat layer (#2) along CL
19m
Layer 2 Clay and Peat
3m
Layer 1 Silty Sand
4m
Fill
2m
Pile
6m
Layer 4 Pleistocene Sand 4m
Fig. 7. Soil profile and CPT results [16].
Layer 3 Silty Sand
7m
Interface
load test performed after 3 weeks of the initial driving is also available. 3.1. The numerical model 4m
In order to benefit from the axisymmetric feature of the 2D numerical model, a 0.3 m diameter pile is assumed
Fig. 8. Axisymmetric pile-soil model for pile driving simulation.
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
CL
5x5cm Pile Elements (Linear Elastic) Uniform Stress=38 kPa (2m fill)
Embedment Depth
Pile-Soil Interface
Pile Toe Viscous Boundary
10x10cm Soil Elements
5x5cm Soil Elements
Viscous Boundary Fig. 9. Finite difference mesh for the axisymmetric pile-soil model.
413
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
the shaft, rest of the layers are mostly sand, but with relatively low qc values from CPT test. The boundary distances in horizontal and vertical directions, from shaft and pile tip, respectively, are 4 m in the axisymmetric mesh. The selected distances are far enough to minimize the stresses resulted and propagated from pile driving, as discussed in the following sections. In accordance with the highest stress gradients near pile shaft and the soil around the pile tip, the mesh sizes are specified 5 · 5 cm within 2 m distance from shaft and tip, and 10 · 10 cm outside this zone to the vertical and horizontal boundaries (Fig. 9). Interface elements between pile shaft and soil are inserted to facilitate the pile-soil slip during driving. The Mohr–Coulomb failure criterion is used for the interface, on the basis of which slip between soil and pile occurs, if the generated shear stress during driving exceeds the shear strength. As a result of hammer blow on the pile head, part of the stress wave enters and propagates within the surrounding soil. The degree and magnitude of the penetrated energy depends on the interface and the pile-soil mechanical parameters [17]. In real pile driving, the penetrated stress waves propagate away from the pile and dissipate gradually, which is referred to as ‘‘radiation damping’’. As the model boundaries are limited, viscous boundaries are assumed at the vertical and bottom horizontal boundaries to absorb such waves and prevent their reflection. Dampers in two perpendicular directions are implemented at each boundary node as proposed by Lysmer [18]. An elasto-perfectly plastic model with Mohr–Coulomb failure criterion is used for the soil material. The required mechanical parameters for each soil layer are estimated on the basis of CPT results with common correlations available in literature [19]. The estimated parameters are presented in Table 1. The cohesion and Poisson’s ratio are assumed 10 kPa and 0.3, respectively, for all the four layers.
(Fig. 9). In the second step, the interface elements with zero shear strength are considered between pile-soil to allow the pile to displace freely under self-weight. For performing the pile driving dynamic analysis, the recorded force–time at pile head by the strain gage (resulted from pile driving) is input on the pile, and then the velocity signal is calculated from the analysis. The calculated Zv at the pile head is compared with the measured Zv, and if not matched, efforts are made in subsequent analyses by trial and error, similar to the CAPWAP type analyses, to match the signals. The measured force, measured Zv, and the fine-tuned computed Zv are shown in Fig. 10. The most important influencing time domain for the shaft resistance adjustments is between peak force (or velocity) to 2L/C (about 9.5 ms in the study case). Therefore, trial and error were made to match the computed and measured velocity signals, by adjusting the interface strength and deformation parameters for different layers. The final parameters are presented in Table 2. In order to adjust the pile tip effects in the signal matching process, the soil strength (angle of friction and cohesion) and deformation (elastic modulus and dilation
1.2 Force - measured
1
Zv - measured Zv - computed
0.8 0.6 F or Zv (MN)
414
0.4 0.2 0 0
4
8
12
16
20
24
28
32
36
40
-0.2
3.2. Signal matching analysis
-0.4
Stage construction option has been used to consider the effect of gravitational stresses during the analysis. In the first step, the soil layers are allowed to reach equilibrium under self-weight stresses, and then the resulted deformations are set to zero. The effect of fill is introduced to the model by a uniform stress of 38 kPa for the 2 m top fill
Soil layer
Friction angle
Cohesion Elastic modulus Poisson’s Density (kPa) (MPa) ratio (kN/m3)
1 2 3 4
36 20 30 36
10 10 10 10
0.3 0.3 0.3 0.3
-0.6 t=0
Time (msec)
Fig. 10. Measured force, and measured and computed Zv wave traces.
Table 2 Interface strength and deformation parameters from signal matching analysis
Table 1 Soil parameters for different layers from CPT correlations
20 20 20 40
t = 2L/C
19 19 19 19
Soil layer
Interface strength parameters Interface deformation parameters Friction angle
Cohesion (kPa)
Shear stiffness, Normal stiffness, Ks (kN/m3) Kn (kN/m3)
1 2 3 4
5 5 25 30
0 0 10 10
50,000 50,000 50,000 50,000
5,000,000 5,000,000 5,000,000 5,000,000
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
angle) parameters were adjusted in a limited zone around the pile tip. Local damping parameters considered for the pile and soil are 3% and 5%, respectively. The final parameters adjusted for soil surrounding the pile tip are presented in Table 3.
0
The peak energy is attained at time 2L/C after which the energy is decreasing and eventually approaching a residual value at about 36 ms. Thereafter, the energy is constant as the force and velocity are almost zero. The reduction of energy between peak to residual is attributed to elastic rebound of the pile tip and shaft, known in practice as soil quake. A relatively good correlation is observed between the two measured and computed curves, indicating a good signal matching effort. The pile head and tip displacements with time for the impact blow resulted from signal matching analyses are
Table 3 Strength and deformation parameters of soil around pile tip from signal matching analysis Soil block dimensions below pile tip
90 · 90 cm
Strength parameters
Deformation parameters
Friction angle
Dilation angle
45
Cohesion (kPa) 20
6
Elastic modulus (MPa) 500
9 8 Enthru Energy (kN.m)
Tables 2 and 3 present the final interface and soil parameters, to achieve a reasonable match between the measured and computed Zv signals shown in Fig. 10. The parameters obtained for around the pile tip from signal matching analysis (Table 3), indicate high magnitudes for deformation parameters (E of 500 MPa). This is attributed to the fact that the pile tip is embedded in a sandy deposit having a high densification potential during the driving efforts. On the basis of an analytical study by Yang [20], the influence zone limit of an axially loaded pile in clean sand and compactable silty sand are 3.5D–5.5D and 1.5D–3D, respectively. The signal matching trial and error efforts indicated that a cube with depth 3D (90 cm) below the pile tip elevation is providing a good match. The cube depth is in accordance with the influence zone ranges proposed by Yang [20] for the silty sand. In order to investigate the other aspects of the system response, further results are presented and discussed below. The applied energy to the pile head due to hammer impact is computed versus time and compared with the measured energy (Enthru) from dynamic testing as shown in Fig. 11. The Enthru was computed using Eq. (1) which is in fact the accumulation of force multiplied by displacement at the elevation of attached gages. The displacement is in fact calculated from v Æ dt Z t W ¼ F v dt: ð1Þ
t = 2L/C 10
7 6 5 4 3 2
Measured Enthru
1
Computed Enthru
0 0
4
8
12
16 20 24 Time (msec)
28
32
36
40
Fig. 11. Measured and computed applied energy (Enthru) versus time.
shown in Fig. 12. The difference between the two curves at peak (about 3.5 mm) is equivalent to the pile elastic compression, which is eliminated when the permanent displacement after 36 ms has reached. The figure also shows that the permanent displacements at both tip and head are 5 mm. The reported average penetration per blow during driving in the last 25 cm is 6 mm (41 blows/25 cm) that shows a very good correlation with the computed 5 mm permanent penetration. The other important observation in Fig. 12 is the elastic pile tip displacement of 2.5 mm (difference between pile tip peak displacement and the permanent displacement). This 2.5 mm displacement is in fact the same as, so-called, toe (tip) quake that has been suggested by Goble et al. [21] to be considered as D/120 that compares very well with the calculated value in which D is the pile diameter. Some results are presented and discussed in this part to briefly touch the soil inertia or the so-called radiation damping effects on the system response. The traveling com15 Pile Head Pile Displacement (mm)
4. Results
415
Pile Tip 10
5
0 0
5
10
15
20
25
30
35
40
Time (msec)
Fig. 12. Pile head and tip displacements versus time.
45
50
416
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
pressive or tensile waves along the pile shaft cause relative displacement between pile and soil, which in turn results in generation of shear wave in the adjacent soil that can propagate radially.pThe ffiffiffiffiffiffiffiffiffi propagation speed of such a wave is equivalent to G=q, in which G and q are shear modulus and mass density of the soil, respectively. Figs. 13 and 14 present the shear stress variations with time at a depth of 12 m and at different distances from the pile skin (0, 0.5, 1, 2, 3.6 m), with and without considering soil viscous damping, respectively. The shear wave propagation speed at such a depth is about 64 m/s using the above equation and the soil layer properties in that depth (silty sand of layer 3) as follows: E 20 MPa ¼ ¼ 7:69 MPa; 2ð1 þ tÞ 2ð1 þ 0:3Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7:69 106 Pa C¼ ¼ 63:6 m=s: 1900 kg=m3
G¼
ð2Þ ð3Þ
This is in accordance with the time required to achieve the peak shear stresses of Figs. 13 and 14 at each specified distance (0.5, 1, 2, 3.6 m) from the pile skin. Both figures represent rapid reduction of shear stress wave amplitudes 70 60 Embedment Depth=12m Shear Stress (kPa)
50 Pile skin 40
0.5m
30
1m
20
2m 3.6m
10 0 0
10
20
30
40
50
60
70
-10 -20
with distance from pile skin, such that near the vertical boundaries it is practically zero. Comparison between the results of Figs. 13 and 14 indicates that the viscous damping reduces the peak shear stress mobilized at each distance. Results of both Figs. 13 and 14 indicate that shear stress wave propagation is important up to a radial distance of about 2 m, which is equivalent to about six pile diameter. Therefore the selected vertical boundaries of 4 m are adequate as was pointed out before. In order to evaluate the compressive and or the so-called p-wave propagations, the results of normal stress variations with time at the pile tip (pile base), 0.5, 1, 2, and 3.6 m below pile tip are depicted in Fig. 15. When the hammer impact wave arrives at the pile tip, a compressive peak (dynamic) normal stress of 5600 kPa is calculated that drops and dissipates rapidly with distance from the pile tip. At a distance of 1 m (about 3 pile diameter) it is reduced to 2600 kPa, which is already much less than that at the pile tip. At distances over 2 m (6-pile diameter), the effect is negligible. It is reminded that the soil elastic modulus within a depth of 90 cm below pile tip was substantially increased to compensate for the increase in soil density due to driving effects. The double increase of elastic modulus and area (block) is in fact a substantial increase in the soil impedance below the pile tip. The last part of the presented results includes comparison between static load test results and static response predictions. Fig. 16 shows the static axial load versus pile head displacement both from the static load test in the field and the simulated results by the numerical model. The analysis was performed using the same strength and deformation parameters obtained for soil and interface, from signal-matching analyses described in the former subsections. It is observed in Fig. 16 that up to a pile head displacement of 30 mm, a reasonably good correlation exists between prediction and static load test. Assuming a failure criterion for the pile at a displacement equivalent to 10%
Time (msec)
Fig. 13. Variation of shear stress in soil at different distances from pile shaft, considering soil viscous damping.
6000 Pile toe 5000
70
1m
Embedment Depth=12m
Normal Stress (kPa)
60 Pileskin
50 Shear Stress (kPa)
0.5m below toe
0.5m
40
1m
30
2m
20
3.6m
10 0 0
10
20
30
40
50
60
70
4000
2m 3.6m
3000
2000
1000
-10 -20 Time (msec)
Fig. 14. Variation of shear stress in soil at different distances from pile shaft without considering soil damping.
0 0
10
20
30 40 Time (msec)
50
60
70
Fig. 15. Variation of normal stress at different distances below pile tip.
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418 Model prediction
Static load test
2000
Load (kN)
1800 1600 1400 1200 1000 800 600 400 200 0 0
10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 Pile head displacement (mm)
Fig. 16. Static load test simulation result compared with the real test result using the parameters obtained from signal matching analysis.
of the pile tip diameter, or in fact 30 mm, the ultimate pile capacity is scaled off as 1200 kN. This is indicating the reasonable accuracy of the numerical model in the static capacity prediction. The distribution of the shear stress on the pile shaft on the basis of the static load test simulation is shown in Fig. 17. The tip and shaft resistances are thus obtained as 650 kN and 550 kN, respectively. The trend of shear stress variations on the pile shaft compares reasonably well with the qc variations of the CPT profile of Fig. 7. 5. Discussion A successful so-called continuum numerical model is developed which proved to perform very well for the wave propagation and dynamic analysis of piles subjected to hammer blows. The developed model is superior to the conventional Smith-type models used in practice for highstrain test result interpretations, in that it more realistically considers the surrounding soil effects. Still there is a long
Shear stress (kPa)
Depth (m)
0
10
20
30
40
50
60
70
80
90
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
Fig. 17. Calculated distribution of shear stress on the pile shaft; compares reasonably well with the CPT profile.
417
way to go on this path to explore and manipulate its merits, however. One important factor that was not considered in this study is the existence of the ground water table (GWT) that has been reported at about 1 m below ground surface for the case pile studied. It is obvious, and of great importance, that the GWT effects should be considered in the analysis. The dominating layers at the pile shaft and at the pile tip are sandy, therefore the effective stress conditions are prevailing, even during the initial driving. But, why this was not considered in the current study, and how its effect has been compensated for, are briefly discussed below. GWT causes the effective stress to reduce to about half, and therefore the shaft resistance reduces to the same degree. According to the b method, the shaft resistance at any depth z is rs ¼ K tgd ðc0 ZÞ:
ð4Þ
It should be noted that in this fundamental equation, K is significantly affected by driving effort of the pile in sandy soils. Around the pile shaft and below the pile tip, the sand is usually densified to some distance from the pile. Densification around the pile shaft contributes to substantial increase in K. In the employed numerical model, however, the pile was embedded within the soil as if it is a bored pile. No attention has been paid on the driving effects yet, in the developed model. In fact, it is not so easy to implement the densification effect of the driving effort in the numerical model. To compensate for this lack of consideration of increase in K due to driving effort, the total stresses are considered in this analysis, meaning application of c instead of c’ (19 kN/m3 instead of 9 in our example case). This is more or less compensating for the increase in K due to driving effects. Otherwise, we would have required using d (interface friction angle) values very higher than the meaningful values, to match the predicted and measured Zv signals. Similar discussion is applicable to the tip resistance. It is reminded that the GWT results in effective stress reduction around the pile tip, and therefore reduction in the pile tip capacity. The reduction in tip capacity and stiffness would cause in major reduction in impedance in the soil zone below the pile tip, and therefore major changes in the F and Zv reflections at 2L/C. However, we had to substantially increase the elastic modulus of a soil block below the pile tip to match the predicted and measured signals. This is attributed to the fact that the sand mass density is significantly increased during the driving efforts of the pile, hence increasing the stiffness and strength parameters of the soil around and below the pile tip. Therefore, the balance between considering versus not considering the effective stress reduction due to GWT is adjusted by the stiffness and strength parameters of the block below the pile tip in the numerical model. Had we considered the effective stress conditions, we would have ended up in higher values for the stiffness and strength
418
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
parameters for the block to match the predicted and measured signals. This shortcoming has to be dealt with in future studies. However, as this study is one of the few attempts towards the application of a so-called continuum numerical model for the signal matching analyses, the effective stress analysis is not considered yet. The main focus has been on the successful operation of the model. It should be pointed out that also in the conventional mass–spring– dashpot models the effective stresses are not distinguished. In fact during the signal matching trial and error process, different parameters are varied such that the best match to be found, irrespective of the effective or total stress issues. 6. Conclusions The main objective of this numerical study was to develop a so-called continuum numerical model to simulate wave propagation and dynamic response of the pile driving problem and using the model for signal matching analyses. The model was successfully developed using a standard finite difference code. The findings of the study are summarized below: 1. If meaningful parameters are considered for the pile material, and the boundary conditions are introduced properly, the axisymmetric solid elements can adequately simulate the wave propagation along the shaft, and the reflections due to the impedance changes at the pile tip. 2. The inertial and/or radiation damping of the surrounding soil mass is automatically considered in this continuum numerical model. This is the most important contribution of the developed numerical model versus the mass–spring–dashpot type models. However, the degree of importance and outcome of such advantages still need to be manipulated and quantified in the future works. 3. The force and velocity signals were matched by varying the strength and deformation parameters of the soil and interface element. The model then proved to be very well in predicting other quantities such as Enthru, pile toe (tip) quake, static load–displacement curve and the load distribution on the pile shaft. 4. The effective stress type analysis, non-linear hardening stress–stress relations for soil, and sensitivity analysis of different parameters are the topics that require further studies to explore the capabilities of the new continuum numerical model.
References [1] Randolph MF. Science and Empiricism in Pile Foundation Design. 43rd Rankine Lecture. Geotechnique 2003;54(1). [2] Salgado R, Abou Jaoude G. Assessment of axially-loaded pile dynamic design methods and review of INDOT axially-loaded pile design procedure. Proposal for Research and Implementation Study, Purdu University, USA; 2003. [3] Smith EAL. Pile-driving analysis by the wave equation. J Soil Mech Found Div 1960;86(EM 4):35–61. [4] Coyle HM, Lowery Jr LL, Hirsch TJ. Wave equation analysis of piling behavior. In: Desai CS, Christian JT, editors. Numerical methods in geotechnical engineering. McGraw-Hill; 1977. [5] Fakharian K, Feizee S. Application of wave equation analysis (WEAP) for design of driven piles. In: Proceedings of 4th Int Conf on Coasts, Ports & Marine Structures (ICOPMAS 2000), Bandar Abbas, Iran; November 2000. [6] Feizee S. Application of wave equation analysis (WEAP) for design of driven piles. MSc Thesis, Department of Civil & Environmental Engineering, Amirkabir University of Technology, Tehran, Iran; 2000. [7] Rausche F. Soil response from dynamic analysis and measurements on piles. Thesis presented to the Case Western Reserve University, at Cleveland, Ohio, in 1970, in partial fulfillment of the requirements for the degree of Doctor of Philosophy; 1970. [8] Rausche F, Goble GG, Likins G. Dynamic determination of pile capacity. J Geotech Eng 1985;111(3):367–83. [9] Likins GE, Rausche F. Correlation of CAPWAP with static load tests. In: Proc 7th Conf Application of stress wave theory to piles. Petaling-Jaya, Selangor, Malaysia; 2004. pp. 153–165. [10] Coutinho ALGA, Costa AM, Alves JLD, Landau L, Ebecken NFF. Pile driving simulation and analysis by the finite element method. In: Proc 3rd Int Conf, Application of stress-wave theory to piles, Ottawa, Canada; 25–27 May 1988. [11] Borja RI. Dynamics of pile driving by the finite element method. Comput Geotech 1988;5(11):39–49. [12] Mabsout M, Tassoulas JL. A finite element model for the simulation of pile driving. Int J Numer Methods Eng 1994;37:257–78. [13] Mabsout ME, Reese LC, Tassoulas JL. Study of pile driving by finiteelement method. J Geotech Geoenviron Eng, ASCE 1995;121(7): 535–43. [14] Kolsky H. Stress waves in solids. New York: Dover Publications; 1963. [15] Geerling J, Smits MThJH. Prediction of load–displacements characteristics of piles from the results of dynamic/kinetic load tests. In: Proc 4th conf application of stress wave theory to piles. Delft University, Netherlands; 12–24 September 1992. [16] Heijnen WJ. General soil conditions of the test sites. In: Proc 4th conf application of stress wave theory to piles. Delft University, Netherlands; 12–24 September 1992. [17] Nath B. A continuum method of pile driving analysis: comparison with the wave equation method. Comput Geotech 1990;10:265–85. [18] Lysmer J, Kuhlemeyer RL. Finite dynamic model for infinite media. J Eng Mech 1969;95(EM4):859–87. [19] Bowles JE. Single piles: dynamic analysis, load tests, foundation analysis and design. McGraw-Hill; 1996. [20] Yang J. Influence zone for end bearing of piles in sand. J Geotech Geoenviron Eng 2006;32(9):1229–37. [21] Goble GG, Rausche F, Likins G. Associates Inc. Introduction to GRLWEAP, the 1997-2 GRLWEAP PROGRAM. Cleveland, Ohio; 1996.
Computers and Geotechnics 35 (2008) 406–418 www.elsevier.com/locate/compgeo
Application of a continuum numerical model for pile driving analysis and comparison with a real case Shahram Feizee Masouleh, Kazem Fakharian
*
Department of Civil and Environmental Engineering, Amirkabir University of Technology, Tehran 15875, Iran Received 4 February 2007; received in revised form 24 July 2007; accepted 28 August 2007 Available online 29 October 2007
Abstract This paper presents the results of a so-called continuum numerical model for wave propagation analysis and soil-pile dynamic response during pile driving. An axisymmetric finite difference numerical model is developed having solid elements for both pile structure and the soil media surrounding and below the pile. Interface elements are used between the pile shaft and the soil to facilitate the sliding between the two media. The performance of the developed model is verified in two stages. First, a simple rod is subjected to a half sinewave force function at the rod head and the corresponding reflections of force and velocity (multiplied by impedance) are presented for different boundary conditions at the rod tip. The model is then used for signal matching analysis of a real driven pile for which complete information of soil layering, dynamic test signals, and static load test results are available. The signal matching analysis was performed successfully and comparison between several other predicted and measured parameters proved the reasonably good performance of the developed continuum model. 2007 Elsevier Ltd. All rights reserved. Keywords: Pile driving; Stress wave analysis; Continuum numerical model; Finite difference method; Signal matching analysis
1. Introduction Despite the fact that pile design approaches have advanced enormously within the past few decades, but yet, the most fundamental aspect of pile design, that of estimating the axial bearing capacity, relies heavily on empirical correlations as discussed by Randolph [1]. The piles are classified in two main groups in terms of their installation method, including non-displacement piles and displacement piles. In non-displacement piles, in case of utilizing casing or drilling mud as well as proper construction conditions, the initial state of stresses are relatively maintained adjacent to the pile shaft and below tip [2]. In displacement piles such as driven piles, on the other hand, the pile is pushed down by hammer impact or jack. The
*
Corresponding author. Tel.: +98 21 22 90 5633; fax: +98 21 22 27 2000. E-mail address: [email protected] (K. Fakharian). 0266-352X/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2007.08.009
changes in stress state such as degree of compaction around shaft or below tip are difficult to measure, thus making it more difficult to estimate the axial pile capacity. Many attempts have been made for developing an adequate dynamic method to accurately estimate the static axial capacity. The main objective of the dynamic methods has been to establish a relationship between the blow counts and the axial capacity of the pile. Many dynamic formulas have been developed and are still used by some contractors during construction, but the most important innovation was the introduction of one-dimensional wave equation analysis in piles (WEAP) by Smith [3], in which pile is simulated by a number of masses attached to each other by elastic springs, and soil is replaced by a number of springs, sliders and linear viscous dampers to simulate the visco, elasto-plastic response of the soil. The WEAP model, however, requires to estimate some key soil parameters such as quake and damping coefficients as well as hammer efficiency. Such factors significantly influence the analysis results [4–6].
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
407
Nomenclature Notations A pile cross section area (m2) C wave speed in pile (m/s) D pile diameter (m) E elastic modulus (N/m2) F force (N)
Rausche [7,8] proposed the signal matching technique and developed an analysis method, CAPWAP, to overcome the shortcomings of WEAP model. The signal matching technique relies on the pile dynamic test results (PDA). The model used in CAPWAP is very similar to WEAP, except that unloading effects are also considered, though such effects on capacity are insignificant [9]. The signal matching analysis is the most reliable dynamic method in practice. This method, however, uses mass–spring–dashpot system to model the surrounding and below tip soil media, therefore imposes some limitations. Pile is a continuum body having infinite number of degrees of freedom surrounded by a continuum media that is soil. In the mass–spring system, the boundary conditions are different than the real system, causing deviations in the system response. Also the effect of soil inertia (also called ‘‘radiation damping’’ or ‘‘geometrical damping’’) can not be easily considered. To overcome the above shortcomings of mass–spring– dashpot models, different attempts have been made to make use of other numerical solutions. Randolph tried to implement the far field effect of soil using additional dashpot (so-called inertial dashpot representing the surrounding soil inertia or radiation damping), as shown in Fig. 1 [1]. Some other studies have tried to use finite element method so that the soil media around the pile shaft and below the pile tip are directly considered in the model. Coutinho et al. [10] conducted a research on the drivability of piles into the Brazilian coast deposits. They considered the effect of some parameters influencing the pile driving including pile-soil interface and damping. Borja [11] adopted a FEM scheme to model the pile driving system, capable of continuously interpolating displacement, velocity and acceleration profiles along the pile length. As opposed to the previous models that were using the impact velocity of hammer as initial condition, Borja utilized the impact force curve of the hammer with respect to time. This allowed a more realistic simulation of the pile–hammer interaction. Some results were presented on the basis of analyses on an assumed example problem, and no case study was used for comparison with the numerical predictions. Mabsout and Tassoulas [12] and Mabsout et al. [13] also used an axisymmetric FEM model (Fig. 2) for pile driving analysis. The surrounding soil was assumed to be saturated normally-consolidated clay and its stress–strain behavior
Kn Ks L v Z q
interface normal stiffness (N/m3) interface shear stiffness (N/m3) pile length (L) velocity (m/s) impedance (EA/C) mass density (kg/m3)
Fig. 1. Dynamic soil-pile shaft model to account for inertial damping. (After Randlolph, 2003).
was simulated by a non-linear bounding-surface-plasticity model. The slide-line algorithm was used for the interface modeling to enable the soil-pile relative displacements. The at-rest lateral stress coefficient was considered equal to unity. Viscous boundaries were used to absorb the radial stress waves, resulted from hammer impact on the pile. The hammer impact was modeled as a force–time function on the pile head. The transient dynamic analysis was performed using the constant-average-acceleration algorithm. A bored pile with embedment depths of 6, 12 and 18 m was assumed for the driving analysis. The tip and shaft resistances were distinguished and the results showed higher resistances with increase in pile embedment depth. It was also observed that resistance to driving increases with increase in number of blows. The excessive pore pressure at pile tip was shown to dramatically decrease with distance from the pile. No comparison was, however, made between the predictions and a real-life pile driving data. It was suggested, however, as their conclusions to verify the predictions of the FEM model against actual pile driving data. The main objective of this paper is to develop a so-called continuum numerical model to overcome some of the limitations of the mass–spring–dashpot system, and present an improvement in comparison with the previous FEM studies, in that a signal matching procedure to be adopted for comparing the predictions of the model with a well-docu-
408
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
Lagrangian calculation which is adequate for large deformation modeling. One of its advantages is possibility of modeling large systems with a limited memory, as no stiffness matrix is required (as opposed to FEM, for example). In addition, large deformations do not significantly increase the run-time as no stiffness matrix updates are required after each load or time increment. The wave propagation analysis in an elastic rod with different end (tip) boundary conditions is conducted first to verify the model performance and to evaluate the sensitivity of some parameters. Three series of analyses are performed on: (1) a rod with fixed-end (tip) condition (2) a rod with free-end (tip) condition (3) a rod with elastic support at the end (tip). No soil or other supports are considered around the rod shaft yet. Fig. 3 shows the axisymmetric mesh and the boundary conditions of the fixed-end, freeend, and elastic support, in 3a, 3b and 3c, respectively, for a 0.5 m diameter, 20 m long rod. More details of the model and results are presented in the following sections. 2.1. Fixed-end rod
Fig. 2. Continuum soil-pile model for simulation of pile driving. (After Mabsout, 1995).
mented actual driven pile. Among the major advantages of the model is the consideration of soil inertia and/or radiation damping effects. An axisymmetric finite difference model is developed using FLAC-2D with its axisymmetric capability to simplify the 3D model to 2D. For preliminary verification of the model performance, an elastic rod with different boundary conditions at the tip (fixed, free, variable elastic modulus) is considered for wave propagation analysis and response from tip. To further investigate the model capability, a real driven pile that was carefully monitored in Netherlands is considered for modeling, for which dynamic and static test results, the soil classification and CPT profiles are available. A signal matching procedure (similar to CAPWAP) was followed by changing the strength and deformation parameters of soil and interface between pile-soil. The fine-tuned parameters were then used to simulate the static axial load–displacement response to be compared with the real test results. 2. Numerical model and preliminary verification A finite difference scheme is adopted to simulate an axisymmetric model for pile driving analysis, using FLAC-2D version 4. This software is developed on the basis of
In the first series of analyses, a fixed-end round rod is assumed having length, diameter, elastic modulus, Poisson’s ratio, and unit weight of 20 m, 0.5 m, 25,000 MPa, 0.15, and 25 kN/m3, respectively. An axisymmetric mesh of 5 · 5 cm is used for modeling the rod half-section, as shown in Fig. 3. The rod gravity is neglected in the analysis. The hammer impact on the rod head is simulated by a half-sine stress wave with an amplitude of 5 MPa and frequency of 320 Hz. Regarding pffiffiffiffiffiffiffiffiffi that the wave speed in the rod is 3162 m/s ðc ¼ E=qÞ, the wave length of the halfsine signal generated in the rod is 4.94 m. Therefore, the force and velocity records are picked up at 5 m below the rod head, to more clearly represent the end (tip) dynamic responses, and prevent the mixing up of the upward stress wave and downward reflection from the rod head. The analysis results for a rod without damping effects are presented in Fig. 4a. The force (F) and velocity (v) variations with time are plotted in which velocity is multiplied by the rod impedance (Z = EA/C) to be equivalent to force (Zv). This is the same procedure used in PDA test output presentations. The horizontal axis is scaled in L/C, which is the time required for stress wave to travel between the pick up records and the rod end (tip). The force sign conventions are positive and negative for compressive and tensile waves, respectively. Downward displacement of a particle on the rod section shows a positive velocity, and upward displacement shows a negative velocity. Fig. 4a shows that after a time equivalent to 2L/C, that is the time required for stress wave to travel between the record point at 5 m below the rod head to rod tip and return to the same point (L = 15 m), the F wave is suddenly shifted up (positive) and Zv wave is shifted down (negative) and their amplitudes are equal to the generated wave amplitude. This observation is exactly in accordance with the 1D longitudinal wave propagation theory in rods
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
CL
CL
409
CL
Mesh size: 5x5cm
Elements with Variable Elastic Modulus
Viscous boundary
(a)
(b)
(c)
Fig. 3. Axisymmetric mesh for a rod with different boundary conditions at tip: (a) fixed-end; (b) free-end; (c) elastic support at the end.
a
b 1200
1200 F
F 800
Zv
400 0 0
2
4
6
-400
8
10
Force (kN)
Force (kN)
800
Zv
400 0 0
2
4
6
-400 -800
-800
-1200
-1200 L/C
L/C
Fig. 4. Dynamic response of the fixed-end rod: (a) without viscous damping; (b) with viscous damping.
8
10
410
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
without damping. The immediate F and Zv waves shifted down (negative) after tip reflections, are the reflections from the rod free head (top) boundary condition. In fact the downward initial compressive wave is reflected compressive at rod fixed-end (tip) and reflected tensile at the rod free-top (head) boundary conditions. Fig. 4b represents similar results for the same rod and loading, but with a 1% assumption of viscous damping ratio for rod. The ‘‘local damping’’ option of the model is activated that is on the basis of increase or decrease in mass of mesh nodes at specific times during oscillations, such that when velocity sign is changing, mass is increased, and at minimum and maximum velocities, the mass is reduced [14]. The F wave and Zv wave amplitudes are attenuated with time due to the damping effects of rod, as expected. 2.2. Free-end rod In the second series of analyses, a free-end (tip) rod is assumed, having the same geometry and material properties of the previous example. The results without and with damping are presented in Fig. 5a and b, respectively. Similar, but opposite in sign, variations are observed as expected, which is in fact representative of the 1D longitudinal wave propagation in a rod with free-end boundary conditions at rod tip and head. 2.3. Elastic support at rod tip In the third series of analyses, an elastic spring was assumed at the pile tip. The spring was simulated by using solid elements similar to the ones used for the rod itself with the same total cross section area of the rod (0.196 m2), but with a length of 1 m (Fig. 3c). The stiffness and impedance of the 1-m rod (acting as a spring) was changed by varying its elastic modulus from 10 MPa to 10,000 GPa, to see the effect on the dynamic response due to the impact on the rod head. A viscous boundary condition was imposed at the bottom of the 1-m rod to absorb the portion of
b
1200
1200
800
800
400
400
0 0
2
4
6
8
10
-400
Force (kN)
Force (kN)
a
the wave not reflected at the tip. The same half-sine stress wave was applied at the rod head to simulate the dynamic impact. The results are presented in Fig. 6a–h, for E values of the 1-m segment varying (decreasing) between 10,000 GPa to 10 MPa, respectively. For the highest value of E (10,000 GPa) in Fig. 6a, the response is identical to the fixed-end rod of Fig. 4a. For the lowest value of E (10 MPa) in Fig. 6h, on the other hand, the response is the same as the free-end rod of Fig. 5a. A gradual transition is observed between compressive to tensile wave responses, when the elastic modulus varies from high to low magnitudes. It is interesting to note that when the elastic modulus is equal to 30 GPa (Fig. 6d), which is very close to the elastic modulus of the rod, both force and velocity reflections are very close to zero, i.e. no significant reflections are observed. This is due to the fact that at such an elastic modulus, the impedance (EA/C) of both the rod and the segment are equivalent to each other, and therefore, no wave is reflected. It is reminded from the basics of 1D wave propagation that the magnitude of the reflected wave is proportional to the impedance variations along its propagation path. For impedances higher than the rod (Fig. 6a–d), the reflections are compressive, whereas it is the other way around for lower impedances (Fig. 6e–h). The observations very well confirm the performance of the numerical model, as a base for more complicated analyses of a real pile embedded within various soil layers below the pile tip and around the pile shaft. The results of the three series of analyses presented above on a 20 m rod with open surroundings, but different boundary conditions at tip (free, fixed, constant impedance), indicates that the developed model is capable of simulating the pile driving analysis and the wave propagation and reflections from tip and head. Both force (F) and velocity (Zv) signals versus time are identifiable and promising to be used in a signal matching type analysis similar to CAPWAP, for example. To further evaluate the developed model capabilities, a real pile driven into a layered soil is considered for which complete site conditions and dynamic and static test results are available.
0 0
2
4
6
Zv
-1200
-800
Zv
- 1200 L/C
10 F
F -800
8
-400
L/C
Fig. 5. Dynamic response of the free-end rod: (a) without viscous damping, (b) with viscous damping.
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
411
E = 1000 GPa
E = 10000 GPa 1200
1200 F
800
F
800
Zv
400 0 0
1
2
3
4
- 400
Force (kN)
Force (kN)
Zv 400 0 0
- 400
- 800
- 800
-1200
-1200
1
2
3
L/C
L/C
E = 30 GPa
E = 100 GPa 1200
1200 F
800
F
800
Zv
400 0 0
- 400
1
2
3
4
Force (kN)
Force (kN)
Zv 400 0 0
- 400
- 800
- 800
-1200
-1200
1
2
3
E = 1 GPa
E = 10 GPa 1200
1200 F
800
F
800
Zv
400 0 0
1
2
3
4
Force (kN)
Force (kN)
Zv 400 0 - 400
0
1
3
4
-1200
-1200
L/C
L/C
E = 10 MPa
E = 100 MPa 1200
1200 F
800
0 1
2
3
4
- 800
Force (kN)
Zv
400
0
F
800
Zv
- 400
2
- 800
- 800
Force (kN)
4
L/C
L/C
- 400
4
400 0 - 400
0
1
2
3
4
- 800
-1200
-1200 L/C
L/C
Fig. 6. Effect of the support elastic modulus at tip on dynamic response of the rod.
3. Case study pile The case pile is selected from an elaborate pile testing research project along with the ‘‘4th International Conference on the Application of Stress-wave Theory to Piles’’ held in Netherlands, 1992 [15,16]. Several precast concrete piles were driven at the test site behind the auditorium of Delft Technological University. Static and dynamic test results are available for the selected test pile used in this study.
The soil profile of the site together with CPT results are shown in Fig. 7. The pile is precast concrete with 0.25 m dimension and elastic modulus of 36,000 MPa. The 19 m long pile was driven by a hydraulic hammer ICE SC-40, to an embedment depth of 18 m. The blow counts for the last 1 m of driving for 25 cm penetrations are 37, 29, 32 and 41, respectively. The force (F) and velocity (Zv) signals recorded by installing strain gages and accelerometers at 1 m below the pile head and related to the last blow at the end of pile driving are available. The results of a static
412
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
having an average correlation with both the cross section area and the circumferential area of the main pile, as shown in Fig. 8. The pile tip has penetrated about 3 m into a medium dense to dense sand (so called Pleistocene sand). Except a 4 m thick clay and peat layer (#2) along CL
19m
Layer 2 Clay and Peat
3m
Layer 1 Silty Sand
4m
Fill
2m
Pile
6m
Layer 4 Pleistocene Sand 4m
Fig. 7. Soil profile and CPT results [16].
Layer 3 Silty Sand
7m
Interface
load test performed after 3 weeks of the initial driving is also available. 3.1. The numerical model 4m
In order to benefit from the axisymmetric feature of the 2D numerical model, a 0.3 m diameter pile is assumed
Fig. 8. Axisymmetric pile-soil model for pile driving simulation.
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
CL
5x5cm Pile Elements (Linear Elastic) Uniform Stress=38 kPa (2m fill)
Embedment Depth
Pile-Soil Interface
Pile Toe Viscous Boundary
10x10cm Soil Elements
5x5cm Soil Elements
Viscous Boundary Fig. 9. Finite difference mesh for the axisymmetric pile-soil model.
413
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
the shaft, rest of the layers are mostly sand, but with relatively low qc values from CPT test. The boundary distances in horizontal and vertical directions, from shaft and pile tip, respectively, are 4 m in the axisymmetric mesh. The selected distances are far enough to minimize the stresses resulted and propagated from pile driving, as discussed in the following sections. In accordance with the highest stress gradients near pile shaft and the soil around the pile tip, the mesh sizes are specified 5 · 5 cm within 2 m distance from shaft and tip, and 10 · 10 cm outside this zone to the vertical and horizontal boundaries (Fig. 9). Interface elements between pile shaft and soil are inserted to facilitate the pile-soil slip during driving. The Mohr–Coulomb failure criterion is used for the interface, on the basis of which slip between soil and pile occurs, if the generated shear stress during driving exceeds the shear strength. As a result of hammer blow on the pile head, part of the stress wave enters and propagates within the surrounding soil. The degree and magnitude of the penetrated energy depends on the interface and the pile-soil mechanical parameters [17]. In real pile driving, the penetrated stress waves propagate away from the pile and dissipate gradually, which is referred to as ‘‘radiation damping’’. As the model boundaries are limited, viscous boundaries are assumed at the vertical and bottom horizontal boundaries to absorb such waves and prevent their reflection. Dampers in two perpendicular directions are implemented at each boundary node as proposed by Lysmer [18]. An elasto-perfectly plastic model with Mohr–Coulomb failure criterion is used for the soil material. The required mechanical parameters for each soil layer are estimated on the basis of CPT results with common correlations available in literature [19]. The estimated parameters are presented in Table 1. The cohesion and Poisson’s ratio are assumed 10 kPa and 0.3, respectively, for all the four layers.
(Fig. 9). In the second step, the interface elements with zero shear strength are considered between pile-soil to allow the pile to displace freely under self-weight. For performing the pile driving dynamic analysis, the recorded force–time at pile head by the strain gage (resulted from pile driving) is input on the pile, and then the velocity signal is calculated from the analysis. The calculated Zv at the pile head is compared with the measured Zv, and if not matched, efforts are made in subsequent analyses by trial and error, similar to the CAPWAP type analyses, to match the signals. The measured force, measured Zv, and the fine-tuned computed Zv are shown in Fig. 10. The most important influencing time domain for the shaft resistance adjustments is between peak force (or velocity) to 2L/C (about 9.5 ms in the study case). Therefore, trial and error were made to match the computed and measured velocity signals, by adjusting the interface strength and deformation parameters for different layers. The final parameters are presented in Table 2. In order to adjust the pile tip effects in the signal matching process, the soil strength (angle of friction and cohesion) and deformation (elastic modulus and dilation
1.2 Force - measured
1
Zv - measured Zv - computed
0.8 0.6 F or Zv (MN)
414
0.4 0.2 0 0
4
8
12
16
20
24
28
32
36
40
-0.2
3.2. Signal matching analysis
-0.4
Stage construction option has been used to consider the effect of gravitational stresses during the analysis. In the first step, the soil layers are allowed to reach equilibrium under self-weight stresses, and then the resulted deformations are set to zero. The effect of fill is introduced to the model by a uniform stress of 38 kPa for the 2 m top fill
Soil layer
Friction angle
Cohesion Elastic modulus Poisson’s Density (kPa) (MPa) ratio (kN/m3)
1 2 3 4
36 20 30 36
10 10 10 10
0.3 0.3 0.3 0.3
-0.6 t=0
Time (msec)
Fig. 10. Measured force, and measured and computed Zv wave traces.
Table 2 Interface strength and deformation parameters from signal matching analysis
Table 1 Soil parameters for different layers from CPT correlations
20 20 20 40
t = 2L/C
19 19 19 19
Soil layer
Interface strength parameters Interface deformation parameters Friction angle
Cohesion (kPa)
Shear stiffness, Normal stiffness, Ks (kN/m3) Kn (kN/m3)
1 2 3 4
5 5 25 30
0 0 10 10
50,000 50,000 50,000 50,000
5,000,000 5,000,000 5,000,000 5,000,000
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
angle) parameters were adjusted in a limited zone around the pile tip. Local damping parameters considered for the pile and soil are 3% and 5%, respectively. The final parameters adjusted for soil surrounding the pile tip are presented in Table 3.
0
The peak energy is attained at time 2L/C after which the energy is decreasing and eventually approaching a residual value at about 36 ms. Thereafter, the energy is constant as the force and velocity are almost zero. The reduction of energy between peak to residual is attributed to elastic rebound of the pile tip and shaft, known in practice as soil quake. A relatively good correlation is observed between the two measured and computed curves, indicating a good signal matching effort. The pile head and tip displacements with time for the impact blow resulted from signal matching analyses are
Table 3 Strength and deformation parameters of soil around pile tip from signal matching analysis Soil block dimensions below pile tip
90 · 90 cm
Strength parameters
Deformation parameters
Friction angle
Dilation angle
45
Cohesion (kPa) 20
6
Elastic modulus (MPa) 500
9 8 Enthru Energy (kN.m)
Tables 2 and 3 present the final interface and soil parameters, to achieve a reasonable match between the measured and computed Zv signals shown in Fig. 10. The parameters obtained for around the pile tip from signal matching analysis (Table 3), indicate high magnitudes for deformation parameters (E of 500 MPa). This is attributed to the fact that the pile tip is embedded in a sandy deposit having a high densification potential during the driving efforts. On the basis of an analytical study by Yang [20], the influence zone limit of an axially loaded pile in clean sand and compactable silty sand are 3.5D–5.5D and 1.5D–3D, respectively. The signal matching trial and error efforts indicated that a cube with depth 3D (90 cm) below the pile tip elevation is providing a good match. The cube depth is in accordance with the influence zone ranges proposed by Yang [20] for the silty sand. In order to investigate the other aspects of the system response, further results are presented and discussed below. The applied energy to the pile head due to hammer impact is computed versus time and compared with the measured energy (Enthru) from dynamic testing as shown in Fig. 11. The Enthru was computed using Eq. (1) which is in fact the accumulation of force multiplied by displacement at the elevation of attached gages. The displacement is in fact calculated from v Æ dt Z t W ¼ F v dt: ð1Þ
t = 2L/C 10
7 6 5 4 3 2
Measured Enthru
1
Computed Enthru
0 0
4
8
12
16 20 24 Time (msec)
28
32
36
40
Fig. 11. Measured and computed applied energy (Enthru) versus time.
shown in Fig. 12. The difference between the two curves at peak (about 3.5 mm) is equivalent to the pile elastic compression, which is eliminated when the permanent displacement after 36 ms has reached. The figure also shows that the permanent displacements at both tip and head are 5 mm. The reported average penetration per blow during driving in the last 25 cm is 6 mm (41 blows/25 cm) that shows a very good correlation with the computed 5 mm permanent penetration. The other important observation in Fig. 12 is the elastic pile tip displacement of 2.5 mm (difference between pile tip peak displacement and the permanent displacement). This 2.5 mm displacement is in fact the same as, so-called, toe (tip) quake that has been suggested by Goble et al. [21] to be considered as D/120 that compares very well with the calculated value in which D is the pile diameter. Some results are presented and discussed in this part to briefly touch the soil inertia or the so-called radiation damping effects on the system response. The traveling com15 Pile Head Pile Displacement (mm)
4. Results
415
Pile Tip 10
5
0 0
5
10
15
20
25
30
35
40
Time (msec)
Fig. 12. Pile head and tip displacements versus time.
45
50
416
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
pressive or tensile waves along the pile shaft cause relative displacement between pile and soil, which in turn results in generation of shear wave in the adjacent soil that can propagate radially.pThe ffiffiffiffiffiffiffiffiffi propagation speed of such a wave is equivalent to G=q, in which G and q are shear modulus and mass density of the soil, respectively. Figs. 13 and 14 present the shear stress variations with time at a depth of 12 m and at different distances from the pile skin (0, 0.5, 1, 2, 3.6 m), with and without considering soil viscous damping, respectively. The shear wave propagation speed at such a depth is about 64 m/s using the above equation and the soil layer properties in that depth (silty sand of layer 3) as follows: E 20 MPa ¼ ¼ 7:69 MPa; 2ð1 þ tÞ 2ð1 þ 0:3Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7:69 106 Pa C¼ ¼ 63:6 m=s: 1900 kg=m3
G¼
ð2Þ ð3Þ
This is in accordance with the time required to achieve the peak shear stresses of Figs. 13 and 14 at each specified distance (0.5, 1, 2, 3.6 m) from the pile skin. Both figures represent rapid reduction of shear stress wave amplitudes 70 60 Embedment Depth=12m Shear Stress (kPa)
50 Pile skin 40
0.5m
30
1m
20
2m 3.6m
10 0 0
10
20
30
40
50
60
70
-10 -20
with distance from pile skin, such that near the vertical boundaries it is practically zero. Comparison between the results of Figs. 13 and 14 indicates that the viscous damping reduces the peak shear stress mobilized at each distance. Results of both Figs. 13 and 14 indicate that shear stress wave propagation is important up to a radial distance of about 2 m, which is equivalent to about six pile diameter. Therefore the selected vertical boundaries of 4 m are adequate as was pointed out before. In order to evaluate the compressive and or the so-called p-wave propagations, the results of normal stress variations with time at the pile tip (pile base), 0.5, 1, 2, and 3.6 m below pile tip are depicted in Fig. 15. When the hammer impact wave arrives at the pile tip, a compressive peak (dynamic) normal stress of 5600 kPa is calculated that drops and dissipates rapidly with distance from the pile tip. At a distance of 1 m (about 3 pile diameter) it is reduced to 2600 kPa, which is already much less than that at the pile tip. At distances over 2 m (6-pile diameter), the effect is negligible. It is reminded that the soil elastic modulus within a depth of 90 cm below pile tip was substantially increased to compensate for the increase in soil density due to driving effects. The double increase of elastic modulus and area (block) is in fact a substantial increase in the soil impedance below the pile tip. The last part of the presented results includes comparison between static load test results and static response predictions. Fig. 16 shows the static axial load versus pile head displacement both from the static load test in the field and the simulated results by the numerical model. The analysis was performed using the same strength and deformation parameters obtained for soil and interface, from signal-matching analyses described in the former subsections. It is observed in Fig. 16 that up to a pile head displacement of 30 mm, a reasonably good correlation exists between prediction and static load test. Assuming a failure criterion for the pile at a displacement equivalent to 10%
Time (msec)
Fig. 13. Variation of shear stress in soil at different distances from pile shaft, considering soil viscous damping.
6000 Pile toe 5000
70
1m
Embedment Depth=12m
Normal Stress (kPa)
60 Pileskin
50 Shear Stress (kPa)
0.5m below toe
0.5m
40
1m
30
2m
20
3.6m
10 0 0
10
20
30
40
50
60
70
4000
2m 3.6m
3000
2000
1000
-10 -20 Time (msec)
Fig. 14. Variation of shear stress in soil at different distances from pile shaft without considering soil damping.
0 0
10
20
30 40 Time (msec)
50
60
70
Fig. 15. Variation of normal stress at different distances below pile tip.
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418 Model prediction
Static load test
2000
Load (kN)
1800 1600 1400 1200 1000 800 600 400 200 0 0
10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 Pile head displacement (mm)
Fig. 16. Static load test simulation result compared with the real test result using the parameters obtained from signal matching analysis.
of the pile tip diameter, or in fact 30 mm, the ultimate pile capacity is scaled off as 1200 kN. This is indicating the reasonable accuracy of the numerical model in the static capacity prediction. The distribution of the shear stress on the pile shaft on the basis of the static load test simulation is shown in Fig. 17. The tip and shaft resistances are thus obtained as 650 kN and 550 kN, respectively. The trend of shear stress variations on the pile shaft compares reasonably well with the qc variations of the CPT profile of Fig. 7. 5. Discussion A successful so-called continuum numerical model is developed which proved to perform very well for the wave propagation and dynamic analysis of piles subjected to hammer blows. The developed model is superior to the conventional Smith-type models used in practice for highstrain test result interpretations, in that it more realistically considers the surrounding soil effects. Still there is a long
Shear stress (kPa)
Depth (m)
0
10
20
30
40
50
60
70
80
90
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
Fig. 17. Calculated distribution of shear stress on the pile shaft; compares reasonably well with the CPT profile.
417
way to go on this path to explore and manipulate its merits, however. One important factor that was not considered in this study is the existence of the ground water table (GWT) that has been reported at about 1 m below ground surface for the case pile studied. It is obvious, and of great importance, that the GWT effects should be considered in the analysis. The dominating layers at the pile shaft and at the pile tip are sandy, therefore the effective stress conditions are prevailing, even during the initial driving. But, why this was not considered in the current study, and how its effect has been compensated for, are briefly discussed below. GWT causes the effective stress to reduce to about half, and therefore the shaft resistance reduces to the same degree. According to the b method, the shaft resistance at any depth z is rs ¼ K tgd ðc0 ZÞ:
ð4Þ
It should be noted that in this fundamental equation, K is significantly affected by driving effort of the pile in sandy soils. Around the pile shaft and below the pile tip, the sand is usually densified to some distance from the pile. Densification around the pile shaft contributes to substantial increase in K. In the employed numerical model, however, the pile was embedded within the soil as if it is a bored pile. No attention has been paid on the driving effects yet, in the developed model. In fact, it is not so easy to implement the densification effect of the driving effort in the numerical model. To compensate for this lack of consideration of increase in K due to driving effort, the total stresses are considered in this analysis, meaning application of c instead of c’ (19 kN/m3 instead of 9 in our example case). This is more or less compensating for the increase in K due to driving effects. Otherwise, we would have required using d (interface friction angle) values very higher than the meaningful values, to match the predicted and measured Zv signals. Similar discussion is applicable to the tip resistance. It is reminded that the GWT results in effective stress reduction around the pile tip, and therefore reduction in the pile tip capacity. The reduction in tip capacity and stiffness would cause in major reduction in impedance in the soil zone below the pile tip, and therefore major changes in the F and Zv reflections at 2L/C. However, we had to substantially increase the elastic modulus of a soil block below the pile tip to match the predicted and measured signals. This is attributed to the fact that the sand mass density is significantly increased during the driving efforts of the pile, hence increasing the stiffness and strength parameters of the soil around and below the pile tip. Therefore, the balance between considering versus not considering the effective stress reduction due to GWT is adjusted by the stiffness and strength parameters of the block below the pile tip in the numerical model. Had we considered the effective stress conditions, we would have ended up in higher values for the stiffness and strength
418
S. Feizee Masouleh, K. Fakharian / Computers and Geotechnics 35 (2008) 406–418
parameters for the block to match the predicted and measured signals. This shortcoming has to be dealt with in future studies. However, as this study is one of the few attempts towards the application of a so-called continuum numerical model for the signal matching analyses, the effective stress analysis is not considered yet. The main focus has been on the successful operation of the model. It should be pointed out that also in the conventional mass–spring– dashpot models the effective stresses are not distinguished. In fact during the signal matching trial and error process, different parameters are varied such that the best match to be found, irrespective of the effective or total stress issues. 6. Conclusions The main objective of this numerical study was to develop a so-called continuum numerical model to simulate wave propagation and dynamic response of the pile driving problem and using the model for signal matching analyses. The model was successfully developed using a standard finite difference code. The findings of the study are summarized below: 1. If meaningful parameters are considered for the pile material, and the boundary conditions are introduced properly, the axisymmetric solid elements can adequately simulate the wave propagation along the shaft, and the reflections due to the impedance changes at the pile tip. 2. The inertial and/or radiation damping of the surrounding soil mass is automatically considered in this continuum numerical model. This is the most important contribution of the developed numerical model versus the mass–spring–dashpot type models. However, the degree of importance and outcome of such advantages still need to be manipulated and quantified in the future works. 3. The force and velocity signals were matched by varying the strength and deformation parameters of the soil and interface element. The model then proved to be very well in predicting other quantities such as Enthru, pile toe (tip) quake, static load–displacement curve and the load distribution on the pile shaft. 4. The effective stress type analysis, non-linear hardening stress–stress relations for soil, and sensitivity analysis of different parameters are the topics that require further studies to explore the capabilities of the new continuum numerical model.
References [1] Randolph MF. Science and Empiricism in Pile Foundation Design. 43rd Rankine Lecture. Geotechnique 2003;54(1). [2] Salgado R, Abou Jaoude G. Assessment of axially-loaded pile dynamic design methods and review of INDOT axially-loaded pile design procedure. Proposal for Research and Implementation Study, Purdu University, USA; 2003. [3] Smith EAL. Pile-driving analysis by the wave equation. J Soil Mech Found Div 1960;86(EM 4):35–61. [4] Coyle HM, Lowery Jr LL, Hirsch TJ. Wave equation analysis of piling behavior. In: Desai CS, Christian JT, editors. Numerical methods in geotechnical engineering. McGraw-Hill; 1977. [5] Fakharian K, Feizee S. Application of wave equation analysis (WEAP) for design of driven piles. In: Proceedings of 4th Int Conf on Coasts, Ports & Marine Structures (ICOPMAS 2000), Bandar Abbas, Iran; November 2000. [6] Feizee S. Application of wave equation analysis (WEAP) for design of driven piles. MSc Thesis, Department of Civil & Environmental Engineering, Amirkabir University of Technology, Tehran, Iran; 2000. [7] Rausche F. Soil response from dynamic analysis and measurements on piles. Thesis presented to the Case Western Reserve University, at Cleveland, Ohio, in 1970, in partial fulfillment of the requirements for the degree of Doctor of Philosophy; 1970. [8] Rausche F, Goble GG, Likins G. Dynamic determination of pile capacity. J Geotech Eng 1985;111(3):367–83. [9] Likins GE, Rausche F. Correlation of CAPWAP with static load tests. In: Proc 7th Conf Application of stress wave theory to piles. Petaling-Jaya, Selangor, Malaysia; 2004. pp. 153–165. [10] Coutinho ALGA, Costa AM, Alves JLD, Landau L, Ebecken NFF. Pile driving simulation and analysis by the finite element method. In: Proc 3rd Int Conf, Application of stress-wave theory to piles, Ottawa, Canada; 25–27 May 1988. [11] Borja RI. Dynamics of pile driving by the finite element method. Comput Geotech 1988;5(11):39–49. [12] Mabsout M, Tassoulas JL. A finite element model for the simulation of pile driving. Int J Numer Methods Eng 1994;37:257–78. [13] Mabsout ME, Reese LC, Tassoulas JL. Study of pile driving by finiteelement method. J Geotech Geoenviron Eng, ASCE 1995;121(7): 535–43. [14] Kolsky H. Stress waves in solids. New York: Dover Publications; 1963. [15] Geerling J, Smits MThJH. Prediction of load–displacements characteristics of piles from the results of dynamic/kinetic load tests. In: Proc 4th conf application of stress wave theory to piles. Delft University, Netherlands; 12–24 September 1992. [16] Heijnen WJ. General soil conditions of the test sites. In: Proc 4th conf application of stress wave theory to piles. Delft University, Netherlands; 12–24 September 1992. [17] Nath B. A continuum method of pile driving analysis: comparison with the wave equation method. Comput Geotech 1990;10:265–85. [18] Lysmer J, Kuhlemeyer RL. Finite dynamic model for infinite media. J Eng Mech 1969;95(EM4):859–87. [19] Bowles JE. Single piles: dynamic analysis, load tests, foundation analysis and design. McGraw-Hill; 1996. [20] Yang J. Influence zone for end bearing of piles in sand. J Geotech Geoenviron Eng 2006;32(9):1229–37. [21] Goble GG, Rausche F, Likins G. Associates Inc. Introduction to GRLWEAP, the 1997-2 GRLWEAP PROGRAM. Cleveland, Ohio; 1996.