Lecture 17: Modeling of Elasto-Plastic and Creep Response I

Related Resources

Study Guide (PDF)

Readings

Sections 6.6.3, 6.6.4

Examples

Problem 6.26

References

Bathe, K. J., et al. “On Some Current Procedures and Difficulties in Finite Element Analysis of Elastic-Plastic Response.” Computers & Structures 12 (October 1980): 607-624.

Snyder, M. D., and K. J. Bathe. “A Solution Procedure for Thermo-Elastic-Plastic and Creep Problems.” Nuclear Engineering and Design 64 (March 1981): 49-80.

Sussman, T., and K. J. Bathe. “A Finite Element Procedure for Nonlinear Incompressible Elastic and Inelastic Analysis.” Computers & Structures 26 (1987): 357-409.

Kojic, M., and K. J. Bathe. “The 'Effective-Stress-Function' Algorithm for Thermo-Elasto-Plasticity and Creep.” International Journal for Numerical Methods in Engineering 24 (August 1987): 1509-1532.

Kojic, M., and K. J. Bathe. “Thermo-Elastic-Plastic and Creep Analysis of Shell Structures.” Computers & Structures 26 (1987): 135-143.

Eterovic, A. L., and K. J. Bathe. “A Hyperelastic-Based Large Strain Elasto-Plastic Constitutive Formulation with Combined Isotropic-Kinematic Hardening Using the Logarithmic Stress and Strain Measures.” International Journal for Numerical Methods in Engineering 30 (October 20, 1990): 1099-1114.

Montans, F. J., and K. J. Bathe. “Computational Issues in Large Strain Elasto-Plasticity: An Algorithm for Mixed Hardening and Plastic Spin.” International Journal for Numerical Methods in Engineering 63 (May 14, 2005): 159-196.

• Download English-US transcript (PDF)

The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high-quality educational resources for free. To make a donation or view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

PROFESSOR: Ladies and gentlemen, welcome to this lecture on non-linear finite element analysis of solids and structures. In this lecture, I'd like to continue with our discussion of the use of constitutive relations in non-linear analysis. In the previous two lectures, we talked about modeling elastic materials in non-linear finite element analysis. We talked about the use of linear stress-strain loss and the use of non-linear stress-strain laws as regard to the T.L. and U.L. Formulations.

I'd like to now discuss with you the modeling of inelastic materials, in particular elastoplasticity and creep in finite element analysis. Let us proceed as follows. Let us first briefly discuss inelastic material behavior as we observe them in the laboratory test results. Let us then briefly discuss how we model these laboratory test results in one-dimension analysis. And then I'd like to generalize these modeling considerations to 2D and 3D stress situations. So let me now walk over here to the view graph projector to share the view grass with you that I've prepared for this lecture.

The most distinguishing feature of an inelastic material is that the total stress is not uniquely related to the current total strain. of course, for an elastic material, the total stress is given uniquely by the total current strain. For the inelastic material, this is not the case. And to calculate the response of the material, we need the response history. We need to work with stress increments and strain increments throughout the total solution.

If we look at the equation that we will be operating upon in the analysis, we have on the left-hand side here a differential stress increment given via the elastic constitutive tensor times the total strain increment, differential total strain increment, I should say, minus the differential inelastic strain increment.

Now, notice that this here is really is the differential elastic strain increment. Total strain minus inelastic strain gives us elastic strain. Of course, always the differential of course. And this, therefore, is the basic relation that we will be using, that we are using, in inelastic response calculations. Notice that I did not include temperature effects here, which would mean I would have to subtract, in addition, the thermal differential strain increment. Later on, we will also consider such relation.

The inelastic response may occur rapidly all slowly in time, depending on the type of problem you're dealing with, depending on the type of problem you want to model in nature. If the inelastic response occurs very rapidly, then we model it using plasticity rules. In other words, in plasticity, the model assumes that the inelastic strain increment occurs instantaneously.

In creep, however, the model assumes that the inelastic strain occurs as a function of time, and in fact, this time scale here might be hours, days, or years. Of course, the actual response is a combination of these two phenomena in nature, and we therefore use plasticity and creep, or alternatively, we may also use a viscoplastic material model to model what happens in actual nature.

In our discussion now, we want to consider materially non-linear only conditions. In other words, we only would focus on material non-linearities. We know that, of course, all the information that we are discussing can directly be applied in an MNO analysis, Materially Non-Linear Only analysis, and it can also directly be applied in a large displacement, large rotation, by small-strain analysis, by using the total Lagrangian formulation.

We talked in an earlier lecture about the fact that using second Piola-Kirchhoff stresses and Green-Lagrange Strains, we can directly use the total Lagrangian formulation with these stress and strain measures to model a large displacement, large rotation response of materials by using the same material relationship that we use to employ in material non-linear only analysis, but with these two new stress and strain measures.

Well, this is once more summarized here, that in other words, in the total Lagrangian information, we use the Green-Lagrange strain component instead of the engineering strain component, and the second Piola-Kirchhoff stress component instead of the engineering stress component, and then we directly can, with the material relationship that we have established for an M&O analysis, directly with these material relationships, we can model large displacement, large rotation, but small strain problems.

