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,