Before we go further, let us summarize what all we did in the last class about isoparametric
formulation. Remember that we talked about isoparametric formulation as a very general
methodology for a variety of elements; the elements may be plane stress, plane strain,
solid elements and so on. What we are going to do is to quickly summarize all the steps
that were involved in the last class and then go ahead to evaluate the stiffness matrix
as well as the right hand side that is the load. If you remember that about two or three
classes or even before, we had derived the stiffness matrix terms. We said that the stiffness
matrix can be written in terms of B transpose DB dv.
We started, in the last class actually we started, there and we introduced the concept
of the natural coordinates and how mapping could be done. Now we transform this equation
into an equation of this form, where psi and eta if you remember, we introduced it as natural
coordinates. Please note that this J we are talking about, is a Jacobian about which we
will see in a minute and that this symbol indicates that it is a determinant of Jacobian.
There should not be any confusion. It is not a matrix and it is just a symbol for determinant
of the Jacobian. We wrote down Jacobian as dow x by dow psi and so on and then we went
ahead to derive the strain matrix. Have a look at the strain matrix. The strain matrix
has u comma x u comma y v comma x and v comma y. In order to transform this into a strain
matrix, we introduced a 1000 0001 0110 and so on
We defined the Jacobian such that u comma x and u comma y can be related to u comma
psi and u comma eta through a Jacobian. Just have a look at that and it will be clear that
we introduced the Jacobian to, I think there is a small mistake.
I will just get back here dow x by dow psi dow y by dow psi and sorry about that; I think
dow x dow eta and dow y by dow eta. From this, we derive that u comma x and u comma y could
be related to u comma psi and u comma eta through an inverse of Jacobian. We did that
in the last class and called that as gamma s.
Ultimately, when we stopped in the last class, we said that u comma x u comma y v comma x
and v comma y could be related to u comma psi u comma eta v comma psi and v comma eta
through a gamma matrix. So, we call this as a gamma matrix. What is the whole idea of
doing this? The whole idea of doing this is to replace the dependent of x in the B matrix
to the dependent of psi and eta. In other words B, a function of x and y should be replaced
in terms of psi and eta. That is what we are trying to do. We have come to a stage now,
where I can write out u comma x u comma y v comma x and v comma y in terms of my psi
and eta. What is that I am going to do here? At the next stage, I am going to plug that
into this particular matrix. You can now see this. So, u comma x, this matrix will be replaced
by this gamma matrix into u comma psi u comma eta matrix; this matrix, this vector. How
do I replace this? This is now replaced in terms of the N’s, because remember u is
equal to sigma Ni ui and so on. u comma psi u comma eta and so on will be replaced by
means of N’s and D’s.
This goes here, so, you will see that let us call this matrix as DN, part of it we called
it as DN, because we are completing it let us call that as DN. Then you would see that
this is gamma; gamma into DN into D is equal to this vector, which goes here. That complete
thing multiplied by this 1’s and zeroes would now give me B.
In other words, epsilon is equal to the 1’s and zeroes matrix multiplied by this gamma
matrix, what we just now saw, multiplied by this DN matrix multiplied by d. Those are
the three matrices that we would have in order to define d. We know that epsilon is expressed
as B into D. So, obviously B now becomes, all these terms put together, all these terms
put together now becomes B. Please note the difference between the DN, which I have defined
here and the same DN, name which I have given in the previous class.
May be if you want to avoid confusion you can call this as DN total, because the difference
is that it has four rows and previous one had two rows. May be, if you are bit confused
about it, you can give another name to it, but any way DN is nothing but dow N’s by
dow psi and dow eta; that is the matrix.
Having done this, what is the next step? Just substitute this B into my first equation here,
first equation here, B transpose DB. What is D? D is nothing but the stress strain matrix;
you know about that already. For plane stress, we know what it is? So, we substitute this
B into this equation. We know all the terms here.
In other words, we can write down that k as phi say psi comma eta t J d psi sorry d eta.
We can write down k like that. Our next step, our next step is to evaluate this particular
integral. This is a very happy situation, because the integral was surface integral,
volume integral which was bit difficult for you to derive, especially when the shape is
quite bad. What is that we have done? We have replaced that by means of psi eta coordinate
system and the coordinate varies now only from minus 1 to plus 1. What is the advantage
now? The advantage is that we can resort to what is called as numerical integration. We
are going to see what numerical integration means. I hope things are clear now. We are
now in a position to calculate stiffness matrix.
Please note again, a reminder that whatever be the shape, whatever be the quadrilateral,
we had just derived it for quadrilateral, whatever be the quadrilateral, even whether
the quadrilateral is regular or whether it is sorry, irregular, whatever be the shape,
the procedure is the same, because we map all of them into a lovely square in the psi
eta coordinate system. The procedure becomes quite well established and it is very easy
now to write down computer program. Let us look at quickly summarize what we are going
to do with respect to this numerical integration. What is numerical integration?
Let us look at an equation of the form, say integral minus 1 to plus 1 phi d psi. We can
evaluate this by different schemes, numerical integration schemes. What are these schemes?
These schemes are called as Gauss quadrature scheme and then Simpson's rule, for example
is one of the numerical integration schemes and so on. Usually in finite element analysis,
it is the Gauss quadrature rule that is extensively followed. There are other rules like Newton-Cotes
and so on. Thus Gauss quadrature rule is what is usually followed.
What does this rule say? The rule states that if you evaluate phi, what is phi? Phi is now
a function of psi. If you now evaluate phi, which is a function of psi at certain designated
positions and multiply it by weight and sum them up and sum them up then that result whatever
you have, would be very accurate depending upon how many positions of this designated
places you include. Let us look at the simplest one dimensional case. Let us say that we are
trying to evaluate that kind of curve. Let us say that we trying to integrate this between
minus 1 to plus 1.
In a one point quadrature rule, the position at which we are going to evaluate this is
given by that position where psi is equal to zero. Obviously, this is our I and that
is the psi; sorry that is the phi and that is the psi. The designated position here is
psi is equal to zero and what is it that you get here
and yeah, that is nothing but phi at psi is equal to zero. Let us call this as say phi
is equal to zero. If I now I ask what is the integral value of this curve by evaluating
that function phi at a position called psi is equal to zero, your immediate answer would
be that look that is only the area under the curve and so you can write this down as 2
phi zero. In other words what does it mean? It means that this two is the, what is called
as the weight function. This is one point quadrature rule. You can have more than one
points as the quadrature rule and let us see what the other rules say. Now, let us look
at the two-point quadrature rule.
Two-point quadrature rule says that you evaluate the same function at two positions, two positions;
no more at this centre, no more at this position, I said two positions. These two positions
are given by plus or minus 1 by root 3. So, you evaluate this function at two positions.
What do you mean by evaluate the function? Phi is a function of psi. What is that you
do? You substitute the value of psi to be plus 1 by root 3 and then again minus 1 by
root 3. So, you will get at two positions, the value of phi. Let us call the positions
as say phi1 and phi2.
The numerical integration rules states that this integral can be written as W1phi1 plus
W2phi2, where W1 and W2 are the weights and they happen to be 1 in this particular problem
or in this particular rule. So, this is called as two point quadrature rule. You have also
three point quadrature rules.
The three points are spread at three different locations and again these locations are given
in terms of psi and eta. The locations happen to be zero, plus and minus root of 0.6 plus
and minus of root of 0.6. In other words, the three point integration
rule states that, one point here, one point here and one point here. These are the three
points at which the value of phi is evaluated and the value of phi at psi is equal to zero
plus root of 0.6 minus root of 0.6. What are the weights here? The weight happens to be
5 by 9 for this position, 8 by 9 this position and of course 5 by 9 for this position. In
other words, now if I use three point quadrature rule, then this, say this integral say I,
now becomes W1 which is say 5 by 9 into phi1. Phi1 is the value of phi at minus root of
0.6 plus 8 by 9 into phi2. Phi2 is the value of phi at psi is equal to zero plus 5 by 9
into 5 by 9 into phi evaluated at root of 0.6. W1 phi1 plus W2 phi2 plus W3 phi3 is
what we get as a three point quadrature rule. You can keep on extending. In fact there are
tables that are available. I am not going into the details of it; you can have four
point quadrature rule and so on. I am not going to go into the details of it. Nevertheless,
it is important that the accuracy of it is well understood in the sense that whenever
you have root of 0.6 and you do not want to calculate it or you want to calculate and
put it, you do not want the computer to calculate it, then this root of 0.6 should be given
to very large number of significant digits, because that will affect our results.
The question now is what is the quadrature rule that you are going to use or you should
use? Before we come to that, let us now extend our concept to 2D aspect, because my K matrix,
which I had and which I transformed is now in 2D; may be we have to extend to it 3D as
well. Let us look at 2D, because if you understand 2D, how to do 2D, then 3D extension is quite
simple. In other words, what is that we are going to look at?
Instead of what i had instead so Now, K is minus 1 to plus 1, minus 1 to plus 1, some
phi t J d psi d eta and phi is a function of psi comma eta. What is phi, psi comma eta
is always confusion for students. Phi psi comma eta is nothing but B transpose DB calculated
at psi and eta.
Let us say that that is the element after the transformation into the natural coordinate
system and so on. What is the rule that we are going to use? Say for example, we use
a 2 by 2 rule. We use two integration point for one-dimensional case. Now, we are going
to use 2 by 2 rule, which means that I will have totally 4 points. How do I get four points?
By extending that psi plus or minus 1 by root 3 and eta plus or minus 1 by root 3. In other
words that, this line, this line and this line are such that, this is equal to plus
1 by root 3 minus 1 by root 3 and minus 1 by root 3 and so on. In other words, the coordinates
of this is minus 1 by root 3, minus 1 by root 3, 1 by root 3 minus 1 by root 3, 1 by root
3 1 by root 3 and so on, so, I get 4 points. Then how do I now evaluate this stiffness
matrix? That is evaluated as two sigma’s; two sigma’s ij both of them vary from 1
to 2. Phi calculated at these positions, four positions, let me call at phi i j w i.
Two sigma’s, let me call Wi Wj and phi calculated at psii etaj, i and j i and j varies from
1 to 2.
I will have now have 4 points or in other words, W1 into W1 phi into, phi calculated
at not into, phi calculated at psi1 eta1, then W1 W2 phi calculated at psi1 eta2; what
is psi1 eta1 and so on. These are points; that is psi1 eta1 which is minus 1 by root
3, minus 1 by root 3. What is psi1 eta2? That is this point and so on. You calculated these
four points, multiplied by the weights then or rather this point and so on. You get or
in other words, to calculate the psis at all the 4 points and the weights happen to be
for this problem 1 and so Wi into Wj into phis and these four points give you the value
of K. In other words, what is that we do? We calculate K at each of these points and
then sum them up. That is what we mean by stiffness matrix. Is that clear? May be you
can have a look at this whole thing again, what all we did, so this stiffness matrix
calculation becomes quite nice and easy.
There are two other things that we have to notice or we have to do. That is quite involved.
Let us see what are the other things that we have to be involved? I hope that I have
followed all the symbols what we did in the last class. Now, what are the other things
that we have to do? We have to calculate the right hand side. In the right hand side, we
have basically two types of loads, apart from the concentrated load; the body force and
the surface force. Body force is not a problem, because the procedure is exactly similar,
but surface forces creates certain troubles. What is the trouble that we are going to have
with the surface forces?
Let us look at this element. Let us say that that is the psi eta coordinate system. Let
us say that this element sits on the boundary of a component, say for example, that is the
boundary of the component and so on. In that component or in that element, we have defined
a surface force or surface traction, pressure whatever we call it. Let us say that we have
given a pressure of say phix and phiy in this particular surface. What is it that I have
to do? I have to now calculate the right hand side.
What is the right hand side? If you remember, the right hand side is N transpose phi ds
at s. In other words, we have to evaluate a surface integral which is N transpose phi
ds. This phi is given as phix and phiy. We know what is N? We have been through it for
quite some time and hence we know what N is. The whole idea now, is to calculate what N
transpose phi ds is. How are we going to calculate? Because, we have already we are into what
are called as natural coordinate system, we have already transformed the element into
a natural coordinate system and hence my job is to now transform this equation as well
to the natural coordinate system.
Or in other words, my first job is to recognize that the phase in which I have given the load
happens to be what? Eta is equal to minus 1 in the natural coordinate system. In other
words, symbolically writing it, if I write this as the natural coordinate system, no,
sorry; psi eta eta is equal to minus 1 and hence eta is equal to minus 1. That is the
phase in which I have given, I have given the load or in other words here is the phase
in which I have mentioned phix and phiy. Now, how am I going to evaluate this? In other
words, my equation now actually becomes, I am sorry, evaluating it as d psi; this implies
I should not say equivalent. But there are lots of things I have to do in order to get
to this.
What are the things that I have left out? I have left out quite a few things. They are
implied; no, they are not equal, because I have left out certain things deliberately.
What are the things that I have left out? I cannot straight away convert it like this.
Though actually I want to have equation of this sort, I cannot straight away do it. First
of all, let me convert this surface into line. Why is that we have to convert the surface
into line, because we are now talking only about a two dimensional situation, a plane
stress example.
My first step is not to write like this, though ultimately I want to do it. First step is
say that look integrate it say along the length of the element. Let me call this as zero to
l N transpose phi dl. That is my first step. What is N transpose phi dl? How do I convert
it? In other words, I can do that as minus 1 1 N transpose phi J into d psi. That is
exactly how I do it. What is J? Determinant of Jacobian or the Jacobian; in this case,
what is J and how do I calculate N? Usually these are questions which are confusing to
the students, the answer is very simple.
How do I calculate N now, because you may wonder N had psi and eta. How do I calculate
N now for one dimensional case, from my previous N’s which I have? How do I calculate that?
The answer is very simple. In the previous N’s, what are the N’s you had? 1 by 4
into 1 minus psi into 1 minus eta plus that is for N1; N2 you had 1 by 4 into 1 plus psi
1 minus eta and so on. In those N’s just substitute eta is equal to minus 1. In other
words, what we mean is this N or N transpose is calculated at eta equal to minus 1. So,
that problem is solved.
When I substitute eta is equal to minus 1 by N, purely become a function of psi. No
problem, so, straight away I can calculate. Phi is a value, phix and phiy that is in other
words, the x and y direction they are given. That is not an issue, because let us for the
moment assume that they are constants or else we may have to do something else. Let us not
confuse the issue. Let us say that they are constants that are acting in the x and the
y direction or else we have to interpolate them as well. Let us not worry about that
or in other words phix is constant through out this phase, phiy is constant through out
that phase. They are just written there in terms of numbers.
The first term is known. Second term is a vector, it has two entries or two rows, one
column and next is J. Once I do that, once I calculate the J, my calculation of this
integral becomes very simple, because I can use numerical integration. Question is, now
what is that J? Let us take a pause and see how to calculate that J matrix? Yeah, so have
a look at this and see for a minute as to how we can now calculate J? In other words,
we have to very carefully define what J is, the Jacobian. After all what is that we have
done here in this equation? We have replaced dl by J d psi. So, what should be J? J should
be dl by d psi.
We have to find out now dl and write it in terms of d psi. That is all we have to do.
Once we do that, in other words, once we calculate dl, write down dl, then it is very simple
for us to calculate J. What is dl or how do you calculate dl, is the question.
dl can be written as, let us look at that figure. Have a look at this figure and see
how we can write down. Have a look that figure and see how we can write down dl. We can write
down dl in terms of dx by d psi whole squared or dow x by dow psi whole squared plus dow
y by dow psi whole squared whole power half d psi. How do I get dow x by dow psi and dow
y by dow psi? That is very simple because, I know that this can be written in terms of
N’s.
When I now substitute dow x by dow psi and dow y by dow psi, where x is equal to Ni xi
and then simplify this, I get that dl is equal to l by 2 d psi. I am not going to do that
exercise; that is very straight forward, just substitute it and you will see that dl and
d psi are related by l by 2 and now what is this l by 2? l by 2 is nothing but the Jacobian.
With that result, we can write down our previous expression for N transpose phi J d psi as
minus 1 to plus 1 N transpose, phi remains the same; now, J becomes l by 2 and we get
d psi.
What is l by 2? What is l? l is the length of that element, length of this particular
edge we are talking about. So that is the l and hence this expression is useful to determine
re, or due to pressure. Please note that phi is actually a vector, phix and phiy. Correspondingly
this is also a vector and it has two entries. I need not tell you what N has? We can very
easily find out, because I have already given you a clue that this N has to be calculated
with eta is equal to minus 1. In other words, what is essentially or what is the crux of
our approach?
Our approach has been that we not only calculated the stiffness matrix by transforming it to
a natural coordinate system and then resorting to what we call as the numerical integration.
We also calculated or we can also calculate the other re’s or the r’s or the right
hand side, to be more precise or in other words the load due to the surface traction
or pressure can also be calculated by means of this natural coordinate system by transforming
the equation. What is this equation, how did we get? Please note this equation we got it
from our virtual work principle. By transforming the equation which we got from the virtual
work principle on to say the natural systems or natural coordinate system, where psi varies
from minus 1 to plus 1. That in a nut shell, shows that all the aspects whether it is the
stiffness matrix or whether it is the force terms, all these things can be calculated
in a very consistent fashion. There is always a tendency to replace this phi by means of
dividing that or multiplying this by total length and then dividing it by 2 and putting
it at 1 and 2.
One of the things we should notice, there are two things I should not say one; two things
that you have to notice. One is that, you should always have a look at how phi is defined.
Please note that in this particular definition we have assumed that phi is defined per unit
length.
If phi is not defined per unit length, then you will have to now look at t, thickness
as well so that, you convert it into unit length. Why is that we are emphasizing on
unit length is because we are now in a plane. So, it is important that you look at how phi
is defined. If phi is defined anywhere else, then your units on the left hand side and
the right hand side should be consistent. Say for example, we have defined a t in the
left hand side for K. Then you should also have a t, if it so happens that you have defined
phi as well as t. In that case, you have to multiply this by t as well. Please note another
This is a very simple thing that people make a mistake; please note how phi is defined,
so that the units what has been given there is consistent and thickness is taken care
of properly by the units that are given for phi. That is a very important thing and the
other one is that this is only an example for eta is equal to minus 1. It so happens
that the same element may be in the other end, at the other side of the boundary and
so that we may have to apply pressure not at eta is equal to minus 1, but say psi is
equal to 1 and so on. You have to be properly accounting for the surface where we give this
pressure. What is the surface that we are interested in, where is that we are giving
pressure and accordingly you have to do this derivation.
Thirdly if the pressure does not happen to be constant, if it varies, say varies, from
one end to the other, then you have to introduce linear interpolation even for phi’s as well.
That makes the problem slightly more complicated, but not difficult. It is very straightforward.
Only thing is that in that case, you have to now introduce interpolation for the phi’s
as well. There is always a tendency to lump these loads. We will see later in the course
that lumping may be very dangerous; it may not be consistent with the result we have
got.
In fact these loads which we have calculated by such a very rigorous procedure of integration
are called as consistent and nodal forces. They are consistent, they are consistent with
our other definitions. So, they are called consistent nodal forces. You have to be very
careful not to miss out this and just to divide say for example length is equal some phi and
if the load say phiy acting is say 5 per unit length, so you cannot multiply 5 into 5 which
is 25 and then divide it by 2 is 12.5, 12.5 and so on. Yes, it will work for linear elements,
but for other higher order elements this kind of procedure, will not work. Even when you
write a program, even if it works for some case, it does not mean that you should not
do this integration. You should go ahead and do this kind of integration and if you want
to write a code, then this is what should be implemented and that is very important.
A very detailed account of how to write the code for this has already been given by Owen;
I think Hinton and Owen, in a book called Finite Element Programming, where you will
see that a very detailed exposition of how this code can be written is given. Anyway,
for our course it may not be right now necessary to go into the details of how programming
is done, but nevertheless, it is nice to know that these things are important. What is the
second issue? We have now looked at what is called as a quadrilateral element. What are
the other types of elements which are very popular and which are important?
Triangular elements; triangular elements become very important. Why is that triangular elements
are important? Because, most of the mesh generation codes or mesh generation packages have triangular
or tetrahedron …. three dimensional counter part as a mesh generation facility. Today
there are lot of mesh generators, which can do quadrilateral elements. So, there is no
problem, but hardly there are any mesh generators which can do what we call as a hexahedron
element and hence it becomes important that we know how to handle, how to handle triangular
elements.
Basically there are two ways in which triangular elements stiffness matrix are determined.
But, by the way, we already noted one thing important about triangular element. If you
remember that, we said that triangular elements are constant strain triangular elements or
in other words we said that the strain happens to be a constant in a triangle. We illustrated
this with a very simple example. We said that, what is that we said? We said that let us
say that we call these three nodes as 1 2 and 3. What did we say? We said displacement
can be a0 plus a1x and a2y can be written like this.
We said, if we want to calculate dow u by dow x, what is that we get? a1. What does
it mean? It means that, when we express the displacement like that as an interpolated
value, then dow u by dow x becomes a constant. It does not matter whether we replace these
a’s by the degree of freedoms at the three nodes 1 2 and 3; does not matter, because
the basis is only a linear interpolation formula like that and hence dow u by dow x becomes
a1 and so on. So, the strain within that element becomes constant.
What is the trouble with this? We have already seen this with an example. We saw that because
of this, we have to have lot of elements, especially when there is bending that is involved.
There are lots of elements that we have to use in order that the bending exercise or
when an element or when the component bends, bending is properly taken care of. I am reemphasizing
it, basically because I have found lot of mistakes during modeling. Why is that there
are lot of mistakes? Because it is very easy to do what we call as free meshing. Take a
three dimensional component, mesh it with a tetrahedron and just go ahead and finish
the problem. But, if you do not do proper mesh sensitivity, as we did on the other day,
then you are going to be in trouble.
With all these drawbacks, still we have to learn about triangular element and it is three
dimensional generalization which is the tetrahedron element. We will concentrate only the triangular
element, because what all we are talking about in the triangular element is very similar
to the tetrahedron elements. There are two procedures, the two ways in which you can
calculate the stiffness matrix. One, procedure number one, is to collapse two nodes in a
quadrilateral element. What do we mean by that?
Let us say that I have an element like this, simple element, just illustrate how we collapse.
Let me call this as 1 2 3 and 4, anticlockwise sets. Fine; now I can develop a triangular
element following exactly the same procedure, whatever I have done so far, exactly the same
procedure, by collapsing 1 and 2 and which results in an element something like this.
Whatever be the shape, it does not matter. What do I do? I collapse x1 to x2.
What do you mean by collapsing? I simply mean that I give the coordinates of 1 and 2 to
be the same, x1 is equal to x2, y1 is equal to y2. Where do I do that? I do that in my
expression for x and y, first step. What is that I mean by collapsing? By substituting
that x1 is equal to x2, y1 is equal to y2 in these two equations; is that all? No, because
we are using isoparametric formulations and hence we have similar values for u and v as
well; for u and v as well. So, u and v are again written and this collapsing exercise
is done. So, I can say that it is 1 comma 2 3 and 4. The collapsing exercise is done
there as well, so, u1 is now the same as u2. In other words, when I do that what essentially
I get as degrees of freedom? Simple; u2 v2 u3 v3 and u4 v4. These are the degrees of
freedom that I get. One way of calculating the stiffness matrix is by doing this collapsing
business. This is very simple and I will do that in the next class by taking an example,
simple example and following all the steps that I have done in this class and the previous
class.
The other way of doing it, there is another way of doing it, is by following what is called
as area coordinate systems. The area coordinate system is some sort of a generalization or
an extension of the natural coordinate system, which we have in place right now. So, the
second way of doing it is by area coordinate system. Whatever we have been talking about
so far has been for linear element. But that is not sufficient. I said that the interpolation
values can be of a higher order polynomial. It can be quadratic, can be cubic and so on.
Whatever we have been taking about, whether it is the quadrilateral element or whether
it is a triangular element, it is possible to raise the polynomial level or polynomial
degree, from a linear exercise to quadratic, cubic and so on.
Our exercise, next exercise is to complete the triangular element by collapsing and get
into this area coordinates and then extend this whole exercise to a higher order element
which is not very difficult. I am not going into each of the element. For example, continuum
elements of 3D, I am not going to the details, because it is exactly similar. I leave that
out as an exercise. We will meet in the next class and see how we tackle these triangular
elements.