Let us now start by considering some observations off the material response in laboratory tests. We only want to look at these observations schematically, and therefore, we do not give any details. We also note that regarding the notation that I want to use now, we will not anymore have the t superscripts on the stress and t superscripts on the strain. We will not have these t superscripts when we look at laboratory test results.

So whenever we have a material law, a laboratory material law that we're looking on, or a material law that is identified by laboratory test results, or, in fact, a material law that we will enter into a finite element analysis, I will not have the t superscript on stress and t superscript on strain anymore. We will, however, use those superscripts again when we take that material law and actually apply it in finite element analysis in an incremental solution.

The material behavior is obtained by performing a tensile test in the laboratory, schematically shown as here. We have a force p up there, a force p here. There is a gauge length. Given here, the initial length is l0. There's cross-sectional area A0, and of course, the engineering strain, notice I now leave out the t here, the upper superscript t. I leave that out, as I just explained.

The engineering strain at any time is given as shown here. The stress is given, shown here. Where l now is the length of this gauge length at that time. And here we have, shown in red, how the specimen actually were [? deformed. ?] Notice A0 here and l right there. So l minus l0 is really the extension of the gauge length. We divide that by the original length, and we get the engineering strain, the small strain, measure that we are so used to seeing in linear analysis. And here the engineering stress that we are so used to seeing in linear analysis, as well.

So this is the test from which we obtain our information regarding the material behavior. And let's look at, systematically, test results. Here we have seen engineering stress plotted vertically, and the engineering strain plotted horizontally. And the test result that we would have obtained are shown systematically in red. Notice fracture occurs here corresponding to an ultimate strain. Notice that the test, of course, is performed in the tensile region in an actual analysis. We also want to apply a constitutive relation in the compressive region. And it is assumed that the same behavior that we have in the tensile region is also applicable in the compressive region. So we simply reflect this curve into the compressive region as shown by the blue dashed line. These are the test resides obtained at a constant temperature.

Now let us look at the effect of the strain rate. If the strain rate increases, in other words, de dt increases, we can see that the test results would change schematically, as shown here. Also, there is an effect of temperature. If the temperature is increasing, but of course at a constant temperature, we perform the test, then as the temperature increases for the different tests, we would see that the behavior of the stress-strain law is shown here by these red curves, schematically.

This is the instantaneous response of the material. Now let's look at the material behavior in time dependent response. In other words, over a very large time. Now at a constant stress, inelastic strains develop. And this is an important effect for materials when the temperatures are high.

A typical response curve, engineering strain, plotted vertically here now, and time plotted horizontally, is shown on this view graph. Notice there is, instantaneously, of course, inelastic or elastoplastic strain that occurs. That's the one we just talked about. And then with time, there is an increase in strain, as shown by the red curve here. We talk about the primary range of the creep curve, a secondary range, and a tertiary range until fracture occurs.

This is a curve at which sigma is constant and the temperature is constant, as well. If we increase the stress, we, of course, get an increase in the instantaneous strain, as shown here, but also the actual creep strain will change. And that is indicated by these red curves. Once again, fracture up here. And these are the results now at a constant temperature, but at different levels of constants stresses. The levels of constant stress means here, one stress, there, another stress, there, another stress, and another stress. But this stress is constant with time.

There's also a very significant effect off the temperature on the creep strength, and that is shown on this view graph. Once again, strain plotted here times there. There's an initial strain, of course, instantaneously reached. That's the one we talked about earlier. And then at a very low temperature, almost no creep effect, as shown here. Basically a straight line. Then as the temperature increases, the material creeps, and of course, doing this for this curve, the temperature is constant.

For this curve, the temperature is constant in time. For this curve, the temperature is also constant, but at a higher level. And for this curve, the temperature is also constant in time, but still at a higher level. And once again, fracture occurring here, and of course, sigma is constant, as well. For each of these curves, we have the same stress. But a different temperature here, then there, then there and there, the temperature being constant in time for each of these curves.

these are some typical laboratory results that we would get. And now we want to look at, how would we model these laboratory test results in a finite element solution? And let's look first at the models we might use for a one-dimensional stress situation.

So we consider a one-dimensional situation as shown schematically here-- a bar under a load p giving a displacement u. Notice now we go back to introduce our left superscript time t, because we want to look now at this [? grade ?] times t. The stress at any time is given by this formula. The area being constant, of course, the cross-section of the area being constant. The strain at any time is given by the displacement measured here, at that time, divided by the original length. Remember, we're talking about infinitesimally small strain conditions.

Let us assume also that the load has now increased monotonically to a final value, p star. And let us assume that the time is long, so that inertia effects are negligible, so we can talk about a static analysis.

If you look at the response that we would see, we see the following. The load increases as shown here, and here we have the time axis. Notice T star is small. It's a very small time interval. And in this time interval, creep effects would not take place because T star is so small. If they do take place, their effect is negligible. That's the assumption.

