Posted in Academic Issues, Geotechnical Engineering

Constitutive Elasticity Equations: Plane Strain and Axisymmetric Cases

In the last post we considered the three-dimensional elastic constitutive equations. Three-dimensional analyses are becoming more common in geotechnical engineering but many problems are considered and analysed using a two-dimensional paradigm. In this post we’ll look at plane strain/axisymmetric problems. The references are the same as the original post.

Let’s begin by defining what we mean by the terms “plane stress,” “plane strain” and “axisymmetric” by using and illustration.

Plane Stress, Plane Strain and Axisymmetric Cases (from Owen and Hinton (1980))

Let’s concentrate on the last two here. With plane strain we assume that

\epsilon_z = 0 (1a)
\tau_{xz} = 0 (1b)
\tau_{yz} = 0 (1c)

This reduces the stress and strain vectors to

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

\epsilon=\left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}0\\ \noalign{\medskip}{\it \gamma}_{{\it xy}} \end{array}\right] (3)

The constitutive matrix is

D_{e}=\left[\begin{array}{cccc} {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & 0\\ \noalign{\medskip}{\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & 0\\ \noalign{\medskip}{\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & 0\\ \noalign{\medskip}0 & 0 & 0 & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(2-2\,\nu\right)}} \end{array}\right] (4)

Applying the Lamé constants the way we did the last time,

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] (5)

Inverting this matrix,

