Posted in Academic Issues, Geotechnical Engineering

Constitutive Elasticity Equations: Plane Stress

In our previous posts we discussed three-dimensional and two-dimensional plane strain/axisymmetric constitutive elasticity equations. In this post we will consider plane stress. This case doesn’t appear very much in geotechnical engineering; it’s more of interest to mechanical engineers and deals with problems such as this.

It’s probably a good idea to look again at the graphic for this:

Plane stress, plane strain and axial symmetry cases, from Owen and Hinton (1980)

The conditions of plane stress are thus:

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

If we compare these to plane strain, we find out that we switch the z-strain for the z-stress. This means that deriving the equations–especially the 3 x 3 case–is the mirror image of the plane strain case, and we will see this clearly below.

Let’s start by noting that, as before, the starting 4 x 4 constitutive matrix is the same as plane strain, thus

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

and its inverse

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

The stress and strain vectors, however, are

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

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

We see Equation (1a) applied to the stress vector to yield Equation (4). Multiplying through as we did before,

\sigma=\left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{{y}}\\ \noalign{\medskip}0\\ \noalign{\medskip}\tau_{{\it xy}} \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}_{xy} \end{array}\right] (6)

\epsilon=\left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}\epsilon_{z}\\ \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)}}\\ \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)}}\\ \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)}}\\ \noalign{\medskip}{\frac{\tau_{{\it xy}}}{G}} \end{array}\right] (7)

If we equation the two third row elements of the vectors in Equation (6), we have

\lambda\,\epsilon_{{x}}+\lambda\,\epsilon_{{y}}+{\frac{\left(1-\nu\right)\lambda\,\epsilon_{{z}}}{\nu}}=0 (8)

Solving for the z-axis strain,

\epsilon_{z}={\frac{\left(\epsilon_{{x}}+\epsilon_{{y}}\right)\nu}{-1+\nu}} (9)

Since the z-axis strain is a function of the x-axis and y-axis strains, the obvious question is this: can we go to a 3 x 3 matrix formulation? The answer is yes if we literally invert the process. We need to do the following:

  1. Consider the basic equation in this form: D_{e}^{-1}\sigma=\epsilon .
  2. Since \sigma_z = 0 , we can remove the third column and row from the constitutive matrix and the third row from the stress (and strain) matrix.
  3. We multiply through to yield

\left[\begin{array}{ccc} -{\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}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\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]=\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)}}\\ \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)}}\\ \noalign{\medskip}{\frac{\tau_{{\it xy}}}{G}} \end{array}\right] (10)

We then invert the inverse to obtain the 3 x 3 constitutive matrix. We use the forward formulation D\epsilon=\sigma and multiply through to obtain

\left[\begin{array}{ccc} {\frac{\lambda\,\left(-1+2\,\nu\right)}{\nu\,\left(-1+\nu\right)}} & {\frac{\lambda\,\left(-1+2\,\nu\right)}{-1+\nu}} & 0\\ \noalign{\medskip}{\frac{\lambda\,\left(-1+2\,\nu\right)}{-1+\nu}} & {\frac{\lambda\,\left(-1+2\,\nu\right)}{\nu\,\left(-1+\nu\right)}} & 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]=\left[\begin{array}{c} {\frac{\lambda\,\left(-1+2\,\nu\right)\epsilon_{{x}}}{\nu\,\left(-1+\nu\right)}}+{\frac{\lambda\,\left(-1+2\,\nu\right)\epsilon_{{y}}}{-1+\nu}}\\ \noalign{\medskip}{\frac{\lambda\,\left(-1+2\,\nu\right)\epsilon_{{x}}}{-1+\nu}}+{\frac{\lambda\,\left(-1+2\,\nu\right)\epsilon_{{y}}}{\nu\,\left(-1+\nu\right)}}\\ \noalign{\medskip}G{\it \gamma}_{{\it xy}} \end{array}\right] (11)

After these things we use Equation (9) to obtain the z-axis strain.

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.