Posted in Geotechnical Engineering, STADYN, TAMWAVE

Relating Hyperbolic and Elasto-Plastic Soil Stress-Strain Models

It is routine in soil mechanics to attempt to use continuum mechanics/theory of elasticity methods to analyse the stresses and strains/deflections in soil.  We always do this with the caveat that soils are really not linear in their response to stress, be that stress axial, shear or a combination of the two.  In the course of the STADYN project, that fact became apparent when attempting to establish the soil modulus of elasticity.  It is easy to find “typical” values of the modulus of elasticity; applying them to a given situation is another matter altogether.  In this post we will examine this problem from a more theoretical/mathematical side, but one that should vividly illustrate the pitfalls of establishing values of the modulus of elasticity for soils.

Although the non-linear response of soils can be modelled in a number of ways, probably the most accepted method of doing so is to use a hyperbolic model of soil response.  This is illustrated (with an elasto-plastic response superimposed in red) below.


The difficulties of relating the two curves is apparent.  The value E1 is referred to as the “tangent” or “small-strain” modulus of elasticity.  (In this diagram axial modulus is shown; similar curves can be constructed for shear modulus G as well.)  This is commonly used for geophysical methods and in seismic analyses.

As strain/deflection increases, the slope of the curve decreases continuously, and the tangent modulus of elasticity thus varies continuously with deflection.  For larger deflections we frequently resort to a “secant” modulus of elasticity, where we basically draw a line between the origin (usually) and whatever point of strain/deflection we are interested in.

Unfortunately, like its tangent counterpart, the secant modulus varies too.  The question now arises: what stress/strain point do we stop at to determine a secant modulus?  Probably a better question to ask is this: how do we construct an elasto-plastic curve that best fits the hyperbolic one?

One solution mentioned in the original study is that of Nath (1990), who used a hardening model instead of an elastic-purely plastic model.  The difference between the two is illustrated below.

Elasto-Plastic Response

Although this has some merit, the elastic-purely plastic model is well entrenched in the literature.  Moreover the asymptotic nature of the hyperbolic model makes such a correspondence “natural.”

Let us begin by making some changes in variables.  Referring to the first figure,



x = \epsilon

Let us also define a few ratios, thus:

A_1 = \frac {E_1}{\sigma_0}

A_2 = \frac {E_2}{\sigma_0}

A = \frac {E_2}{E_1}

Substituting these into the hyperbolic equation shown above, and doing some algebra, yields

y=\frac{x A_1}{1+x A_1}

One way of making the two models “close” to each other is to use a least-squares (2-norm) difference, or at least minimising the 1-norm difference.  To do the latter with equally spaced data points is essentially to minimise the difference (or equate if possible) the integrals of the two, which also equates the strain energy.  This is the approach we will take here.

It is easier to equate the areas between the two curves and the \sigma_0 line than to the x-axis.  To do this we need first to rewrite the previous equation as

y'=1-\frac{x A_1}{1+x A_1}

Integrating this with respect to x from 0 to some value x_1 yields

A_{hyp} = \frac{ln\left( 1+x_1 A_1 \right)}{A_1}

Turning to the elastic-plastic model, the area between this “curve” and the maximum stress is simply the triangle area above the elastic region.  Noting that

E_2 = \frac {\sigma_0} {x_2} ,

employing the dimensionless variables defined above and doing some additional algebra yields the area between the elastic line and the maximum stress, which is

A_{ep} = \frac {1}{2 A A_1}

Equation the two areas, we have

ln\left( 1+x_1 A_1 \right) - \frac {1}{2 A} = 0

With this equation, we have good news and bad news.

The good news is that we can (or at least think we can) solve explicitly for A , the ratio between the elastic modulus needed by elasto-plastic theory and the small-deflections modulus from the hyperbolic model.  The bad news is that we need to know A_1 , which is the ratio of the small deflections modulus to the limiting stress.  This implies that the limiting stress will be a factor in our ultimate result.  Even worse is that x_1 is an input variable, which means that the result will depend upon how far we go with the deflection.

