Tip:
Highlight text to annotate it
X
We have been talking about the vibration problem till now. So let us continue our discussion
further. We had said that the certain properties of the Eigen value problem that we formed
in the free vibration case.
One is that the finite element that the Eigen values or the natural frequencies obtained
from the finite element solutions are higher than those obtained by the exact solution
or the actual modes. This is sometimes referred to as the finite element system being stiffer
than the exact one. This is some feature of the finite element solution that one has to
keep in mind. Why do we have to keep it in mind? Because, if I am writing a program and
to do the computation of the natural frequencies of a system, which is the very simple extension
of what we had been developing till now. Then as a benchmark, I should test the program
for specific problems for which we can construct the exact solution. For example, the problem
that I have been considering, here in this case, when I go and use the finite element
approximation and get the natural frequencies, I should get the natural frequencies to be
higher than that obtained from the exact one and as the mesh gets refined, they should
converge to the exact one from above. These are certain features which have to be kept
in mind. If it does not happen then something is wrong in the way I have done the programming
or in the way I have developed the finite element approximation.
Let us also look at some implementation issues. We had said that here I have taken this example
of a 5 node system and we had let us say that we are taking a piece wise linear approximation.
I can again take higher order approximations just like we have done the standard finite
element computation. Here this is the element 1. This is the element 2, element 3 and element
4. So just like we did the procedure for static analysis, here we have to obtain the global
matrices K and M. We have to obtain these global matrices from our finite element computation.
Again I will say that look what is this K was as I had done ‘0’ to ‘L’ integral
of EA phii prime phij prime dx was Kij and this can now be written as we have done earlier.
That is L is equal to 1, 2, in this case 4, the number of elements ‘xl’ to ‘xl+1’
EA phii prime phij prime. I can break the integral into sum of the integral over each
of these elements.
Which of the phii’s are active in the particular element? The only phi which is in the active
in this element, let us say element 3 is phi3 and phi4. Exactly what we have doing till
now will now be done and we will compute these matrices at the element level, assemble them
in the global matrices and that, we have done. But we have to do for the mass matrix also,
where the integral will be integral of rho phii phij.
So we have to now compute the Element Level Stiffness and Mass Matrices. How do I go about
doing it? Let us see. So we have our generic element with node l, l+1. This is the element
Il. Here the basis function from the figure I have drawn which is non zero. The basis
function will be phil and phil+1. These will correspond to the element level N1 of the
element and this will be N2 of the element. If I want to write the U(x), the approximation
of it in the element Il will be equal to alphal phil plus alphal+1 phil+1. These things will
not change. These things remains the same as what we have been doing till now. We are
going to essentially find the effect of these in the integrals.
We will have integral at the element level xl to xl+1 EA Ni for the element l prime Nj
for the element l prime dx, prime means the derivative with respect to x. This is equal
to Kij for the element l. i is equal to 1 to 2. In fact, if we use higher order approximation,
it goes up to P+1 and j is equal to 1 to 2. So this will give me 2 by 2 matrix in the
case of the element. This will give me K for the element l. Similarly, I will have integral
from xl to xl+1 rho Ni for the element, Nj for the element, dx. This is equal to Mij
for the element l. What is the assembly procedure? Assembly procedure will be this i locally
at the element level corresponds to some global I and similarly, the j at the element level
corresponds to some global J. So i = 1 corresponds to l, i = 2 corresponds to l+1 globally. j
= 1 corresponds to l globally and j = 2 corresponds to l+1 globally. That is I again use the local to global enumeration
routine that we had made long ago in the 1 dimensional code to give me this information
and then I will go and assemble in the following way.
Kll will be equal to, these are the global level, Kll plus K11 from the element. Kl,
l+1 will be equal to Kl, l+1 plus K12 for the element. Kl+1, l will be Kl+1, l plus
K21 for the element. Kl+1, l+1 will be equal to Kl+1, l+1 plus K22 for the element. This
is the simple assembly procedure for the K matrix is exactly the same as we what we had
done and similarly, we will do for the mass matrix. The assembly that is Mll is equal
to Mll plus M11 l and Ml, l+1 is equal to Ml, l+1 plus M12 l and Ml+1, l is equal to
Ml+1, l plus M21 l and finally, Ml+1, l+1 is equal to Ml+1, l+1 plus M22 for the element.
We do this assemble procedure. We loop our all the element and we construct the global
K and the M matrices.
Now the question is that the K and the M, if I take again our problem that we had considered,
the mesh that we have taken, this will be a 5 by 5 system, but from the boundary condition
at this end, I have to put alpha1 equal to 0 because the U(x) has to be 0 at the point
x equal to 0. So how do I enforce alpha1 equal to 0?
One way is to actually take the global systems and eliminate alpha1. Set alpha1 equal to
0 from boundary condition. We can do that reduce this 5 by 5 system to 4 by 4 system
that is that remove the variable alpha1 and the first equation. We actually solving for
the system which will have K22, K23, K24, K25 all the way down to K52, K53, K54, K55,
this into alpha2, alpha3, alpha4, alpha5. This will be equal to lambda into M22 M23
and so on, M25 M52 M53 and so on up to M55 into alpha. So this will be the reduce system
or we can use the tool that we have developed for doing the elimination.
We go to the 5 by 5 system. How did we set alpha equal to 0 on both sides? We simply
eliminate. We will do this, set this term to 1, first entry, then make it 0 then I will
have K22 K23 K24 K25 and so on, into the alpha array, alpha vector is equal to lambda into,
I will put again here 1 0, 0, 0, 0, 0, 0, 0, 0 and then M22 and so on. First row gets
completely decoupled and we get essentially alpha1 is equal to lambda alpha1. This will
give me lambda is equal to 1. So we will get 1 Eigen value as 1, throw it out. We will
not consider this Eigen value in our system of natural frequencies. We will throw all
these Eigen values out and the remaining 4 Eigen values will be the ones, we are going
to take. We will take lambda2, lambda3, lambda4, and lambda5. This will become our modes. This
become the first mode, this becomes second mode and so on. This will actually lead to
our omega1, this will lead to omega2, this will lead to omega3 and this will be omega4.
We should check that indeed the Eigen values are the same as what we would have obtained
by doing this. These are the next 1 would be equivalent, but throw this Eigen value
out on a certain equations. With this, we should be in a position to set up the Eigen
value problem and solve it.
One more thing that we would like to talk of is again while doing this integration,
we go to the master element, use all the two set we have developed as for as integration
over the master element is concerned. Do the numerical integrations. Create the matrices
and we will get a solution.
Let us go back and see what we had done as far as the exact solution is concerned. We
had these functions as our modes corresponding to the natural frequencies omega1, omega2,
omega3 and omega4 and we set these modes are oscillatory in nature and we set keep a note
of that because this was needed in what we do.
The first mode shape is like this. Second one is like this. Third one is like this.
Imagine, let us go back and recreate this picture in what we are doing.
Let us come here. This is x is equal to l and here is the exact U1 as the function of
x corresponding to omega1. What will be the finite element approximation do if I took
the one element solution? That is if I took a mesh with only one element. This is l. This
is 0. I simply put one element. In that case, we will see that we will get a 1 by 1 system
because one element will have the two nodes. This node will be knocked off will be left
only in terms of this one and the solution may look some thing like this. This will be
we get as U1FE from the finite element solution. It will be a straight line. This curve is
quite a nice curve. So this will not be such a bad approximation and we can see that if
we solve the one element problem, that is a single degree of problem, we will get an
Eigen value which is a very close to the exact point.
In fact, we can do it here. Let us do it here. Here if this is 0, this is l the first shape
function will be essentially this one and this one. This shape function will be x by
L and this one will be 1-(x by L). So we are only interested in this. uFE as the function
of x will be equal to alpha x over L. Integral 0 to L, EA, the derivative of x by L. So 1
by L into 1 by L, dx is equal to K11. This will be equal to EA by L. Similarly, the mass
matrix will be 0 to L, rho into x squared by L squared dx. This will be equal to rho
by L squared into x cubed by 3. So it will be (1 by 3) L3. So this will be equal to rho
L by 3. We get this is M11. If I do this problem, I will get K11 into alpha is equal to lambda
M11 into alpha. So lambda will be equal to K11 by M11 which will be equal to EA by L
divided by rho L by 3. This will be equal to 3 EA by rho L squared.
Lambda will be 3 EA by rho L squared, lambda from the finite element solution. So this
is what we will get for that problem. Once we have this, now what about the alpha? The
alpha can now be computed quite easily but let us look at this lambda. What was the lambda
exact?
If I go back to the earlier result, lambda exact will be equal to, it will be essentially
the square of this, because the way we have done omega squared was our lambda. So let
us go to omega and so it will be 2 into n minus 1, into pi by 2L into EA by rho, whole
squared pi by 2L into EA by rho whole squared.
That will be pi by 2L into EA by rho whole thing squared. This will be equal to pi squared
by four into EA by rho L squared. If we see that pi squared will be approximately 10,
so 10 by 4 and this will be approximately 2.5 and this is approximately 2.5 EA by rho
L2. This is the lambda exact. The lambda exact is less than the lambda, we have obtained
and the natural frequency will be square root of this. This is quite close as an approximation.
Not bad with one element, we could create a solution which is quite close to the one
which we actually needed.
Now let us go back and look at the two term solution, if I did. What will happen, if I
did the two term solution? Put this, what will happen if I do the two term solution
that is when I put two elements? This approximation of the mode shape will become better. As this
becomes better, I will see that the lambda1 FE will come much closer to lambda1 exact
and further, if I took this 1 term solution, I just cannot expect to do a good job for
the second mode. There is no second mode which I can capture. With the two term solution
now, I will have two Eigen values. Two frequencies can seemingly be approximated with the two
term solution. What was the second mode? The second mode would do this. In the two term
solution, for the one mode, very good, for the second mode, I will do most probably,
I will do this. This is certainly not a very good approximation of the second mode, but
nevertheless, it is an approximation and we will see that when I took two element mesh,
the lambda1FE will be very close to the lambda1 that I exact but the lambda2 FE will not be
as close as to the lambda2 exact but it will again be higher, not so bad, but badness will
be certainly a mode as compared to the lambda1. As we refined the mesh further and further,
the quality of the approximation of the natural frequency or the Eigen value corresponding
to this problem becomes better and better.
With mesh size h, the lambdaFE or in that sense omegaFE tends to the omega exact. There
is a very interesting way of measuring it and that is given by if we say in this case,
where we have distinct Eigen values. Lambdai minus lambdaiFE due to the finite element
solution for a given mesh, I could have hundred elements in the mesh or I could have five
elements in the mesh but as the mesh is refined that is has the h changes, I am using the
h method of the finite element method, this will be less than equal to some constant into
the strain energy of Ui exact minus UiFE. This is the very important result. This U
is the strain energy which means that is given by integral 0 to L, EA into Ui exact prime
minus UiFE prime whole squared dx. So this depends on the strain energy of the error,
this is called the strain energy of the error in the Eigen mode or in the Eigen function.
How well can finite element method approximate the corresponding mode also reflects on how
nicely can we obtain the corresponding frequency?
This number is dependent on this squared and here, in the case of Eigen value problem,
essentially for such a case, lambdai exact minus lambdai , actually I should have it
other way now but since I have taken the absolute value. It is something C h to the power of
2P. C will change for the corresponding mode that is as we go to the higher and higher
frequencies, this constant will degenerate further. We will get bigger and bigger as
the weight of the convergence will be this. We get twice the rate of convergence of the
normal finite element solution. We had talked about the energy norm that is square root
of this strain energy. Here, the error in the Eigen value converges as the strain energy
itself. Eigen value problems are much more forgiving as compared to getting information
about the point by stresses point by strains. We can obtain our Eigen modes with the much
coarse mesh, pretty good values of the Eigen values or the natural frequency can come with
the much coarse mesh as compared to what will be needed to get point wise estimation of
the stresses and so on and other issue.
This is one thing that one should keep in mind that depending on the particular problem,
we have to design the approximation in such a way that we get accurate results and we
should know something more about how the solution behaves, in order to be confident of what
we are computing and what we are giving to the particular user. Incidentally, this same
feature what we had discussed here also works in the two dimensional problem. Let us say,
I am talking of a domain with the crack, a crack domain. Here also I may be interested
in solving the elasticity problem. I may be interested in getting the natural frequencies.
In that case, what we will get is essentially the nature of solution and how would we do
the approximation in the vicinity of the crack tip is going to effect the quality of the
natural frequencies that we compute for the crack rate. So the meshing in the vicinity
of the crack tip becomes very important. We can not ignore it and in fact, in this case,
lambdai exact minus lambdaiFE will go essentially as some constant h to the power of minimum
of 2 alpha, P.
It depends on whether alpha is smaller P is small and in general, the crack leads to where
alpha is the exponent of the leading singular part of the solution. In the case of the crack,
it is essentially half. Any P take will be certainly greater than half. So alpha will
dominate that is the nature of the singularity at the crack tip is going to dominate the
rate of conversion and then we have to do all the things that we did for the stress
computation to get the good solution but nevertheless, the rate of convergence has the two in the
power. It is certainly better than the point wise data. Here, we are talking lambda as
a global data.
Similarly, Eigen value problems can also be there in the case of buckling, buckling of
plates, beams etc. In this case, as we take the P against deflection for a beam, let us
say the tip deflection then where I have how is the load applied. I have a compressive
load being applied to the structure this center. This compressive load should lead to only
an actual deformation but they comes a value if I keep on increasing the load, that comes
the value of the load at which the beam or this bar does not deflect inclined. It wants
to do this. It tries to have a transverse deflection. How do you take care of that?
For that we will look at equilibrium in the deform configuration and there is a whole
theory behind how to do the analysis. Nevertheless, as the P is increases, this is the transverse
deflection. Let us say it is of the tip. This increases that comes the critical value then
it falls. This is essentially what we call as a buckling load of the critical and we
try to approximate it through this is solving a non-linear problem. We try to approximate
it using a linearized problem. In the non-linear problem, we use concepts from finite deformation
elasticity and all these things, finite elasticity but we are not going to use it here.
The idealization is based on assuming that up to the critical load, the member essentially
behaves in a linear elastic way that is it is only going to give actual deformations,
negative deformation because of the compression.
At the critical load, it develops, transverse deflections also and without going into the
details of it, we will get that let us say that I am applying a load P is just less than
P is equal to P critical minus a very small number, in that case equilibrium under the
action of the load for this member is governed by standard linear elasticity that is I will get integral 0 to
L in the energy sense EA du divided by dx. I am just writing the weak form of the barred
problem. dv divided by dx into dx is equal to P into v at L or minus P into v at L. This
is what drives the problem.
At the critical load, at this point, increases the load little bit so that I get P critical.
At that point, this strain we had assumed is given by du divided by dx, not any more,
at this point, the strain will be governed by the non-linear or the finite strain or
the green strain. It will be given as essentially du divided by dx plus as an approximation
half of dw divided by dx whole squared where w is the transverse deflection. This is essentially
by throwing out some higher order terms. This is called Von Karman theory approach. Cut
anywhere the actual forces P and then we can show that we will have when we increase the
load little bit. This part is assumed that I am in equilibrium corresponding to this
load. From here, I am applying the extra bit that is I am applying the extra load to go
to the new equilibrium position. The new equilibrium position given this will have an additional
displacement in the plane and the actual direction and the transverse direction. This is generally
given in terms of what we call as perturbation displacements. From the current equilibrium
configuration corresponding to this and then it can be shown that this whole thing can
be written in terms of
the following weak form. This will be equal to P integral 0 to L into dw P divided by
dx into dw bar P dx.
This essentially gives if I take now w P as the incremental or the perturbational transverse
deflection, write it using the finite element representation using the c 1 elements, this
one will be approximated using the finite element way and this will be corresponding
test function or the rate function and this is again on the right hand side, I have the
unknown w P and here is the test function. This will lead again to an Eigen value problem
of the following type.
The K, this is for the bending into alpha is equal to P into what we call as K G into
alpha, where is this coming from, this is called geometric stiffness matrix or geometric
part, in the sense, that it is coming from this part of the strain that is by representing
the large strain components. The entries of K, Kij is equal to integral 0 to L, EI phii
double prime, phij double prime and Kij geometric is equal to integral of 0 to L. I will go
back and see what do I have phii prime, phij prime dx. I am not writing these things in
detail. I am really interested in finding this load P actually this happen at which
the transverse deflection is present. There again I setup this. Now phii or c 1 functions,
this is very important that these will be c 1, because we want phii double prime phij
double prime to be defined.
We do the same thing. Again this is the where the w P as the function of x is given as sum
alphai phii . i is equal to 1 to some number of unknowns N. This gives me the Eigen value
problem where lambda is equal to P that is the buckling and that is the load unknown
load and then the critical load out of this is essentially equal to the minimum Eigen
value that we have. That is what I am looking at. I would like to obtain the minimum Eigen
value corresponding to this problem. Here again I get P critical is equal to the minimum
Eigen value and how many Eigen values will get for this system will get N by N depending
on the number of unknowns.
Similarly, if I went and solve this in the continual form that is we will get, if we
do the formulation of this in an exact form P.
This would be the exact differential equation that we have and again if we solve using what
we had already used for the free vibration problem, I will actually get infinite values
of P because this is an Eigen value problem at the continual level, infinite values of
P. They will be distinct value of P and corresponding to each P, I will get a buckled mode shape.
That is exactly what we do in our buckling analysis and if I will take a rod like this
and apply to this, a compressive load, it deforms. It will deform like this. So the
rod will become this at P critical. When the load becomes P critical, the rod will become
this.
As we had talked about earlier here, we have lambda1, lambda2 up to lambda N such that
lambda1 is less than equal to lambda2 is less than equal to lambda3 and so on. The first
Eigen value is the Eigen value I am interesting in. This is very important also from an implementation
point of view because in this analysis, we are only interested most of only interest
in the critical load. I am interested only the first Eigen value. If I am interested
in the first Eigen value, I really do not need to compute for the remaining Eigen values.
If I do not need to compute for the remaining Eigen values, I can use faster solvers which
give me the first Eigen value very nicely. The cost of computation comes down tremendously,
if I am interested in only the first few Eigen value and I am not demanding all the Eigen
values.
So as an implementational issue computation cost, for the Eigen value problems, this is
the big issue. When I use I want to get information about the first few Eigen values or one Eigen
value and the cost comes down tremendously. What will be the accuracy of the P critical
that we obtained? Again, it will be determined by essentially, this bending strain energy
of the perturbation transverse displacement w P. How well are we approximating w P that
is how well is the finite element solution going to represent this mode? With one element,
I will get this, so exactly the same theory. Certainly, there are certain differences but
the same ideas can be extended to the Eigen value problem under computation.
There are lots of Eigen values solvers available in the open source. We can download it from
the net or numerical recipes in Fortran and C then so many books on those they also have
Eigen value solvers which are generally a little bit more difficult to make program
and finally make it useful for computation as compared to the solvers, we had for the
inversion of the matrices. From the open source, we can have many such algorithms which are
available. We can go to a site called netlib.org that is an open source site, we can download
lots of bookings from there. There is something called ARPACK which is very powerful, generalize
Eigen value solvers package available again as a free open source package. It can be used
to compute Eigen values of large systems. Here we did not really say what the size of
the system is.
Depending on the problem, 1-D problem, meshes are small. Approximation is small. When I
want to do this Eigen value problem in 2-D or 3-D for example, in the 2D, very important
thing is when I am interested in plate buckling, in this case, the plate is subjected to a
compressive load or even to a shear load. Here is the plate subjected to some shear
loads. In that case, I am interested in the critical value of this load. In that case,
if the plate has cut outs everywhere, crack sitting some where at various sites. I have
to resolve these geometrical details properly so that the finite element solution corresponding
to the buckled mode shape is not bad. The finite element solution has to be good in
order to have a good value of the P critical. We have to do the meshing which is detailed
enough to take care of these details. This leads to a very large sparse Eigen value problem.
The mass matrix and the K and the stiffness matrix could be 20000, 30000 or 40,000 by
40,000 size. Then getting the Eigen values becomes the tough job and some of these algorithms
have to be used. So there is the lot of literature available on how to simply solve a large sparse
Eigen value problem. We can use straight away from the literature.
In this lecture, we are trying to cover some classes of Eigen value problems that arise
in mechanics. Some typical example that we have taken is the free vibration problem or
the buckling problem, both of which give rise to Eigen value problems which have to be solve
first and formulated in the finite element sense and then they have to be solved and
then check whether those numbers are reliable or not. Benchmark problems have to be solved
and for the students it will not be very difficult to program this in the existing static analysis
codes because we already have the structure with which we have the shape function routines,
we have our element calculation routines, where all we have to do is add this extra
parts which may be for computation of the mass matrix or computation of the geometric
stiffness matrix add them. Do follow the same procedure that we are following for the assembly.
Assemble these matrices. Obtain, download or write whichever way we like to do. Get
an Eigen value solver. Put it and check. This way, we can now expand our codes that we have
developed the 1D and the 2D code, to also take care of Eigen value problems and so we
can create much more information which is useful to a designer or to an analyst.
In the next lecture, we are going to look at some problems which fall in the category
of non-linear problems. Till now we have only dealt with linear problems. We will add non
linearity in to the system which is not forcibly done but this is where the real situation
is actually, most of the problems in real life are non-linear in nature. Those nonlinearities,
if we add to the system then what happens to the response of the system in terms of
the finite element solutions, that is how do I go ahead and do a finite element computation,
which will give me the information that I need correctly? So we are going to look at
finite element formulation for nonlinear problems and from there we will start. A very simple
example is this. I take a bar again a stick to a bar that we had. Here is P. Here is f(x)
and I will now make it lie on elastic supports on either side such that this can be represented by distributed
actual spring which can be given with the spring constant k(x). We had taken it earlier
as k0. I will add to it to a part which is this.
The spring now has the nonlinear part in the definition of the spring constant. How do
I take care of the response of the bar with this kind of a support spring? Example, this
could be the effect of the rubbery support, where the rubber behaves in the nonlinear
way. We should try to find its information will be respect to the force which is applied.
This could be a model of rubbery spring, something like that. We could add more complexities
but primarily, we would take this problem as our problem for which we will develop the
non linear analysis tools. We will solve the problem and we see, how well or how badly
do we do with respect to, what is the exact solution of this problem.