D_{e}^{-1}=\left[\begin{array}{cccc} -{\frac{\nu}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & 0\\ \noalign{\medskip}{\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & -{\frac{\nu}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & 0\\ \noalign{\medskip}{\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & -{\frac{\nu}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & 0\\ \noalign{\medskip}0 & 0 & 0 & {G}^{-1} \end{array}\right] (6)

Again using the same method as before, we determine the stress and strain vectors thus:

\sigma = \left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{{y}}\\ \noalign{\medskip}\sigma_{{z}}\\ \noalign{\medskip}\tau_{{\it xy}} \end{array}\right]=\left[\begin{array}{c} {\frac{\left(1-\nu\right)\lambda\,\epsilon_{{x}}}{\nu}}+\lambda\,\epsilon_{{y}}\\ \noalign{\medskip}\lambda\,\epsilon_{{x}}+{\frac{\left(1-\nu\right)\lambda\,\epsilon_{{y}}}{\nu}}\\ \noalign{\medskip}\lambda\,\epsilon_{{x}}+\lambda\,\epsilon_{{y}}\\ \noalign{\medskip}G{\it \gamma}_{{\it xy}} \end{array}\right] (7)

\epsilon=\left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}0\\ \noalign{\medskip}{\it \gamma}_{{\it xy}} \end{array}\right]=\left[\begin{array}{c} -{\frac{\nu\,\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{y}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{z}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}-{\frac{\nu\,\sigma_{{y}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{z}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{y}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}-{\frac{\nu\,\sigma_{{z}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{\tau_{{\it xy}}}{G}} \end{array}\right] (8)

We can write the third row of Equation (8) as follows:

{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{y}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}-{\frac{\nu\,\sigma_{{z}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}=0 (9)

From which we can determine \sigma_z :

\sigma_{z}=\nu\,\left(\sigma_{{x}}+\sigma_{{y}}\right) (10)

This suggests that we can reduce the matrix from 4 x 4 to 3 x 3 and the vectors to 3-row. The elasticity equation D \epsilon = \sigma would look thus:

\left[\begin{array}{ccc} {\frac{\left(1-\nu\right)\lambda}{\nu}} & \lambda & 0\\ \noalign{\medskip}\lambda & {\frac{\left(1-\nu\right)\lambda}{\nu}} & 0\\ \noalign{\medskip}0 & 0 & G \end{array}\right]  \left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}{\it \gamma}_{{\it xy}} \end{array}\right] = \left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{{y}}\\ \noalign{\medskip}\tau_{{\it xy}} \end{array}\right] (11)

and inverting the matrix, the equation D^{-1}_e \sigma = \epsilon is

\left[\begin{array}{ccc} {\frac{\left(-1+\nu\right)\nu}{\lambda\,\left(-1+2\,\nu\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+2\,\nu\right)}} & 0\\ \noalign{\medskip}{\frac{{\nu}^{2}}{\lambda\,\left(-1+2\,\nu\right)}} & {\frac{\left(-1+\nu\right)\nu}{\lambda\,\left(-1+2\,\nu\right)}} & 0\\ \noalign{\medskip}0 & 0 & {G}^{-1} \end{array}\right]  \left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{{y}}\\ \noalign{\medskip}\tau_{{\it xy}} \end{array}\right]  = \left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}{\it \gamma}_{{\it xy}} \end{array}\right]  (12)

Multiplying Equations (11) and (12) through,

\sigma=\left[\begin{array}{c} {\frac{\left(1-\nu\right)\lambda\,\epsilon_{{x}}}{\nu}}+\lambda\,\epsilon_{{y}}\\ \noalign{\medskip}\lambda\,\epsilon_{{x}}+{\frac{\left(1-\nu\right)\lambda\,\epsilon_{{y}}}{\nu}}\\ \noalign{\medskip}G{\it \gamma}_{{\it xy}} \end{array}\right]  (13)

\epsilon= \left[\begin{array}{c} {\frac{\left(-1+\nu\right)\nu\,\sigma_{{x}}}{\lambda\,\left(-1+2\,\nu\right)}}+{\frac{{\nu}^{2}\sigma_{{y}}}{\lambda\,\left(-1+2\,\nu\right)}}\\ \noalign{\medskip}{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+2\,\nu\right)}}+{\frac{\left(-1+\nu\right)\nu\,\sigma_{{y}}}{\lambda\,\left(-1+2\,\nu\right)}}\\ \noalign{\medskip}{\frac{\tau_{{\it xy}}}{G}} \end{array}\right]  (14)

with Equation (10) rounding out for \sigma_z .

Axisymmetric Case

The axisymmetric case is achieved if we replace Equation (2) with

\sigma=\left[\begin{array}{c} \sigma_{{r}}\\ \noalign{\medskip}\sigma_{{z}}\\ \noalign{\medskip}\sigma_{{\theta}}\\ \noalign{\medskip}\tau_{{\it rz}} \end{array}\right] (15)

and Equation (3) with

\epsilon=\left[\begin{array}{c} \epsilon_{{r}}\\ \noalign{\medskip}\epsilon_{{z}}\\ \noalign{\medskip}\epsilon_{{\theta}}\\ \noalign{\medskip}{\it \gamma}_{{\it rz}} \end{array}\right] (16)

The subscript swap for polar coordinates is as follows:

  • x becomes r
  • y becomes z
  • z becomes θ
  • xy becomes rz

The math after that is identical, including the reduction for \sigma_{\theta} .

Note: after this was first posted, some errors were discovered. These were corrected with the assistance of Jaeger, J.C. and Cook, N.G.W. (1979) Fundamentals of Rock Mechanics. London: Chapman and Hall. This book has an excellent treatment of basic theoretical solid mechanics in general and elasticity in particular. Although they do not put the equations in matrix form, they do use the Lamé constants in their formulation, albeit in a different way than is done here. We apologise for any inconvenience or confusion.

Posted in Academic Issues, Geotechnical Engineering

Constitutive Elasticity Equations: Three-Dimensional Formulation

This is the beginning of what hopefully will be an enlightening series on elasticity and plasticity, primarily as it appears in finite element code. It is taken from a number of sources, including (but not limited to) Owen and Hinton (1980), Smith and Griffiths (1988) and Cook, Malkus and Plesha (1989.) Although finite element code has certainly advanced, the basics are still the same, especially with elasticity (plasticity is another matter.) More information on finite element analysis in geotechnics (with a link to Owen and Hinton) can be found here.

The basic elasticity relationship is familiar to engineers and can be stated as follows:

\sigma = E \epsilon (1)

where

  • \sigma = stress
  • E = Young’s modulus, or modulus of elasticity
  • \epsilon = strain
Element with normal and shear stresses (from Belyaev (1979))

It can be used in one-dimensional analyses such as tension and compression testing; we’ll come back to that later. But real problems are always in three dimensions in one way or another. Let’s start by considering an infinitesimal element (these are discussed in Verruijt) that looks like the one in Figure 69. There are many ways of looking at this element. For one thing, it is in static equilibrium, which means that the stresses can be changed along with the coordinate system and the results still be in static equilibrium. Related to this is the fact that there is one coordinate system at which point the shear stresses vanish and the three normal stresses that are left are the principal stresses; this is discussed here.

These, however, are not our main interest here. Our main interest is how the element deforms under stress. For now we will restrict our discussion to elastic, path-independent stresses. We will rewrite Equation (1) is a slightly different form:

D \epsilon = \sigma (2)

Now both \sigma and \epsilon are vectors; you can see the multiple stresses in Figure 69. We use the notation D because the modulus of elasticity E is a material property and thus is a part of the elasticity matrix D .

In any case, when the actual vector and matrices are substituted into Equation (2), for three-dimensional problems such as shown in Figure 69 the equations look like this:

\left[\begin{array}{cccccc} {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & 0 & 0 & 0\\ \noalign{\medskip}{\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & 0 & 0 & 0\\ \noalign{\medskip}{\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\nu}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(1-2\,\nu\right)}} & 0 & 0 & 0\\ \noalign{\medskip}0 & 0 & 0 & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(2-2\,\nu\right)}} & 0 & 0\\ \noalign{\medskip}0 & 0 & 0 & 0 & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(2-2\,\nu\right)}} & 0\\ \noalign{\medskip}0 & 0 & 0 & 0 & 0 & {\frac{{\it E}\,\left(1-\nu\right)}{\left(1+\nu\right)\left(2-2\,\nu\right)}} \end{array}\right]\left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}\epsilon_{{z}}\\ \noalign{\medskip}\gamma{}_{{\it xy}}\\ \noalign{\medskip}{\it \gamma}_{{\it yz}}\\ \noalign{\medskip}{\it \gamma}_{{\it zx}} \end{array}\right]=\left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{{y}}\\ \noalign{\medskip}\sigma_{{z}}\\ \noalign{\medskip}\tau_{{\it xy}}\\ \noalign{\medskip}\tau_{{\it yz}}\\ \noalign{\medskip}\tau_{{\it zx}} \end{array}\right] (3)

