Posted in Uncategorized

NAVFAC DM 7.1 Gets 2022 Refresh —

Nearly every geotechnical engineer has referred to the venerable 1986 NAVFAC DM 7.1 Manual on Soil Mechanics. There are numerous useful charts, correlations, tables, and equations for a variety of geotechnical challenges, which is probably […]

NAVFAC DM 7.1 Gets 2022 Refresh —

That didn’t take long…after our post The “Before and After” of NAVFAC DM 7, I was made aware of this. This has been posted to our Soil Mechanics page. As far as a print version is concerned, we are considering this, the old version continues to be available.

Posted in Uncategorized

Elasto-Perfect Plasticity with Mohr-Coulomb Failure

This is a continuation of the previous post on Mohr-Coulomb Failure. The references and nomenclature there apply to this post as well.

With the Mohr-Coulomb failure criterion defined, the elasto-plastic constitutive matrix can be developed. Since pure plasticity is assumed, thus once failure is reached and F=0 , the stresses can “rearrange” themselves to remain at the failure surface but cannot move from it unless the strains are reduced to the point where the stress state is within the failure surface and F<0. This is easier to see if we look at the diagram above. We should also note that, for this presentation, the plane strain/axisymmetric case will be considered; the derivation would proceed a little differently for other cases such as plane stress or three-dimensional stress.

To do this, both Equations 8 and 10 from the previous post are differentiated. Considering that

\sigma=\left[\begin{array}{c} \sigma_{x}\\ \sigma_{y}\\ \sigma_{z}\\ \tau_{xy} \end{array}\right] (1)

for the two-dimensional model, then

\frac{\partial F}{\partial\sigma}=a_{F}=C_{1}a_{1}+C_{2}a_{2}+C_{3}a_{3} (2)


