Geotechnical Engineering

My Thoughts on “DFI’s Rumble: Numerical Modeler vs Skeptic”

Posted on 13 July 2022 by Don Warrington

I found myself interested in the “rumble” between Augusto Lucarelli and Dr. Brian Anderson of Auburn University (not so far from where I teach) about numerical modelling, and the pros and cons of same. For this website, the timing is interesting, coming in the middle of my series on constitutive modelling (there’s more to come, Lord willing.)

You’d think that someone who a) has his terminal degree in Computational Engineering and b) teaches the geotechnical component of the civil engineering undergraduate program at the University of Tennessee at Chattanooga would be an enthusiastic proponent of the use of numerical models for the work. And they’re certainly valuable, but let me just lay out some observations about them in light of the discussion above:

At the start mention was made of using things like SPT tests and Atterberg limits for various determinations. The truth of the matter is that use of lab and field data alike is frequently empirical in its application. If the empirical correlation is not applicable to a certain situation for whatever reason, or the soil data has problems, going from “closed form” (more about that later) solutions to numerical models isn’t going to fix problems of this kind, although it may be better at hiding their deficiencies.
A distinction needs to be made between many of the models we developed in the past and those which are more commonly characterised as “numerical models” now. The USACOE program CBEAR was mentioned (although programs such as CWALSHT, MAGSETII, CSETT, CSANDSET and really Pile Buck’s SPW911 are similar in nature.) These are really computational aids for “closed form” solutions or similar methods. Today a numerical solution is generally finite element or finite difference, such as the COM624/LPILE family, WEAP/GRLWEAP or its inverse CAPWAP, in addition to the numerous finite element codes in use. The former are really automations of things engineers used to do completely by hand; the latter are beyond hand calculation. The advantage of the former is that, although it’s certainly possible to mess up the results with bad input, the possibility of computational error in the process is reduced. With the latter direct replication of the results with hand calculations or closed form solutions is impossible except for very simple cases, and it is here that uncritical acceptance of results is especially dangerous.
It is certainly possible, and desirable during model development, to make a reasonable comparison between a “closed form” solution and one generated by a numerical model provided that the case furnished to the latter is simple enough for reasonable comparison. There are two problems here. The first is that many “modellers” do not make the effort to develop a model in this way by first starting with a simple case and then proceeding to a more complex one. The second is that many of the “closed form” solutions we have in geotechnical engineering are in reality empirical correlations with variable validity. An example of this are the static pile capacity equations, for which there is a proliferation of solutions. This indicates that the science is certainly “not settled,” and in any case runs into the strength vs. service problem (more about that shortly.) The program director for my PhD program said that a good analytical solution was better than a modelled one all other things equal (an example of such a comparison in a simpler case is here.) But finding such a “good analytical solution” in geotechincal engineering isn’t always easy.
I think that the way we’ve taught civil engineering on the undergraduate level does not prepare students for a meaningful understanding of the way finite element and finite difference programs actually work. In some ways it’s more serious than that: the geotechnical component is the first place where students have to face two-dimensional stress distribution in a meaningful way, which has been a major challenge for me teaching these courses over the years. In this respect the way it’s done across the pond (you can see this in texts such as Verruijt and Tsytovich) is better, which may explain some of Dr. Lucarelli’s optimism about this topic. Compounding this problem is the fact that students aren’t required to write code as they used to (Python may change that or may not) which increases the “black box” aspect of packaged models. (In recent years, I find some of my students struggling with engineering spreadsheets…)
Traditionally we’ve divided analytical solutions in geotechnical engineering into three parts: those based on elasticity, those based on plasticity theory and those based on consolidation theory. The advantage of models such as finite element ones is that they can handle more than one of these at a time. The disadvantage is that how these are handled depends upon the model used by the program, something the practicing engineer may or may not understand. For example, those programs which use Mohr-Coulomb theory implemented elasto-plastically use what is admittedly a crude model, but one which fits with a wider range of soils and is better integrated with our current testing scheme. A model used with parameters which have to be “massaged” to fit it is a model asking for trouble.
The entrance of plasticity into models–unavoidable in geotechinical engineering–brings up the whole issue of path dependence, which in turn brings up uniqueness issues for the solution. These are unavoidable; attempts by engineers to produce the single “right” answer are doomed to failure. The happiest result we can hope for is the “best” solution. This pushes things toward a more probabilistic way. Since the advent of LRFD and other statistical methods there is a greater appreciation of this among engineers which may not be shared by clients or those in the legal profession.
The topic of the overemphasis on strength loading (and the corresponding de-emphasis on service loading) is a critical one in geotechnical engineering. Although certain problems (such as slope stability and bearing capacity failure) are catastrophic in nature as they push the soil to the upper bound failure state, most geotechnical problems are with settlement. This is especially true with deep foundations, where service load failure is the greater danger. I know that Bengt Fellenius has beat on this topic for a long time and I attempt to emphasise this with my students as well.
The tendency to throw the solution into the hands of “big data” or “machine learning” without using these tools to sharpen our predictive models is a scientific disaster waiting to happen, and that’s not only true of geotechnical engineering.

