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.
The title of this lecture is "Basic Considerations in
Non-Linear Analysis." In this lecture I would like to
discuss with you first the application of the principle
of virtual work.
We will spend quite a bit of time on the basic principle of
virtual work, because it's such an important ingredient
for non-linear finite element analysis.
In fact, the displacement-based finite
element method, which is widely used for non-linear
analysis, is based on an application of the principle
of virtual work.
I also would like to discuss with you and particularly
emphasize some basic requirements of mechanics.
There are a few very important requirements of mechanics that
it is important to always keep in mind, and I'd like to share
some experiences with you regarding these requirements.
And finally, to emphasize all the points that I have tried
to make in the lecture, I like to present and share with you
some experiences regarding some actual analyses.
And I've chosen two analyses.
First we will look at the analysis of a plate with a
hole, and then we will look at the analysis of a
plate with a crack.
Let me now go over to my view graph and start with a
discussion of the principle of virtual work.
The principle of virtual work, in essence, says that this
integral here on the left-hand side is equal to what is on
the right-hand side.
And we should, of course, discuss now in detail what we
have on the left- and on the right-hand side.
On the right-hand side, we have this tR--
this is a script R, an unusual symbol, but some time ago,
some long time ago, I decided to select this
symbol for this quantity.
This script t, superscript t R, contains two integrals.
The first integral is an integral over the volume tv.
Notice it's the current volume of the body of tfiB, the body
forces at time t multiplied by the virtual
displacement, del ui.
This product here is integrated over the
volume at time t.
We add to this integral another integral, which is an
integral over the surface area of the body.
And here we have surface forces tfiS per unit area, and
these are virtual displacements, delta uiS.
We integrate this product over the total area of the body.
So this defines the right-hand side.
On the left-hand side, we have one quantity here, which is t
tau ij, t once again referring, of course, to time.
tau is a stress, and ij being the components of the stress.
This stress is called the Cauchy stress.
It's a force per unit area at time t.
And on the left-hand side, we have also a quantity which is
a virtual strain.
Delta here meaning virtual, teij meaning the strain.
And it's defined as shown here.
Now let us look briefly at what these quantities are.
We have a partial of del ui with respect to the current
coordinate txj, and we have a partial of del uj with respect
to the current coordinate txi.
Here I have prepared a view graph which separates out,
once more, these important quantities.
And on the top here in red, we see the virtual displacement
and the virtual strain.
Notice that the virtual strain here carries a subscript t, a
left subscript t.
This left subscript t means that it's a strain referred to
the current configuration.
del always meaning virtual, of course.
And it's important to recognize the following.
It's important to recognize that these are the virtual
strains corresponding to these virtual displacements.
That's the most important point to recognize.
In green here, we have some other quantities--
the current volume, the current
surface area of the body.
Current meaning always at time t.
And here in blue, we see the externally applied forces.
tfiB being the body forces, the volume forces, and tfiS
being the surface forces per unit volume, per unit surface.
These are externally applied forces.
So we have here quite different quantities--
namely, virtual quantities and real quantities.
So let me go back to this view graph to discuss with you in
more detail what the principle really tells us.
It tells that the internal virtual work consisting of the
product of the actual stresses--
t tau ij, which in general analysis
are, of course, unknown--
times the virtual strains, which we impose, this internal
virtual work--
there's an integration over the whole volume--
must be equal to the external virtual work.
And the external virtual work, once again, is given here.
This equation must hold--
and this is an important point--
for arbitrary variations in displacement, or arbitrary
virtual displacements that satisfy the geometric boundary
conditions, that satisfy the essential boundary conditions,
or you might say, the displacement boundary
conditions.
In other words, this is an equation that must hold for
any virtual displacements and corresponding virtual strains.
Remember, I just said that for virtual displacement, we
always have corresponding virtual strains.
This equation must always hold, provided the
displacements satisfy the displacement boundary
conditions.
Here we have, once again, the unknown forces, or stresses, I
should say, forces per unit area.
Notice these are Cauchy stresses.
Here we have the word Cauchy.
In an infinitesimal displacement analysis, we
really only talk about one kind of stress--
the engineering stress.
In large displacement analysis, large strain
analysis, we have different kind of stress measures.
And the stress that I now talk about here is the Cauchy
stress, which is the force per unit area at time t, which is,
in other words, an actual physical stress.
Later on, we will introduce another kind of stress, a very
important stress measure, which, however,
is not quite physical.
This one here is a physical stress, the
force per unit area.
And that is the stress that you, for example, in a
computer program, would like to get printed out.
This is the force per unit area that you would design
your structure with.
This here is, once again, the virtual strain.
And notice that this strain is referred to the current
geometry, the coordinate at time t.
Therefore, stress at time t times virtual strain at time
t-- that's what we're looking at here, and that gives us the
total internal virtual work when we
integrate over the original--
over the current volume, current, I should say, very
important--
the current volume of the body.
So if we have a body that moves through space and
undergoes large displacements, large rotations, large
strains, then this here is the quantity of the internal
virtual work with the current stresses at time t in the
configuration at time t.
Let us look at a schematic example, because there are so
many important considerations.
Here we have, in a three-dimensional Cartesian
coordinate frame, denoted with x1, x2, x3 as coordinate axes,
a body in which we are showing two material particles.
Notice one of them in the body itself, in the volume of the
body, rather, and one on the surface.
The body is here supported, and this is the
configuration at time 0.
Now this body moves through a certain amount.
It undergoes large displacements, large
rotations, large strains, and takes on the red
configuration here.
Notice that in that red configuration, we have certain
body forces, which we already talked about earlier, and
surface forces, that we also talked about earlier.
Notice that the particles have now moved to new
configurations, of course, to new positions.
The surface particle has moved there, and the particle within
the volume has moved there.
We denote this as tui, the displacement at time t, or the
displacement from configuration 02,
configuration t, for this particular particle.
Similarly, here we denote this displacement as tuiS.
Now the principle of virtual work says that we will be
imposing on to this red configuration the
configuration at time t, a variation in displacement,
which is the virtual displacement.
And here we see now one such virtual displacement.
Notice that it is shown by the blue line for the whole body,
and it's a virtual displacement from the
configuration at time t.
Notice that the particle, which was originally here, or
which was at time t here, I should more precisely say, has
moved to this point in the virtual displacement.
And similarly, the particle that was here at time t has
moved to this configuration, to this point here.
Now, this is one variation.
However, the principle of virtual work, of course,
states that the left-hand side must be equal to the
right-hand side--
I think you know now what I mean by left-
and right-hand sides--
for any variation.
So we have here another variation.
All such variations, however, must satisfy the displacement
boundary conditions, and those displacement boundary
conditions are given right down here.
Now notice that the particle, which I looked at earlier
already, has now moved to this configuration, or in this
variation of displacement.
Similarly, this particle has moved here.
Now the principle of virtual work states once again that
for any variation that is satisfying the displacement
boundary conditions, the left-hand side, the internal
virtual work, must be equal to the right-hand side, the
external virtual work.
And these are two variations that we looked at.
Of course, there are many, many more.
In fact, we could think of millions, billions such
variations.
And we will, in finite element analysis, however, discretize
the whole structure using finite elements, and then only
allow certain sets of variations.
Namely, those contained in the interpolation functions of the
finite elements.
Well, let us now look at one important point.
Namely that if we were to integrate the principle of
virtual work by parts, we would obtain the governing
differential equations of motion plus the force or
natural boundary conditions, just like an infinitesimal
displacement analysis.
However, the difference, of course, is that we have
referred all our static and kinematic variables now to
time t, to the current configuration of the body.
Let's look now at an example.
It's a simple example, in some sense, but it still
illustrates many of the basic points that I just mentioned.
The example is a truss stretching
under its own weight.
Here we have the truss in its original
configuration shown in black.
It has an original length L0.
It's subjected, of course, to gravity, g.
A typical cross-sectional area is given here, 0A.
Notice that this bar, under its own weight, moves down to
the configuration at time t here.
tA is this cross-sectional area now.
In other words, this area has now moved to that area.
The length of the truss, of course, has changed, and in
fact has changed to tL.
In order to analyze this problem, we have to make
certain assumptions.
The first assumption is, plane cross-sections remain plane.
Second assumption--
constant uni-axial stress on each cross-section.
And then if we make these two assumptions, then we have a
one-dimensional analysis.
Using these assumptions, we then obtain directly, from the
general principle of virtual work, the following.
Here we have, on the left-hand side, the general left-hand
side of the principle of virtual work, Cauchy stress,
the virtual strains, integrated over the current
volume at time t, which then becomes, on the right-hand
side here, the one stress only, because we have only one
stress in this bar, the vertical stress, times the
virtual strain corresponding to that motion of the bar,
namely, vertically downwards.
And this product is integrated over the current area at time
t, and of course, the longitudinal coordinate tdx.
Notice that here we are integrating over the current
volume, because we have a tA and a tL here.
The right-hand side becomes, in this particular case, an
integration over the current length of the gravity force,
which involves the mass density at time t, in other
words, the current mass density.
Notice the mass density, of course,
changes as the bar extends.
g is the gravity constant.
And here we have the virtual displacement, area at time t,
and of course the differential increment tdx.
Once again, we here also integrate over the
length Lt at time t.
Hence, if we now equate these two quantities, we directly
obtain what's shown here.
Now it's important to recognize that in this
equation, these here are the virtual strains corresponding
to these virtual displacements.
And once again, we apply the principle of
course, at time t.
So stresses at time t, meaning forces per unit area at time
t, area at time t, and so on.
This virtual strain, also used to describe this quantity, is
given right here.
Notice that here we have the partial of del u
with respect to tx.
tx, the coordinate at time t.
And this del u is exactly the del u that we see right there.
Now this principle must hold for any arbitrary variation in
displacement, any arbitrary virtual displacement, that
satisfy the displacement boundary conditions.
And the displacement boundary condition here, of course, is
simply that the displacement is 0 at the top of the bar.
If we use, now, integration by parts on this equation that I
just discussed with you, we obtain this equation here.
And that equation contains really two terms,
two parts to it.
First, an integral that is taken over the length L of the
bar, and then a quantity that corresponds to the length of
the position 0 and L, tL.
And all of the sum of these two must be equal to 0.
Since the variation del u is arbitrary, except, as I
mentioned, right at the top of the bar, we directly obtain
that this equation must hold, which is really extracting
this term here from the integral.
And this, of course, is the governing
differential equation.
And also, we must have that at the lower part of the bar,
this force or natural boundary condition must be satisfied.
Therefore, this is an example which shows, really, how the
principle of virtual work contains the basic
differential equations of equilibrium.
Here it is just one, because we have only one direction for
the motion.
And also the stress boundary conditions, the natural or
force stress boundary conditions.
The finite element application of the principle of virtual
work proceeds as follows.
Here, once again, virtual strains, virtual displacements
in the volume of the body and on the surface of the body.
These are the virtual quantities.
The actual applied forces, externally applied forces, are
these blue quantities.
We deliberately chose a different color to emphasize
that these are real forces, blue, applied to the body, and
these are virtual
displacements and virtual strains.
The stresses here at time t, of course, are unknown, and
those we want to calculate.
In the finite element method, what's shown in this big box,
we interpolate displacements, virtual displacements.
We calculate from those interpolations, strains, real
strains, virtual strains.
Via the stress-strain law, we got stresses, and we will see
in the later lectures that in fact, what we are arriving at
is a force vector tF and an external load vector.
Now, this tF here, shown in green, is a result of this
stress t tau ij.
And it's really the nodal point force vector that
corresponds to these internal element stresses, t tau ij.
This R here is a nodal point force vector that corresponds
to these externally applied forces.
This del UT, on left-hand side and on the right-hand side, is
a vector of virtual displacements, nodal point
displacements, and which, in the finite element analysis,
we then set equal to the identity matrix, invoking the
principle of virtual work in turn for each degree of
freedom, just like we do in linear elastic analysis.
This becomes the identity matrix, that becomes the
identity matrix, and the result is, of course, simply
that tF must be equal to tR.
The procedure that we are following in this box, of
course, is of much interest in non-linear analysis.
And what's happening in this box here, we will be talking
about in some considerable detail in the later lectures.
But for the moment, let's now assume that the solution at
time t is known.
Hence the stress at time t, the volume at time t, et
cetera, are known, and we want to obtain the solution at time
t plus delta t.
In other words, we want to move ahead 1 increment in time
delta t and establish the solution at that time.
The principle of virtual work applied now at time t plus
delta t in exactly the same way as we just applied it at
time t results in this set of equations.
A vector of nodal point forces that is equivalent to the
current internal element stresses at time t
plus delta t now.
And these are the nodal point externally applied forces at
time t plus delta t.
Well, of course, we don't know t plus delta tF, and for the
solution, we now say that t plus delta tF is equal to tF,
which we do know, plus a matrix tK--
we call this as tangent stiffness matrix--
times del U. If we set this left-hand side equal to t plus
delta tR, the externally applied loads, we directly
obtain this equation.
We are putting the unknowns on the left-hand side, and the
unknowns, of course, are the incremental displacements.
The tangent stiffness matrix is known.
It's the matrix corresponding to the
configuration at time t.
The load vector is known, and tF is assumed to be known.
Now, this is here an out-of-balance load vector,
and this set of equations then result into the solution of an
increment in nodal point displacements delta U. This
increment in nodal point displacement is added to the
displacement at time t to obtain the displacement at
time t plus delta t.
But--
and here is the big but--
by writing this equation up here, we had to
linearize the problem.
You see, we are dealing with a large amount of non-linearity
that can come from the kinematics, from the material
relationships and so on.
And we had to linearize.
We linearized here about the configuration at time t to
obtain this term, and this means that down here in this
equation, we only can have an approximation sign.
In other words, what we have calculated are not the exact
increments in nodal point displacement, but an
approximation to it.
More generally, then, we repeat that process.
We apply the same process in an iterative way to obtain a
very accurate solution.
And this is done in the following way.
On the right-hand side now, we have a load vector, the one
that we discussed already before, minus the nodal point
forces corresponding to time t plus delta t in
iteration i minus 1.
I should say at the beginning of iteration i, or end of
iteration i minus 1.
This gives us an out-of-balance load vector.
We're calculating an increment in displacements, delta Ui.
And that increment in displacement is added to the
displacements which were known already.
At the end of iteration i minus 1 corresponding to time
t plus delta t, we knew these displacements.
We add the increment that we calculated, and we obtain a
better guess, so to say, for the current displacement at
time t plus delta t.
Now of course in this iteration here, we need some
initial conditions, and the initial conditions
are given down here.
The t plus delta t F 0, in other words, when i is equal
to 1, then we need t plus delta t F 0 is given by tF.
The displacement t plus delta t U0, when i is equal to 1, is
given by tU.
These initial conditions are used in the iteration here,
and we notice that this set of equations in the first
iteration, when i is equal to 1, reduces to the equations
that I already discussed with you with the
previous view graph.
Well, if we look at this set of equations, we
recognize the following.
Nodal point equilibrium is satisfied when the equation t
plus delta tR minus t plus delta tF, i minus
1, is equal to 0.
When the externally applied loads are equal to the nodal
point forces corresponding to the stresses at time t plus
delta t, when these two vectors are equal, or rather,
when R minus F is equal to 0, then of course we have nodal
point equilibrium.
Compatibility is satisfied provided a compatible element
mesh is used.
The stress-strain law enters in the calculation of the
tangent stiffness matrix and the nodal point force vector
F, which corresponds to the current stresses at time t
plus delta t, iteration i minus 1.
Well, we will get back to these considerations just now.
But most important is to recognize that we have to
calculate F t plus delta ti minus 1 very accurately from
the displacements t plus delta t Ui minus 1.
The general procedure is depicted on this view graph.
Here we have the displacement at time t plus delta t ui
minus 1, which give us strains.
We can directly, via differentiations of the
displacements, get strains, which gives us stresses if we
enter here with the appropriate
constitutive relation.
Of course, we have to use the constitutive relation
corresponding to the problem that we want to solve.
We get stresses, and the stresses, then, give us this
nodal point force vector that I talked about.
Notice that in this relationship here, we have to
perform an integration, and that integration is expressed
here, where we obtain the stresses at time t plus delta
t in iteration i minus 1 by taking the stresses at time t
and integrating the stress-strain law times the
incremental strains from time t to time t plus delta t,
where we use, of course, the strains at time t, and the
strains at time t plus delta t
corresponding to that iteration.
This step is most important, because you have to calculate
the F vector correctly so that when you say that this
equation is satisfied, that indeed you have obtained the
right configuration for the finite element mesh that
you're considering.
So it is most important to calculate the F correctly.
Of course, the tangent stiffness matrix should also
be calculated appropriately, and we will talk about that as
well later on.
In all this discussion, I should point out we assumed
that the nodal point loads are independent of the structural
deformations.
In other words, they vary only as a function of time.
Here's a small problem that schematically
shows what we mean.
The load R is always acting vertically, never mind how
much this beam has displaced.
And that load variation is shown here on the graph as a
function of time.
We satisfy the basic requirements of mechanics if
we perform this finite element analysis appropriately.
We satisfy the stress-strain law because we are evaluating
the stresses correctly corresponding
from the given strains.
We also satisfy the compatibility requirements
provided we use compatible elements, and of course we
satisfy the displacement boundary conditions.
Equilibrium we also satisfy corresponding to the finite
element nodal point degrees of freedom.
But there is one important point.
This means that we only satisfy global
equilibrium in general.
We satisfy equilibrium locally only if we use a fine enough
finite element mesh.
In other words, I'd like to distinguish now between
satisfying global equilibrium and local equilibrium.
Local equilibrium means that we want to satisfy the
differential equations of equilibrium at every
point in the mesh.
And that is only achieved if we have a large enough finite
elements for the given problem--
in other words, if our discretization is fine enough.
A check on whether we are satisfying local equilibrium
well enough is to look at the stresses along the boundary
and see whether the stress boundary
conditions are satisfied.
That is one check.
Another check is to test whether along boundaries of
the elements, there are small stress jumps only.
In other words, whether the stresses from one element to
another element do not jump too much.
We call that the stress jumping because as we will
just talk further, when we calculate the stress at a node
for an element, we get, of course, one value.
And if at that same node, the stresses are calculated using
another element, you would get, in
general, another value.
This stress difference should not be very large, and of
course, if we have a homogeneous material, if we
have other conditions satisfied as well, in other
words, if typically the exact solution would have a
continuous stress variation, then the stress jump should be
very small for an accurate solution.
And this is a good check, an indication of whether we
indeed have obtained a valid solution.
Well, I like to now demonstrate to you some of
these considerations by showing you the results of two
example analyses.
But it is best if we show these examples on a different
reel, on a second reel of this lecture.
Thank you very much.
We just discussed some issues regarding how equilibrium is
satisfied in a displacement-based finite
element analysis.
We discussed briefly that global equilibrium is always
satisfied for any finite element mesh
that you have selected.
By global equilibrium, we mean equilibrium at the nodal
points of the finite element mesh.
However, locally, equilibrium is only satisfied if the mesh
is fine enough.
And an indication to see whether a mesh is fine enough
is given, for example, by whether these stressed
boundary conditions are well satisfied, and whether there
are no stress jumps between elements.
I would like now to share with you experiences relating to
two problems.
These are actually linear elastic analysis problems that
exemplify how we are looking at stress jumps in order to
choose an appropriate finite element mesh.
The first problem is the analysis of a
plate with a whole.
The plate has certain geometric dimensions, of
course, and it has also material data.
The material data are given right there.
Notice that we also give material data for an
elastoplastic analysis.
We will not perform now the elastoplastic analysis.
We will consider it later.
For the moment, we just need these two material data.
The| geometric data for the plate are given here.
Now, this plate with a hole has a
doubly symmetric structure.
In other words, there are two symmetry lines, as shown here.
So we only need to look at 1/4 of the plate in the analysis.
Notice that the plate is subjected to a pulling force
up here and similarly down here.
The purpose of the analysis is to accurately determine the
stresses in the plate, and of course, once again, assuming
that the loads are small enough to warrant a linear
elastic analysis.
The considerations that are valid in linear elastic
analysis regarding the choice of a mesh are of course also
valid in a non-linear analysis, and we more easily
can discuss them now, just using a linear elastic
analysis as an example.
Using symmetry, as I pointed out already, we can just focus
our attention on 1/4 of the plate.
Here is that quarter of the plate.
Here is the hole.
Notice we have roller boundary conditions here, roller
boundary conditions there, and the load is applied up there.
Accuracy considerations in displacement-based finite
element analysis are that first of all, of course,
compatibility must be satisfied.
And that is achieved by using a compatible
finite element mesh.
The material law must also be satisfied, and that is
achieved by just using the appropriate material data.
Of course, in non-linear analysis, this item needs a
lot more attention, but we are performing now a linear
elastic analysis, so this is a very simple condition, and
it's satisfied by simply selecting the proper Young's
modulus, Poisson's ratio for the model.
Equilibrium, once again, is locally only satisfied if we
choose a fine enough mesh.
Otherwise we only approximate equilibrium locally.
And as I pointed out earlier, we can observe the equilibrium
error by plotting stress discontinuities.
Here we have the results of a two-element
mesh for this problem.
Notice this quarter of the plate, idealized using just
two 8-node elements.
This is the line z equal to 0, and we will later on look at
stresses along this line.
This is a line y equal to z.
We will also look at stresses along this line.
But we will plot the stresses along this line here as a
function of y.
You should keep that in mind.
In other words, the first stress quality that we will
see will actually start somewhere here, because we are
plotting this stress as a function of y.
The displacements, somewhat amplified after load
application, are shown here.
Notice the maximum displacement here, uz, is
0.0285 millimeters.
The maximum stress is at this point here, 281 MPa,
megapascal.
Well, if we plot the stresses along the line z equal to 0,
and we're evaluating the stresses, I should also say,
at the nodal points of the elements, then we obtain these
values here.
One value there, one value there, one value there.
And once again, we are plotting the stresses along
the y-axis, measured along here, distance, and the stress
magnitude is measured along here.
Tau zz.
These are nodal point stresses.
And we simply have drawn a smooth curve through these
nodal point stresses.
The applied stress, far away from the hole, is 100 MPa.
If we plot the stresses along the line y equals z, along
this line here, we notice that we have a choice of choosing
either the nodal point stresses from the element
below it or from the element above it.
Well, the nodal point stresses from the element below it are
given by the crosses here.
The nodal point stresses from the element above it--
and notice, these are the nodal point stresses at these
nodes here are given by this point.
I should say, actually, that we are plotting here the
principal stress, the maximum principal stress.
We only look at one stress.
Of course, you could plot such curves for all three stress
components, but here we're looking just at sigma 1 being
the maximum principal stress.
Now there we can see the stress discontinuity.
We see a stress discontinuity right at that point.
Let's see once, what does this stress continuity result from?
Notice that there is a node here--
and I just like to put it in one color--
there's a node here coming from this element.
And there's another node coming in from that element.
Of course, this red point corresponds to this mark here.
So this node here gives us a stress in the stress graph
given by that symbol.
This blue point there is a node corresponding to the
element with the cross, and there we get also a stress
value, and that stress value corresponds to the cross mark
on the graph.
So let's go back to the graph once more.
Once again, same location, because this blue node lies
really on top of the red node.
They are put together.
And yet although we're looking at the same material particle,
we have this stress discontinuity.
And this stress discontinuity tells us that the mesh is
really not fine enough.
So we go to a finer mesh.
And here now we have selected a 64-element mesh.
Once again, 8-node elements, as shown here, for that
quarter of the plate with the hole.
Here is the typical 8-node element.
This is the undeformed mesh, and when we apply the loads,
we obtain this deformed shape.
The maximum stress is now 345 MPa, and the displacement
right at that point here is 0.0296 millimeter.
Well, we can now once again plot stresses.
We can plot stresses along this line--
or I should rather show it here, along this line, and
along that line, same way as we have done before.
And let's look at the stresses along the line z equal to 0.
Of course, y is the coordinate along which we are plotting
the stresses, and you can see now here a stress jump, a
stress discontinuity.
This is where two elements adjoin each other at one node.
Otherwise, we see very little stress discontinuity.
In fact, too small to be visible by eye.
If we plot stresses along the line y equal to z--
and here I should be more specific.
In fact, we only plot, once again, the
maximum principal stress.
We find that we get these curves here.
Now notice that at that location, we get four distinct
stress values.
The reason being that at this node, four elements couple
into that node, and that gives us four stress
values at that node.
If you go further away, we find that there are very few
or very small stress discontinuities.
So however, since we have still one stress discontinuity
there at that one location that I pointed out, we went
ahead and analyzed an even finer mesh, 288-element mesh.
This mesh is shown now here, 8-node elements, once again,
in the undeformed shape, and here in the deformed shape.
Notice that the maximum stress is now 337 megapascals and the
maximum displacement here is 0.0296 millimeters.
We can again plot stresses along this line and along that
diagonal line as a function of y.
And if we do so, we find that along the line z equal to 0,
we find for tau zz virtually no more stress discontinuity.
A nice, smooth stress curve.
And along the line y equal to z, for the maximum principal
stress, also almost no discontinuity.
You can see here a slight discontinuity, but that is
really a very small one that in engineering analysis, we
generally can live with.
Of course, these stress jump plots are very useful and very
helpful in identifying whether a mesh is
actually a good mesh.
However, we should realize that in order to identify
whether a mesh is really good, we would have to plot stress
jumps along many lines in the mesh, and that can, of course,
provide quite some cost.
It also means that we have to interpret a lot of data.
And one perhaps more convenient way to look at
whether a mesh is a good mesh is to plot pressure bands.
We have been lately using this procedure to identify whether
meshes are good, and I'd like to just share that latest
experience with you.
You plot bands of constant pressure, where the pressure
is defined by taking the mean of the stresses tau xx, tau
yy, and tau zz.
Of course, pressure meaning a negative sign in front here.
The two-element mesh gives us this pressure band plot.
Notice that these are bands of 5 MPa difference.
In other words, here you have a jump from this point to that
point of 5 MPa, and once again, from
here to there again.
5 MPa.
These are the pressure band plots for
the two-element mesh.
Here we have it as the same kind of plot for the
64-element mesh.
Let's look at this pressure band plot a
little bit more carefully.
Here we have, as you can see, again, differences of 5 MPa.
This is a non-dark area, unshaded area,
provides a band of 5 MPa.
That shaded area provides a band of 5 MPa.
Now, if we have a good mesh, then the pressure band should
be continuous.
In other words, this band here should be continuous, and
there should be no pressure band break.
Let's look at this one here.
It's nice and continuous, however, there's a break.
And so there is no continuity in stresses here, which is
quite simply displayed by the pressure band
picture shown here.
For a continuous stress field, we would have no breaks in the
pressure bands.
We have selected here 5 MPa.
Of course, you could also select 1 MPa, and that would
be then a much tighter criteria.
You could also select 10 MPa.
Of course, then that would be a very loose criterion, and 5
MPa for this particular problem, since you can see,
we're getting quite a number of bands along here, is a
reasonable way to proceed.
If you go to the 288-element mesh, the pressure band plots
look like this.
Beautifully smooth.
You can see here the smoothness in the pressure
bands, no discontinuity at all.
And this shows really that this mesh certainly satisfies
the criterion of smooth pressure bands, smooth
pressure variations, over the whole mesh.
And that is a strong indication that the stresses
are smooth, and therefore that the mesh is really an
acceptable mesh.
We see, therefore, that stress discontinuities are
represented by breaks in the pressure bands.
As the mesh is refined, of course, these breaks in the
pressure bands disappear, or certainly become smoother
until finally they disappear.
The stress state everywhere in the mesh is represented by one
picture, and that is, of course, the beauty of this
whole procedure.
If you plot stress jumps, stress component jumps, then
you have to plot for each component along all the lines
where elements join each other, these
components of stresses.
And there is a voluminous amount of information that is
being generated.
Here we have just one picture that basically gives us a lot
of information about the smoothness
of the stress situation.
Of course, the pressure band is plotted, maybe plotted by a
computer program.
Actual magnitudes of pressure are not given, but we are not
really interested in seeing those from these pictures.
What we would like to see are simply whether the pressure
bands are smooth, and whether, therefore, the mesh is an
acceptable mesh.
Let's summarize the results.
Here we now performed three analyses of this plate.
One two-element mesh was used, a 64-element mesh was used,
and a 288-element mesh was used.
The number of degrees of freedom are given here.
Notice that the number of degrees of freedom goes up
very rapidly as the number of elements goes up.
The relative cost is given in this column, and you can see
here that we have normalized the cost to the cost of the
analysis using the 64-element mesh.
So the two-element mesh is very cheap, 8% of the cost,
and the 288-element mesh solution is quite expensive.
The displacement at the top of the plate is given in this
column, and we notice that with a 64-element mesh, we
really obtain a very good displacement prediction.
In fact, there is no change going to a
larger number of elements.
And the stress concentration factor, interestingly enough,
increases first and then slightly decreases.
Notice this is not an infinite plate--
it's a finite plate, and therefore you do not get the
value of three, for example, that you would see in an
infinite plate, for a hole in an infinite plate.
So the two-element mesh cannot be used for stress predictions
is certainly a conclusion.
It's not accurate enough.
The 64-element mesh gives reasonable results for
stresses, and certainly very accurate results for the
displacements.
The 288-element mesh is probably overrefined for
linear elastic stress analysis.
However, for other types of analysis, for example, an
elastoplastic analysis, if you want to have very accurate
results, we might want to use, actually,
this fineness of mesh.
One question that one, of course might have is, why did
he not use 9-node elements?
We used 8-node elements.
Here we show now a mesh of 9-node elements.
Notice each of these elements now is a 9-node element.
Would the solution significantly improve if we
had used 9-node elements?
In other words, one 9-node element for each 8-node
element that was used before?
Well, the answer is no.
The solution does not improve very much.
If you look at this column here, let's look at the 64
8-node elements solution, in which case we had 416 degrees
of freedom.
With a 9-node element, we will have 544 degrees of freedom.
The displacements at the top are unchanged between the
8-node and 9-node element results.
The stress concentration factor is very
little effected here.
Only in the fourth digit, as a matter of fact.
And the stress jump and pressure band plot do not
change significantly, either.
Well, the second example that I'd then like to discuss with
you is the analysis of a plate with a crack.
We will once again look at different meshes for this
problem, and we will also discuss
one additional aspect--
namely that for certain types of answers that you're looking
for, you actually can use sometimes very coarse meshes.
Whereas if you're looking for the answer of a much more
difficult question, of course you will have to use a much
finer mesh.
Here we have a plate with a crack.
Notice the plate is subjected to a pull of 100 MPa right,
and at the left-hand side, same pull.
The crack is lying right there.
Notice it's not a a symmetric crack.
This is a crack tip A, that's the crack tip B. To model this
plate, we have to really model half of it.
We have to model half of it, because the crack is not a
symmetric crack.
The material data for this problem are given down here.
E, Young's modulus, nu, Poisson's ratio, Kc is the
fracture toughness.
And here we have the fact that the thickness of the plate
shall be just 1 centimeter.
It's a plane stress condition that we're looking at.
And the question that we want to address here is, will the
crack propagate?
Well, let us look at some background information.
Since we want to perform a linear elastic analysis and we
want to apply linear elastic fracture mechanic theory when
considering whether the crack will propagate in the plate,
we calculate a stress intensity factor, and we have
a mode 1 stress condition, so we calculate K1.
K1 determines the strength of the 1 over square root r
stress singularity at the crack tip, r of course being
the radius from the crack tip.
And as very well established in linear elastic fracture
mechanics, when K1 is greater than Kc,
the crack will propagate.
Kc is a property of the material.
The computation of K1 is performed by this here.
It's equal to the square root out of E times G. G is
obtained as the partial differentiation of pi with
respect to a, where pi is a total potential energy of the
structure, and a is the area of the crack surface.
G is known as the energy release rate for the crack.
This is quite well established, and what we want
to do, of course, in the finite element analysis, is to
calculate this G value.
The question is, how do we do that?
Well, in this finite element analysis, each crack is
represented by a node.
And if you look at a cut through the plate where this
is here the crack, we have this here as the node
at the crack tip.
This node in the finite element analysis is, so to
say, moved forward a differential amount, so that
the change in crack area is given via
this little bit here.
Now that the original crack tip location is here, the new
crack tip location is right there.
Then we can, in this process, evaluate del pi over del a as
del pi del l times 1/t, where t is a
thickness of that plate.
Of course, this here is only conceptually being done, that
we're actually moving the crack.
What we're really doing is the differentiation of del pi with
respect to l, where l is the length of the crack and del l
is the differentiation of l into this direction here.
The quantities del pi del l may be efficiently calculated
by an algorithm that we have just recently been publishing.
Of course, there are different algorithms available as well,
but we would like to use this procedure to calculate
data pi over l.
And you can look up how it's being done in
this particular paper.
Let's look at the results.
Finite element analyses have been performed using a
17-element mesh, and this is the mesh here.
We call this a coarse mesh.
Notice the crack is right there.
Notice that we have only two elements above the crack, to
the boundary of the plate, and we have here once again also
only two elements below it.
Notice that these nodes around the crack tip B have been
moved to the quarter point.
And similarly, they've been moved here to the quarter
point, meaning that this node is not at the midpoint of this
element side, but at the quarter point of
that element side.
And similarly for the other nodes around the crack tips.
If we look at the results, the stress results, maybe plus the
stresses on the line of symmetry, we find that tau yy
varies as shown here, then of course increases very rapidly
near the crack.
This is where the crack is.
Increases very rapidly, then drops down to very close to 0.
Of course, it should be 0 right in here.
And here we have, coming from the other side to crack tip A,
this very rapid variation in stresses.
These are the stresses tau yy.
If you look at the pressure band plot--
I introduced you to those kinds of plots
in the earlier example--
we find that the pressure jumps are really
larger than five MPa.
For example, you can see it right here, you can see it
right there, you can see it right there.
So it's a coarse mesh, and certainly the stresses are not
very well represented, not very well picked
up with this mesh.
However, if we look at the stress intensity factors that
we wanted to estimate in this analysis, because they will
tell us whether the crack will propagate or not, we find that
for the stress intensity factors, for example, KA, we
get an excellent solution.
For KB, the solution is not as good, but perhaps still
acceptable.
We now go to a finer mesh, a much finer mesh.
128 elements are used in this mesh.
And we can see once again here, 1/2 the plate.
Here is point A, crack tip A. Here is point B, crack tip B.
Notice the crack lies in between.
And we will be plotting stresses along this line.
This is a line of symmetry here, and we will plot tau yy,
which is the stress across here, as a function of z.
But let's first look a little bit closer at the mesh.
If we look at one detail closer to the crack tip A, we
see this is the mesh that we use, and we can go still
further and show the elements just around the crack tip,
right A. Here is A. Notice that we use triangular
elements, these are isoparametric degenerate
elements, and that we have moved, once again, the nodal
points from the midpoint to the quarter point.
Because we know that if we do so, we are generating in the
finite element the 1 over square root r singularity that
we would like to have in the analysis, because we know that
the stresses vary as a function of 1
over square root r.
So it's of benefit to shift the nodes to the quarter
points, as was done in this example.
Now we look here at the stresses.
Stresses tau yy, as I pointed out, is a function of z.
And you can see here a stress that very rapidly grows, then
on the other side comes down to practically 0 in the
cracked area, and on the other side, coming again now from
the right here, a stress that varies very rapidly up close
to tip A.
We also can see here that the stress variations
are nice and smooth.
Of course, that indicates to us that the mesh is really a
good mesh, but we can also look at our
pressure bands again.
We find that the pressure jumps are quite small if we go
not too close to the crack tip.
In other words, here we have some pressure band variation,
some pressure band breaks, but we have still quite smooth
bands that tell us that the mesh is quite fine.
If we go very close to the crack tip A, then we find that
the pressure jumps are, of course, large right here.
You can see here breaks in the pressure bend.
Right here, for example, a typical example, if you look
very closely here, there's a black band coming
into a white band.
So there's clearly a break in the pressure band.
So here, then, in this area, the stresses are not
continuous, and would also not be very accurately predicted.
However, we know that our scheme of calculating the
stress intensity factor works for even a coarser mesh, and
for this fine mesh, we've got indeed excellent results, as
you can see for KA and for KB.
Both are very well predicted.
The point of this analysis is really that the degree off
refinement needed in a mesh depends on what kind of
question you have, what kind of problem you're solving, and
what kind of question you're asking.
If you're only interested in predicting displacements
accurately, a coarse mesh may do.
If you want to calculate stress intensity factors, also
a coarse mesh may do.
If you want to calculate lowest natural frequencies,
associate mode shapes in a dynamic analysis, once again,
a coarse mesh might very well do.
However, if you want to calculate stresses very
accurately, then you may need a fine mesh.
And in general non-linear analysis, we need accurate
stresses, because the stresses influence the material
relationship that has to be used as we go through the
incremental solution.
And therefore, frequently we need, in non-linear analysis,
a finer mesh than what would do well in a linear analysis.
So we should keep that in mind, that in fact what might
be a sufficiently fine mesh in linear elastic analysis may
not be sufficiently fine in a non-linear elastoplastic
geometrically non-linear analysis.
Of course, how fine the mesh has to be in a general
non-linear analysis depends very much on many criteria,
criteria that we will still discuss in the next lectures
of this course.
But I hope this gave you a bit of an overview of what are
some important considerations in linear as well as
non-linear analysis.
Thank you for your attention.