Finite Element Method

Prof. C.S. Upadhyay

Department of Mechanical Engineering
Indian Institute of Technology, Kanpur

Module - 3 Lecture - 3
In the last lecture, we had stopped at the derivation of the bilinear form for the finite element
solution and we had talked about the Orthogonality of the Error; what is orthogonality condition?
Let me reiterate, what the orthogonality condition was.
(Refer slide time 00:29)

We got B (u minus u
FE
For all such v(s) this is true, here we said that one candidate v is u
,v) is equal to 0 for all v given by this representation. All vs which could
be represented as a linear combination of our global basis functions that we have used for the
approximation and such that v also satisfies a geometric condition. For example, for the problem
we have at x equal to 0, u equal to 0, so v at 0 has to be also 0.
FE
itself because u
FE
satisfies
zero condition at the point x equal to 0 for our problem and it can also be represented in terms of
the shape function as basis function as we have done, so from this condition we reach the
conclusion that B (u minus u
FE
) is equal to B (u
FE
, u
FE
). We are going to use it for determining
what is B (u minus u
FE
, u minus u
FE
(Refer Slide Time: 02:59)
), first of all physically what is this equal to? Let us go back
to our definition of B.

So B of, if I put some function u, u is equal to integral over x is equal to 0 to L, EA du dx whole
square dx. What is this quantity equal to? This is as we have done earlier, is twice the strain
energy due to the function u. If I have a deflection u, due to the defluxion u or displacement u,
we have the strain in the material, strain needs to stress, the strain energy of the body is equal to
my calligraphic u and it turns out that d by the definition is nothing but twice the strain energy of
u, so when I am talking of B (u minus u
FE
). First of all (u minus u
FE
) as we had defined is equal
to error in the solution. This is the exact, this is the finite element the difference of the two is the
error. I could write B (u minus u
FE
) in the short form as this is equal to B (e minus e); this is
nothing but, by our definition twice the strain energy of the error. What does it mean? If I give
the defluxion or displacement which is equal to (u minus u
FE
) to my body of this bar then the
strain energy due to the displacement is going to be u (e). If I have the strain energy then I can
expand it.
(Refer slide time 05:00)

This will be equal to B (u,u) this is a very simple job to do because these are nothing but integral.
We expand the integral is now and we can see the bilinear form coming out to this type minus
2B (u,u
FE
) plus B (u
FE
,u
FE
). This quantity that we have here is already known. How is it known?
We have written this quantity in terms of the equivalent definition as 2B (u
FE
minus u
FE
). So B
(e,e) becomes minus B (u
FE,
u
FE
). If I look at B (e,e) by the definition of the integral that we have
given here any expression instead of u I put e then this integral, the quantity under the integral
sign is always greater than or equal to 0. Why so? Because the integral if we see is the square
term, if it is square term that is the area under the curve is always the positive area so this has to
be greater than equal to 0. This is the true for the u if I put e the same thing to happen. What I
know is that here, in our expression here this quantity greater than equal to 0. It is equal to 0,
when e is equal to 0.
(Refer slide time 08:00)

What do I have then? Because of this I will have B (u,u) minus B(u
FE
,u
FE
) is greater than equal
to 0. What is this expression equal to? This is equal to two of the strain energy due to u minus the
strain energy due to the u
FE
The strain energy due to the exact solution is greater than equal to strain energy due to the finite
element solution. How does it translate in to the parlance that is u in to the finite element? This
means that since, strain energy due to the exact solution is the greater than that of the finite
element solution we say that the finite element solution is stiffer than u. One thing we should
keep in mind that all this could be done provided u
is greater than equal to 0.
FE
We know that the finite element solution is stiffer in what sense that is the strain energy due to
the finite element solution is always lesser than that of the exact solution. This could be a very
good check for the finite element program that one writes. If I have a problem with geometric
constraints which are homogeneous in nature then this is true irrespective of bilinear formulation,
this is the generic representation.
satisfies homogeneous boundary condition
that is the only to the problem where the geometric constrains are homogeneous in nature. If they
are not, we cannot say anything about that problem.
If I have that problem I solve that I should see that as the solution is refined that is I take more
element in the mesh; that is N is increased, N is the number of element, increase N or I change
the P. P is increased in both the cases, I would expect that the strain energy due to the finite
element solution should approach the strain energy due to the exact solution from below. It
should be a monotonic conversion to the strain energy of the exact solution. If strain energy, the
computed strain energy goes like this that is it is sometimes increasing, sometimes decreasing
something is wrong with finite element program or the computation that I have done. Obviously,
one question that will come to mind is how do I compute the strain energy of the finite element
solution?
(Refer Slide Time: 11:32)