I am glad that DFI set this dialogue up, it was informative and timely. I am also glad the Equipment Corporation of America, who was Vulcan Iron Works‘ mid-Atlantic dealer for many years, has sponsored this effort.

Constitutive Elasticity Equations: Uniaxial Cases

Posted on 11 July 202220 October 2022 by Don Warrington

Our last cases to consider are the cases of uniaxial loading. Although these cases may seem trivial, they are important. As it was with the two-dimensional cases of plane stress and plane strain, there are two variations of the uniaxial loading to consider: uniaxial stress and uniaxial strain.

Uniaxial Stress

We start by considering the constitutive equation of the three-dimensional case, which is

$\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]$ (1)

For uniaxial stress, the following parameters are zero:

$\sigma_y = \sigma_z = 0$ (2a)
$\tau_{xy} = \tau_{xz} = \tau_{yz} = 0$ (2b)

As a result, stress and strain vectors are

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

It’s worth noting that we’ve zeroed out the shear stresses and strains. That’s because, with a uniaxial loading, the x, y and z stresses and strains are guaranteed to be principal stresses and strains, and thus the shear stresses are zero. We will see that the shear stresses and strains will return with axis rotation.

In any case, doing the usual multiplications and inverse multiplications, we have the following:

$\epsilon=\left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}\epsilon_{{y}}\\ \noalign{\medskip}\epsilon_{{z}}\\ \noalign{\medskip}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0 \end{array}\right]=\left[\begin{array}{c} -{\frac{\nu\,\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}{\frac{{\nu}^{2}\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}\\ \noalign{\medskip}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0 \end{array}\right]$ (5)

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

Considering the first rows of Equation (5) gives us

$\epsilon_{{x}}=-{\frac{\nu\,\sigma_{{x}}}{\lambda\,\left(-1+\nu+2\,{\nu}^{2}\right)}}$ (7)

and from the next two rows

$\epsilon_{y}=\epsilon_{z}=-\nu\epsilon_{x}$ (8)

It’s possible to formulate the computation of the stresses based on the strains, but in uniaxial loading it is usually easier to know the stresses first.

Although uniaxial loading is common in tensile and compressive stresses of materials such as metals and plastics, in geotechnical engineering it is mainly concentrated in two types of tests: Triaxial tests and Unconfined Compression (UC) tests. Since in Mohr-Coulomb theory we use a maximum shear stress to determine the point of failure, and the shear stresses and strains don’t appear in the above equations, an explanation as to how shear stresses and strains apply is reasonable.

Rather than go into a math-intensive derivation of the problem, recourse to Mohr’s Circle is more instructive.

Various conditions for triaxial and UC tests (from NAVFAC DM 7.01)

Let’s look at the top diagram from a UC test. The main uniaxial stress $\sigma_1 = \sigma_x$ . If we go 90 degrees anti-clockwise around Mohr’s Circle from $\sigma_1$ , we arrive at the point of maximum shear. This means that the plane of maximum shear is 45 degrees (half Mohr’s Circle angle) from the centre axis of the test specimen. Thus we would expect failure along that plane for a soil that follows Mohr-Coulomb elastic-perfectly plastic behaviour.

In the case of materials with internal friction, things are a little more complicated, because the angle from the principal stress to the point where Mohr’s Circle meets the failure envelope is greater than 90 degrees, thus the angle between the centre axis and the failure plane is greater than 45 degrees, and depends upon the interaction of the cohesion and the internal friction angle.

Of course people who actually run these tests will tell you that these specimens don’t always fail along these planes, or simply bulge at failure. That’s because soils don’t always obey Mohr-Coulomb elastic-perfectly plastic failure theory, which is one of those things that makes geotechnical engineering challenging.

Uniaxial Strain

Although this case doesn’t seem at first glance to be very useful, it’s really an important one for geotechnical engineers. We’ll start again with Equation (1) and zero out the following parameters:

$\epsilon_y = \epsilon_z = 0$ (9)

Of course Equation (2b) applies as well. Doing this results in the following vectors, and I’ll go ahead and multiply and inversely multiply through as before:

$\epsilon=\left[\begin{array}{c} \epsilon_{{x}}\\ \noalign{\medskip}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0 \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}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0 \end{array}\right]$ (10)

$\sigma=\left[\begin{array}{c} \sigma_{{x}}\\ \noalign{\medskip}\sigma_{y}\\ \noalign{\medskip}\sigma_{z}\\ \noalign{\medskip}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0 \end{array}\right]=\left[\begin{array}{c} {\frac{\left(1-\nu\right)\lambda\,\epsilon_{{x}}}{\nu}}\\ \noalign{\medskip}\lambda\,\epsilon_{{x}}\\ \noalign{\medskip}\lambda\,\epsilon_{{x}}\\ \noalign{\medskip}0\\ \noalign{\medskip}0\\ \noalign{\medskip}0 \end{array}\right]$ (11)

Extracting the important results,

$\sigma_{{x}}={\frac{\left(1-\nu\right)\lambda\,\epsilon_{{x}}}{\nu}}$ (12)
$\epsilon_{x}={\frac{\sigma_{{x}}\nu}{\left(1-\nu\right)\lambda}}$ (13)
$\sigma_{y}=\sigma_{z}=\frac{\nu}{1-\nu}\sigma_{x}$ (14)

Although it may seem esoteric, this condition is actually important because it simulates confined specimens under uniaxial load. The best known of these, of course, is consolidation testing and all of the variations of that which are employed. It may seem strange to employ elastic theory in this situation, but that speaks to Jean-Louis Briaud’s pet peeve regarding consolidation testing and the partial solution to the problem he poses. Some discussion of this relationship can also be found in Verruijt.

The whole concept of confined specimens is based on the fact that the entire ground upon which we build is a confined specimen, confined by the semi-infinite soil surrounding the loaded portion. Thus this unlikely looking case is very relevant for real geotechnical engineering.

Academic Issues, Geotechnical Engineering

Constitutive Elasticity Equations: Plane Stress

Posted on 9 July 20229 July 2022 by Don Warrington

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:

Consider the basic equation in this form: $D_{e}^{-1}\sigma=\epsilon$ .
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.
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.

Academic Issues, Geotechnical Engineering

Constitutive Elasticity Equations: Plane Strain and Axisymmetric Cases

Posted on 6 July 20229 July 2022 by Don Warrington

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,

Inverting this matrix,

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.

Academic Issues, Geotechnical Engineering

Constitutive Elasticity Equations: Three-Dimensional Formulation

Posted on 5 July 20226 July 2022 by Don Warrington

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.