So plasticity effects predominate here, and those we take into account. And then from here onwards, there will be no more plastic deformations in the material. Now the material will creep with time. And so all the inelastic strains are modeled from here onwards as creep strains, and from here to there as plastic strains.

Now let's see how we would model these strains. Here we see how we would model the plastic strains in this uniaxial situation. Notice we assume a bilinear material description, bilinear meaning a Young's modulus until the yield stress is reached, and then another constant value of a new Young's modulus, if you want to. It's not really Young's modulus. We call it ET. It's a strain-hardening model. Notice that ET is generally considerably smaller than E. Typically it might be 1% of E.

So we're coming up here. After the yield stress, we have a constant Young's modulus, and beyond the yield stress, we have a strain-hardening modulus, ET, to describe the change in stress as a function of the strain. Notice now that if we had moved up to here and we were to unload from this level of stress, we would follow the original Young's modulus in the unloading that we had in the loading. In other words, here we had the same E in unloading, indicated by the dashed line, as we have for loading here.

In our particular modeling situation, we are coming up here to the yield stress, then we proceed, going along ET, up to that point, and that point corresponds to our time t star, where P star has been reached. The maximum load has been reached.

Notice that this part here is a plastic strain corresponding to t star, and this part here is the elastic strain, corresponding to t star. And that is obtained by simply dropping a vertical line down from that point t star onto the horizontal axes. The stress is then given by the Young's modulus E, that's E, or the slope is E here, times te at any particular time. In other words, we could have here, instead of the t, the t star.

At any particular time this holds, and the total inelastic strain is given by simply the plastic strain, because remember in our model, we assume that t star is so small that all the inelastic strains can be modeled using the plasticity assumption. This assumption. Of course, once again, these relationships hold for any particular point on this curve here, in particular, the point time t star, where these are the values that we're using here.

Now from time t star onwards, we have creep effects, and those would be typically modeled as shown on this view graph. Here is t star, which is very small, almost equal to 0. And here we have the time axis, the creep strain axis, shown here, and here we have a typical creep law that is widely used in engineering practice, the power law, shown here. A0, A1, A2 are constant. That has to be determined from the laboratory tests, and they will change for different materials.

Notice that the stress is still generated by the Young's modulus times the elastic strain, and the total inelastic strain now is made up of the creep strain obtained by this formula plus the plastic strain, the plastic strain already identified on the earlier view graph.

Another way of modeling the instantaneous inelastic strains, which we call plastic strains, and the time-dependent inelastic strains, which we call it creep strains, is to use a viscoplastic material model. And in this case, we have this basic relationship here where we say that the total strain times a time derivative of the total strain, is given by the time derivative of the stress divided by the Young's modulus plus a fluidity parameter times this bracket here. And this bracket here needs to be read as shown down here. Notice it is equal to 0. If sigma is smaller than sigma yield, the yield stress, and is equal to sigma minus sigma yield, sigma is greater than sigma yield.

This part here is the time derivative of the viscoplastic strain. Of course, gamma depends on the material-- on the particular material that you're looking at. And also, sigma yield will change possibly with time, because there might be strain hardening in the material, the same kind of strain hardening that we have observed in our plasticity description.

Typical solutions for the 1-D specimen would be, if you have a non-hardening material, in other words, the sigma yield, the value does not increase, then we would have here in elastic strain, and with gamma, determined by gamma, a viscoplastic strain that, of course, varies with time. And you can see the total strain here is made up of this elastic strain plus a viscoplastic strain, which also depends on gamma. If we have a hardening material, then the situation for the total strain looks as shown on this side of the view graph. Total strain plotted, again, vertically, time here. We have an elastic strain. And with time, a viscoplastic strain that, however, will reach a steady state value, an ultimate value.

This amount here depends on the increase of the yield stress as a function of the viscoplastic strain. Notice that also, this curve, this red curve, depends on the value of gamma. As gamma increases, you are getting a curve that goes more sharply up to this green line, this steady-state solution.

Let us now look at plasticity in more detail. In other words, how we are more dealing plasticity relationships. plastic materials, phenomena of plasticity in finite element analysis. So far, we've only considered loading conditions. However, we of course, in an actual analysis, might have that the material is loading, unloading, again reloading and so on. So we have to look briefly at, how do we take such situations into account? Well, on this view graph, systematically, we see what kind of load history a material might be subjected to. Here we have a load increasing. Say, up to here, we have elastic conditions, then we go into plasticity. At this point, we unload elastically, and right here, we go back into plasticity. From here to there, plasticity, and then elastic behavior again.

So how would we model such a load situation on the material level? That is what we like to discuss now. Here we show a material assumption, or a material relationship, that would be typically used. Stress plotted up, strain plotted horizontally. Here we load the material, going to plasticity in the loading cycle, in the loading regime, and now we unload elastically. Up to this point, we are elastic. We now are plastic, and we unload again elastically, from the compression region.