It should be immediately apparent that there are many repetitious quantities. There are many ways of dealing with the problem. For this analysis we will substitute Lamé’s constants, which are

\lambda = \frac{\nu E}{(1+\nu)(1-2\nu)} (4)

G = \frac{E}{2(1+\nu)} (5)

where \nu is Poisson’s Ratio and G is the shear modulus of elasticity. The latter is very useful in soil mechanics (which can be seen, for example, here) and in fact Equation (4) can be restated as follows:

\lambda = \frac{2\nu G}{1-2\nu} (6)

At this point we need to observe one important thing: no matter whether you use Equation (4) or (6), the Lamé constant \lambda blows up (or becomes singular, as we say in polite circles) when \nu = \frac{1}{2} . Materials where this is the case (such as very soft clays) are in reality fluids, the normal stresses become equal and the shear stresses only appear when the fluid moves. Avoiding this condition is one of the challenges of geotechnical finite element analysis.

It’s also worth noting that this is the same Lamé who came up with the theory behind shrink fits that is important with Vulcan hammers.

In any case, if we apply Equations (5) and (4) or (6), we obtain the following:

\left[\begin{array}{cccccc} {\frac{\left(1-\nu\right)\lambda}{\nu}} & \lambda & \lambda & 0 & 0 & 0\\ \noalign{\medskip}\lambda & {\frac{\left(1-\nu\right)\lambda}{\nu}} & \lambda & 0 & 0 & 0\\ \noalign{\medskip}\lambda & \lambda & {\frac{\left(1-\nu\right)\lambda}{\nu}} & 0 & 0 & 0\\ \noalign{\medskip}0 & 0 & 0 & G & 0 & 0\\ \noalign{\medskip}0 & 0 & 0 & 0 & G & 0\\ \noalign{\medskip}0 & 0 & 0 & 0 & 0 & G \end{array}\right]  \left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}\epsilon_{{z}}\\ \noalign{\medskip}\gamma{}_{{\it xy}}\\ \noalign{\medskip}{\it \gamma}_{{\it yz}}\\ \noalign{\medskip}{\it \gamma}_{{\it zx}} \end{array}\right]=\left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{{y}}\\ \noalign{\medskip}\sigma_{{z}}\\ \noalign{\medskip}\tau_{{\it xy}}\\ \noalign{\medskip}\tau_{{\it yz}}\\ \noalign{\medskip}\tau_{{\it zx}} \end{array}\right] (7)