C_{2}=cos\theta\left[\left(1+tan\theta tan3\theta\right)+\frac{sin\phi\left(tan3\theta-tan\theta\right)}{\sqrt{3}}\right] (3)
C_{3}=\frac{\sqrt{3}sin\theta+cos\theta sin\phi}{2J'_{2}cos3\theta}


a_{1}=\left[\begin{array}{c} 1\\ 1\\ 1\\ 0 \end{array}\right]
a_{2}=\frac{1}{2\sqrt{J'_{2}}}\left[\begin{array}{c} \sigma'_{x}\\ \sigma'_{y}\\ \sigma'_{z}\\ 2\tau_{xy} \end{array}\right] (4)
a_{3}=\left[\begin{array}{c} \sigma'_{y}\sigma'_{z}+\frac{J'_{2}}{3}\\ \sigma'_{x}\sigma'_{z}+\frac{J'_{2}}{3}\\ \sigma'_{y}\sigma'_{x}-\tau_{xy}^{2}+\frac{J'_{2}}{3}\\ -2\sigma'_{z}\tau_{xy} \end{array}\right]

Inspection of Equation 3 will show that C2 and C3 are undefined for values of θ =30°. This is the well known “corner” problem which has occupied the literature for many years, as summarized by Abbo (2011). For this study the “corner cutting” of Owen and Hinton (1980) is used, and

C_{2}=\frac{1}{2}\left[\sqrt{3}-\frac{sin\phi}{\sqrt{3}}\right],\:\theta\approx\pm30^{\circ} (5)

Although the accuracy of this “corner cutting” can be improved with methods such as those of Abbo (2011), they come at more computational expense.

In like fashion for the plastic potential function,

\frac{\partial Q}{\partial\sigma}=a_{Q}=C_{4}a_{1}+C_{5}a_{2}+C_{6}a_{3} (6)


C_{4}=\frac{sin\psi}{3}\\C_{5}=cos\theta\left[\left(1+tan\theta tan3\theta\right)+\frac{sin\psi\left(tan3\theta-tan\theta\right)}{\sqrt{3}}\right]\\C_{6}=\frac{\sqrt{3}sin\theta+cos\theta sin\psi}{2J'_{2}cos3\theta} (7)

and by extension

C_{4}=\frac{sin\psi}{3}\\C_{5}=\frac{1}{2}\left[\sqrt{3}-\frac{sin\psi}{\sqrt{3}}\right],\:\theta\approx\pm30^{\circ}\\C_{6}=0 (8)

and the vectors a1 , a2 , a3 are the same as with Equation 2.

At this point we need to pause and consider how we plan to present our result. It’s certain possible, for example, to multiply through Equations 2, 3, and 4 and use this result in subsequent derivations. This approach suffers from two problems: the element expressions in the matrices and vectors become very complicated, and the computation cost is increased by the repetitive calculations. For some good advice on this subject–for this and other problems–we would suggest that you review this.

Now the elasto-plastic constitutive matrix can be considered. Here the difference between path indepedence and path dependence becomes significant. In the elastic case, we can use the elastic consitutive matrix De and, knowing the strains, compute the stresses, and perform the inverse as well. With plasticity, it is necessary to compute the strains and stresses incrementally, moving from one stress and strain state to the next and adding the results to what has gone before. We do this with both elastic, plastic and combined strains and stresses for computational consistency. The elastic results will be the same except for the additional accumulation of numerical error due to the larger number of computations.

For any strain increment partially or totally beyond the yield surface, that strain increment will contain both elastic and plastic portions (Griffiths and Willson (1986)), or

d\epsilon=d\epsilon^{e}+d\epsilon^{p} (9)

Strain increments occur normal to the plastic potential surface, thus

d\epsilon^{p}=\lambda_{p}a_{Q} (10)

For elastic materials, the incremental stress-strain relationship is simply

d\epsilon^{p}=\lambda_{p}a_{Q} (11)

where, for plane strain and axisymmetric problems, as before,

D^{e}=\left[\begin{array}{cccc} {\frac{\left(1-\nu\right)\lambda}{\nu}} & \lambda & \lambda & 0\\ \noalign{\medskip}\lambda & {\frac{\left(1-\nu\right)\lambda}{\nu}} & \lambda & 0\\ \noalign{\medskip}\lambda & \lambda & {\frac{\left(1-\nu\right)\lambda}{\nu}} & 0\\ \noalign{\medskip}0 & 0 & 0 & G \end{array}\right] (12)

and also, for axisymmetric problems only,

\epsilon^{e}=\left[\begin{array}{c} \epsilon_{r}\\ \epsilon_{z}\\ \epsilon_{\theta}\\ \epsilon_{rz} \end{array}\right] (13)

For a perfectly elasto-plastic material, i.e., one without hardening or softening, stress changes will take place only during elastic action. Thus in these cases Equations 9, 10 and 11 can be rearranged to yield

d\sigma=D^{e}\left(d\epsilon-\lambda_{p}a_{Q}\right) (14)

Since stresses on the failure surface can “rearrange” themselves without leaving same surface,

a_{F}^{t}d\sigma=0 (15)

Combining and rearranging Equations 9 through 15 yields the following relationships:

\lambda_{p}=\frac{a_{F}^{t}D^{e}d\epsilon}{a_{F}^{t}D^{e}a_{Q}} (16)

d\sigma=\left(D^{e}-\frac{D^{e}a_{Q}a_{F}^{t}D^{e}}{a_{F}^{t}D^{e}a_{Q}}\right)d\epsilon (17)


d\sigma=\left(D^{e}-D^{p}\right)d\epsilon (18)


D^{p}=\frac{D^{e}a_{Q}a_{F}^{t}D^{e}}{a_{F}^{t}D^{e}a_{Q}} (19)

and finally

D^{ep}=D^{e}-D^{p} (20)

We need to break this down for greater understanding. Let us first define

a_{\beta}=\frac{1}{a_{F}^{t}D^{e}a_{Q}} (21)

This is a scalar quantity. (Admittedly sometimes it is difficult to follow matrix equations because they combine scalar and vector/matrix quantities.) We can thus rewrite Equation 19

D^{p}=a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)D^{e} (22)

and Equation 20 can be likewise rewritten as

D^{ep}=\left(I-a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)\right)D^{e} (23)

where I = the identity matrix.

If plasticity has not yet been initiated, the quantity a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right) (which is a 4 x 4 matrix) is null, and the elasticity matrix governs. Once it is initiated the magnitude of a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right) begins to increase and Dep begins to decrease, which indicates a ‘softening’ of the material. This is what we would expect with plasticity. The whole topic of matrix magnitude is not a straightforward one and is usually viewed in relationship to matrix norms, which are discussed in Gourdin and Boumharat (1989).

In any case

