We had started our discussions on an alternate formulation for stiffness matrix in a triangular
element. Remember, we had already emphasized the importance of triangular element and hence
it is important that we understand what we are in for, when we do triangular elements.
Unfortunately in the industry, again I am just emphasizing the importance; unfortunately
in the industry if you really look at how people do finite element analysis, most of
them do not have the time to put a very proper mesh. They do not have the time to put, for
example, a hexahedron element. What they usually tend to do is to rely on this kind of what
they call as automatic mesh generation. In fact, I have seen from my experience about
70 to 80% of the people, why it may be more than that also, use only automatic mesh generation.
That is what is given in the software, but there are lots of pitfalls in it. It is nice,
it is good to use it, because you can finish your job very fast, but there are lots of
pitfalls to it and that is the reason why we are going through this kind of derivations.
Let us now look at, concentrate on this triangular element. We are now going to develop or rather
look at what we call as a simple coordinate system. What is the coordinate system that
we are going to give or on alternate coordinate system?
Let us say that I am looking at a point P. Let me join. P can be at any point whether
here, here, here or whatever it is. Let me join the three vertices of the triangle and
call as A, the area of this triangle and A1 and A1 to be the area of that piece that is
area of the triangle P23 and so on. This is called as A1, this is called as A2 and this
is called as A3. Is that clear? So, A1, A2 and A3. Let me say that any point P, at this
point or any point P, can be expressed by means of a coordinate system, say, psii; psii,
which is equal to Ai by A. Is that clear? That means that this point P can expressed
by means of coordinate systems psi1 psi2 and psi3.
Yes, I can see in your face. Immediately there is a question that we are in a plane and we
need only two coordinate systems, independent; the basis all of you know some linear algebra.
The bases in a plane are only two, what we call as i and j are what are required. Now,
you are talking about three coordinates. Where is the lacuna? Yeah, I can understand. Immediately
that is the question that comes to your mind. Exactly, it is correct. In other words, psi1
psi2 and psi3 which are A1 by A, A2 by A, A3 by A are not independent or in other words,
psi1 plus psi2 plus psi3 is equal to 1.
In other words, this brings us to an important concept that at psi1 psi2 psi3 are related
and that they can be used as a first set of equation to solve or to express x and y in
terms of psi1 psi2 and psi3, to express in terms of psi1 psi2 and psi3. In other words,
I can write down, say for example, x to be psi1 into x1 plus psi2 x2 plus psi3 x3. My
whole idea immediately, remember what we did? We tried to develop a mapping between a new
coordinate system which we had put few days back and x and y. I am going to develop the
same thing. I am going to say that, let us say that x is written as psi1x1 plus psi2x2
and plus psi3x3. Similarly, I can write down y to be psi1y1 plus psi2y2 plus psi3y3, so
that, taking these three equations, I can write them in the matrix form and say 1 x
y, 1 x y, is equal to psi1 psi2 psi3.
What are these things? Simple. 1 1 1; that is the first equation. Second equation is
x1 x2 and x3. Third is y1 y2 and y3. Let me call this matrix, say A matrix. If you like
another name you can give it, but be careful, you should not confuse with other matrices
that you know.
So, that psi1 psi2 psi3 now can be expressed as A inverse. A inverse can
actually be expressed in a very explicit fashion. A inverse can be expressed in terms of an
explicit fashion. Let us see how we express A inverse, let me do that.
I mean, just for, is equal to 1 by 2A into x2y3 minus x3y2, x3y2, y23 and x32. That is
the first row. x3y1 minus x1y3; what is y23? y23 is y2 minus y3; y23 is y2 minus y3. People
who know it by heart, this is not a problem for, if you know how to do an inverse of a
2 by a 3 by 3 matrix. y31 and x13; this is available in any standard text book. Next
is x1y2 minus x2y1 y12 and x21. Not very difficult to remember, because 23 32 31 13 12 21 and
similarly 23 32 31 13 12 21; 12 21 12 21. Now, what is that we have achieved? The first
step is that we have achieved or we have determined or we have put down, one what we call as a
coordinate system. I am not happy with only this, because I want to convert this into
a natural coordinate system and follow the procedures which I know. In other words, I
do not what to stick to 1 2 3, but can I convert it into a form which is very similar to the
form which I had for my quadrilateral? How do I do? How am I going to do it?
Simple; what I am going to do is, I am going to say that let me call psi1 is equal to a
psi coordinate system; let me call psi1 is equal to psi. Just wait a minute, till I write
that and then psi2 is equal to eta, so that psi3 becomes 1 minus psi minus eta. Let me
write like this, so that I still stick to my psi eta coordinate system. I am defining,
what is that I am doing? I am defining psi and eta to be like this. Psi1 is equal to
psi, psi2 is equal to eta, so that psi3 automatically becomes 1 minus psi minus eta. That is my
first step, I have defined. What all I have done, I am very simple; whatever I have done
is only took a what is called as a area coordinate system and defined my natural coordinates
psi and eta in terms of this area coordinate system. But, my job actually begins there.
What is my job?
Let us trace, let us keep always in mind what we are trying to do, so that you do not get
confused. What is that we are doing in isoparametric element? I am recapitulating for the nth time,
but does not matter, because this has to get into your mind when you follow this derivation.
First is, I put one coordinate system, one natural coordinate system. This is the first
step I am doing. The second step is what? I am going to do, I am going to write, your
determinant, it comes later you know. I think, what I should do is, at the end of the class,
what I am going to do is I am going to give or I am going to ask you to give me a framework
to code this whole thing. That is why I told you in one of the classes that coding is important,
because you will understand all the steps very clearly.
The next step is determination or writing down the shape functions. That is important.
I cannot use obviously my 1 by 4 into 1 minus eta into 1minus psi and those shape functions,
which I had used in the earlier classes, I cannot use it. Now, I have to write down shape
functions for this. Once I write shape functions, then I determine what is called as the Jacobian
through the determination of DN. Is it clear? Then, I determine the DN value, then I determine
Jacobian, then I determine B matrix and then I determine stiffness matrix. So, that is
how we proceed. This is common whatever be or however I do it, these steps proceed. In
fact, what I am going to do is to illustrate this for a six noded triangular element, so
that you will understand that the procedure is exactly the same.
For a linear triangular element of what we have given, I can straight away use this to
be N1, this to be N2 and this to be N3. N1 N2 and N3, shape functions, but let us now
move; that is not very difficult and you can work it out yourself. But, let us move to
a slightly more difficult or not difficult actually, but a different type of element.
So far, we have been looking at what are called as linear element or first order element.
Whether I put these three nodes in this element or the four nodes that I had put for a quadrilateral
element, the quadrilateral element becomes what is called as a bilinear element. Here
it is a linear element. However I put it, it is very obvious that we are looking at
only linear interpolations. On the other hand, let us look at an element which looks like
this.
A 6 noded triangular element, 6 noded triangular element, whose nodes, note this carefully,
how I am defining 1 2 3 4 5 and 6, whose nodes I am defining it as 1 2 3 4 5 6. What is the
peculiarity of this element? Though I have written it as straight sided, in the x y coordinate
system can someone guess, first of all how many nodes are there? 6 nodes are there. So,
I have to necessarily go above my linear interpolation. What is that we are doing in isoparametric
element? Shapes are also linearly interpolated for that previous element, because x and y,
remember I had written it in terms of psi1 and psi2. In fact, u and v also I am going
to write it in the same fashion. x is equal to N1x1 plus N2x2 plus N3x3. That is what
I will write.
In this case also, both x as well as u, or both x and y as well as u and v, I will now
write down by a polynomial which is more than, the orders of which is more than, linear.
It would become in this case, quadratic. What does this physically mean? What does this
mean? It means that in my x y coordinate system, in my x y coordinate system, this triangle,
this triangle what I had put can take what shape? It takes a shape something like this;
something like this, because these chaps are defined by quadratic interpolation. What is
the quadratic interpolation that we are going to use or in other words, see, please note
that I am going to use the same shape functions, sorry, same coordinate system psi and eta
or psi1 psi2 and psi3, either way. You can use psi1 psi2 and psi3 or psi and eta; the
relationships are already given. How am I going to express the, what I call as the shape
functions, because I am interested in the shape functions to express my shape of the
element as well as u and v.
I write this as, say for example, how many shape functions are going to be there? Obviously
6; so, I am going to write this as N1 is equal to say psi1 into 2psi1 minus 1. Just look
at that for a moment and tell me whether it is okay with you. Just check whether it is
okay with you, in the sense that does it, I mean I have just put it down by inspection.
By experience you can write it down like that. You can see that, what is that I should get?
I should have N1 to be equal to 1 at node 1 and at all other nodes, it should be equal
to zero. It happens quite nicely there and N2 is now given by psi2 into 2psi2 minus 1.
N3 is equal to psi3 into 2psi3 minus 1 and N4 is equal to 4psi1 psi2 and so on. N5 is
similarly 4psi2 psi3 and N6 is equal to 4psi3 psi1 and so on. What is that I have done?
I have essentially put down, what I would call as a natural coordinate system by means
of areas and express them either as psi eta or you can keep it as psi1 psi2 psi3. Either
way you can substitute, in terms of psi and eta in this as well and written down for a
six noded quadrilateral a shape function, which is, sorry, six noded triangle, a shape
function which happens to be quadratic. Is that clear, a shape function which is quadratic.
What does this mean in terms of or what is the improvement I have done? What is the improvement
I have done here? Pardon? Strain; beautiful. One of my major drawbacks in my previous,
say, derivation or our previous discussion, one of the things which we said is that that
triangle which I had put down had a problem basically because it was a constant strain
triangle. What will happen to this? What will happen to this? Now, it is quadratic. You
need not even workout equations; you can have a feel and just tell me what will happen?
The displacements are now expressed in terms of psi squared and so on, psii squared and
so on. So, it will become linear; beautiful.
If we call the previous triangle which we had as constant strain triangles, CST as people
call it, then this triangle, we call this as a linear strain triangle or LST. Is it
clear? What is my next step? What is my next step? Once I have written down shape functions,
what is that I am going to do with these shape functions? Yeah, DN; correct. Why not you
try and write down DN for this? Next step is, because I am moving towards Jacobian,
moving towards Jacobian. I have to calculate or I have to write down DN. By the way what
would be the size of DN? I mean, I am talking about the complete DN. 6 by, 6 by 6; 6 by
2, you can write down the coordinates as x1y1 x2y2 x3y3 x4y4 or if you want to write down
as x1y1 x2y2 completely, either way you can write it down, but you have to be careful
as to how you write down. After all Jacobian is given by dow x by dow psi and so on. That
is the way you write down for this as well.
We can calculate DN, put down a DN and then you can calculate. You can write it down in
terms of psi and eta as well, then you can calculate DN. I am not going to repeat all
that, because I am sure all of you know it and once I calculate DN, then what I am going
to do is to calculate J, J inverse and then substitute it into my stiffness matrix expression,
from which I can calculate stiffness matrix. Till that point, I leave it to you. I am sure
you can all work it out. But, I am going to stop there at that point and look at the stiffness
matrix. How am I going to estimate the stiffness matrix?
What is the peculiarity in estimating the stiffness matrix? I cannot use just like that
minus 1 to plus 1, minus 1 to plus 1 I have to keep them as area coordinates. The type
of integration that I am going to use is different there. It is to be area coordinates. I am
not going to again go into the details of it, lot of text books have this information.
That is in other words, the exact say Gauss point quadrature and other things that are
used, I am not going into the details of this for the triangular element, because it is
a big chart. Any standard text books, like for example R.D.Cook, he has listed down what
are the Gauss points and the corresponding weightages and so on. I would like you to
look at that book for weightages and so on.
But, the only thing I want to mention is that if it is a phi dA, then this can be written
as sigma Wi Ji phii, where the values of the W’s, weights and the corresponding positions
what are called as the Gauss point positions, the point positions, can be obtained from
these tables. They are not just multipli….; it is not 2 by 2 so on. They are, I mean,
they are quite involved. I am going to leave that because that is going to consume lot
of time to look at them. But, it is anyway available in standard text books or standard
even tables that are available for it; you can work it out from that.
In a nut shell, what it really means is that it is possible to express or possible to calculate
the stiffness matrices of any type of element through the same procedure. What does that
mean? It means that, whether I am going to, what I looked at as a constant, sorry, linear
strain triangle or higher order elements, then also it will be possible for me to express
or to calculate stiffness matrix in the same fashion. Just as a small exercise, let us
now go into what we call as a higher order element for a four sided element. Let us look
at this for a four sided element.
In other words, let us look at an element which looks something like this. Let us now
look at a higher order element of this shape. How many nodes will there be? Why 8? I specifically
asked this question, basically because all of us have a tendency to say initially, that
there are nodes only on the edges. Pardon? Yeah, you are prompting what I am going to
teach before hand. But nevertheless, this is 1 2 3 4, which we have, which we know.
Why I ask this question though he prompted before hand a type of element that I am going
to teach, but what is important or what should have come from you is a different question.
What is it? I am going to give you a minute to think about it.
How did we arrive at the shape functions for this element, 4 noded elements? Forget about
these curves, just 4 noded elements. Just look at that element and how did we arrive
at the shape functions for this? There are two ways in which if you remember, we arrived
at that. There are two, in fact three ways. One of the very mathematical way or mathematically
correct way, I do not mean others are wrong, but rigorous way which we arrived at is by
Lagrange's interpolation. In fact, we interpolated it from a one dimensional case we interpolated
it in two dimensions. Remember that, we interpolated it in 1D and 2D. This is for a linear; go
back and look at your Lagrange's interpolation. Suppose I want to go to quadratic, next step.
How many nodes should be there? How many nodes should be there?
Three nodes. Now, three nodes, if I put three nodes like this on one side, is this complete
and can I write down Lagrange's interpolation from this? Because, what did we do? What did
we do in the linear case?
We wrote, in this case we wrote interpolation in one side, we wrote interpolation in the
other side and then we combined them. No, no, no, no, wait a minute, wait a minute.
I want, see, I want linear in this direction and linear in this direction. That is how
we got in this fashion. Here, if I have to do it, it has to be quadratic throughout if
you say from a Lagrangian interpolation point of view, then actually I should have one more
node in the center; one more node in the center, then only there will be, you know, quadratic
in all the lines that I take.
That is why we call this as a Lagrange element, a 9 noded Lagrange element, whose interpolation
function can be determined from Lagrangian formula, Lagrangian interpolation formula.
The first thing that you should keep or I mean take away from your mind is that there
should be elements only on the side, there can be elements also, sorry, nodes only on
the sides of the elements. There can be nodes also inside and such an element becomes Lagrange
element. It is a very good element, you know. It behaves very nicely and so on but only
thing is cost becomes very expensive when you keep using more nodes and so on. So, there
is a variation to Lagrange element, which are called as serendipity elements.
The difference is that I do not have one node at the center. There is no node at the center
of a serendipity element and many times the shape function for a serendipity element is
written by inspection; by inspection, by intuition. Is it clear? What is the difference between
a Lagrange element and a serendipity element? A Lagrange element has a node at the center
and its interpolation functions can be determined straight away from the Lagrangian interpolation
function in a similar fashion as we did for the one dimensional case. Serendipity elements
do not have a node at the center and the shape functions are written by inspection. Is that
clear? That is the first thing.
I mean, of course, it is not very difficult to write it. May be for a minute you can see
how do you write it, or I will give you a minute to think as to how whether you can
give a procedure. First of all can you give me a procedure as to how I can write the shape
functions for such an element or in other words, the clue is, can we proceed from the
shape functions of the 4 noded elements? Can we proceed from the shape functions of the
4 noded elements? Let me give you, why this serendipity elements become important later,
but just look at it and see whether it is possible to proceed from the shape functions
of the 4 noded element? Just look at that element, I will give you a second to see how
this can be done. I just want a procedure from you.
Yeah, the first thing is that let us see how this shape function looks like, looks like
for a 4 noded element. See, that can be done something like this. Let us say that that
is an element. I just make it lie down and just plot the shape functions. How will it
look like? Something like this; this is a, do you understand? It is a N1 would be a linear
surface, such that what is this value? That is this value. 1 and shape functions are zero
at all other places. Now, if I want to adopt this, what should I do? There are two things
I should do. I have put down nodes at these places. I hope you understand that this is
a surface like this, sloping surface. In fact, I think it would be better if I put it like
this may be, the picture becomes better, sorry, something like that. What is that I should
do now? I should remember that there are nodes at these places and make these shape function
go to zero at these places. That is my first thing. How do I do it? Usually how it is done,
is something like this. I will put down a shape function, say for example, 5 6 7 and
8. Let me call them as 5 6 7 and 8. Let me call them as 5 6 7 and 8. How does my shape
function N5 look like in the same mode?
It looks something like this. That is the one; that is the shape function for this and
it would be something like this at the center and goes to zero at other nodes. It is quadratic
in one direction for this and linear and goes like that in the other direction or in other
words, N5 can be written as half into 1 minus psi squared into 1 minus eta. Half into 1
minus psi squared into 1 minus eta is a good candidate for shape function for N5. Is it
clear? When I put down this candidate, just check whether it satisfies all my requirements.
N5, see that they are zeros at all positions. Do they become zeros at all other nodes? Obviously.
What is for N8? What is for N8? Similar thing for N8, how do I write down for N8? This is
the psi eta coordinate system. I am going to come to that again. I think, I forgot one
of the things that I said psi eta coordinate system and how do I write this?
1 minus eta squared into 1 minus psi; 1 minus psi. Before I proceed, I think that one of
the comments which I left out is that please note psi and eta coordinate itself is a curved
coordinate system. They are not straight lines, obviously they are curved coordinates, because
after all psi and eta we chose as one of the curves in the element and hence psi and eta
itself are curved coordinate system. Is that clear? What I did was, yes, you can say how
did I write? I wrote it out of a considerable experience; looking at it, saying that this
can be quadratic, so 1 minus eta psi squared. So, by lot of things, may be you would have
taken one hour to write it, I just like that wrote, because I knew what the thing is. But,
this was written purely by inspection, by looking at lot of candidates and saying look,
this looks like one candidate which can satisfy whatever I want. N8 can also be written like
this.
What is that I want to do? I want to now put down one new N1 in such a fashion that this
N1 can force the values of N’s at these nodes here go to zero. Is that clear? What
does it mean? It means that if I can make it go something like this, the shape go to
zero, then, I can get a new N1.
or in other words what it means is, my previous N for a quadrilateral element should go to
zero at this point and that point; is that clear, should go to zero at these two points
and should not, two things; it should go to zero at these two points and should not now
get to non zero value at all other nodes. Wherever those fellows are zero, they should
remain the same. I mean, look at that very carefully, I have written down N5 and N8 by
inspection, by saying that 5 and 8, it should be, you know N should be equal to 1 at 5 and
8 respectively; should go to zero at all other positions. By inspection, by considerable
experience, I have written N5 and N8 like that.
What I want to now do is to force this N1. I will tell you, why? What is important later,
why we are writing it like this later? I want now this N’s to go to zero at these two
places and retain zeroes wherever it is going, same level. What does it mean or is there
any clue? You understand what I am saying? Is there any clue from what I have done or
what I have said? In other words, can you now play around with N1 N5 and N8, because
N5 and N8 are not zeroes at the places where I want it to go to zero, 5 and 8. Is that
clear?
Can I now play with N5 and N8, so that I will make this N1 go to zero at 5 and 8 as well
and not make it raise above zero at 2 3 4 and so on. No, no, no, wait, wait, wait; have
a look at it. Do not rush, do not rush. Do you understand? Now, I have this much value.
I have to make this fellow, multiply it or by something subtract it. How do I do? Is
that clear? I have to, I have to make this fellow, this length here, go to zero. The
way I can do it is to do something with the shape function. Yeah, how do you or what do
I do? So, that is equal to 1 according to this, N5. Say for example, take only 5, N5;
so, one place we will make it go to zero. What is it that
I do? What is the value here? 0.5. 2N5 if I put, what will happen? If I put 2N1, then
what will happen? It is good, you know, I am happy with these answers. If I put 2N1,
what will happen to N1 at 1? It will become 2. It will become 2, so, I cannot use 2N1.
Yeah, come on, someone said that answer. N5 by, if we put N5 by 2 alone, correct; that
is it. N1, my linear fellow, minus half N5. I am spending time, because purely by inspection,
by intuition; that is why it is called serendipity element. Serendipity is what? What is serendipity?
By accident, by chance; so, this element came about by chance and shape functions are written
purely by inspection. Yes, there is one more way of writing it, by our what? By straight
lines. We will do that later, because this gives an insight as to how the shape functions
look like. That is why I am going through this. Is it complete? Is it complete? You
have not completed it. Yeah, but you have made one fellow go to zero. This fellow goes
to zero, but you have not touched this fellow. What? What do I do? Because, what happens
to N5 at this position? It is zero. Is it plus? Why, why should I have a plus there?
Minus N8 by 2, period; so minus N8 by 2.
N1 is equal to that linear N1(l) minus N1, sorry, N5 by 2 minus N8 by 2 will be my new
shape function. Similarly, look at this. This is quadratic one side and linear on the other
side. Similarly, I can write down. I am not going to do it now, but you can say that similarly
you can write down for N2 N3 and N4 N6 and N7. Is it clear? What is the advantage of
writing like this? First of all coding becomes very simple. See, for example you may decide,
you may decide to use higher order elements at some places and stick to your linear elements
at some other places. So, there will be a transition element. What will be the transition
element or how would it look like? How would it look like, a transition element?
Say, I am going to have 8 noded elements here, all these 8 loaded elements. Then, I want
to now go to a linear element in the next step. So, I will have an element like this.
Look at this element, look at that element. That element has 5 nodes, 5 nodes; 3 one side,
other sides have 2 each; 5 elements. If I develop my shape function like this, by inspection,
by writing down N5, knocking off from N1 such that making it zero at other positions and
so on. What will happen if I have an element like this? Is it first valid or not? Sorry,
is it valid? Can you write down now shape function for this? Say for example, 1 2 3
4; 5 is absent, 6 is absent, 7 is absent. So, I have only 8. Can you handle this element?
Of course, you can handle it. How do you now write the shape functions for this?
You write down first the shape function for N8 in the same fashion as we did right now.
Half into 1 minus eta squared into 1 minus psi and then what do you do? Knock off from
N1, half N8, period; knock off from N1 half N8 and from what else do you knock off? N4;
N4 knock them off, so that you will now get a new element itself, 5 noded element. Shape
functions can be very nicely written. Is that clear? What we have done essentially right
now, is to write shape functions for quadratic elements. I have written down by inspection.
If you take some time to sink, my suggestion is go and think about this, this particular
thing and now we will extend this concept. By the way, what do I do now with J, I mean,
next step? Write DN, write J B transpose DB T J the psi d eta. Steps are the same; they
are no more different. We will comment, is this clear, any question on this? Only thing
is, if you know how to write the shape function, period; the whole thing follows. I am not
going to emphasize more on this. The only thing, major emphasis that is required is
on our numerical integration. There are lots of issues on numerical integration. We will
talk about numerical integration, the issues involved in numerical integration in the next
class.