Now notice in this particular material description, we assume isotropic hardening. And this assumption means that once we have loaded up to this point, the new yield stress is given by this value here. And before we would yield again, in compression, typically, or if we unload here and come back, of course, intention, we would have to have a gain as stress that is equal to this yield stress. So notice in compression, we have this value being here, the same as that one. And that means we have to go all the way down here until we will get back into plasticity, and at this point, we go into plasticity, as shown here.

The total plastic strains is a sum of this part plus that part. All of these are plastic strains that have been accumulated in the loading cycle. Notice, once we have loaded in compression further, we have a new yield stress, sigma yield 2, which is larger than sigma yield 1. So whenever we go into the plasticity regime and load further, our yield stress increases, and that means we assume isotropic hardening.

I should also point out at this point that we make an assumption of the bilinear material behavior, namely, bilinear in having one Young's modulus here and one strain-hardening modulus here. If you have a complex stress-strain relationship measured in laboratory, you may want to use one et up to a certain stress level, and then another et from that stress level on to the next stress level, and maybe another et beyond that. So you have a multilinear description of the material behavior. We don't have time to go into such details. It's an extension of this material description.

Instead of using isotropic hardening assumption, we may want to use a kinematic hardening assumption. And that assumption is depicted on this view graph. Notice stress, again, upwards, strain, horizontally, total strain horizontally noted, and here we're coming up in our loading cycle to this point, and the unload now, same way as before in isotropic hardening-- the difference is right now occurring in this regime, namely our yield stress that will make us yield in compression is not the same yield stress that we talked about in isotropic hardening. Instead, this difference in stress is simply twice the initial yield stress.

So we reach yield in compression much earlier than in isotopically hardening. Namely, when the stress has decreased by twice this amount here, we come into the regime of yielding again. And now we would be yielding all the way along here, up to this point, and now we unload elastically again. Notice, this is now the plastic strain that we accumulate in the compression region. This is the plastic strain that we accumulate in the tensile region.

These were a few words regarding how we model one-dimensional elastoplastic behavior, and let us now look at the multi-axial plasticity behavior. How we model, in other words, plastic behavior in multi-axial stress conditions. In that case, we need a yield condition, a flow rule, and a hardening rule, and I'd like to consider now isothermal, in other words, constant temperature conditions. The yield condition flow rule and hardening rule are expressed using a stress function. And very widely used stress functions are the Von Metus function and the Drucker-Prager function.

The Von Mises function is given on this view graph here. tF is equal to 1/2 tsij tsij minus t kappa. Of course, we're summing here over all stress components sij where the stress component sij is called the deviatoric stress, and is given as t sigma ij minus t sigma mm over 3. This is summing all the normal stress components, in other words, t sigma 11 plus t sigma 2 2 plus t sigma 3 3, and dividing by 3. This is the Kronecker delta that we encountered already earlier. The t kappa value is given here as 1/3 t sigma yield squared. Notice the yield stress here changes as a function of loading, as a function of time. And in fact, we are here using isotropic hardening conditions.

In the Drucker-Prager stress function, assumption we have this equation here. tF is equal to 3 alpha t sigma m plus t sigma bar minus k. Notice t sigma m is defined as here. Notice the sum here over all the stress is ii, in other words, t sigma 1 1 plus t sigma 2 2 plus t sigma 3 3. That is t sigma i i. And here we have t sigma bar defined, again in terms of the deviatoric stresses.

We will talk in this lecture more about the Von Mises yield condition, which is very widely used for modeling metal behavior. The Drucker-Prager yield condition is widely used to model rock and solid mechanics types of problems. We'll not talk in this lecture about this condition any further. Please refer to the textbook, where there is a bit more information given about the Drucker-Prager stress function.

We will also use, in our discussion, both matrix notation and index notation. In matrix notation, for example, the vector DEP of plastic strains is defined as shown here. Notice we are summing these two components here for the shear strains. The stress vector is shown here, and if we use index notation, we imply these components. dep 1 2 and dep 2 1. Notice, once again, in these vector or matrix notations, these two components are added to get the total plastic shear strain, 1, 2.

The basic equations are, then, for the Von Mises yield condition, as shown here. This is a yield condition. tF, t sigma ij, t kappa. t kappa is a function of the plastic strains. These are the current stresses. And [UNINTELLIGIBLE] yielding tF is equal to 0. That's what we're saying here.

The equivalent for the one-dimensional response, in other words, would be this equation. Notice once again that tS is equal to 0 throughout the total plastic response, throughout all of the yielding that we observe. The flow rule, and here we use the associative flow rule, because tF is here the same tF that we also use in the yield condition. The flow rule is given by this equation. Notice we have on the left-hand side the plastic strain increments. Here we have a scalar, t lambda. And the partial of tF with respect to the current stresses. t lambda is a positive scalar. In one-dimension analysis, we would have that the plastic strains 1 1, 2 2, and 3 3 are as shown here. And notice, of course, they only depend on the one-dimensional stress, or the one stress that is only active. And notice that the sum of these plastic strains is 0.

