# Mod-11 Lec-3 Global Finite Element Assembly and Imposition of Boundary Conditions

The last class we studied the finite element
analysis of transverse vibration of beam and we developed elemental equation for the mass,
difference and force vector and this particular development of the elemental equation. For that we used Galerkin method and there
is another method which is called Rayleigh-ritz method that method can be used for development
of the assembler elemental equation. So far gravity we are not covering that method. Now, in the present lecture we will try to
see the application of the finite element which we have developed earlier or we can
able to use that for different end conditions or for a system. And for that initially we will be taking 1
simple example of simply supported beam and we will try to find out the natural frequency
of simply supported beam. And even corresponding mode shapes will plot
and we will show the detail procedure how it can be obtained. Even for a rotor system we will formulate
using infinite element method the procedure how to get the unbalanced response if a particular
disk is having some unbalanced in the rotor system. And then we will see some other techniques
specially some kind of a condensation schemes if the elements numbers are very high, then
how we can able to reduce the system equation of motion. So, that the computation times time is less. So, let us take on example of a simply supported
beam. So, we are considering 1 simply supported
beam .this the support of the beam now for the property of the beam is diameter
of the shaft is and say ten mm the length total length of the beam is 3 meter Young’s
modulus of the beam material is this much density is of steel and so with the help of
this we can able to calculate the second moment of area that comes out to be this much area
of the cross section of the beam would be this much. And now we want to a discretize the beam in
to various number of elements and for that for simplicity of illustration we will divide
this element into 3 number. So, let us see equal number of the size of
the length of the element was same 1 meter 1 meter, so there are 3 element 1 2 and 3,
I can able to give them node numbers also 1 2 3 and this is the forth node .So, now
we will be writing the elemental equation for each of these elements 1 by 1. This is the finite element equation for element
1. In this we have first the mass matrix which
is 4 by 4 size on this instead of variable we are keeping all the value of the variable,
because this matrix is symmetric. So, I am not writing the full matrix and here
we have the acceleration vector which is as given as w 1 double dot that is a time derivative
then w 1 double dot prime. Prime represent the derivative with respect
to spatial that is x variable and this is for w 2 and this is w 2 prime double dot prime
represent the slope. So, this is the mass matrix similarly the
stiffness matrix will take this form this also four by four matrix, this matrix is also
symmetric and here the displacement and slope vector are there for node 1 and node two because
element 1 is having node 1 and node two. So, corresponding displacement and slopes
are having subscript 1 and two, then in the right hand side we have the shear force and
bending moment there is the reaction forces and moments. So, this is the elemental equation for the
first element now we will be writing the elemental equation for second element, element two. So, in this again I am giving the detail steps. So, that the procedure is clear. This is the mass matrix this is because the
element size is same and dimensions of the element are also same. So, we are having this mass matrix identical
to the previous 1, if there is a step change in the diameter of the shaft then may be will
be having changes specially in this terms; because they contain the dimensions of the
shaft. Here, again we will be having some changes
because now the element two contains node two and node three in left and right side
of the element two .So, we have w 2 w 2 prime w 3 w 3 prime and these are accelerations. So, double dot terms; will be there then comes
the stiffness matrix. Here is similar to the previous element, element
1 because this element is of the same geometry, same material property is the previous 1. So, they remain the same this again both are
symmetric here, we have displacement and slope vectors at node two and node three, then in
the right hand side, because we are considering only the free vibration. So, there are no external forces and moments
only internal forces and moments are there. So, here correspondingly we will be having
subscript two and three in the shear force and moments. So, now shown how to the two element equations
1 element is still remaining. So, for completeness I am showing that also. So, this is the element three in this I am
just showing the mass matrix which is exactly same as the previous 1 whatever the changes
are there I am mentioning here So, here since element two, three is having node a three
and four. So, the subscripts will be accordingly changing
to three and four. This the mass matrix stiffness matrix is also
same because the dimension and the material property of the element is same as the previous
once only these subscript will change to three and four corresponding to the displacement
and the slope. And in the right hand side again there will
be shear force corresponding to node three and node four. So, these are the three elemental equation
which we have illustrated, because now we want to get the system equation. So, we need to assemble all these three elements. So, now I will be showing the assembly procedure
how we will be assembling these three element to get the system equation. Because this is a big matrix. So, I will be showing 1 matrix in this particular
slide the mass matrix and similarly; I will be showing the stiffness in the next slide. So, here we have 1.46 into 10 raise to minus
3 coefficient which is common then we have 156, 22 these entries are coming from the
first element, if you see previous mass matrix for element 1 these are the entries in that. So, this coming from the element 1. Now, for element two we can use in different
color 156 plus 22 54 minus thirteen plus 4 thirteen minus 3 this is red color entries
are coming from the element two, or if we want to see how they are coming here, let
us first write the vectors. So, here w 1, w prime, double dot is there
then we have w 1 double dot prime ,w 2 double dot prime is not there here, then w 2 double
dot prime, w 3 double dot, w 3 double dot prime, w 4 double dot ,w 4 double dot prime. So, we can see that the first four entries
in this vector corresponding to element 1, and the element two contains the node two
also. So, there will be overlapping when we will
assemble the element 1 and two there will be overlapping of the element 1 and element
1 mass matrix specially, at node two which is common to both element 1 and two. So, these are the overlapping terms; we can
able to see next we will be entering the mass matrix for the third element. So, third element is having common node three
which is there for node two also and node three also if you see here or this element
three node two and three this is element two. So, element two this is three and this is
four. So, element two and three is having three
as common node. So, we are getting all these overlapping of
the terms of the mass matrix. So, here it will be fifty-four and then last
term will be thirteen, then here it is 4 13 minus 3 here 156 minus 22 and last 1 is 4. So, we can see that now the size of this matrix
is 8 into 8.Earlier it was it was 4 into 4 now it has become 8 into 8. Similar; assembly will be doing for the stiffness
matrix. Now, I am showing the stiffness matrix. So, these are the entries from the first element
and we have this vector w 1 w 1 prime, w 2 prime, w 3
,w 3 prime, w 4 and w 4 prime. For second element we can use different color. So, in this we will be adding 12 plus 12 plus
minus 6 plus 6 plus 4 plus 4 minus 12 these are coming from the second element these. Now, for the third element again we 12 plus
12 plus 6 there is 1 term missing here of the second element this. So, here we have two more entries 1 is minus
12 another is six. So, this plus 4 minus 6 and 2 then below this
twelve is there minus 6 and 4.So, we can see that these black entries are coming from there
element three and rest of the terms here they are 0. In previous mass matrix also rest of the terms
are 0. So, all these are 0 here also we have all
the terms 0 here also. So, this is the Stiffness matrix and previous
was the mass matrix. This is equal to the internal forces that
is minus S 1, M 1 then S 2, M 2 this is coming from the first element from second element
these are the entries this is S 3, M 3 then from the fourth element again we have this
is minus S 3 this is minus M 3, S 4, M 4. So, you can see that these terms are becoming
0. So, whatever the reaction forces we are acting
at nodes in various elements once we join them; obviously, they will cancel each other. So, because of that they are becoming 0 at
these two locations. And here which are the shaft ends there, because
it is supported through that is simply supports are there. So, this will be having shear force and bending
moment it depends upon the end condition because while development of the element we do not
consider the Boundary condition. So, now we assemble the all the three elements. Now, need is to find out what are the boundary
condition and apply suitably. So, let us see, now the Boundary conditions This was the node 1, 2, 3 and four. So, displacement, because this is a simply
supported end. So, displacement here it is 0 and bending
moment is 0 this is also simply supported. So, here w 4 displacement is 0 and bending
moment is also 0, because this is simply supported and. So, will not be having the shear force 0 and
the slope. Slope will also be having some finite value,
but for other kind of end conditions we will be having different Boundary conditions. So, these are the Boundary conditions we need
to apply to the, our equation assemble equation. So, let us see where will be substituting
these values. So, this is the mass matrix. So, here we can see that w 1 is 0 and w 4
is 0 similarly in the stiffness matrix we have these conditions. And in the internal force this M 1 is 0 and
M 4 is 0. So, in this particular equation in the right
hand side the equation is having this form I will just write in the compact form M w
double dot plus kw is equal to the internal forces. So, in this particular right hand side you
can see that this term is there in which S 4 is non 0 and S 1 is non zero. So, and all other entries are 0. And if we go back to the previous slide you
can see that w 1 is 0. So, when we multiply this vector with the
matrix k, whole column will be getting multiplied by the 0.So, it will not contribute anything. Similarly, this 1 w 4 is 0. So, last, but 1 column will not contribute
anything because they are getting multiplied by 0. So, we need to eliminate them, even the first
row also we need to eliminate because in this S 1 is unknown. So, we are eliminating them also last but,
1 row we need to eliminate. Similarly, in the mass matrix on the same
lines we can eliminate these 2 columns, the first column and the seventh column; and the
first row and the seventh row. So, whatever remaining terms are there that
will be writing separately that is nothing, but the system equation. after the application of the boundary conditions. So, that will be having reduced form and it
will take this form this is the mass matrix. So, we eliminated the first row and seventh
row and similarly these columns. So, we got these numbers also whatever the
contribution from various element was coming at the common node that we are adding up and
showing here is the resultant numbers. So, you can see that now the equation has
reduced from 8 by 8 to 6 by 6 and here we have vector corresponding to w double prime
1.So, w 1 is not here, because we eliminated them. Because, its values is 0 then w 2 double dot
is there, w 2 prime double dot is there, w 3 double dot is there, w 3 double dot prime
is there, W 4 is not there and w 4 double dot prime is there. So, this prime represents slope its value
is told. Rest of the terms here they are 0, these are
also 0 similarly this side also. Here also all terms are 0. So, this is the reduced form of the mass matrix. Similarly, we will be having Stiffness matrix. That is, this is the Stiffness reduced stiffness
matrix. So, I am writing the elimination of the first
row and seventh row and we are adding all the terms, which had from various elements. So, this is also; obviously, 6 by 6 matrix
remaining terms are 0.And here also as in the mass matrix we have w 1 prime , w 2 prime,
w 3 prime, and w 4 prime w 4 is not there because its values 0. And toward the right side we have a vector
contains 6 terms all 0. They are corresponding to the shear force
and bending moment which are already got cancelled at the intermediate nodes. So, this is a transpose of the vector. So, this particular matrix which we obtained
is of this form. This is a standard form of the governing equation
of a dynamic system. So, once we have the system equation in this
form now, we can able to analyze for the free vibration because in this case we are not
consider any external force. So, now we will be solving the free vibration
analysis of this. So, for that if, we resuming the displacement
vector as some amplitude and let us say Omega is the natural frequency of the system. If we differentiate this twice we will get
j Omega square ej Omega t, because this is independent of time. And if we substitute in this equation coining
equation, we will get minus Omega square M plus k and w bar which is amplitude will be
common this equal to 0. So, this is standard Eigen value problem
and we can get the Eigen values of this system in 2 forms, let us see, what are those 2 forms? So, if we in the previous equation, if we
multiply the whole terms by k inverse will get K inverse M minus just yes k inverse. So, taken Omega this side omega is nothing
but, natural frequency this is square of this. So, if we obtain the Eigen value of this let
us say, this is a A matrix. So, inverse of that and square root will give
us the natural frequency. So, whatever the Eigen value will be getting
from here we have to do the inverse of that and then take the square root that will be
the natural frequency of the system. Sometimes this Stiffness matrix is singular
then what happens we cannot able invert this the Stiffness matrix then other form of this
is we in the previous equation we can see that instead of multiplying by the inverse
of k we can do the inverse of M. So, that case we will be having the following;
form of the this is M inverse k this is also a standard Eigen value problem. And if this is a B matrix Eigen values of
whatever we will get if we take the square root of that that will be the natural frequency
of the system. Now, once we obtain the natural frequency
of the system there are corresponding mode shapes are also present. So, what we need to do once we obtained let
us say, n number of natural frequency of a particular system. So, 1 of the natural frequency will chose
first we will be going back to the previous equations. These equations in which yeah this particular
equation we will be substituting the natural frequency which we obtain here and with the
help of this, we will get the relative displacement this is the displacement vector. So, relative displacement of various nodes
will be getting it and that is nothing but, the mode shape or the Eigen vector of the
system. And for each natural frequency will be having
a unique or I should say a 1 set of Eigen vector all are relative because may be 1 of
the 1 of the Eigen vector value we have to choose and remaining will be normalize accordingly. So, for every natural frequency will be having
a mode shape.And now we will see for the problem which we are tackling how we can able to get
the various natural frequency and also we will see the convergence, as well increase
the number of element how the value of the natural frequency it convergence to the actual
value. So, now I am because these, Eigen value problem
can be solve using MATLAB very easily. So, now I am directly giving the Eigen value
for various number of elements, because here we have shown for 3 elements, but we can able
to easily increase the number of element with the procedure we have shown. I am showing the natural frequency comparison
of for let us say, mode 1, mode 2 and mode three. The exact solution, that I will be giving
the expression of that. So, let us say first is the finite element
method with 3 element we got 14 point 1 nine radian per second all unit is radian per second. First mode natural frequency, second mode
is 57.point 3, 9 and third is 141 point 6. And if, we increase the element to 6 it gives
us fourteen point 1, 8, 56 point 7, 7 128 point 1. And for 100 element I am increasing suddenly,
because once we have the program for the assembly procedure number of element is not a problem. 14 point 1,8 . So, hardly any improvement
with the increase in the element for the first 1 Second 1 slightly improved third 1 there
is some improvement, but now it is quite close to the previous 1. And if we see the exact formulation natural frequency for simply supported beam is given as n square pi square EI by m l 4. Where n is the mode number 1, 2, 3 and we have m as the mass per unit length row A and here E and the other things we already mentioned previously. So, if we substitute this value here we will get a 14 point 1 eight n square n is the mode number. So, for n is equal to 1 that is: with 6 element we are getting with quite close to the close form solution. So, close form solution, I am writing here that is exact is 14 point 1,28 for node 1 mode 1 and for mode 2 it is 50 for mode 2 it is 56.12 and for third mode it is 127 point six. So, you see that when we took 3 element the first natural frequency was quite close to the exact 1, Second 1 is marginally, but third 1 there was larger large difference, but with 6 element we could able to reach quite close to this and there is marginal improvement with number of elements. So, even the 6 element is quite accurate for this particular case. Now, we will see how extract the Eigen vector from Eigen vector or the mode shape from the Eigen vector because for if we are using a MATLAB corresponding to each natural frequency or the Eigen value there will be a Eigen vector. And I am just showing for 3 element case, how the for corresponding to the first natural frequency Eigen vector look like and then how to extract the mode shape from that.We had the for 3 element these were the displacement and slopes after applying the boundary conditions, that was the 6 by 6 matrix and the the vector was 6 into 1. So, the first column of the Eigen vector matrix which we get from the standard Eigen values routines, subroutines you will find these kind of numbers. First column will be having 6 entries if you are using 3 elements and these entries are in the number form, but we have to interpret them what are they this is positive. So, you can see that whatever the vector we had. So, there indicating what are, the values And for mode shapes we know that we need the displacement and w 1 we already know that because of this safety supported end condition they are 0.And w 4 is also 0. So, basically;
we extracted the mode shape which is the displacements as minus point 4 2 minus point 4, 2 and 0. So, if we want to plot them this node 1, 2, 3,4. So, you can see that because negative
sign if we ignore these 2 are 1 4 2, 1 4 2 these are 0.So, this is something if join
them by quadratic or the cubic curve the shape function will be getting a this kind of curve
which is nothing but, half sin wave. And we know that is simply supported beam the first
mode is half sin wave. So, similarly we can go for the second column and we can extract
the mode shape corresponding to the second natural frequency and similarly for the third
and fourth and so on. So, these extraction of the mode shape up to various modes already
showing in the a platform. So, this is the a mode shape of a particular
different beam having think looking repeated this yeah.So, now I will be showing the mode
shape for simply supported beam a for various number of modes. So, this is a mode shape
part which shows up to fifth mode for simply supported beam. So, you can see that with
a line with circle is this representing the first mode. And there are various other modes
second mode is by this star. So, this is the second mode in which the full sin wave is
there first it was the half sin wave but, this is the full sin wave. Similarly, we can
able to see the third fourth and fifth. So, we have already seen how we can able to obtain
the natural frequency and mode shape for a particular beam when it is undergoing transverse
vibration. Now, we will see another example in which the similar simply supported a shaft
is there is carrying a disk and it is having some unbalanced. How we can able to obtain
the unbalanced response and because the response or displacement of a rotor is is a vector
quantity. So, it will be having magnitude and phase. So, this is the aim this particular
shaft which is simply supported at ends. The dimension of this is exactly same as the previous
example in which where 3 meter beam, only difference here is here also I am dividing
this into 3 element, but we are keeping a disk at the third element. And this disk is
having mass of 1 point 5 kg it is having unbalanced mass of point 0 0 5 kg at the radius of point
0 0 5 point 0 5 meter and the phase of that with respect to some reference point on the
shaft is 30 degree. And here we have the phase something like this in which let say this
is the y axis, this is the z axis unbalance of the shaft is here. So, this is the phase
with respect to. So, we are measuring the phase with respect
to horizontal axis y. So, now we have described the problem. Now, we need to develop the elemental
equation for 3 element. Slight difference will be there in this present case as compared
to the previous 1, because one thing is there that. Now, because of unbalanced force the
motion of the beam or the shaft in horizontal plane and vertical planes are now coupled;
that means, we need to consider the equation of motion of horizontal plane and vertical
plane together in the previous free-vibration analysis we consider the 1 plane motion only.
So, for that let us see the elemental equations and how we can able to get the unbalanced
response from that. See this is the elemental equation of element
1 and here we have consider both the plane motions we can able to see that here we have
the displacement slope in 1 plane then in the same plane of the same node because node
1. Now, has 2 displacement 1 in horizontal direction another in the vertical direction
similarly slope. And similarly at the node 2 will be having 4 coordinates 2 for linear
displacement and 2 for the angular displacement and this is the rearrangement when we assemble
the 2 plane motion they are coming from the elemental equation. This is for the mass matrix
and this is Stiffness matrix and here this is the external force which is 0, because
on element 1 there is no eccentricity, we consider the eccentricity in the disk only
which is we are assuming to be at element 3 these are the internal forces.
And for the second element, because we are considering the disk in the second, third
element. So, second element will remain same as the element 1 only thing in place of node
1 and 2 in element 2 the node numbers will be 2 and 3. But, to the form of the mass matrix
and this Stiffness matrix will remains same. But, there will be change in the third element.
So, let us see the how the third element take the form. So, this is the element finite element equation
for the third element and you can able to see that there is a contribution of the disk
mass which is coming here and here corresponding to the displacement in node 3 and in the horizontal
and vertical directions. So, they are corresponding to these. So, 1 point 5 was the mass of the
disk and because this is common here this term: . So, we need to divide that by that
amount. So, this is the contribution of the disk Stiffness matrix is not getting change
because disk is not affecting the Stiffness external force is getting a change because
now this represent the unbalanced mass So, this is and where Omega is the spin speed of shaft
and this particular quantity is j theta. Theta is phase. So, this term comes from the phase
term. And you can see that the force in y direction which is this is 1 and force in
the z direction which is here if rotor is rotating counter clockwise direction. So,
z will be lagging by ninety degree with respect to y. So, we will be having relationship between
the Fz is equal to minus j Fy because minus j represent the phase difference between the
y Fy and Fz is 90 degree and Fz is lagging Fy by 90 degree. So, the same thing is appearing
here we are multiplying by minus j because this is the Fz force which is lagging Fy by
90 degree. This is the internal forces corresponding to node 3 and 4. So, now, we have obtain the
we assemble the we have obtained, the elemental again. Now we have obtain the elemental equation
of all the 3 element. Now we need to assemble them and apply the boundary condition so that
we can get the system equation from where we can obtained, the unbalanced response. So, this is. So, this is a mass matrix of
the assemble equation in which we have all the 3 element equations combined including
the disk also. So, next is the Stiffness matrix this is stiffness
matrix is also of all the 3 elements. And this is the total external force vector and
these are the internal forces you can see that the internal forces at the common nodes
will be becoming 0. Now, we will be applying the a boundary conditions
we know that what are the boundary conditions we have that is at node 1 we have w 1 and
v 1,0. And at last node we have both the displacement in the horizontal and vertical direction 0.Similarly,
moments in the 2 directions at node 1 is 0 and in node 4 is also 0. So, if we apply this
boundary conditions here, this is the reduced form of the equation. So, once we apply the
boundary condition, the system equation will get reduced. I am showing only this Stiffness
matrix how it looks like. So, this is this Stiffness matrix in the reduced form similarly
we can have the reduction in the mass matrix. This is the external force of all the internal
forces we are keeping 0 here all other non-zero internal forces we those equation we eliminated
as in the previous case. So, once we apply the boundary condition the
assembled matrix will reduce to a smaller matrix In this particular case because we
have these boundary conditions. we can able to see that correspondingly we will be having
only these second, 4, 5, 6 and 7, 8 up to 12,14 and 16 columns and rows will be retaining
others will eliminate. Because, if you see this particular Stiffness matrix is having
16 into 16 size. And we are eliminating 4. So, there will be total 12 into 12 matrix
of the mass and Stiffness will remain and others will be getting eliminated. So, what
we can do I am just showing the one of the stiff matrix that is the Stiffness matrix.
The size of this matrix is 12 into 12 we are already eliminated those rows and columns
which are giving us no contribution. Because of there are getting multiplied by 0.And in
this case here you can see these forces are acting at node 3 and both in the horizontal
and vertical direction. And now this particular equation is in this
form compact form which is standard dynamic equation and it contains the force also. In
previous case only free-vibration was there. So, this force term was there. Now, again
we can take the solution of this as amplitude and a harmonic function which is Omega is
nothing, but spin speed of the shaft. So, if we differentiate this twice we will get
minus Omega square w ej Omega t. And if you substitute this 2 in the equations of motion.
Here and here we will minus Omega square w. So, we will get minus omega square M plus
k matrix w bar is common which is amplitude? And we have the force again we are writing
as amplitude because the force can be written as this can be written as j Omega t.
Because this is coming from the unbalance and correspondingly we are choosing the response
of the same frequency. And what over phase information because for the present case there
is not damping. So, the phase between the force and the displacement will be same. But,
in case there is a damping then phase information will be contained in this and it will be then
complex in nature, but at present or because there is no damping. So, it it will remain
real. Now, this particular matrix I can able to call this as A matrix. So, we can write this is as A which function
of speed w is bar is equal to F bar. Now, to get the response we need to invert the
A matrix A matrix. So, once we got this here you can see that Omega is the variable. So,
Omega can be change from 0 up to some Omega max. So, 1 by 1 we need change Omega if we
take Omega is equal let us say, 100 rpm. We substitutes on this and then we know the force
we can get the response for that Omega and let us say, 100 rpm. Then we can change this
again we have to obtain, the response corresponding to this separately by its different run, and
so on, we can get the these w bars for various value of omega the spin speed.
Today we have seen how to the application of the finite element method for finding the
natural frequency and mode shape of a simply supported beam. In this example: we have taken
only 3 element of the beam and we show how to write the elemental equation and then how
to assemble those equations and then after getting the assemble equation the application
of the boundary condition, and how to get the reduced form of the system equation which
can be used for obtaining the natural frequency and mode shape. Subsequently we apply the
apply the finite element method for finding the unbalanced response of a rotor system
in which there was 1 disk having some eccentricity and through detailed steps we showed how the
unbalanced response can be obtained for the a rotor system.