This last point makes sense if we consider the two integrals.  The integral for the elasto-plastic model is bounded; that for the hyperbolic model is not because the stress predicted by the hyperbolic model is asymptotic to the limiting stress, i.e., it never reaches it.  This is a key difference between the two models and illustrates the limitations of both.

Some additional simplification of the equation is possible, however, if we make the substitution

x_1 = n x_2

In this case we make the maximum strain/deflection a multiple of the elastic limit strain/deformation of the elasto-plastic model.  Since

x_2 = \frac {sigma_0}{E_2} = \frac {1}{A_2} = \frac {1}{A A_1}

we can substitute to yield

ln\left( 1+\frac{n}{A} \right) - \frac{1}{2A} = 0

At this point we have a useful expression which is only a function of n and A .  The explicit solution to this is difficult; the easier way to do this is numerically.  In this case we skipped making an explicit derivative and use regula falsi to solve for the roots for various cases of n .  Although this method is slow, the computational time is really trivial, even for many different values of n .  The larger value of n , the more deflection we are expecting in the system.

The results of this survey are shown in the graph below.


The lowest values we obtained results for were about n = \frac{x_1}{x_2} = 0.75 .  When n = \frac{x_1}{x_2} = 1 , it is the case when the anticipated deflection is approximately equal to the “yield point.”  For this case the ratio between the elasto-plastic modulus and the small-strain hyperbolic modulus is approximately 0.4.  As one would expect, as n increases the elasto-plastic system becomes “softer” and the ratio A = \frac {E_2}{E_1} likewise decreases.  However, as the deflection increases this ratio’s increase is not as great.

To use an illustration, consider pile toe resistance in a typical wave equation analysis.  Consider a pile where the quake (x_2 ) is 0.1″.  Most “traditional” wave equation programs estimate the permanent set per blow to be the maximum movement of the pile toe less the quake.  In the case of 120 BPF–a typical refusal–the set is 0.1″, which when added to the quake yields a total deflection of 0.2″ of a value of n = 2.  This implies a value of A = 0.2139950 .  On the other hand, for 60 BPF, the permanent set is 0.2″, the total movement is 0.3″, and n = 3, which implies a value of A = 0.1713409.  Cutting the blow count in half again to 30 BPF yields n = 5 or A = 0.1383195.  Thus, during driving, not only does the plastic deformation increase, the effective stiffness of the toe likewise decreases as well.

Based on all this, we can draw the following conclusions:

  1. The ratio between the equivalent elasto-plastic modulus and the small-strain modulus decreases with increasing deflection, as we would expect.
  2. As deflections increase, the effect on the the equivalent modulus decreases.
  3. Any attempt to estimate the shear or elastic modulus of soils must take into consideration the amount of plastic deformation anticipated during loading.  Use of “typical” values must be tempered by the actual application in question; such values cannot be accepted blindly.
  4. The equivalence here is with hyperbolic soil models.  Although the hyperbolic soil model is probably the most accurate model currently in use, it is not universal with all soils.  Some soils exhibit a more definite “yield” point than others; this should be taken into consideration.
Posted in Deep Foundations, TAMWAVE

TAMWAVE: Pile Toe Resistance, and Some More on Pile Shaft Resistance

With this post we begin to discuss our “other” project: the TAMWAVE project.  It’s been around a long time but is now being revised.  The concept is to afford students a method of getting acquainted with several aspects of computer-aided driven pile design, including the following:

  • Estimating axial capacity of the pile;
  • Estimating the axial load-settlement of the pile;
  • Estimating the lateral load-settlement of the pile; and
  • Determining the drivability of the pile with a given hammer.

The current version of the online software is here.  One thing we’re doing is to designate the entire project as “TAMWAVE,” even though much of the routine isn’t really part of the wave equation program.