This stress-strain relationship is given here on this view graph-- stress increment on the left-hand side, elastic stress-strain, material [UNINTELLIGIBLE] go into the CE. Here's the total strain increment minus the plastic strain increment. The 1-D equivalent, in other words, for the one-dimensional stress situation, we would have this equation here.

Now our goal is to determine a CEP matrix such that we have on the left-hand side the stress increment, and on the right-hand side, that CEP matrix times the total strain increment. Of course, the kind of stress conditions, yielding, et cetera, conditions must all go into the determination off CEP.

Let us now look at the general derivation of this CEP matrix. This is quite a general derivation. It's also given in the textbook. You may want to refer to the textbook as well while we discuss this deviation or afterwards.

We define these quantities here, tqij being the partial of F with respect to the current stresses with the plastic strains fixed. We also define tpij equal to minus partial F with respect to the plastic strains with the stress fixed. These are two convenient definitions for the derivation. If we use matrix notation, we assemble in this vector at tq-- we are writing actually, the transposed here, just for ease of writing-- these components are assembled in this vector. Notice we have 2 times tq12, 2 times tq 2 0 3, and 2 times tq 3 1. And these tools here result from our definition of the plastic strain and stress increment vectors.

Remember in the plastic strain vector, we had the sum of Dep 1 2 plus Dep 2 1 components, we got the total plastic share strain assembled in that plastic strain vector. And because we have defined our plastic strain vector and our stress increment vectors as shown earlier on the view graph, we now carry these two along here. tp, transposed again, is written out here.

We know determine t lambda in terms of de. That is our first step. Because if we have this relationship, then we can say, we have the plastic strain increment in terms of a total strain increment. Using tf equal to 0, throughout the plastic deformations, because we know that our yield condition has to be satisfied throughout all of the plastic deformations, we can directly write down this equation here.

And with our definition just given, we rewrite this right-hand side as shown here. And we also substitute for dep, t lambda times tq. Notice, once again, we have now in the fourth fifth, and sixth component of tq this factor two, because of the definition of dep shown earlier. Well, there's nothing else than a matrix rewriting of this relationship here, of this right-hand side here. And we now use this relationship in our further developments.

But before doing so, let us see how we can rewrite this left-hand side here. tq times [UNINTELLIGIBLE] t sigma. We can write this left-hand side also as shown on the right-hand side here, because we can substitute for t sigma simply the stress-strain law. The flow rule, of course, gives us this relationship here for the plastic strain increment. Hence, if we use what we have up here and substitute for the plastic strain increments, we directly obtain what I'm circling here.

Next we use what we derived on the previous view graph, namely, tq times post t sigma is equal to what's shown here on the right-hand side, what I'm circling now. And this together then gives us an equation that we can solve for t lambda. Having solved for t lambda, and I will show you the result on the next few graph charts now, we can substitute 40 lambda back here, and we get our plastic strains directly.

Well, let's look at the result. Here we have t lambda, notice, given in terms of the differential total strain increment. And now substituting, as I said just now, 40 lambda, we get the plastic strains in terms of the total strain increment.

This relationship is now used in the stress-strain law again. Notice the stress-strain law is shown here. And we directly obtain, by expressing dep in terms of de by the relation that we just have obtained on the previous view graph, we directly obtain this final expression for CEP. Of course the current stresses go in there, and the current material conditions go in there.

And we will apply this general expression to an example just now, namely the example of the Von Mises flow rule and yield condition. However, before doing so, we know need to break, because we have to continue our discussion on a second reel.

So let us now apply this general relationship here which we just derived to the case of the Von Mises yield condition with isotropic hardening. In that case, yielding is reached when the principal stress is here on the right-hand side, substituted into this right-hand side, give us the yield stress, t sigma yield. Or another way of saying it, our tf function is given right down here, where we have here the deviatoric stresses squared, t kappa is shown here.

We had in fact this relationship already on an earlier view graph. Notice that this relationship on TF is really the same as what you see up here. But of course, here we are talking about the use of principal stresses, whereas here we are talking about the use of deviatoric stresses.

Pictorally, we can represent the yield surface as shown here. In plain stress analysis, this would be the yield surface. In other words, t sigma 2 and t sigma 1 are the principle of stresses that are generally non-zero, but t sigma 3 is 0. Here we have the yield stress, and on the right-hand side, we see the [? general ?] yield condition in three-dimensional stress space.

Notice initially we have sigma yield, [UNINTELLIGIBLE] radius of this cylinder, and that radius expands with time, and it expands uniformly, as shown, with the green circle here to t sigma yield. In other words, the cylinder grows uniformly, the radius grows with isotopically hardening. Notice that for the plane stress case, the hardening would be stress by having this ellipse grow, as briefly indicated, by this green line, which of course goes all the way around it.

Let us now compute the derivative that we need to have. First we consider tpij. tpij was defined as shown right here. And simply substituting for tf in terms of the deviatoric stress components and the current yield stress, and differentiating, we obtain this part here.