We see from this that we have simplified assembling the constitutive matrix considerably. From a coding standpoint, doing this eliminates not only the possibility of coding error but also computational cost by eliminating the sheer number of operations. This is especially useful if these matrices are reconstituted during the analysis (as one would expect if strain-softening is applied, for example.)

If we multiply the left hand side through, we have

\left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{{y}}\\ \noalign{\medskip}\sigma_{{z}}\\ \noalign{\medskip}\tau_{{\it xy}}\\ \noalign{\medskip}\tau_{{\it yz}}\\ \noalign{\medskip}\tau_{{\it zx}} \end{array}\right] = \left[\begin{array}{c} {\frac{\left(1-\nu\right)\lambda\,\epsilon_{{x}}}{\nu}}+\lambda\,\epsilon_{{y}}+\lambda\,\epsilon_{{z}}\\ \noalign{\medskip}\lambda\,\epsilon_{{x}}+{\frac{\left(1-\nu\right)\lambda\,\epsilon_{{y}}}{\nu}}+\lambda\,\epsilon_{{z}}\\ \noalign{\medskip}\lambda\,\epsilon_{{x}}+\lambda\,\epsilon_{{y}}+{\frac{\left(1-\nu\right)\lambda\,\epsilon_{{z}}}{\nu}}\\ \noalign{\medskip}G{\it \gamma}_{{\it xy}}\\ \noalign{\medskip}G{\it \gamma}_{{\it yz}}\\ \noalign{\medskip}G{\it \gamma}_{{\it zx}} \end{array}\right] (8)

Inversely, if we want to know the strains, we would compute these by the relationship

\epsilon = D^{-1} \sigma (9)

Inverting,