If I look at the strain energy for this simple problem it is equal to half integral of x is equal to 0
to L. I am writing it in a loose way because we should have taken summation over the elements it
should be EA (du
FE
by dx) square dx. What is u
FE
equal to? u
FE
is equal to sigma of u
i
phi
i
. This
will be equal to I can rewrite it integral x is equal to 0 to L EA in to sigma i is equal to one to NP
plus 1. P is the order of approximation, n is the number of elements so total number of unknown
is NP plus 1 in to u
i
d phi
i
by dx whole square dx. What is this will be equal to? This will be
equal to if I expand it out through the summation i is equal to one to NP plus 1 summation j is
equal to one to NP plus 1, u
i
, u
j
integral x is equal to 0 to L EA phi
i
prime, phi
j
prime.
What is this quantity? This is nothing but our global stiffness entry K
ij
.

(Refer Slide Time: 14:38)
So this summation will
be equal to half of U transpose in to K in to U. When I solved the problem I obtain this vector U.
I already have stiffness matrix, global stiffness matrix all I have to do is take this product. This is
the strain energy of the finite element solution. This is very easy to compute and it is certainly
for any one is writing the code it is healthy practice to check it. Create problem for which I know
what the finite element solution should do and the check whether my code is really doing that or
not. We have said that finite element solution is preferred?

Then the question is that this is equal to u of the error. What does this signifies? If we see that
when the finite element solution is equal or is the same as the exact solution, the displacement
due to the difference of these two solution that is (u minus u
FE
) becomes smaller and smaller. As
u
FE
goes closer u, (u minus u
FE
) becomes smaller and smaller because of the strain due to the
error which is (u minus u
FE
I have a finite element solution converges to the exact solution. The strain energy due to the error
should fall down. At what rate the u(e) comes down? Let us take a simple problem. I take a
loading f is equal to x to the power 2 here i put P is equal to 1 and I put the material EA is equal
to 1 and there is no distributed string.
) becomes smaller, when smaller and the strain becomes smaller, the
stress becomes smaller so this also becomes smaller.
Length is one, this is x is equal to 0, this is x is equal to 1. The length of the bar one and taken
the loading as x

to the power 2. Obviously, I know that this problem becomes the differential
equation is given by this and u
0
(Refer Slide Time: 17:00)
is 0 du dx at x equal to one is equal to one.

For this problem since this is the second derivative of u is equal to x

to the power 2. The exact
solution is fourth one so u exact will be equal to x

For these meshes I will find the finite element solution. As the mesh is becoming finer and finer
the error that is the (u-u
of 4 by half plus Ax plus B where the A and
B are easily obtained from the boundary condition. This is going to be my exact solution, for this
I would like to take a mesh this is a uniform element that is all element has the same size h. The
mesh size is h and the number of element N is equal to 1 by h. I will use an order of
approximation over the full mesh of either order 1 or 2 or 3. I will Fix P that is Fix P and then
refine the mesh that is I am going to take N is equal to 2, 4, 8, 16, 32, 64 and so on. For this
various measure that is I can taken P equal to either 1 or 2 or 3 and for this P, I take the number
of elements in the mesh 2, 4, 8, 16, 32, 64 so it is doubling every time.
FE
) should become smaller and smaller that is u
FE
should approach u. If
we expect as we do this e should decrease. How we would like to know, how fast does the e
decrease with respect to the h? That is mesh size and the fixed order of the approximation P
where the fixed order of the approximation P and h with respect to this how fast does it decrease.
(Refer slide time 19:54)