Notice t sigma ij is fixed. Therefore, these are constant here, and of course drop out. We get this result, and all we have to evaluate the partial of the yield stress with respect to the plastic strain. That is achieved by the concept of the plastic work. The plastic work per unit volume is the amount of energy that is unrecoverable when the material is unloaded. This energy, in other words, has been used up to create a plastic deformation. What we will do is use this concept for the one-dimensional case, and then generalize it to the multidimensional situation.

Pictorally, in the 1D example, in other words, only one stress component being non-zero, one normal stress component to be non-zero, we have this stress-strain law. You're looking at the bilinear case, which we discussed earlier. And the plastic work is simply obtained by this integral here. Stress times plastic strain. Integrate it over the whole path. So this shaded area is equal to the plastic work.

In general, of course, we have the current stress times the plastic strain and in component forms. And we are summing over all the components. If we look at the one-dimensional case, we can directly express twp as shown here, in terms of the Young's modulus and the strain-hardening modulus, and the yield stress at time t, this one right there, and the initial yield stress, that one right there. This expression here gives the shaded area. We can now differentiate and obtain the yield stress differentiated with respect to the plastic work, this expression down here.

Using this expression in all our general multidimensional dimensional case, we directly obtain tpij. In terms of what's shown here on the right-hand side, notice we take sigma yield stress at time t with respect to the plastic work corresponding to time t, and then partial plastic work with respect to the plastic strains. Substituting here, we get this bracket for the sigma yield with respect to the plastic work, and plastic work with respect to the plastic strains is equal to the current stress. And immediately we arrive at this final relationship here for tpij.

Alternatively, we could have also used the effective stress in effective plastic strain. The effective stress is defined as shown here, using the deviatoric stresses, and the effect of plastic strain is defined here. Notice this is a differential increment in the effect of plastic strain. And if we use these two quantities, we can also obtain tpij via this relationship, and the result is of course the same that I just showed you on the previous view graph.

Let's now consider tqij, because that is the other quantity we need. And here we recognize we have to take the partial of f with respect to the current stresses with the plastic strains fixed. We substitute for f, and we differentiate, recognizing now with the plastic strains fixed, the yield stress is a constant. Therefore, all we need to do is differentiate the product of the deviatoric stresses, and that is achieved here.

The substitute for the deviatoric stress is shown on the right-hand side here, and recognize that the t sigma kl differentiates with respect to t sigma ij is nothing else than the product of these two Kronecker deltas. Delta ik, of course being equal to 1 when i is equal to k and equal to 0 when i is not equal to k. And we also recognize that this relationship here reduces directly to that relationship, and the whole right-hand side then becomes simply the deviatoric stress at time t.

With these relationships for tqij and tpij substituting into our general derivation of CEP, we directly obtain this CEP matrix for the three-dimensional case of the Von Mises flow condition which isotropic strain-hardening. Let's look at some of these terms. In front, we have the Young's modulus and the Poisson's ratio down here. Young's modulus here, Poisson's ratio here. Notice here, again, Poisson's ratio, the same way as we have it in elastic analysis, here we have now the deviatoric stress component entering. Similarly here, the deviatoric stress component entering. And so on.

There's also a factor beta that goes in there, and beta is defined down here in terms of the yield stress, Young's modulus and strain-hardening modulus, and Poisson's ratio.

So this is a general expression, and let's now see how good we use that general expression. Well, it would enter in the evaluation of the stresses from time t to time t plus delta t. We have a stress increment which is added to the current stress at time t to obtain the stress at time t plus delta t. This relationship here is expressed as shown here with the CEP matrix. In other words, d sigma is given as CEP times de. Now, this integration has to be performed at every integration point used in the finite element. Here we're showing 3 by 3 integration used in an 8-node, two-dimensional element. And notice once again that this integration has to be performed at every one of these integration points, which of course can be quite expensive. So we need to perform the integration as effective as just possible.

One simple scheme to use is Euler forward integration. And without subimplementation-- I will explain subimplementation just now-- without subincrementation, we evaluate the integral on the left-hand side here by simply saying, let us evaluate CP corresponding to time t. Since we know the stress at time t, we know the deviatoric stress at time t, we can directly evaluate CEP, entering into that matrix equation that I gave you on one of the previous view graphs. You can get CEP at time t, since we know the stress at time t, and simply multiply it by the total increment in strain.

Now, this would be OK if the strain increment is very small, and consequently also, say it leads to a very small stress increment. But in general, this linearization involved in using the Euler forward method can introduce a too-large error that we don't want to see in the analysis, and therefore, we use what we call the subincrementation.

Within subincrements, the subincrementation proceeds as shown this view graph for the evaluation of this integral. Notice CEP, times t, times delta e the over n, where n is the number of subincrements, once again. This means that we now can evaluate the stress at the end of the first subincrement. We use that stress, which corresponds to t plus delta tau-- tau is delta t over n-- to evaluate this CEP matrix and multiply the result again by delta e/n.