When most of the methods we use today were developed back around forty years ago and earlier, there wasn’t a really good way to distribute them away from mainframe computers.  The advent of DOS changed that, but with the shift towards Windows software most of these packages’ successors became proprietary.  Today we have DOSBOX to run these programs but current students, glued as they have been the last decade to their smartphones, find these hard to use. And, although geotechnical engineering isn’t the fastest moving branch of civil engineering, newer methods have been developed to analyse driven piles.

Overview of Toe Resistance

We’ve discussed in detail some newer methods of estimating the shaft resistance of driven piles, for sands and clays.  Although the original idea was to use them to enhance STADYN, they’re certainly applicable here, albeit with a few modifications.

With toe resistance, one of the advantages of 3D FEA code like STADYN is that it obviates (in theory at least) the need to estimate the toe resistance of the pile, let alone its progressive mobilisation.  That’s illustrated for drilled shafts in Han, Salgado, Prezzi and Lim (2016).  That’s not the case with a 1D routine like TAMWAVE, and so consider we must the toe resistance.  Han, Salgado, Prezzi and Zaheer (2016) (to whom we had recourse earlier) have a convenient listing of the toe methods that “go” with the shaft methods we discussed earlier, along with many others.

Starting with the toe resistance in sand, we have the following:

q_b = \left( 1-0.0058D_r \right)q_c

In this case q_b is the unit toe resistance of the pile, D_r is the relative density in percent, and q_c is the uncorrected cone resistance.  For toe resistance there are several schemes for averaging q_c around the toe, dating back to Schmertmann’s research, which is discussed in Fellenius.

For clays, the corresponding formula to Kolk and van der Velde (1996) is this:

q_b = 0.7 \left( q_t - \sigma_{vo\,toe}\right)

q_t is the corrected cone resistance at the toe; correction of q_c is also discussed in Fellenius\sigma_{vo\,toe} is the vertical total stress at the toe.

Soil Property Input and CPT Implementation

The original routine used the method of Dennis and Olson which really requires choosing whether the soil is cohesive or cohesionless and then answering some additional questions which are specific to the method.  The bifurcation of methods between the two soil types for driven piles is common but misleading; soils are seldom entirely one or the other but exhibit characteristics of both.  We plan to address this issue later for STADYN but for now we will stick with it for TAMWAVE.

One of TAMWAVE’s features which is carried over before is that there is only one soil type allowed for the entire length of the pile.  This is largely to preserve the academic nature of the software and discourage commercial use (which is prohibited anyway.)  That simplifies the writing of the code considerably but we must still choose how we should input the soil properties.

For the new version of TAMWAVE we opted to input the soil properties using two parameters.  The first is the two-letter unified code (SM, ML, etc.) for the characteristic soil type for the pile under consideration.  The second is the consistency or density of the soils, which is given using the verbal designations (“loose,” “hard,” etc.) which are customary in geotechnical engineering.  These are translated into actual properties using the “typical” correlations found in the Soils and Foundations Manual and are shown at the top of the page.  This isn’t a very exact method of proceeding but for the purpose of the routine it is adequate.

Use of these correlations gives us the following information:

  • Unit weight of the soil
  • Internal friction angle (cohesionless soils) or undrained shear strength/unconfined compression strength (cohesive soils)
  • SPT blow counts, generally corrected to N_{60} .

Conspicuously absent from this list are CPT results.  The general trend in pile capacity formulae in recent years is to correlate them to CPT results.  While the advantages of CPT testing are undeniable (and it’s certainly more consistent than SPT testing) the fact is that many of the soil borings that practitioners deal with feature SPT data, as do the typical values that TAMWAVE adopts.  Fortunately we have the correlations developed by Robertson and Campanella which relate the two.  Since the relationship between the two is based upon soil type, and we have that already, it is possible to automate the process and estimate equivalent CPT data from the typical SPT data we already have.  This relationship (and its limitations) is discussed in detail in Fellenius.

