Hello, welcome to basics of finite element
analysis course, today is the last day of
this week and in today's lecture we will discuss
a particular Eigen value problem and actually
develop its finite element formulation and
solve it, so this problem relates to vibrations
in a bar and the governing differential equation
of the bar is.
? a times ?u with respect to time - ?e, a
??u ??x with respect to x equals Q and to
explain 
suppose this is the point where, which I am
considering u is the displacement in next
direction and it can vary with time and position
this is my x-coordinate okay, ? is the density
of the bar in this case we are assuming that
? is homogeneous or it is same on that particular
small element of the bar.
A is the radius of cross section, e is the
Young’s modulus and Q is the traction per
unit length, it is traction per unit length
so it is units of Newton's per meter. Our
goal is to find Eigen value problems, Eigen
values, and Eigen vectors using FEA okay.
Now our first step is so this is our second
order partial differential equation it is
second-order both in time and in space, so
the first thing we do is we try to convert
this into a second order ordinary differential
equation and we do that by realizing that
if the bar is vibrating like this, this motion
will be and if it is vibrating on its own
it will have some sort of a harmonic motion
with some natural frequency ? right.
So we say that u(x, t) equals a function u(
x) times e j?t oh excuse me it is e i?t, so
this is equation 1, this is equation 2, so
the amplitude of this motion can vary from
place to place but at every point the point
will vibrate in a harmonic way but each point
may have a different amplitude which is u(x)
but the motion of each point will be harmonic,
that is what it means. So now what I do is
put two in one, so what I get is ?A and I
am going to just write u rather than u(x)
for purposes of privacy so ?Au and because
when it is differentiated by time two times.
So I get –?2- ?? over ??x EA ??u ??x and
this entire thing ei?t right and because we
are interested in finding out the Eigen values
of the problem the right side of the system
will be 0. Now we note that ei?t is common
to the entire expression so I can eliminate
this right so I will erase this, the other
thing I note so this entire equation is now
only in x because this U.
U is nothing but a function of x only which
means that ?? over ??x equals d over dx right
so I will replace these partial differential
operators with regular differential, operators
total derivative operators, so it becomes
d over dx. So now I see that I can express
this entire expression as EAU´ is equal to
? AU?2 or I can write it as ?2 times ?AU,
and there should be a negative sign here okay,
so this if I look at this it is in the same
form as A(U) equals ?B (U) which is the standard
form for Eigen value problem.
So I have developed this converted this form
into an Eigen value problem statement okay,
now at this stage I am going to solve these
equations using the finite element method
and all of you already know how to solve this
problem using the finite element method. The
first step in the finite element method if
we want to solve this problem is we will break
the domain into small elements, the second
step will be we will assume interpolation
function for U, we will plug those things
here, we will develop and find the residue
of the overall thing over the entire element
multiplied by a weight function, equate that
residue weighted into the thing 0, weaken
the differential form, weaken the statement.
And develop a system of equations and then
we will solve those equations to get our answer,
so that is the overall process and that is
what we are going to do.
So our overall equation is EAU´ and I remove
the negative sign on this side equals ?2 and
what I will do is just for purposes of convenience
I will replace ?2 by ?, so it is ? times ?AU
okay, so the weak form, so now we are going
to get the weak form but before that we will
express it as a weighted residual form, so
at this stage so suppose my bar is this long
and this is one, so I have broken it into
several elements, this is node 1, node2 let
us say the element is he and so for the eth
element I have length is he and there are
only two nodes so I am assuming that this
is a linear element but before I get to the
element order I will develop a weighted residual
form, so I integrate it from zero to he, I
have a weight function w multiplied by x-EAU´
the entire thing differentiated with respect
to x- ? A?U dx equals 0 dx? because I am in
local coordinates.
From this I develop a weak form, so my weak
form is 0 to he and I reduce the differentiability
on the first term and shift it to w so I get
w´ and when I do that the negative sign before
EA goes away right -? wA?Udx? equals on the
right side I have, so I will get some extra
boundary terms and I am shifting them to the
right side and what I get is w evaluated 0
times Q1e + w evaluated1 no he times Q2e where
Q1e equals - AU´ evaluated at 0 and Q2e equals
–A this entire thing is evaluated.
AU´ evaluated at he okay. At this stage my
next thing is that I assume an interpolation
function for u, so u is equal to ? of Uj and
Øjx and this is for the eth element and j
equals 1 to m and in this case I am using
M is equal to 2 because I assume that it is
a linear element, also my weight function
equals Øi for the eth element which is the
variation using principles of variational
mechanics this is nothing but it represents
a variation in U.
So then I get 0 to he dØi over dx? EA dØj
over dx? all for eth equation –? A? over
here by the way a equals e times A so coming
back to my weak form this is my first term
–? A ? and w I have assumed as Øi for the
eth term and u will be Øj eth term, and because
I have to some, you is ? of uj Øj so I have
to multiply it by Uj and sum it up j is equal
to 1 to m and integrate it over the entire
domain or, because this is may not appear
on your screen clearly so I will again j is
equal to 1 to m and I am going to multiply
it by Uj e dx? okay, and this equals w evaluate
at 0 Q1e+ w evaluated he Q2e. Now we look
at these two terms, the one in green and the
one which is underlined as red okay, the green
term has EA embedded in it A is the area of
cross section and e is the Young’s modulus
of the material.
And so it represents kind of the stiffness
of the system, if you have a stiffer system
e will be high if you have a very compliant
system e will be less right so this represents.
Stiffness of the system, on the other hand
the red term the one which is underlined red
if we do not worry about ? then it is a times
? times dx is, so suppose you have a small
element which is dx long, cross-sectional
area is a times ?, ? is the density so it
represents mass of the system.
Okay so what we get is two matrixes k – ? m
and both these matrices are for the eth element
and this entire thing is multiplied by an
unknown vector U, and on the right side I
get a Q vector a force vector, but these forces
are internal to the system okay, these are
internal to the system. Here ki j for the
eth element equals integral of dØi over dx?,
dØj over dx? times EA dx?, also the mass
matrix so the stiffness matrix is defined
earlier and then mass matrix is defined as
integral of 
A ?ØiØj dx? and of course everything is
with the super script d so that is my stiffness
matrix and that is my mass matrix.
So in our earlier problems we had only one
metrics which was multiplied by the unknown
vector u, here we have two matrices one is
the stiffness matrix the other one is the
mass matrix right and they are and the mass
matrix is multiplied by that ? which is the
Eigen value of the system, and on the right
side you have the Q vector which is basically
the internal reaction forces at the two ends
of the element okay, couple of things.
So this is stiffness matrix if you look at
it, it is symmetric because if you replace
jby i and i by j you get the same answer,
so k ij is equal to k ji. The mass matrix
is also symmetric in nature because when you
replace i by j so M i,j is equal to Mji okay
so this is the, so these are the two matrices
k and M, so now at the element level we have
developed the Eigen value formulation for
the problem.
So our next step is to assemble all these
things okay so.
That is what we are going to do, so before
I do that actually I will write down the values
of these terms so K for the eth element is
equal to EA over he and this is for the case
when m is equal to two, it is for a linear
element, for a quadratic element it will be
a three-by-three matrix okay for a linear
element it will be a two-by-two matrix and
the terms are 1 -1 -1 and 1, and then the
mass matrix for the eth element is ? Ahe over
6 2 1 1 2, there is one thing interesting
usually you will note at the mass matrix,
the overall mass of the system which is an
element he long will be what, its cross-sectional
area times density times he right.
So you have in mass matrix at least in this
case the numerator ? Ahe is the mass it is
divided by 6 and then you look at all the
terms in the matrix, you have 2 2 1 and 1
if you add them up two plus, two plus, two
plus, four plus one five and one 6 that divided
by six is so, if you add up all the elements
of the mass matrix they end up with the total
mass and the mass matrix can never be negative
because mass is never negative it can never
be negative okay
So there is something interesting about your
mass matrix, so with this understanding what
we will do is we will develop assembly for
our system, so our system is like this, so
it is a long bar and at this end of the bae
it is rigidly fixed, at the other end of the
bar I have a spring whose stiffness is K okay,
now what I am going to do is I am going to
split this bar into four elements.
So this is element number one, element number
two, element number three, element number
four, so I have a total nodes number of nodes
is one, two, three, four and five okay, and
I am interested in finding out the Eigen values
of this problem okay. One thing you will immediately
notice that suppose this spring was not there
then the bar would be would have less problem
vibrating back and forth right, once the spring
is there then the overall stiff system becomes
stiffer.
So bar will have a higher angular frequency
because its over all stiffness that is k it
goes up right so we have to figure out how
to incorporate the
Influence of this k in the overall assembly
level equations and to do that first we have
to develop the assembly level equations, so
let us write down these assembly level equations
for element 1. Our assembly level equation
is EA over he 1 -1 - 1 1 – ? ?Ahe over 6
2 1 1 2, the entire thing is multiplied by
U1 U2 e and in this case e is this the first
element so e is one and this equalsQ11Q21
okay, so this is the these are the two equations
for the first element. The equations for the
second element are identical except that the
superscript 1becomes two okay because we have
assumed that the length of the element is
same in all the cases.
So for element 2 EA over he 1- 1 -1 1 –? ?Ahe
over 6 2 1 1 2, this entire thing is multiplied
by U1 U2 for the second element and that equals
Q vector when Q1 Q2 for the second element.
Likewise I will have another set of equations
for element 3 and here it will be the U vector
will be U1 U23 = Q1 Q23 as the superscript
and for element four again similar set of
equations, U1 U244 equals Q1 Q244 okay, so
we have developed four sets of element level
element equations total number of equations
her is eight right and we have five nodes.
1, 2, 3 and 4 and 5, we have five nodes, so
the first thing we do 
is we established the continuity conditions
for U, what are those continuity conditions
at node 2 U21 is equal to U 12 is equal to
U2, so this becomes U2, this becomes U1 same
using same approach so this is U2 I am replacing
all these local degrees of freedom with their
global values, this is U3, this is U4, this
is U4, and this is U5 okay.
The second thing I have to note is that the
total force let us say node 2 so this is continuity
and they are 3, 3 more relations for that
and then the second one is force balance,
and for force balance the total force at node
2 is Q22 + Q13 right, similarly the total
force at node yeah total force at node 4 is
Q32 + Q14 , and the total force at node two
is the sum of these two correct, so if the
total force is that much then I have to add
these equations so I have to add this equation
with this equation, I have to add these two
equations, and I have to add these two equations,
and then in that way I will get total of five
equations, right now I have eight so I have
five equations and I have five degrees of
freedom.
So the total number of, so the final set of
equation becomes if I do that is EA over h0
1 -1 0 0 0o -1 2 -1 0 0 0 0 0 -1 2 -1 0 0
-1 2 -1 0 0 -1 1, this is my K matrix global
K matrix - ? ?Ahe over 6 and 
I get global mass matrix. What is the mass
matrix, 21000, 14100, 01410, 0 this is 00141,
00014 okay, and this entire thing would be
2, so this entire thing is multiplied by u1,
u2, u3 u4, u5, and this equals the Q vector
and the sum of Q12 this guy and this guy is
0 so what I will get is Q1 000 and Q5 okay.
Now let us look at Q5, Q5 is the force which
is being applied on the system due to the
spring, so if the node 5 it moves by a distance
u the spring will apply a force in the opposite
direction and its value will be K times u5.
So Q5 is equal to -K times u5okay, so this
guy becomes -K times u5 and what I, and because
this is unknown u5 I ship this, I move this
u5 on the left-hand side. So here what it
gets is I get another term here and this is
plus, I do not have space to write here, this
is negative and it will be K times h0 by EA.
Because there is a common term out here okay,
so this is coming this extra term is coming
because of this term, and once I have done
this then on the right side this term goes
away so this entire the Q5 the term the constant
on the right side in the fifth equation becomes
zero okay.
So this is the first boundary condition which
is this one, the second boundary condition
is that at x is equal to 0 that is at node
1, u is 0 so I specify this to be 0 okay.
So now I have to worry only about finding
u2 u3 u4 and u5 okay, which means yeah. So
I do not have to worry about my first equation,
I do not have to worry about my first equation
because u1 I have to only worry, I have to
only solve for u2, u3, u4 and u5.
So I have to worry about the four, last four
equations and also in the last four equations
when I multiply this minus 1, oh excuse me
by u1 it will be 0 so contribution of these
terms will also be 0 so I will consider only
this block in the K matrix, and for the same
logic I will consider only this block in the
m matrix. So first I have assembled, now I
have reduced the number of equations.
So ultimately what I get is, [[KRED] – ?[MRED]]
u2 u3 u4 u5 is equal to 0000 okay, now these,
these are four equations and four unknowns,
these equations will be true in two situations.
Situation one case 1, u2 = u3 = u4 = u5 = 0
this is one solution. But we do not like this
solution it is a trivial solution for this
I did not have to do all this finite element
analysis, the second solution will be case
two when the determinant of |[Kred]-? [Mred]|
equals 0, if this determinant is zero and
I know all the values in K reduced, I know
all the values in M reduce the only thing
which I do not know is ?, that is the only
thing I do not know.
So because these are four equations, I will
get a fourth-order algebraic equation in ?, fourth
order equation, so it will have ?4 times something
plus ? Q times something and so on and so
forth. So I will get four values, four values
of ?s, each of these values is called an Eigen
value because it is a characteristic of the
system, that is what it literally means okay.
Now once I get that ? so, I so I get four
?s, ?1, ?2, ?3, ?4. But what is our goal we
have to find u right, we have to find u that
was our goal so we have found these ?s.
Next what I do is I consider ?1 and 
I plug that ?1 in [Kred - ?1 Mred] times this
reduced vector u2 u3 u4 u5 equals 0. Now I
know everything in this matrix, I also know
?1, so I can find the relationships between
u2, u3, u4 and u5. When I get those relationships
that solution for u is called Eigen vector
for system corresponding to ?1, corresponding
to the first Eigen value. Similarly I consider
now the second Eigen value so I get a second
Eigen vector, then I take the third Eigen
value I get the third Eigen vector, and if
their system is having n number of equations
reduced number of equations then I will get
n Eigen values, n Eigen vectors and each Eigen
vector will be corresponding to one particular
Eigen values okay. So that is the, how I figured
out my Eigen vectors okay.
And finally our original goal was what to
find u, our original goal was to find u.
And we had said that u(x, d) is u(x)ei?t right,
so u(x) so in ? we know that ? equals ?2 that
is what we had assumed right. So I know what
is ?, so I take first ?1 come corresponding
to ?1, I put ?1 in this equation, and I also
put u1 which is the first Eigen vector. So
that is my first solution, then I take the
second Eigen value, I find out its Eigen vector,
put both these Eigen vectors and Eigen values
in this solution, and I get the second solution
for the system, that is the that is the second
solution and each Eigen vector in this case
is also called a mode shape of the system.
Mode shape of the system means so you have
four nodes suppose two three four and five,
first node does not move because why? Because
it is fixed, so what does a mode shape mean
that when you are, when the system gets excited
by first Eigen frequency maybe u1 is like
this, u3 is this, this and this, so this is
your first mode shape, which means this is
the shape or the boundaries between which
the system is going to vibrate you know, maybe
the second mode shape will be this.
So it all depends on the system you know,
so you will have four different shapes corresponding
to four different Eigen values and the overall
solution for each Eigen value is this, so
this completes our treatment for the Eigen
value problem. Next week we will, which will
be our final week we will extend this discussion
into two important topics. The first topic
will be time dependent problems, how do we
solve time-dependent problems, and the second
topic we will cover will be numerical integration
that when we develop all these K matrix relations
how do we numerically integrate them using
our computers. So that is pretty much it and
look forward to seeing you tomorrow. Bye.