Now this then, together, gives us this stress at the end of the second subincrement. And like that, we proceed to forward march the stress over each of the subincrements up to time t plus delta t. This is nothing else really but the application of the previous formula, the one on the previous view graph, over each subincrement, always updating, of course, the stress that goes into the evaluation of the CEP matrix.

Pictorally, we have the following situation. Along here, we are plotting time t we have here, time t here, time t plus delta t there at the end. Here we have the subincrements 1, 2, 3, et cetera. You know that by the dot. And notice at the end of the first subincrement, we have time t plus delta tau. At the end of the second subincrement, we have time t plus 2 delta tau, and so on. The subimplementation schematically proceeds by calculating this stress increment so that we have t plus delta tau sigma, then evaluating a new stress-strain law, based on this stress, multiplying by the subincrement in strain, total strain, to get a stress increment. Thus you get this stress, and like that, you proceed, forward marching, again and again, until you have the final stress at time t plus delta t.

Let us now look at the summary of the procedure the way that you would use it in a computer program, and that summary is given on the following view graph. We assume that the strain corresponding to time t plus delta t is known. We have calculated that. The stress corresponding to time t is also known, and of course the strain corresponding to time t is also known. We then, as a first, meta-calculate a strain increment. DELEPS is equal to STRAIN minus EPSILON. This is a strain increment from time t to time t plus delta t.

We use that strain increment to calculate a stress increment, assuming elastic behavior. CE denotes a matrix of elastic material behavior. In other words, in this matrix go the Young's modulus and Poisson's ratio, and of course this one varies depending on what kind of conditions, plane-stress, plane-strain, axis-symmetric, or three-dimensional analysis you're looking at. But I think you know what I mean by that.

Having delsigma, we calculate tau by updating sigma for that stress increment. With this stress now, we check whether the strain was actually taken elastically. In other words, whether our assumption that we have made up here was correct. And that is done by entering into the yield condition with that tau. If f of tau is smaller or equals 0, then the strain increment is elastic, and in this case, tau was correct, and we return If that was not the case, then we have to continue in our solution procedure, and we do so as follows. If the previous state of stress was plastic, we set this ratio to 0 and go to g, because the whole strain, then, was taken up plastically, also in this step. Otherwise, if this condition is not true, we have to find that stress value at which yielding started.

And we do also by entering into the yield function with sigma plus ratio times delsigma, where ratio is now unknown. Delsigma is known, sigma is known, so we can calculate from this equation here the value of ratios. And ratio gives us that stress increment, which brings us to initiation of yielding.

Having ratio, we update the stress from sigma too tau, a new value of tau, where ratio goes in here, and we calculate the elastic-plastic strain increment, using also ratio. And that is a total elastic-plastic strain increment here. This value depth now is used integrate out the stress occurring during the plastic deformations.

And that is done by the subincrementation process that we discussed earlier. We divide depth into subincrements DDEPS, and we calculate this tau value by cycling through this subincrementation here, which we described on an earlier view graph. And this then completes the calculation of the stress value corresponding to time t plus t. So this is a very compact summary of how the stress integration from time t to time t plus delta t is performed. Of course, there are many important details here that we haven't really looked at in depth, but I think it gives you a summary that is quite valuable, because it tells, overall, how the integration is being performed.

Let us now look at two example solutions. First, an example solution for which we have prepared some slides, and after that, an example solution for which we have prepared an animation. Let me walk over here to share with you the information on the slides. And here we consider the following problem. You have a rigid punch here subjected to a load, lying on an infinite half-space. As the load increases, plasticity develops down here, and there is an ultimate load, a limit load, that we want to calculate. In other words, there's a load here at which there will be very large deformations down here, and the material can take no higher load. This plane strain punch problem was considered by many investigators, but Hill actually predicted the plastic zones, and also the ultimate limit load. And we want to compare our finite element solution with the solution given by Hill.

On the next slide now, we see the finite element mesh used. Notice we are only idealizing 1/2 of the domain, as shown by this boundary here. Notice also that we have moved this boundary here and that boundary here far enough away from the load application. Notice that we have a fine mesh here, a coarse mesh there. In the fine mesh, we use 8-node elements, in the coarse mesh, we use 4-node elements. And there is a transition region here, there, from the 8-node over to 7-node, into 6-node, and then finally into 4-node elements.

It's an interesting small point here. There is an incompatibility here between this element and these two elements. An incompatibility because this element here, being a 6-node element, has a parabolic displacement distribution along this line, whereas this element here, being a 7-node element and a linear assumption along this line here for displacement, only admits a linear displacement distribution. So there is small incompatibility here, and similarly here. However, that incompatibility we deem to be small, that it would not affect the results very considerably.

Now, in an elastoplastic plastic, analysis like in any non-linear analysis that you would like to perform, I pointed out in an earlier lecture, it is always important to start with the linear solutions, just to see how well your mesh performs and whether you have any data input errors, possibly, even. It is good to perform first a linear solution where you can compare with analytical or other results.

