Category: Academic Issues

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.

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

General Remarks About Plasticity

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

(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

(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

(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

(4)

and the first invariant is

(5)

The deviator stresses are defined as

(6)

The second and third deviatoric invariants are

(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:

(8)

where

(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 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

(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

• Strain
• Elastic Strain
• Plastic Strain
• Strains in Polar Directions
• Lamé’s Constant, Pa
• Elasto-Plastic Normality or Plasticity Constant
• Internal Friction Angle, Degrees or Radians
• Dilitancy Angle, Degrees or Radians
• Deviator Stress in Cartesian Directions, Pa
• Principal Stresses, Pa
• Normal Stressees in Cartesian Directions, Pa
• Shear Stress, Pa
• Shear Stressees in Cartesian Directions, Pa
• Lode’s Angle, Radians or Degrees
• Plastic Function Vector
• Failure Function Subvectors
• Constant for Substepping Algorithm
• Failure Function Vector
• Cohesion, Pa
• Failure Function Subconstants
• Elastic Constitutive Matrix
• Plastic Constitutive Matrix
• Elasto-Plastic Constitutive Matrix
• Mohr-Coulomb Failure Criterion, Pa
• Shear Modulus of Elasticity, Pa
• Second Deviatoric Stress Invariant, Pa
• Third Deviatoric Stress Invariant, Pa
• First Stress Invariant, Pa
• Tangent Stiffness Matrix
• Time or Newton Step Number
• Plasticity Potential Function, Pa

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

(1)

For uniaxial stress, the following parameters are zero:

(2a)
(2b)

As a result, stress and strain vectors are

(3)

(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:

(5)

(6)

Considering the first rows of Equation (5) gives us

(7)

and from the next two rows

(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.

Let’s look at the top diagram from a UC test. The main uniaxial stress . If we go 90 degrees anti-clockwise around Mohr’s Circle from , 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:

(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:

(10)

(11)

Extracting the important results,

(12)
(13)
(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.

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:

The conditions of plane stress are thus:

(1a)
(1b)
(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

(2)

and its inverse

(3)

The stress and strain vectors, however, are

(4)

(5)

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

(6)

(7)

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

(8)

Solving for the z-axis strain,

(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: .
2. Since , 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

(10)

We then invert the inverse to obtain the 3 x 3 constitutive matrix. We use the forward formulation and multiply through to obtain

(11)

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