We do not look at the error directly that is we look at the strain energy of the error. Let us look at
the plot of the decay of this error because we know here what is exactly solution error? We can
find the decay of the error as h is reduced that is the number of element in the mesh is increased.
For the fixed e we will do this and what we are going to plot is log of square root of the strain
energy of the error verses log of h. These quantities we are going to split to plot. Let us look at
the plot from the plot we see that first line which is the blue line.
(Refer Slide Time: 20:07)

Blue line corresponds to P is equal to 1 approximation and we see that the error decays the log of
the square root of the strain energy of the error decays with the slope of 1. The slope is
approximately 1.9743 but approximately 1.
(Refer Slide Time: 20:56)

If I am going to fix P is equal to 2, then when I plot this.
(Refer Slide Time: 21:03)

The slope of the line is approximately 1.980 which is 2.
(Refer Slide Time: 21:10)

For P is equal to 3 that they decay the error, the strain energy of the error approximately 3. So I
can reach the general conclusion out of this, the answer is here.
(Refer Slide Time: 21:24)

Though this is one particular example we have taken it to demonstrate the nature of decade of the
error.
(Refer Slide Time: 21:39)

It turns out the square root of strain energy of the error is less than or equal to some constant C in
to h to the power of P. What should the C define? Do not bother out it now what it says? If I fix
the P and I change the h that is refine the mesh then the strain energy of the error decades or falls
down at this rate. Why the strain energy is measured? Because the strain energy is never
negative. When the strain energy becomes 0 to we know that our exact solution and the finite
element solution are the same. This is called the A priori, error estimate for h version of FEM.
What do we mean by the h version? There are two or three various approaches we can take to
reduce the error. Our job is to get the finite element solution as close as possible to the exact
solution. We can achieve the goal to different approach. One approach that we have taken is fix e
we only remain the mesh so this is called the h version. Only h is changed, the parameter h
another approach that we can have is fixed h, fixed mesh and increase the P. If we increase the P
then we get the P version of the finite element method.
We are not refining the mesh we are only increasing the P. We hope that the error comes down
and indeed for such smooth problem this can be exponential, so the decaying error is
exponential. When this is exponential why dont we use this? The answer is that this is defined
on the smoothness of the solution. This C plays the row, how good is the exact solution is done?
In problem where we have domains in corner such that the solution is not smooth. In that case, P
metal will not do good job as we have shown here. Similarly, the h metal will not do good job as
we shown there in that case we have to have interplay between the h and the P.
(Refer Slide Time: 24:38)

This gives rise to hP in a form, so in the hP method both the h and the P are suitably changed in
order to get good conversion of the solution. This is about the conversions of the error in the
strain energy. How does the strain energy of error converts? Obviously, there is other question
that come to mind is how do the point wise value of the error we get?
(Refer Slide Time: 25:27)

If I have P is equal to 1 let us look at the next figure. From the next figure, we see that the error
at the corner node is 0 and because it is the same problem that we have taken the exact solution
of x

power 4 defines that the error at the node value is 0.
(Refer Slide Time: 25:44)

For this problem, we get the nodal values are exact, we are interested in the derivative
information much more than that the value of the displacement itself.
(Refer Slide Time: 26:04)

In the second part of the figure we will note that the derivative of the error vanishes at the
midpoint of the element.
(Refer Slide Time: 26:23)

So midpoint of the element derivative of the error is 0 so the derivative of the finite element
solution is very good at the midpoint of the element. What else we note?
(Refer Slide Time: 26:48)

We see that at the end points of the element, the derivative of the finite element solution is very
far from the exact form. If we look at the error plot, this error is quite high at the nodes of the
elements. This is very important feature of the finite element solution. In fact this is the
drawback of it. We know that in the physical situation we continue the derivative of the solution
has to be continued. Unless we have the material interface here we have the no material interface
which is only two elements, start the interface of the element we should have continuity of the
derivative but it is not. In fact, that is very bad and from the two sides the value need not be the
same, value of the difference that is the error of the derivative need not be the same it is also
called the scattering phenomena.
(Refer Slide Time: 27:51)

Let us now and see whether the situation improves when I take P is equal to 2 that is the next
figure.
(Refer Slide Time: 28:03)