This is what we have done here, and let me show you on the next slide what we have considered. We have applied a point load here to the half space, used the same finite element mesh that I just showed you, and looked at the stresses. At this depth A, below the surface, it's 0.697 inches. At that depth, we're looking at the stresses along this line here on a horizontal line. These stresses here are plotted on this graph for different components. Sigma zz, Sigma yy, tau yz. And they are compared with the theoretical solution, the [UNINTELLIGIBLE] theoretical solution.

Notice excellent comparison, obtained really for the stress components, and this means that our mesh is quite reasonable, and we can proceed with the elastoplastic analysis. Notice please that we have used here two-point integration, 2-by-2 integration, for all the elements, including the 8-node element. We then performed the same analysis, also one, using three-point integration, and these results are shown on this slide. Same picture here. Here stress is as a function of the distance from the center line in inches. Notice with three-point integration, also excellent comparison for the stress components.

On the next slide now, you'll see the elastoplastic results. And here we see the applied load, normalized by a factor plotted versus the displacement of the punch, also normalized by factor. And we see two curves here, one corresponding to three-point integration and one corresponding to two-point integration.

Now notice that the three-point integration results are higher than those of the two-point integration results. And that's a very interesting phenomenon. It is due to the fact that the three-point integration puts a higher incompressibility constraint into the mesh. In fact, a too-high incompressibility constraint.. This is a very interesting phenomenon, but it is best discussed in another lecture.

Notice that however our solution results all give us a horizontal line here for the two-point integration and three-point integration, and certainly, the two-point integration results are quite close to the theoretical value. If we were to use the finer mesh and still use three-point integration, we would get also closer to the theoretical value. Of course, we have here, I should have pointed out, an elastic, perfectly plastic material assumption for the material under the footing. Please refer to the details of this example solution given in the paper in which we have reported this sample solution. And the reference to that paper is given in the study guide.

The second problem that I'd like to discuss with you concerns the analysis of a plate with a hole that we are subjecting to very large tensile forces. And there is one interesting feature that we observe in that analysis which we have not seen here on the slides, namely, the spread of plasticity through the domain. So let me introduce you to this problem using just two view graphs, and then show you the animation which gives the results to the problem.

Here we have the plate with a hole. It's subjected to tensile forces, so you can see. And we want to calculate, really, the ultimate limit load of the plate, and we want to see now the spread of plasticity. Of course, that spread of plasticity also occurs in the analysis of the rigid punch on an elastoplastic footing, but I did not show you that spread of placiticity there. I'd like to show you the spread the plasticity now here.

And you will see that the plasticity actually spread through the elements, integration point by integration point. And of course initially, as a load increase is not so rapid, but then all of a sudden you see a very rapid spread of the plasticity, right through this domain here, that domain there, and so on, by symmetry conditions.

Now, because of symmetry conditions, we also will look at only one quarter of the plate-- this quarter right here. The material assumption that is used in the analysis is shown on this view graph. Stress here, strain there. Notice Young's modulus, as shown here, Poisson ratio there, and the strain-hardening modulus is shown here. We assume isotropic hardening. The same problem would also be considered further in a later lecture, when we actually use the ADINA program to show an application. So you will encounter this problem there again.

Notice that in this particular analysis that you will be seeing now, in the animation, we assumed a materially non-linear only condition. In other words, we assumed infinitesimally small displacement and strains. However, in the very last load steps, just before the plate ruptures, we encounter, actually, larger strains, but only those last load steps. And in order to capture this very last behavior of the plate very accurately, one would have to perform a large strain analysis of the plate, which we did not do here, and which really would be the object off another lecture, as well.

So let's look now at the response of the plate. Yes, it is subjected to the tensile forces and how the plasticity spread through the plate until final rupture. Let's look at the animation. We will look first at a fast animation, where the response occurs very fast, and then we'll look at the same animation again but slower, so that I can better talk about it while you're watching it, and you can better understand what I'm saying.

Here we see one quarter of the plate discretized using 288 eight-node isoparametric plane stress elements. The time course gives the load step. We perform the static analysis and the load corresponding to each solution step is also shown. Here you see now the load increasing, and in the first load step, the plate remains elastic.

Plasticity then developed as a whole and spreads rapidly through the plate as the load increases.

The deformations of the plate are given here to a magnified scale. This was a rather fast look at the plate behavior, so let's look at it again and over a somewhat longer time scale. Here we have again the plate at time 0-- that is at zero-load application. Now the load application starts, and initially, the plate remains elastic.

Be sure that an element integration point has gone plastic by darkening the regions around it.

Since we use 3-by-3 Gauss integration, considering an element on an average 1/9 of the element [? area ?] corresponds to an integration point.

Here you can see how the plasticity spreads from element to element, and indeed, from integration point to integration point. We will discuss the analysis of the plate again in lecture 22. We will also perform analyses and study the results, including large displacement and large strain effects.

Once again, the deformations of the plate are shown here in the animation to a magnified scale. The actual displacements are discussed in lecture 22.

So this brings me then to the end of this lecture. I discussed everything I wanted to discuss with you in this lecture. Thank you very much for your attention.