D_{e}^{-1}=\left[\begin{array}{cccccc} -{\frac{\nu}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & 0 & 0 & 0\\ \noalign{\medskip}{\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & -{\frac{\nu}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & 0 & 0 & 0\\ \noalign{\medskip}{\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & {\frac{{\nu}^{2}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & -{\frac{\nu}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}} & 0 & 0 & 0\\ \noalign{\medskip}0 & 0 & 0 & {G}^{-1} & 0 & 0\\ \noalign{\medskip}0 & 0 & 0 & 0 & {G}^{-1} & 0\\ \noalign{\medskip}0 & 0 & 0 & 0 & 0 & {G}^{-1} \end{array}\right] (10)

and the strains can be computed as follows

\left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}\epsilon_{{z}}\\ \noalign{\medskip}\gamma{}_{{\it xy}}\\ \noalign{\medskip}{\it \gamma}_{{\it yz}}\\ \noalign{\medskip}{\it \gamma}_{{\it zx}} \end{array}\right] = \left[\begin{array}{c} -{\frac{\nu\,\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{y}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{z}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}-{\frac{\nu\,\sigma_{{y}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{z}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}+{\frac{{\nu}^{2}\sigma_{{y}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}-{\frac{\nu\,\sigma_{{z}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{\tau_{{\it xy}}}{G}}\\ \noalign{\medskip}{\frac{\tau_{{\it yz}}}{G}}\\ \noalign{\medskip}{\frac{\tau_{{\it zx}}}{G}} \end{array}\right]  (11)

Relationships such as Equations (8) and (11) will be more useful when we consider two-dimensional problems.

References

  • Cook, R.D., Malkus, D.S. and Plesha, M.E. (1989). Concepts and Applications of Finite Element Analysis, John Wiley & Sons, Inc., New York, NY
  • Owen, D.R.J., and Hinton, E. (1980). Finite Elements in Plasticity: Theory and Practice. Pineridge Press, Swansea, Wales.
  • Smith, I.M., and Griffiths, D.V. (1988). Programming the Finite Element Method. Second Edition. John Wiley & Sons, Chichester, England.

Note: after this was first posted, some errors were discovered. These were corrected with the assistance of Jaeger, J.C. and Cook, N.G.W. (1979) Fundamentals of Rock Mechanics. London: Chapman and Hall. This book has an excellent treatment of basic theoretical solid mechanics in general and elasticity in particular. Although they do not put the equations in matrix form, they do use the Lamé constants in their formulation, albeit in a different way than is done here. We apologise for any inconvenience or confusion.

Posted in Civil Engineering

A Quarter Century of vulcanhammer.net

It comes around every year, but this time it’s very special: today is the twenty-fifth anniversary of the start of this site. I wrote pieces of this for the tenth anniversary and the twentieth anniversary and one in the wake of COVID, which did as much as anything to move it in a more educational direction and not just a document download site.

You’d think that a site this old would have run its course, but it hasn’t: it continues to grow in traffic and popularity. The shift to a more educational emphasis has paid off, especially in traffic from outside of the U.S.. This is the fulfillment of a dream.

Want to support the site? Buy the books listed in the sidebar (or at the bottom, for those of you on mobile devices.) That’s all we ask. Again, as with my class videos, thanks for visiting and God bless.

Posted in Deep Foundations, Geotechnical Engineering

It’s Coming: Geotechnical Analysis One Grain at a Time

In going through some papers, I noticed this one, Friction anlysis of large diameter steel cylinder penetration process using 3D-DEM, whose abstract is as follows:

Large open-ended cylinder piles have been widely used for engineering foundation of port. The penetration process of the large-diameter steel cylinder exhibit complex behaviors, which is difficult to be measured by test and reproduced in numerical models. This study presents a friction analysis of large diameter steel penetration process by using the discrete element method (DEM), which can simulate large deformation and nonlinearity well. Centrifugal model and full-scale model were developed to analyze the sliding friction of the cylinder during installation and the contact force chain of soil particles. The validity of the DEM model was examined by comparing with theoretical values and published studies. Parametric studies were carried out to study the effects of contact parameters on side friction. Simulation results showed that, unlike pile penetration, there is no obvious soil-plug effect during the penetration process of large-diameter steel cylinder. Besides, the inside friction is smaller than the outside friction for large-diameter steel cylinder. What’s more, the computational cost of full-scale model based on the upscale theory was less than the centrifugal model. There is a close relationship between the side friction and micro contact parameters, which provides a reference for the follow-up study of cylinder or pile penetration using DEM.

Friction anlysis of large diameter steel cylinder penetration process using 3D-DEM

Put into simple terms, DEM models the soil by a grain-by-grain analysis of its response to load, which in this case came from 22 m O.D. open-ended pipe piles for the Hong Kong-Zhuhai-Macau Bridge. Believe it or not, I seriously considered using something like this for the STADYN project, but was dissuaded from doing so because of complexity and stability issues. That’s probably for the best, because the program I was in encouraged a “roll your own” approach to computer code, and to be honest at my stage in life and background I wasn’t up to such a development.

One thing that helped the authors of this paper was the fact that the soils they were driving into was cohesionless, without the chemical bonding that comes with cohesive soils. It’s also interesting to note that these piles were vibrated into place; the whole subject of vibratory driving, its performance prediction and static capacity determination, is another long-running subject in this business.

Although the issues of cohesion and verification (always the fun part of geotechincal engineering) need further resolution, DEM is, in my opinion, the ultimate solution of the soil interaction question, and needs to be disseminated further in our industry.

Posted in Academic Issues, Geotechnical Engineering

“Locating and quantifying necking in piles through numerical simulation of PIT” cites my paper “Comparison of Numerical Methods to Closed Form Solution for Wave Equation Analysis of Piling”

This recent paper, authored by Tarek Salem, Atef Eraky and Abdallah Almosallamy, was published in Frattura ed Integrità Strutturale. You can obtain the paper here but the gist of it can be found in the preprint, which can be found here.

My original paper was a derivative of my master’s thesis Closed Form Solution of the Wave Equation for Piles. The paper cited can be downloaded here.

The paper states that I used three methods of prediction for wave propagation in piles (WEAP, ANSYS and the closed form solution via MAPLE) to obtain the solution. Although all three of these and more cannot be applied to all solutions, the idea that we can use a variety of software to solve this problem was an important point of the original research, and I hope that this idea will continue to be pursued.

Note on the graphic at the top: it is one of the figures in my original paper, reconstructed in colour from the original CAD drawing. Many of the graphics in my thesis and my paper were done in CAD; we’ve come a long way in computer graphics.