If I look at P is equal to 2. We observe that the nodal values are exact and at some point in the
middle most probably at the center of the element which is also in the lagrangian element ,the
solution node at these 3 point we find as the solution is very good.
(Refer Slide Time: 28:27)

We can only say in the general case node for this one we can observe this thing that the nodal
values are exact. It comes out that here also the mid point is good but we cannot generally make
that statement.
(Refer Slide Time: 28:51)

What about the derivatives? How we see that there are accurately two points in the element?
There the derivative is very good.
(Refer Slide Time: 29:00)

So du
FE
(Refer Slide Time: 29:24)
dx is good or exact in this case had two points in the element. Obviously, it is quite clear
from the figure that this is the three element mesh.

We have three element of the same size for which we find that the each of the element that the
derivative of the error vanishes at two points when P is equal to 2. This is the very curious thing
which is the feature of the finite element solution.
(Refer Slide Time: 29:48)

Also when I have P is equal to 2, I have two points where the derivative of the error is 0. Let us
go to the next figure which is the P is equal to 3 figure.
(Refer Slide Time: 30:03)

Here we look at the behavior of the error. The error has the change at these end points of the
element and the derivative again jumps a lot that is the error at the end points of the error in the
derivative is quite high. Within the element there are three points.
(Refer Slide Time: 30:29)

As far as the nodal values are concerned, nodal values are exact but there are three points in the
elements. The du
FE
(Refer Slide Time: 31:16)
is very good at three points in the elements. In each of the elements it is a
trend.

There is generally for the element we should get P points of good du
FE
dx. These points in the
general case are called points of super convergence. Remember that, instead of getting answer
they are more questions in the sense at the inter element interface, if I want to get the derivative
of the finite element solution accurately at this interface of the two elements at the corner nodes.
So track something that we are not going to discuss now. Good values of derivative they can be
obtained by suitable averaging. Obviously, very simple thing that one can do if I look at the
linear case one averaging that we can certainly do, if I have these three elements.
(Refer slide time 32:53)

We find that at the mid point of the element the derivative is very good. We can pass a straight
line between the values of the finite element derivative of the finite element solution so that
straight line this is u
FE
prime at this point u
FE
There may be the value here is this and between these two I can pass the straight line so I can do
various kind of fitting a greater value of the derivative that is by itself a whole field and we are
not going to indulge in that. We have seen that how the solution converges in the global mark,
global measure define as the strain energy and as well as the point wise time. This is very
specific especially, the point wise results to problems with EA equal to constant and K is equal to
0. If EA is not equal to constant if I have this kind of a situation. I have a bar with access string
distributed this way. In both the cases we can see that the nodal values is not exact, the u
prime at this point.
FE
Similarly, the derivative will not be the way it has shown here. It is very good at this points that
we have shown with our examples. We have talked about the error purpose problem and how it
converges? One can go ahead and write a finite element program and use it to solve the problems
and based on kind of accuracy we want and the kind of computing power one have.
is not
exact at the node that is the nodal values are not exact however these nodal values are very good.
(Refer Slide Time: 35:59)

Again remember that we think that we have look at are accuracy and economy computational
part is very important. Why economy? Because here if we have Pth order element then we have
NP plus 1 unknown here the rate of conversion in the strain energy is also type Ch
The cost will be given in terms of the number of operation will be type (NP plus 1)
power p. We
will have interplay between these two things because this if this is so then the number of
computation or number of floating point operation that have to be done in order to find the
inverse of the k matrix which is also of this type is NP plus 1 that is K.
3
Once we have determined these things. Let us now again reiterate or go over what are the various
pieces in standard finite element program? What we have done is? If we remember that we first
started of making the nodes.
if I use the
standard gross format. We have to look at the cost of the computation as compare to the accuracy
we want and accordingly we will have to decide on which method to which approach to use that
is whether i should use a higher P or smaller h or whatever the combination of these two things.
For the example problem, we have taken here obviously taken higher P is more economical that
is fix a mesh take the two elements mesh and increase the P then the cost of the computation will
be small.
(Refer Slide Time: 38:03)

