## Determining the Degree of Consolidation

This is the last (hopefully) post in a series on consolidation settlement. We need to start by a brief summary of what has gone before. Note: the material for this derivation and those that preceded it have come from Tsytovich with some assistance from Verruijt.

## Review

In the post From Elasticity to Consolidation Settlement: Resolving the Issue of Jean-Louis Briaud’s “Pet Peeve”, we discussed the issue of how much soils (especially cohesive ones) settle through the rearrangement of particles. We were able to start with the theory of elasticity and, considering the effects of lateral confinement, define the coefficient of volume compression $m_v$ by

$m_v = \frac{\beta}{E}$ (1)

where E is the modulus of elasticity and $\beta$ is a factor based on Poisson’s Ratio and includes the effects of confinement, be that in an odeometer or in a semi-infinite soil mass. We also showed that, for a homogeneous layer,

$\delta_p = m_v H_o \sigma_x$ (2)

where $\delta_p$ is the settlement of the layer, $H_o$ is the thickness of the layer and $\sigma_x$ is the uniaxial stress on the layer. The problem is that $m_v$ is not constant, and the settlement more accurately obeys the law

$\delta_p = \frac{C_c H_o}{1+e_o} \log{\frac{\Delta p + \sigma_o}{\sigma_o}}$ (3)

where $C_c$ is the compression index, $e_o$ is the initial void ratio of the layer, $\Delta p$ is the change in pressure induced from the surface, and $\sigma_o$ is the average effective stress in the layer.

Turning to the post Deriving and Solving the Equations of Consolidation, we first determined that the change in porosity $\Delta n$ could, for small deflections, be equated to the change in strain $\epsilon$. From this we could say that

$\Delta n = m_v \Delta \sigma_x$ (4)

The change in porosity, for a saturated soil whose voids are filled with an incompressible fluid (hopefully water) induces water flow,

${\frac {\partial }{\partial x}}q(x,t)=-{\frac {\partial }{\partial t}} {\it n}(x,t)$ (5)

where $q(x,t)$ is the flow of water out of the pores and $n(x,t)$ is the porosity as a function of position and time. The flow of water is regulated by the overall permeability of the soil, and all of this can be combined to yield

${\frac {k{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)}{{\it \gamma_w }}}=m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (6)

where $k$ is the permeability of the soil and $\gamma_w$ is the unit weight of water. Defining

$c_v = \frac{k}{m_v \gamma_w}$ (7)

and making some assumptions about the physics, we can determine the equation for consolidation as

$c_{{v}}{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)={\frac {\partial }{\partial t}}u(x,t)$ (8)

where \$latex u(x,t) is the pore water pressure. If we invoke the effective stress equation and solve this for the boundary and initial conditions described, we have a solution