Randolph’s Lateral Earth Pressure Coefficient in Sand using CPT Data

In this post we discussed Randolph’s lateral earth pressure coefficient for sands.  The value for K_{max} can also be determined using CPT data as follows:

K_{max} = 0.02 \frac{q_c}{\sigma'_{vo}}

The rest of the formula is the same.


We have developed a new method of inputting soil data into this routine, along with outlining new methods of estimating the ultimate capacity of piles.  It is now necessary to implement these, which we will outline in a subsequent post.


  • Fei Han, Rodrigo Salgado, Monica Prezzi, Jeehee Lim. (2016) “Shaft and base resistance of non-displacement piles in sand.” Computers and Geotechnics, Volume 83, 2017, Pages 184-197, ISSN 0266-352X,
Posted in Deep Foundations, STADYN

Shaft Friction for Driven Piles in Clay: Alpha or Beta Methods?

In a previous post we discussed beta methods for driven pile shaft friction in sands, which are pretty much accepted, although (as always) the values for \beta can vary from one formulation to the next.  With clays, also as always, things are more complicated.

Since the researches of Tomlinson in the 1950’s, the shaft friction of piles in clays has been thought to be a function of the undrained shear strength of the clay multiplied by an adhesion factor \alpha , thus

f_s = \alpha c_u

This was seriously challenged by Burland (1973) who noted the following:

Whereas the use of undrained shear strength for calculating the end bearing capacity of a pile appears justified there seems little fundamental justification for relating shaft adhesion to undrained strength for the following reasons:

  1. the major shear distortion is confined to a relatively thin zone around the pile shaft (Cooke and Price (1973)).  Drainage either to or from this narrow zone will therefore take place rapidly during loading;

  2. the installation of a pile, whether driven or cast-in situ, inevitably must disturb and remould the ground adjacent to the pile shaft;

  3. quite apart from the disturbance caused by the pile there is no simple relationship between the undrained strength and drained strength of the ground.

Burland buttressed his case by noting that

\beta = K tan \phi

and presenting a graph similar to the following:

Beta Image 1

where, as seen earlier,

  • K_o = 1 - sin \phi is in red.
  • tan \phi is in blue.
  • \beta is in green.

Since, for the ranges of drained friction angles for clay (20-25 deg.) the value for \beta was relatively constant, value of \beta were relatively invariant with friction angle, and thus could be estimated with relative accuracy.  His empirical correlation was very successful with soft clays, not as much with stiff ones.

The year after Burland made his proposal, McClelland (1974) noted the following:

It is not surprising that there is a growing dissatisfaction with attempts to solve this problem through correlations of \alpha with c_u .  This is accompanied by a growing conviction that pile support in clay is frictional in character–that load transfer is dependent upon the effective lateral pressure acting against the side of the pile after it is driven.

However, \beta methods–which would embody McClelland’s preferred idea–have never been universally accepted for pile shaft friction in clays.  A large part of the problem, as noted by Randolph, Carter and Wroth (1979) is that the lateral pressure itself is dependent upon the undrained shear strength of the soils.  It is thus impossible to completely discount the effect of undrained shear strength on the shaft friction, even with the remoulding Burland and others have noted.

This has led to the “hybrid” approach of considering both undrained shear strength and effective stress.  This is embodied in the American Petroleum Institute (2002) specification.  A more advanced version of this is given in Kolk and van der Velde (1996).  They give the \alpha factor as

\alpha = 0.9\left( \frac {L-z} {d} \right)^{-0.2} \left( \frac {c_u} {\sigma'_{vo}} \right)^{-0.3} \leq 1

The notation is the same as in this post except that we add c_u , which is the undrained shear strength.

In this case the unit shaft friction is given by the equation

f_s = 0.9\left( \frac {L-z} {d} \right)^{-0.2} \left( \frac {c_u} {\sigma'_{vo}} \right)^{-0.3} c_u

There are a couple of things worth noting about this.

The first is that we can transform this into a \beta method of the form

f_s = \beta \sigma'_{vo}

with the following multiplication

f_s = 0.9\left( \frac {L-z} {d} \right)^{-0.2} \left( \frac {c_u} {\sigma'_{vo}} \right)^{0.7} \sigma'_{vo}

(A similar operation appears in Randolph (2005).)

in which case

\beta = 0.9\left( \frac {L-z} {d} \right)^{-0.2} \left( \frac {c_u} {\sigma'_{vo}} \right)^{0.7}

The only thing we would have to do is to find a way to incorporate the limiting condition for \alpha , which we will discuss shortly.

The second thing is that the term \left( \frac {L-z} {d} \right) appears in both this formulation and that for sands in this post.  The difference is that, while Kolk and van der Velde (1996) use the term in a power relationship, Randolph (2005) uses it in an exponential way.  The basic concept in both is the same: the term is at a maximum at the pile toe and decays toward the mudline.

The two are compared in the figure below.


Here the quantity \left( \frac {L-z} {d} \right) is at the x-axis and the following is at the y-axis:

  • Kolk and van der Velde Method for Clays, \left( \frac {L-z} {d} \right)^{-0.2} in red.
  • Randolph Method for Sands, e^{-\mu \left( \frac {L-z} {d} \right)} in blue, where \mu = 0.05 .
  • e^{-\mu \left( \frac {L-z} {d} \right)} in green, where \mu = 0.02 .

The graph illustrates the problem (from a computational standpoint) with the Kolk and van der Velde method: there is a singularity in their coefficient using the power relationship at the pile toe, while the exponential relationship yields a value of unity at this point.  The last correlation in green is approximately the best fit of the exponential relationship with the power relationship of Kolk and van der Velde, using either 1-norm or 2-norm methods.  It is not very good; it would be interesting, however, to see what kind of value for \mu might result if this had been in Kolk and van der Velde’s original statistical correlation equation.

In view of all this, perhaps the best way to enforce the limit is to do so as follows:

\left( \frac {L-z} {d} \right)\geq1

From all this, we can say that it is certainly possible to compute shaft friction for driven piles with a \beta method provided we include the effects of the undrained shear strength.


In addition to the original study and previous posts, the following references are noted:

Kolk, A.J., and van der Velde, A. (1996) “A Reliable Method to Determine Friction
Capacity of Piles Driven into Clays.” Proceedings of the 28th Offshore Technology Conference, Houston, TX, 6-9 May.  OTC 7993.

McClelland, B. (1974) “Design of Deep Penetration Piles for Ocean Structures.”  Journal of the Geotechnical Engineering Division, ASCE, Vol. 111, July.

Posted in Deep Foundations

Deep Foundations for Transportation Facilities: A Historic Perspective

This is a presentation slide show given by Mr. John G. Delphia, Texas Department of Transportation, Bridge Division, Geotechnical Branch Manager.  It’s a nice overview of deep foundations for transportation structures, including both drilled shafts (which TxDOT has excelled at since the days of O’Neill and Reese) and driven piles.

This slideshow requires JavaScript.

In addition to our own terms and conditions, please note the terms and conditions of the slide show, which are contained in the last slide and which we agree with.  We should also note that this slide show contains content from our companion site, especially from our pages on differential acting hammers, leaders and onshore hammers.

Posted in Uncategorized

An Important Announcement About

After a summer of eclipses and hurricanes, we’re pleased to announce that has finally moved to its new platform as of yesterday. Click here and check out what we have to offer. Most of the content has gone with the site; we’ve added many photographs and used the transition to correct many of the […]

via We’ve Moved At Last, But Took the Online Pile Routines Somewhere Else —