d\sigma=D^{ep}d\epsilon=\left(D^{e}-D^{p}\right)d\epsilon=\left(I-a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)\right)D^{e}d\epsilon (24)

The first thing to note about Equation 20 is that, if same is used to reconstruct a tangent stiffness matrix KT in a true Newton stepping scheme, same KT will not be symmetric if φ ≠ ψ. For many engineering materials, and especially those where both of these quantities are zero, this is not an issue; the flow rule is associated, Dep is symmetric and KT will be also. “Regular” engineering materials are used for, say, the steel or concrete portions of a system, but for many geotechincal applications interest in plastic deformation of these components is limited. It is also not an issue with purely cohesive soils. For situations with cohesionless soils, this is not the case; generally φ ≫ ψ and a non-associated flow rule is necessary. This aspect will be important in many decisions regarding the structure and types of schemes used in the model. An overview of the values of the dilitancy angle ψ can be found in Warrington (2019).

Turning to the inclusion of plasticity itself, it is certainly possible to compute the plastic stresses from Equation 20, and also possible to explicitly derive Equation 19 (Griffiths and Willson (1986)). However, it is not always optimal to do so, either from the standpoint of a workable algorithm or from a computational efficiency standpoint. To compute the final stress state in a load or time step where failure takes place, some type of iteration or multiple steps are required. Potts and Zdravkovic (1999) state that there are two basic types of algorithms to accomplish this: substepping (such as Sloan (1987)) and return (Ortiz and Simo (1986)). For STADYN a return algorithm was chosen, and implemented as follows:

  1. For a load or time step, the estimated incremental strain was computed.
  2. The estimated incremental stress was computed, based on the assumption that the strains were still in the elastic region.
  3. The resulting incremental stresses were added to the stresses at the beginning of the step. If these were elastic, then plasticity was not considered and the incremental stresses were added to the original ones. If they were not, then the plasticity routine was invoked.
  4. The plasticity routine began by computing aF , aQ , F , θ and aβ.
  5. The plasticity constant λp was determined by a modification to Equation 16, namely
    \lambda_p = F \alpha_\beta (25)
  6. The incremental strains in the return step were computed by the equation (see Equation 14)
    d\epsilon_{r}=\lambda a_{Q} (26)
  7. The return incremental stresses were thus
    d\sigma_{r}=D^{e}d\epsilon_{r} (27)
  8. Both of these were subtracted from the current estimated strain and stress, thus
    d\epsilon_{ep}^{n+1}=d\epsilon_{ep}^{n}-d\epsilon_{r}\\d\sigma_{ep}^{n+1}=d\sigma_{ep}^{n}-d\sigma_{r} (28)
  9. The current stress state was checked against the previous stress state. If the norm of the vector difference of the two stress states was within the convergence tolerance, the iteration was stopped and the computed elasto-plastic stresses and strains were accepted. If not, the cycle was repeated.

Potts and Zdravkovic (1999) criticize the return algorithms because they obtain a result based on information and computations in illegal stress space. However, overall the experience in this study is that the return algorithm worked well, with most of the convergence to the new failure surface stress state taking place in the first iteration and the rest refinement steps. A fair way of differentiating between the two is that, with substepping, one starts with the existing stress state and works one’s way to the failure surface, while the return algorithm starts by overshooting the failure surface and then coming back to it. The goal in both cases is to return to the failure surface, and in principle the result should be the same. Another way of differentiating between the two is that substepping is, with its avoidance of illegal stresses, more of an engineering type of approach to the problem, while the return algorithm is a more strictly mathematical method of arriving at a solution.

Posted in Uncategorized

Superposition, and Using Point Loads in Place of Distributed Ones

Elastic solutions to stresses induced by loads at the surface have their limitations, but they allow the use of the principle of superposition. The principle of superposition states that you can add the effects of different loads on a single point. Superposition requires that the stress state of the point be path independent, which is the case with elastic conditions. No matter how you load and unload a point in a system, if elasticity is maintained the result will be the same for a given set of loads.

This is illustrated for point loads in the graphic below, but it applies to distributed loads (such as this and this) as well.

Action of a Number of Concentrated Forces (from Tsytovich (1976))

The stresses that result from each load affect the total stress at the point of interest. They can be computed and added together. So, since “point” loads are physically impossible, is their computation be useful? The answer to this question is “yes” but it takes some judgement, like so many things in geotechnical engineering.