I deliberately, first make the nodes, the n nodes of the elements. If I have to make the nodes how
will I make the node? First I am going to choose my P, rate from the P now we know P and the N
which is the number of elements given N we know h is equal to L by N if I using the uniform
mesh. Let us take a uniform mesh, h is equal to Lby N. Given the P, I have to put these P minus
1 equally spaced points in each element. How do I find those points? I will call something as the
h bar is equal to h by P this is the spacing between two consecutive solution node. How do we
now go ahead and make the mesh? We have to make the mesh solution node so I start with x
1
is
equal to 0 or given value and then we say x
1

plus i
equal to x
1
So if I have this thing then I know that it will start with plus one to P in to k minus one plus P
plus one. So these are the nodes in the elements this is the set of nodes just sitting there in the
element. However we know that the first and the P plus 1 nodes of the extremity of the elements,
we need x
plus i h bar. This is the given mesh.
Once I have the mesh then what do we do? We have to construct the shape functions of the mesh
generally falls Preprocessor. What other information we should create here? After making the
mesh we should also make the element information that is in the element k which is the nodes
this node that we have created these points, very important because we would like to know the
extra of the element as well as we need to point to construct the shape function.
p(k minus 1)
plus 1 and x
p(k minus 1)
plus (P plus 1). The given these two nodes we can do
the integration. All the nodes are now used to define the shape function. The shape function is
defined by taking the product rule we can use available symbolic computation software like
mathematica, Apple or Matlab. To get the derivatives of these expressions there is the shape
function because of the linear quadratic easy to do, for the higher order elements it may be little
bit cumbersome. Better let the software do it, could do it through hand also but still take time.
Once I have this information we see that the other things which is there I really lead these
numbers. Once I have these numbers it is also local to global numbering that we were talking
from the assembly point of view.
This number give all that information that from the assembly point where does the local entry go
to globally as well as the what is coordinate of the points corresponding to this location. Thats
why generally in the one way approximation where we use this kind of a lagrangian shape
functions. In these nodes the solutions points they are indistinguishable from the geometric
points. They are both called nodes for the solution value is for the degree of nodes that is every
nodes corresponding degree of freedom. This is problem has P plus NP plus 1 degrees of
freedom other nodes. Once I have the mesh information available now for the elements I have all
the all the relevant information.
(Refer Slide Time: 43:06)

Then I will do the element calculation which we have done in the some detail which will give us
the element stiffness matrix and the element load vector then we will do the assembly process
then we apply the boundary conditions and then solve, solve means use gauss elimination or any
solver which we must already studied. The whole part is called processor in a finite element code
which is called a processor. We obtain the solution as we have seen that direct finite element data
especially the derivative information may not be too good even it is acceptable to us then we will
have to display that. So display becomes very important point or we have to output the values of
this finite element solution as well as the derivatives at some point.
(Refer slide time 44:22)

Output response info here we may do direct information is outputted or some processed
information. This whole part of outputting in the form which this is called the post processor. For
the finite element program we will have these three big components and these 3 big components
we have several smaller jobs. We can see very clearly that which is already broken in to modular
form and module can be separately programmed, integrated and used. However here there is
certain problem that we have to face.
Till now we have dealt with very simple problem where the EA was constant, F was relatively
simple and still when we go to higher Ps the whole business of integration becomes very
cumbersome. When the material is varying, the cross section may be varying or the material is
varying or in case the load is not a typical load which is the given really thick load. All the cases
doing this job in close form way or the computation till now in the close form. In the close form
may not be so easy. Let the computer do the integration process for which many of us must be
aware of it, we use numerical integration.
In numerical integration that is instead of writing integral we would like to use summation. To
do that we will have to change the way we have done think a little bit, the whole process remains
the same but how we are making it more amenable to computer information. Till now, most of
the things done by the hand or via symbolic computing platform. Here we would like to write a
program which we can do the numerical calculation given this stiffness matrix and the load
vector matrix and go ahead and do the assembly and everything else. What we see in the next
sort of lecture is how to do this job or transferring what we have done till now to a computer?
We will introduce the concept of master element. How we define everything with respect to the
master elements? How do we transform of physical domain to the master element and do the
integration converted to numerical integration and so on.
Thank you.