$\sigma_{x}(x,t)=p\left(1-\frac{4}{\pi}\left(\sin(1/2\,{\frac{\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/3\,\sin(3/2\,{\frac{\pi\,x}{h}}){e^{-9/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/5\,\sin(5/2\,{\frac{\pi\,x}{h}}){e^{-{\frac{25}{4}}\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}\cdots\right)\right)$ (9)

## The Degree of Consolidation

One thing that our theory presentation demonstrated was the interrelationship between pore pressure, stress and deflection. We know what the ultimate deflection will be based on Equation (3) above (or more complicated equations when preconsolidation is taken into consideration.) But how does the settlement progress in time?

We start by defining the degree of consolidation thus:

$U = \frac{\delta(t)}{\delta_p}$ (10)

where $\delta(t)$ is the settlement at any time before complete settlement. For the specific case (governing equations, initial equations and boundary conditions) at hand, the degree of consolidation–the ratio of settlement at a given point in time to total settlement–can be determined as follows:

$U_{o}=\intop_{0}^{h}\frac{\sigma_{x}(x,t)}{ph}dx$ (11)

In this case the result is divided by the uniform pressure p and the height h. Let us further define the dimensionless time constant

$T_{v}=\frac{c_{v}t}{h^{2}}$ (12)

That being the case, if we integration Equation (9) with Equation (11), we obtain

$U_{o}=1-\sum_{n=1}^{\infty}4\,{\frac {{e^{-1/4\,{\it Tv}\,{n}^{2}{\pi }^{2}}}\left (\cos(n\pi )\cos(1/2\,n\pi )-\cos(n\pi )-\cos(1/2\,n\pi )+1\right )}{{n}^{2}{\pi }^{2}}}$ (13)

otherwise put

$U_{o}=1-8\,{\frac {{e^{-1/4\,{\it Tv}\,{\pi }^{2}}}}{{\pi }^{2}}}-{\frac {8}{9}}\,{\frac {{e^{-9/4\,{\it Tv}\,{\pi }^{2}}}}{{\pi }^{2}}}-{\frac {8}{25}}\,{e^{-{\frac {25}{4}}\,{\it Tv}\,{\pi }^{2}}}{\pi }^{-2}\cdots$ (14)

As was the case with Equation (9), only the odd values of n are considered; the even ones result in zero terms.

It is regrettable that, in defining $T_v$, the value $\frac{\pi^2}{4}$ was not included, as using Equation (14) would be much simpler. For certain cases, it is possible to use the first two or three terms. In any case the usual method for determining $T_v$–and by extension the degree of consolidation–is generally done either using a graph or a table, as is shown in the graph at the start of the post (repeated below:)

The notation is a little different. We use the variable $U_o$ to emphasise that we are dealing with the “standard” case. The above graph also gives approximating equations; it is easy to see that, for $T_v > 0.2$, the equation given is simply the first two terms of Equation (14). The distinction between the drainage length h ($H_{dr}$ in the graph above) and the layer thickness H is clear.

## Conclusion

We have covered the basic, classic case of consolidation settlement in this post and its predecessors From Elasticity to Consolidation Settlement: Resolving the Issue of Jean-Louis Briaud’s “Pet Peeve” and Deriving and Solving the Equations of Consolidation. We trust that this presentation has been enlightening and informative.

## Deriving and Solving the Equations of Consolidation

In an earlier post we discussed consolidation settlement. For situations where soil a) decreases its volume due to rearrangement of the particles and b) does so over a relatively long period of time due to difficulties in expelling pore water pressures, we need to know how long it takes to reach maximum settlement, in addition to know what that settlement is and what deflections we might achieve along the way.

This is generally a two part process: a) determining the dissipation of pore water pressure and b) determining the amount of settlement at a given time. This post will focus on the first part of the process; we will deal with the second in a later post. This is a well-worn path in geotechnical engineering but hopefully this derivation is a little simpler than others.

## Developing the Differential Equation

### Relating Porosity to Strain

Let us begin by considering the system above (from Tsytovich (1976)), a uniformly loaded soil layer which is saturated (always) and clay (usually.) The first thing we need to note is that, while the diagram uses z as the variable of length in the vertical direction, from this post we will use the variable x (sorry for the confusion.) The second thing is that we assume the water and the solids to be incompressible. The third thing is that all the changes that take place do so because of changes in the voids where the water is resident. We first note the definition of porosity as

$n={\frac {V_{{v}}}{V_{{t}}}}$ (1)

where

• n = porosity
• Vv = volume of voids
• Vt = total volume

That being the case, the relationship between the porosity of the soil (due to changes in the volume of the voids) and the flow rate of the water can be expressed as

${\frac {\partial }{\partial x}}q(x,t)=-{\frac {\partial }{\partial t}} {\it n}(x,t)$ (2)

We can envision a differential volume having a height x and an area A. Since the problem is one-dimensional, the areas cancel out and the ratio of the void height to the total height is

$n={\frac {x_{{v}}}{x_{{t}}}}$ (3)

where

• xv = height of the voids
• xt = height of the solids

We want to determine the change in porosity from some state 0 to some state 1, just as we did with void ratio in this post. That change can be expressed as follows:

${\frac {x_{{{\it v0}}}}{x_{{{\it t0}}}}}-{\frac {x_{{{\it v1}}}}{x_{{{\it t0}}}-x_{{{\it v0}}}+x_{{{\it v1}}}}}$ (4)

We can assume that the change in void volume xv0 – xv1 << xt0, in which case Equation (4) can be simplified to

$\Delta{{{\it n}}}={\frac {x_{{{\it v0}}}-x_{{{\it v1}}}}{x_{{{\it t0}}}}} = \frac{\Delta x}{x_t}$ (5)

Now we can say, for the small increments we are dealing with here, that the change in strain is equal to the change in porosity,

$\Delta n = \Delta \epsilon$ (6)

From this and our previous post, we can thus use the change in strain to come to the following:

$\Delta n = m_v \Delta \sigma_x$ (7)

Stating this differentially,

${\frac {\partial }{\partial t}}{\it n}(x,t)=m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (8)

In this way we emphasise that both porosity and uniaxial stress are functions of depth and time. We now combine Equations (2) and (8) to yield

${\frac {\partial }{\partial x}}q(x,t)=-m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (9)

### Including Permeability

Darcy’s Law (or more properly d’Arcy’s Law) states that

$q(x,t) = -k{\frac {\partial }{\partial x}}H(x,t)$ (10)

where k is the coefficient of permeability and H is the hydraulic head. We then combine Equations (9) and (10) to obtain

$k{\frac {\partial ^{2}}{\partial {x}^{2}}}H(x,t)=m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (11)

Noting that

$H(x,t) = \frac{u(x,t)}{\gamma_w}$ (12)

where $\gamma_w$ is the unit weight of water and $u(x,t)$ is the excess pore water pressure generated by the decrease in porosity, we substitute this with the result of

${\frac {k{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)}{{\it \gamma_w }}}=m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (13)

### Putting It All Together

We now define the parameter

$c_v = \frac{k}{m_v \gamma_w}$ (14)

At this point we need to stop and make an observation: this whole process wouldn’t be worth too much if $c_v$ wasn’t constant (or reasonably so.) The reason it is is that the coefficient of permeability $k$ and the coefficient of volume compressibility $m_v$ both decrease as the void ratio/porosity decrease, and do so at roughly the same rate. That being the case, we substitute Equation (14) into Equation (13) to have at last

$c_{{v}}{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)={\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (15)

## Tidying up the Physics, Governing Equation, Boundary and Initial Conditions, and the Solution

It should be evident that getting to Equation (15) was a major triumph for geotechnical theory. But it’s also evident that there’s one glaring problem: the dependent variable on the left hand side is not the same as the one on the right. It is here that we need to clarify some assumptions behind our equations.

A basic assumption in consolidation theory is that, when the load at the surface is applied, all of this additional load is initially borne by the pore water pressure. Because of the aforementioned permeability, the water will want to “head for the exits,” i.e., the permeable boundaries of the layer being compressed. When the particles have rearranged themselves and the excess pore water has been squeezed out, the settlement should stop (until secondary compression kicks in.) During this process the load on the pore water is being progressively transferred to the soil particles until consolidation has stopped (which, in theory, it never does, as we will see) and the load is completely handed off to the particles.

That being the case, the consolidation equation (15) should be rewritten as

$c_{{v}}{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)={\frac {\partial }{\partial t}}u(x,t)$ (16)

At this point we should invoke the effective stress equation

$\sigma_x(x,t) = p - u(x,t)$ (17)

where p is the applied pressure, and do the following to determine the time and distance history of the vertical stress $\sigma_x(x,t)$:

• Determine the solution of Equation (16), which will be our governing equation.
• Determine the initial and boundary conditions for the problem.
• Solve the problem for $\sigma_x(x,t)$ using Equation (17).

There are several ways to solve the governing equation for this problem. Verruijt employed Laplace Transforms to accomplish this; these are very useful, as was demonstrated by Warrington (1997) for the wave equation. For this analysis, we will use separation of variables, a method which is not as fundamental as one might like (this is a demonstration of a method that is) but which is easier to follow than most of the others.

The separation of variables begins by assuming the solution is as follows:

$u(x,t) = X(x)T(t)$ (18)

This means that the solution is a product of two functions, one of time and the other of distance, and that those functions are separate one from another. Substituting this into Equation (16) and rearranging a bit yields

${\frac {{\frac {d^{2}}{d{x}^{2}}}X(x)}{X(x)}}={\frac {{\frac {d}{dt}}T(t)}{{\it c_v}\,T(t)}}$ (19)

Since the left and right hand sides are equal to each other, they are equal to a third variable, which we will designate as $\beta^2$. In reality they are equivalent to the eigenvalue $\lambda = \beta^2$, but the utility of the squared term will become evident. In any case,

${\frac {{\frac {d^{2}}{d{x}^{2}}}X(x)}{X(x)}}=-{\beta}^{2}$ (20a)

${\frac {{\frac {d}{dt}}T(t)}{{\it c_v}\,T(t)}}=-{\beta}^{2}$ (20b)

The solutions to these equations are, respectively,

$X(x)={\it C_1}\,\cos(\beta\,x)+{\it C_2}\,\sin(\beta\,x)$ (21a)

$T(t)={e^{-{\it c_v}\,{\beta}^{2}t}}{\it C_3}$ (21b)

We need to pause and consider the boundary conditions. As the problem is shown at the top of the article, the boundary conditions are as follows:

$X(0) = 0$ (22a)

$X'(h) = 0$ (22b)

Equation (22a) represents a Dirichelet boundary condition and Equation (22b) represents a Neumann boundary condition. While the equation for $X(x)$ can certainly be solved for this boundary condition, a simpler way would be to do the following:

• Mirror the problem shown above at x = h so that you have two permeable boundaries.
• You now have a layer of 2h thickness with Dirichelet boundary conditions on both sides. At the centre of this new layer, there is no flow; the water above it flows upward, and the water below flows downward.
• Problems in the field can have either one or two permeable boundaries. The distance h is NOT the thickness of the layer but the longest distance the trapped pore water must travel to escape. Confusing h with the thickness of the layer is a common mistake and should be avoided at all costs.

That said, the second boundary condition is now

$X(2h) = 0$ (22c)

Returning to Equation (21a), because of Equation (22a) $C_1 = 0$ as the cosine is by definition unity at this point. If we then set x = 2h and X(2h) = 0, for real, non-zero values of $C_2$ the boundary condition can be satisfied if and only if

$\beta = 1/2\,{\frac {n\pi }{h}}$ (23)

This is the square root of the eigenvalues. Any integer value of n > 0 is valid for this, and we will have recourse to them all to produce a complete orthogonal set (Fourier Series) to solve the problem.

This change will affect both Equations (21a) and (21b). Combining constants into the coefficient $B_n$, substituting the results into Equation (18) and making it an infinite sum for the Fourier series yields

$u(x,t) = {\it B_{n}}\,\sin(1/2\,{\frac{n\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{n}^{2}{\pi}^{2}t}{{h}^{2}}}}}$ (24)

At this point we need to consider our initial conditions $f(x)$. Although many different distributions of initial pore pressure are possible, the simplest one–and the one most commonly used–is a uniform pressure p, which is the same as the surface pressure compressing the layer. For a complete orthogonal set, the coefficients $B_n$ can be computed by (Kreyszig (1988)) as

${\it B_n}=\int _{0}^{2\,h}f(x)\sin(1/2\,{\frac {n\pi \,x}{h}}){dx}{h}^{-1}$ (25)

Substituting $f(x) = p$ and performing the integration,

$B_n = 2\,{\frac{p\left(1-\cos(n\pi)\right)}{n\pi}}$ (26)

Substituting this into Equation (24) and taking the complete sum, we have at last

$u(x,t)=\sum_{n=1}^{\infty}2\,p\left(1-\cos(n\pi)\right)\sin(1/2\,{\frac{n\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{n}^{2}{\pi}^{2}t}{{h}^{2}}}}}{n}^{-1}{\pi}^{-1}$ (27)

More simply we can say that

$u(x,t)=\frac{4p}{\pi}\left(\sin(1/2\,{\frac{\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/3\,\sin(3/2\,{\frac{\pi\,x}{h}}){e^{-9/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/5\,\sin(5/2\,{\frac{\pi\,x}{h}}){e^{-{\frac{25}{4}}\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}\cdots\right)$ (28)

The vertical pressure on the soil skeleton is determined by combining Equations (17) and (28) to yield

$\sigma_{x}(x,t)=p\left(1-\frac{4}{\pi}\left(\sin(1/2\,{\frac{\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/3\,\sin(3/2\,{\frac{\pi\,x}{h}}){e^{-9/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/5\,\sin(5/2\,{\frac{\pi\,x}{h}}){e^{-{\frac{25}{4}}\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}\cdots\right)\right)$ (29)

It is interesting to note that, although Equation (27) includes all positive non-zero values of n, only the odd ones end up in Equations (28) and (29). For even values of n, $B_n = 0$.

## Conclusion

At this point we have derived the equations of pressure dissapation for consolidation settlement. In a future post we will deal with the second part of the problem, namely how much of the total anticipated settlement has taken place at any given time from initial loading.

## Reference

Kreyszig, E. (1988) Advanced Engineering Mathematics. Sixth Edition. New York: John Wiley and Sons

I posted this ten years ago today on another site. I had no idea that my PhD program would take some of the twists and turns it did before I walked in December 2016. With the challenges the U.S. faces these days, it would do well to consider this problem–and the technologies that were core to the program.

Next month, Lord willing, I start the second year of my PhD adventure in Computational Engineering.  “Adventure” is a good way to describe it, especially in my superannuated state.  In the midst of getting the coursework started last fall (yes, Europeans, sad to say we have coursework with PhD programmes) we had a “freshman” orientation.

One of the things my advisor oriented us about were the substantial computer capabilities that the SimCentre has.  “Substantial” is a relative term; at one point we were in the top 500 computers in the world for power, but he chronicled our falling ranking which eventually left us “off the charts”.  In spite of that the SimCentre continues to tout its capabilities and the rigour of its curriculum, and I can attest to the latter.

In my desperate attempt to re-activate brain cells long dormant (or activate a few that had never seen front-line service) I dug back and discovered that most of the “basics” of the trade were in place when I was an undergraduate in the 1970’s.  Part of that process was discovering, for example, that one of the books I found most useful for projects like this had as a co-author one of my undergraduate professors!  The biggest real change–and the one that drove just about everything else–was the rapid expansion of computer power from then until now.  What has happened to the SimCentre was not that the computer cluster there had deteriorated; it has just not kept up with the ever-larger supercomputers coming on-line out there.  Such is the challenge the institution faces.

One of my advisor’s favourite expressions is “in all its glory”, but in this case what we have is a case of “fading glory”.  The SimCentre isn’t the only person or institution to face this problem:

If the system of religion which involved Death, embodied in a written Law and engraved on stones, began amid such glory, that the Israelites were unable to gaze at the face of Moses on account of its glory, though it was but a passing glory, Will not the religion that confers the Spirit have still greater glory? For, if there was a glory in the religion that involved condemnation, far greater is the glory of the religion that confers righteousness! Indeed, that which then had glory has lost its glory, because of the glory which surpasses it. And, if that which was to pass away was attended with glory, far more will that which is to endure be surrounded with glory! With such a hope as this, we speak with all plainness; Unlike Moses, who covered his face with a veil, to prevent the Israelites from gazing at the disappearance of what was passing away.  (2 Corinthians 3:7-13)

When Moses came down from Mt. Sinai, having received the law directly from God (a truth that Muslims, along with Christians and Jews, affirm) he wore a veil:

Moses came down from Mount Sinai, carrying the two tablets with God’s words on them. His face was shining from speaking with the LORD, but he didn’t know it. When Aaron and all the Israelites looked at Moses and saw his face shining, they were afraid to come near him. Moses called to them, so Aaron and all the leaders of the community came back to him. Then Moses spoke to them. After that, all the other Israelites came near him, and he commanded them to do everything the LORD told him on Mount Sinai. When Moses finished speaking to them, he put a veil over his face. (Exodus 34:29-33)

Moses’ fellow Israelites were afraid to even come near him on account of the glory of God on his face, thus the veil.  And Paul dutifully replicates this reasoning in 2 Corinthians 3:7.  But Paul throws in another reason Moses needed to wear a veil: to hide the fact that, like the SimCentre’s computers,  the glory was fading!  That is what I call a “2 Corinthians 3” problem, and although Paul applies it to the Jewish law, it has broader application as well.

Today, in spite of our “flattened” society, we see people and institutions lifted up and presented as glorious.  This is especially clear in the “messianic” streak that has entered our political life, but it’s also clear in the adulation given to celebrities, corporations and their products.  It would behove us, however, to take a critical look at such things with one question: are we looking at real, enduring glory, or is this just another example of fading glory which is being covered up with hype?

As Paul and other New Testament authors note, the purpose of Jesus Christ coming to the earth was to solve the problem of fading glory, and specifically of a system that could not produce an entire solution to the problem of our sins and imperfections.   Once that happens within ourselves, we can stop living behind the hype of whatever fading glory we’ve tried to hide behind and live in the truth.

‘Yet, whenever a man turns to the Lord, the veil is removed.’ And the ‘Lord’ is the Spirit, and, where the Spirit of the Lord is, there is freedom. And all of us, with faces from which the veil is lifted, seeing, as if reflected in a mirror, the glory of the Lord, are being transformed into his likeness, from glory to glory, as it is given by the Lord, the Spirit. (2 Corinthians 3:16-18)

## If you want glory that lasts, take a look at this

Posted in Academic Issues, Geotechnical Engineering

## Explicit and Implicit Methods and Plasticity: A Warning for Code Writers and Users

It’s time to put a wrap on our series on constitutive equations, both elastic and plastic. We’ll do this in a more qualitative way by pointing out a trap in code writing that can have an effect on the results.

We generally divide problems into two kinds: static and dynamic. With static elastic problems, we can solve the equations in one shot, although we have to deal sometimes with geometric non-linearity (especially when elastic buckling is involved.) With dynamic problems time stepping is involved by definition. When we include plasticity, because of the path dependence issue we’re always using some kind of stepping in the solution. That stepping takes at least two levels: the stepping like we discussed in our last post and load and/or displacement stepping as the load or displacement is increased. Therefore, in a problem involving plasticity, we are stepping the problem one way or the other.

But how large do the steps need to be? In the case of dynamic stepping, that problem is wrapped up with the whole question of implicit vs. explicit methods. To make things simple, explicit methods are those which take information from the past and use that in each step to predict the state at the end of the step, while implicit methods use information from both start and finish of the step to predict the step end state. Implicit methods have the advantage of allowing the model to have larger time steps, which is one reason why they’ve gained currency with, say, computational fluid dynamics models.

Probably the most popular dynamic method used for analyses such as this is Newmark’s Method, which admits both static and dynamic implementations. An example of Newmark’s Method for a simple system (a ram/cushion/pile system) can be found here. In that analysis, except for inextensibility issues (the cap not sticking to the pile, for example) the analysis was entirely elastic.

It would seem that using an implicit method would be the optimal solution for a dynamic system. However, that runs into a serious problem, as discussed by Warrington (2016):

Consider the elasto-plastic model as depicted (above). At low values of strain, elasticity applies and the relationship between stress and strain is determined by the slope of the line, the modulus of elasticity. In the elasto-plastic models considered, the reality is that the relationship between stress and strain is always linear; the key difference between the elastic and plastic regions is that, upon entrance into the plastic region, there are irrecoverable strains which take place. Cook, Malkus and Plesha (1989) observe that, in the plastic region, there is a plastic modulus, which is less than the elastic modulus and, in the case of softening materials, actually negative. They also observe that, in this region, the acoustic speed is lower than that in the elastic region…This is a similar phenomenon to that of the variations in elastic modulus and acoustic speed based on strain, which complicated the determination of the applicable soil properties.

With a purely elasto-plastic model, the plastic modulus is zero, and thus the acoustic speed is also zero. This effectively decouples the mass from the elasticity in the purely plastic region. This result is more pronounced as the time (and thus the distance) step is increased; the model tends to “skip over” the elastic region and the inertial effects in that region. Thus with larger time steps inertial effects are significantly reduced, and their ability to resist pile movement is likewise reduced.

For problems such as this, the best solution is to use an explicit method with very small time steps to “catch” all of the effects of the plasticity. One major advantage of that approach is that explicit methods allow us to dispense with assembling the global stiffness matrix, giving rise to “matrix-free” methods you hear about. That significantly reduces both memory requirements and computational costs in a model. It has also led to changing even problems we consider as “static” (like static load testing) to dynamic problems with small time steps, no stiffness matrices and long run times, an acknowledgement that, strictly speaking, there are no real “static” problems, a fact that makes many civil engineers freak out.

Chances are that the effects of this have been built into whatever code you might be using, and that you’re not aware of how your code is handing it. The wise user of numerical modelling solutions would do well to familiarise him or herself of how the code being used actually handles problems such as this and becomes a more informed user of the tools at hand.

Posted in Academic Issues, Geotechnical Engineering

## An Overview of Mohr-Coulomb Failure Theory

We now move from elastic considerations to plastic ones. The material that follows is adapted from this dissertation, the development of the STADYN code.

As we transition from elastic to elasto-plastic consideration, there are two things that we need to consider.

The first is that there is more than one way to implement plasticity in finite element and finite difference codes. The one we’ll be considering here–elastic-perfectly plastic behaviour using Mohr-Coulomb failure theory–is just one of them. It appears in many FEA codes, and although is justifiably described as ‘crude’ (Massarsch (1983)) it is the basis (explicit or implicit) of our current testing scheme and many of our ‘closed form’ solutions.

Alternatives not only include models such as Cam Clay but also models such as hyperbolic models. The latter have been discussed extensively on this site and will not be considered further here.

The second is that, because it is plasticity, the response of the soil to load becomes path-dependent. With elastic solutions, the response is path independent; you apply a certain load and the soil responds in a certain way irrespective of how the load is applied. With path dependence how the load is applied affects how the soil responds. This is a large part of another characteristic of plastic solutions: uniqueness issues. As a practical matter true uniqueness, which means a linear transformation which is one-on-one and unto, is impossible. We can only attempt to come up with the best solution. Whether that ‘best solution’ is acceptable depends upon the methodology used and the application. In many cases it is, in some it is not.

### Failure Theory

The following is not meant to be a comprehensive treatment of the subject. Much of it is drawn from Nayak and Zienkiewicz (1972) and Owen and Hinton (1980).

According to Mohr-Coulomb theory, failure occurs when the combined stresses find themselves outside of the failure envelope defined by the equation

$\tau=c+\sigma sin\phi$ (1)

This is illustrated in Figure 1.

It should be noted that both Equation 1 and Figure 1 are based on the typical geotechnical sign convention of positive compression.

Failure takes place when Mohr’s Circle for the stress state either intersects or goes above the failure line. The state of intersection, in terms of the principal stresses, is expressed by Verruijt as

$\sigma_{1}-\sigma_{{3}}-2c\cos(\phi)+\sin(\phi)\left(\sigma_{{1}}+\sigma_{3}\right)=0$ (2)

For many applications, such as triaxial testing, this is sufficient, as the goal is to determine parameters c and φ from the tests. But what if it is necessary to analyze stress states other than those on the failure line, as is certainly the case with finite elements? For these cases the failure criterion can be defined as

$\sigma_{1}-\sigma_{{3}}-2c\cos(\phi)+\sin(\phi)\left(\sigma_{{1}}+\sigma_{3}\right)=F$ (3)

This is illustrated in Figure 2.

There are three possibilities for the right hand side:

1. F<0 , failure has not been achieved.
2. F=0 , failure has been achieved.
3. F>0 , the stress state is beyond failure.

How the model responds to each of these states depends upon how the response is modeled. In any case, once F ≥ 0, the effects of further stress are irrecoverable.

For use in finite element code, it is frequently more convenient to express these using invariants. For the case of plane strain/axisymmetry in this problem (and the general case for problems of this kind) it is assumed that

$\tau_{yz}=\tau_{xz}=0$ (4)

and the first invariant is

$J_{1}=\sigma_{x}+\sigma_{y}+\sigma_{z}$ (5)

The deviator stresses are defined as

$\sigma'_{x}=\sigma_{x}-\frac{J_{1}}{3}$
$\sigma'_{y}=\sigma_{y}-\frac{J_{1}}{3}$ (6)
$\sigma'_{z}=\sigma_{z}-\frac{J_{1}}{3}$

The second and third deviatoric invariants are

$J'_{2}=\frac{\sigma{}_{x}^{'2}+\sigma{}_{y}^{'2}+\sigma{}_{z}^{'2}}{2}+\tau_{xy}^{2}$
$J'_{3}=\sigma'_{z}\left(\sigma{}_{z}^{'2}-J'_{2}\right)$ (7)

The whole business of invariants, standard and deviatoric, is discussed in further detail here.

With all of this, Equation 3 can be restated thus:

$\sqrt{J'_{2}}\left(cos\theta-\frac{sin\theta sin\phi}{\sqrt{3}}\right)-c\cos(\phi)+\frac{J_{1}}{3}\sin(\phi)=F$ (8)

where

$\theta=-\frac{1}{3}arcsin\left(\frac{3\sqrt{3}}{2}\frac{J'_{3}}{J_{2}^{'\frac{3}{2}}}\right)$ (9)

The Mohr-Coulomb failure criterion is frequently depicted using a three-dimensional representation on the principal stress axes as is shown (along with the Drucker-Prager criterion) in Figure 3. On the left is the failure surface in true three-dimensional representation, and on the right is same in the octahedral plane. The significance of Lode’s Angle $\theta$ can be clearly seen.

Inspection of Equation 8 shows that the failure function F is not the result of a unique combination of stresses. Additional information is available in the plastic potential function, which in turn is a function of the dilitancy of the material. The simplest way of determining this is by substituting the dilitancy angle ψ for the friction angle φ in Equation 8 (Griffiths and Willson (1986)), or

$\sqrt{J'_{2}}\left(cos\theta-\frac{sin\theta sin\psi}{\sqrt{3}}\right)-c\cos(\psi)+\frac{J_{1}}{3}\sin(\psi)=Q$ (10)

## References

Abbo, A.J., Lyamin, A.V., Sloan, S.W., and Hambleton, J.P., “A C2 continuous approximation to the Mohr-Coulomb yield surface”, International Journal of Solids and Structures 48, 21 (2011), pp. 3001–3010.

Gourdin, A. and Boumahrat, M., Méthodes Numériques Appliquées (Paris, France: Téchnique et Documentation-Lavoisier, 1989).

Griffiths, D.V. and Willson, S.M., “An explicit form of the plastic matrix for a Mohr-Coulomb material”, Communications in Applied Numerical Methods 2, 5 (1986), pp. 523–529.

Massarsch, K. Rainer, “Vibration Problems in Soft Soils”, Proceedings of the Symposium on Recent Developments in Laboratory and Field Tests and Analysis of Geotechnical Problems (1983), pp. 539-549.

Nayak, G.C. and Zienkiewicz, O.C., “Elasto/plastic stress analysis. An generalisation for various constitutive relationships.”, International Journal of Numerical Methods in Engineering 5, 1 (1972), pp. 113–135.

Ortiz, M. and Simo, J.C., “An Analysis of a New Class of Integration Algorithms for Elastoplastic Constitutive Relations”, International Journal for Numerical Methods in Engineering 23, 3 (1986), pp. 353–366.

Owen, D.R.J. and Hinton, E., Finite Elements in Plasticity: Theory and Practice (Swansea, Wales: Pineridge Press, 1980).

Potts, D.M. and Zdravković, L., Finite Element Analysis in Geotechnical Engineering: Theory (London, UK: Thomas Telford Publishing, 1999).

Sloan, S.W., “Substepping Schemes for the Numerical Integration of Elastoplastic Stress-Strain Relations”, International Journal for Numerical Methods in Engineering 24, 5 (1987), pp. 893–911.

Verruijt, Arnold and van Baars, S., Soil Mechanics (The Netherlands: VSSD, 2007).

Warrington, Don C., “Application of the STADYN Program to Analyze Piles Driven Into Sand“, vulcanhammer.net (2019).

## Nomenclature

• $\epsilon$ Strain
• $\epsilon^e$ Elastic Strain
• $\epsilon^p$ Plastic Strain
• $\epsilon_{r,z,rz,\theta}$ Strains in Polar Directions
• $\lambda$ Lamé’s Constant, Pa
• $\lambda_p$ Elasto-Plastic Normality or Plasticity Constant
• $\phi$ Internal Friction Angle, Degrees or Radians
• $\psi$ Dilitancy Angle, Degrees or Radians
• $\sigma'_{x,y,z}$ Deviator Stress in Cartesian Directions, Pa
• $\sigma_{1,3}$ Principal Stresses, Pa
• $\sigma_{x,y,z}$ Normal Stressees in Cartesian Directions, Pa
• $\tau$ Shear Stress, Pa
• $\tau_{xy,yz,xz}$ Shear Stressees in Cartesian Directions, Pa
• $\theta$Lode’s Angle, Radians or Degrees
• $a_Q$ Plastic Function Vector
• $a_{1,2,3}$ Failure Function Subvectors
• $a_{\beta}$ Constant for Substepping Algorithm
• $a_F$ Failure Function Vector
• $c$ Cohesion, Pa
• $C_{1,2,3}$ Failure Function Subconstants
• $D^e$ Elastic Constitutive Matrix
• $D^p$ Plastic Constitutive Matrix
• $D^{ep}$ Elasto-Plastic Constitutive Matrix
• $F$ Mohr-Coulomb Failure Criterion, Pa
• $G$ Shear Modulus of Elasticity, Pa
• $J'_2$ Second Deviatoric Stress Invariant, Pa
• $J'_3$ Third Deviatoric Stress Invariant, Pa
• $J_1$ First Stress Invariant, Pa
• $K_t$ Tangent Stiffness Matrix
• $n$ Time or Newton Step Number
• $Q$ Plasticity Potential Function, Pa