Let us consider the following case, from NAVFAC DM 7.01.

Separate Column Footings Problem, from NAVFAC DM 7.01

To solve this problem, DM 7 converted the square footings to circular ones and then used the chart shown in the post Going Around in Circles for Rigid and Flexible Foundations. This chart, like others, is hard to read. Is it possible to use point loads as a substitute?

There are three things we need to note here. The first is that the load on each column is 27 tons, and is the same for each column.

The second is that the “r” shown in the table above is not the same as it is for the point loads. The variable “r” for the point loads is the horizontal distance from the load to the point of interest.

The third is that there are three column positions shown in the diagram, with three corresponding values of r:

  • The column on top of the load, Column B2
  • The columns in the mid-point of the edges, Columns A2, B1, B3, and C2. These have an value of r of 15′ from the point of interest.
  • The columns in the corners of the square, Columns A1, A3, C1, and C3. These have a value of r of 21.2′ from the point of interest.

Now, instead of the chart, we apply the formula derived earlier for the influence coefficient, which is

K=3/2\,{\frac {1}{\pi \,\left (1+(\frac{r}{z})^{2}\right )^{5/2}}}

from which the stress is computed by the equation

\sigma_z = K \frac{P}{z^2}

Using this formula, we can construction the table below for this problem.

(1) Z, ft(2) r/Z for B2(3) r/Z for A2, B1, B3, C2(4) r/z for A1, A3, C1, C3(5) K for B2(6) K for A2, B1, B3, C2(7) K for A1, A3, C1, C3(8) Stress for B2, tsf(9) Stress for each of A2, B1, B3, C2, tsf(10) Stress for each of A1, A3, C1, C3, tst(11) Total Stress at Z, tsf
Results from the Separate Column Footings problem in NAVFC DM 7.01, using point loads. Note that the stresses in Columns 9 and 10 are due to a single column, and not all four of them. The total stress in column (11) is the sum of the stress in Column 8 plus 4 times the result in Column 9 plus four times the result in Column 10.

We could have opted to add the influence coefficients and then compute the stresses since, for each elevation Z, both the elevation and the column load were the same. We did not for clarity; it is certainly possible to have columns of different loads.

The results are conservative, they tend to be higher, especially at the lower value elevations. It’s worth noting that the total stress at Z = 2′ is higher than the distributed loads on the footings. One way to make this better is to use the center formulae for stress under circles for Column B2 and point loads for the rest. Given, however, the limitations of the method in general, and the considerably lower effort in obtaining the influence coefficients, the method is a reasonable one to use.

Superposition doesn’t only apply to point loads; the example given in the post Analytical Boussinesq Solutions for Strip, Square and Rectangular Loads uses it for rectangular and square loads. Nevertheless, with appropriate engineering judgement, using point loads in place of distributed loads can be a viable option.

Posted in Uncategorized

Estimating Load-Deflection Characteristics for the Shaft Resistance of Piles Using Hyperbolic Strain Softening

The latest research from this site is an offshoot of the STADYN project. The abstract is as follows:

Hyperbolic soil modelling and strain softening have become better understood in recent years; however, their application to modelling load-deflection characteristics for the shaft of both bored and driven piles is in the early stages of development. This paper proposes a method that is based strictly on the strain-softening changes in shear modulus that is experienced to varying degrees around the pile shaft. A dimensionless method is proposed which can be transformed to a physical estimate using simple soil parameters.

You can download the paper here or find it on ResearchGate.

Posted in Uncategorized

Derivation of Flamant’s Equations for the Elastic Response Under a Line Load

In the last post I showed the derivation of the elastic response of a soil to a point load using Boussinesq theory. This is a common part of elementary Soil Mechanics courses but it is uncommon to see the derivation.

The same situation exists with a line load, the elastic solution usually attributed to Flamant. The solution itself also presented in NAVFAC DM 7.01 and Arnold Verruijt’s Soil Mechanics. The derivation is presented below, and comes from the Manual of the Theory of Elasticity, by V.G. Rekach. Some discussion of similar solutions is in my post Analytical Boussinesq Solutions for Strip, Square and Rectangular Loads and my class presentation is at Soil Mechanics: Elastic Solutions to Soil Deflections and Stresses.