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.