Tip:
Highlight text to annotate it
X
Today we will see the Galerkin's approach to the one dimensional finite element problems
and in the Galerkin's method, we have already seen that this is the equation that we get
for the virtual work done. That is this part is the virtual work done against the internal
forces and this is the external virtual work done and this two will be equal in magnitude
and the total virtual work done will be 0. This we have already seen.
Now if we talk of one dimensional problems, we formulated one dimensional problem like
this that if we take a body like this, we are splitting it into a set of 1 D elements
like this. And if you take any individual element of uniform cross section having two
nodes 1 and 2, we said that this has the displacement given by q1 q2 or the element displacement
vector we said was q1 q2 which we are representing as the vector q or the matrix q.
Now in this Galerkin's approach we are giving it a virtual displacement. So if we give this
a virtual displacement, let's say the virtual displacement at this end is phi1 and the virtual
displacement at this end is phi2. And we will make the virtual displacement vector just
like this which will be phi1, phi2 transpose and we will call that let's say the matrix
psi. If you use this and we want to find out the virtual displacement at any point inside this element, if I take
any arbitrary point here, if I want to find out the actual displacement we have written
the equations earlier we said u will be equal to N times q where N is nothing but N1 N2
where N1 and N2 are the shape functions.
So we have said u is equal to N times q. Similarly, if you want to find out the virtual displacement
at any point inside this element, let's say the virtual displacement here is phi. We will
again say that this phi will be a linear combination of phi1 and phi2. So we will say phi will
be equal to N times psi or this is same as N1 phi1 plus N2 phi2. This u equal to Nq that
is the same as N1 q1 plus N2 q2. So we will say phi will be equal to N time's psi. The
next step if you notice, the first term that we have here has a term of epsilon of phi
and we said that this epsilon, this strain is the strain that will be caused by the virtual
displacement that we are giving. So this strain caused by the virtual displacement will be
nothing but d phi by dx. That is if you look at the actual strain, this I am saying is
epsilon of phi.
If I look at the actual strain epsilon that is equal to du by dx and this we have derived
last time, if I differentiate this equation I will say this will be equal to B times q.
Similarly if I differentiate this phi equal to N time psi, I will get this will be equal
to B times psi and this derivation is straight forward. All that I have to do is differentiate
N1 and N2 with respect to x and this matrix B is the matrix given by
this expression. This matrix B is the same as the matrix here that we had and the derivation
for this is also absolutely the same. So epsilon of phi will be equal to B time's psi and epsilon
is equal to B time's q and if we say epsilon is equal to B times q, we can say sigma will
be equal to E times B times q where E is the Young's modulus, B is the element strain displacement
matrix and q is the element displacement matrix. So we have an expression for sigma and we
have an expression for epsilon of phi. If I take these two expressions, I can put them
into the first term that is this term which is which corresponds to the internal work
done.
And if I expand this term for the internal work done, my first term here which is the
volume integral of sigma transpose epsilon of phi dv, this will be equal to sigma is equal to Ebq. So sigma transpose
will be equal to q transpose B transpose multiplied by E which is a constant anyway, it's not
a matrix. So I will put that over here, so this will be equal to the integral of E q
transpose B transpose multiplied by epsilon of phi. Epsilon of phi is B time's psi, so
I will put that over here. B times psi, dv which is the differential volume element that
will be equal to area time's dx, so it will be area of the element multiplied by dx.
Now this integral will be integral over the length of the element. And now in this you
will notice that E is a constant, q is a constant B transpose B and psi and A they are all constants.
They don't change with x and dx, the integral of dx is nothing but the length of the element.
So if I put all these values here, this will be equal to E into A multiplied by integral
of dx which will give me le multiplied by q transpose B transpose B psi. Now again B
transpose, B is this expression, so B transpose will be 1 over le into minus 1 1, so I will
put these values over here. This will be equal to, I will get q transpose into EA le into
1 over le multiplied by minus 1 1 into 1 over le multiplied by minus 1 1 and this whole
thing will be multiplied by psi. This will be equal to q transpose into EA over le, I
have le square in the denominator, le in the numerator and the product of these two matrices
will giving 1 minus 1 minus 1 1. So this thing will be again multiplied by psi. And now this
expression is the same as the expression for ke which is the stiffness matrix. So this
will be equal to q transpose into ke into psi. Sir, what is q and what is psi. What
is q and what is psi? Student: This psi we have defined as phi1 phi2 virtual display.
So q is, q is the actual displacement.
See this is my element, at this point one the actual displacement is q1, at this point
two the actual displacement is q2. In addition to that I am giving a virtual displacement
this phi1 and phi2. q1 q2 give me a vector which is q, phi1 phi2 give me a vector which
I am calling a psi. Student: but the strains we have found out using psi. See using psi
I followed what, I am calling epsilon of phi, the strains which will be caused only due
to the virtual displacement. The actual strains because of the loading, their epsilon which
is nothing but du by dx. The strains caused due to the virtual displacement are d phi
by dx.
So when I give it a virtual displacement, the effort involved in giving the virtual
displacement will be doing some work against the internal stresses. Though internal stresses
are actually the stresses caused by the loading. So that internal work done is this term. That
is why I have a term of the actual stresses multiplied by the strain caused due to the
virtual displacement and the product of the two integrated over the volume will give me
the work done. So this internal work done is equal to q transpose ke into psi. This
is the internal work done within the element. Any question up to this point?
Now if you look at this expression for the Galerkin's method, the first term we have
simplified. Similarly you will simplify the other terms and then we will try to get a
set of equations for solving for q using this general equation. So let's take up the second
term now.
The second term is integral of, volume integral of phi transpose into f into dv. Now if I
take up this, again phi is equal to N time's psi, I am writing it directly. So phi transpose
will be psi transpose multiplied by N transpose, f is the body force within the element, dv
is A times dx where A is the area of cross section in the element. And again now this
integral would be over the length of the element. So this will be equal to, psi transpose is
a constant, f is a constant, A is a constant. This integral will give me psi transpose into
integral of N transpose dx and this whole thing will be multiplied by f and A. This
integral is over the length. And the integral of N transpose dx, we have derived this integral
earlier, I will quickly repeat it here. Psi transpose, N transpose will given me a column
vector of N1 N2 integral of that with dx and this column vector, we have multiplied by
fA. And again if you remember integral of N1 dx, we had said earlier that N1 is the
shape function which along the length of the element is changing from 1 to 0 like this.
At the point 1, it is equal to 1, at this point it is equal to z, at the point 2 it
is equal to 0. And the integral N1 dx will be the area under this curve and this integral
is equal to le by 2 base into height into half.
If we put that over here, this will be equal to psi transpose into le by 2 into 1 1 into
fA. So this is again equal to psi transpose into fA le by 2 into 1 1. Again we get a result
which is very similar to what we had last time. That is when you are using the Rayleigh-Ritz
method, that fA into le is the total body force acting on the element, f is the body
force per unit volume and A into le is the volume of the element. So this is the total
body force acting on the element, so we are saying half of that is acting on the first
node, the other half of that is acting on the second node and this we can write this
as psi transpose into let's say a body force vector Fe and effectively it means that if
you have an element like this. The total body force that is acting is splitted equally into
two parts, one acting on the first node, the other acting on the second node. And the virtual
work done will be the virtual displacement of the first node into this force plus the
virtual displacement of the second node into this force. That is the virtual work done
against the body force acting in this element. So this term can be simplified into psi transpose
times Fe.
So again if you look at this expression, this term again we have simplified this now. Similarly
we can simplify this term also. And I will leave it for you to derive this that the the
this term which is integral of phi transpose Tds, the surface integral of this. This will
be equal to psi transpose times Te where Te will be equal to T times le which is the total
tractive force acting on the element divided by 2 and that would be split equally between
the two nodes.
So effectively the total tractive force will also be split equally between the two nodes.
So the third term, this term will evaluate to this expression. If you notice the terms
that we are getting are quite similar to the terms that we had in the Rayleigh Ritz method
except that their instead of psi we were having q, that's the only difference. And now if
we take up this complete expression what I will get will be, I will just rewrite the
same thing first. Sigma transpose epsilon of phi dv minus volume integral of phi transpose
fdv minus surface integral of phi transpose Tds minus sigma phi transpose p is equal to
0.
Now this is the integral over the complete volume. So integral over the complete volume
will be the same as summation over the integral of each of the elements. So this would be,
this expression will give us, if I do a summation over each of the elements of the integral
for each element that will be the integral over the element sigma transpose epsilon of
phi dv minus summation over all the elements into the integral over the element of phi
transpose fdv minus summation over all the elements integral over the element surface
area phi transpose Tds minus summation of phi transpose p will be equal to 0. Now this
term, this term is what we have just shown to be q transpose ke times psi.
So this will, this will give me summation of q transpose ke time psi minus summation
of this term is what we just derived as psi transpose times Fe. We will say minus summation
of psi transpose Fe, this term is what we have just mentioned over here minus summation
of psi transpose times Te minus summation of phi transpose times P is equal to 0. And
again ke that is the stiffness matrix is the symmetric matrix. So just for the sake of
convenience, what we will do is this term q transpose ke psi, I will write that as what
I am saying is q transpose ke times psi.
This is the same as psi transpose ke q because ke is symmetric and I can take a transpose
without changing nothing. The product of these three is just a single number so I take the
transpose, I will get this and so this is the same as this. So I will use this expression
from now onwards and what we will get, this thing will now become summation of psi transpose
ke times q minus summation of psi transpose Fe minus summation of psi transpose Te minus
summation of sorry phi transpose times P.
Now if I look at these three terms, in each of the three terms I have a force matrix.
This is the force matrix multiplied by in the mutual displacement, this is also a force
matrix multiplied by the virtual displacement and here I have got point forces multiplied
by the virtual displacement at that point. And just as we did earlier, I can combine
these three terms into one single force matrix. So this would be, this is equal to 0. I can
rewrite this as summation of psi transpose ke times q minus summation of psi transpose
F where this F is the element body force matrix and I will say this will be equal to 0. And
if I follow the same method as I did earlier for converting this element matrices into
the global matrices, I can say that this is the same as the writing. I will convert this
psi into a global matrix phi, this ke into a global matrix k and this q into a global
matrix q.
And this again will be minus, this psi I am changing to global matrix phi and this F is
changed I am changing that into a global force matrix and I will say this is equal to 0 where
this global force matrix phi will be equal to phi1 phi2 till phin if I have n elements.
The definition of K, Q and F will be the same as what we had earlier. Global matrix K will
be assembled from the individual stiffness matrices just like we did earlier. The global
force matrix will be assembled from the individual force matrices and the point loads just as
we did that earlier. So we will get this equation for the Galerkin's approach. And this I can
write that as phi into K Q minus F will be equal to 0 where phi is the global, is the
global matrix of virtual displacements, Q is the global matrix, Q will be equal to q1
q2 till qn transpose as a global matrix for deformations, F is the global force matrix,
K is the global stiffness matrix. So we will get phi into K Q minus F will be equal to
0. And this is the equation that we will be using for the Galerkin's approach.
And we have already said that as for the Galerkin's method, this equation has to be satisfied
for all values of phi that means whatever be the virtual displacement I give to this
system, this equation should be satisfied that means any virtual displacement which
is consistent with the boundary conditions. We have if you remember earlier we said that
in the Galerkin's approach the first thing we were saying was that phi has to be has
to take the same form as q that is the virtual displacement and the actual displacement should
have a similar description.
And that we have ensured by saying that q is equal to q1 q2 transpose and psi is equal
to phi1 phi2 transpose and by saying that phi is equal to N time's psi and u is equal
to n time's q. So we are taking a similar formulation for phi as well as for u and the
other constant is that this phi has to be kinematically admissible. That means this
has to satisfy the boundary conditions. We can only give it a virtual displacement which
is consistent with the boundary conditions. So if we keep, if we keep this constant in
mind then this equation has to be satisfied for all values of phi. So this equation, we
can use for solving for Q and if it has to be satisfied for all values of phi, we can
essentially say kQ is equal to F and then attach the boundary conditions to that.
We attach the boundary conditions so that we will come back to the same equations as
we had for the Rayleigh Ritz method. So essentially we will see that whether we use the Rayleigh
Ritz method or the Galerkin's approach, the final equation that we are getting, they will
be the same sorry. The final expression for the stiffness matrix, force matrix and the
resulting values that we will get for Q they will be the same irrespective of whether we
use the Rayleigh Ritz's method or the Galerkin's approach. Any question on the Galerkin's approach?
So essentially if you have if you look at the method that we are using for the one dimensional
problems, let's say in the Galerkin's approach we started with this equation.
Then we simplified this first term and we got an expression like this, an expression
like this.
Sorry, this is the second term. For the first term we got an expression like, we got this
expression.
And for the second term we got this expression and a similar expression we got for the third
term that is for this term. Now whether you have one dimensional problem or two dimensional
problems or three dimensional problems, our aim will always be the same. We will take
expressions like this, simplify them using a set of shape functions. What we have done
is we took a set of shape functions as we had an element like this which are the one
dimensional element and we said that we will take two nodes on this element 1 and 2.
If you specify the displacements at these two q1 and q2, we will find out the displacement
at any arbitrary point by taking two shape functions N1 and N2. So we say u will be equal
to N time's q when N1 and N2 are two shape functions. Those shape functions are basically
deciding the contribution of q1 and q2 at this point. And then we are using this formulation
to simplify this expression, this expression as well as this expression. And we finally
simplify these expressions to simple matrix products and those matrix products are again
in terms of q's.
Now whether we have one dimensional element, 2 D elements or 3 D elements, my approach
is going to be the same. Let's say in 2 D, we will try to take higher triangular elements
or quadrilateral elements maybe with 3 nodes, 4 nodes or higher number of nodes. I have
shape functions for each of these nodes and using that we will simplify this. And finally
the expression that we will get will be, expression will be the same expressions, psi transpose
times Fe or this expression or there will be expressions of this type. The final expression
will always be the same irrespective of what is the shape of the element we are choosing.
Our aim will be that the final expression should be the same.
So basically when we say that we will take a particular type of element that is triangular
or quadrilateral with 3 nodes, 4 nodes or 6 nodes or 10 nodes or whatever we will finally
come to the same set of equations. The only thing that will be different will be how to
get K and how to get F. We will take different shape functions with every element and we
will get different methods of assembling the K matrix and assembling the F matrix but once
this K and the F matrices have been assembled, our method of solution is going to remain
the same. And the method of solution is one method that we have seen that is elimination
approach and we have one more approach that is called the penalty approach. So let's briefly
see the penalty approach and then we will see an example of the penalty approach may
be in the next class.
Now now when we were talking of the penalty approach, this is for solving the system of
equations. If you look at the potential energy expression that we have had earlier in the
Rayleigh Ritz method, we had said that the potential energy pi will be equal to half
of Q transpose kQ minus Q transpose F. And for this we had a let's say some body, around
this body we were considering some constraint that some node over here is fixed or has a
fixed displacement. Let's say this is node number 1 and again we consider the constraint
that Q1 is equal to a1. So what we will do is we will approximate this body by a body
like this. We will add a spring at this location. This spring node has a spring constant of
C and this node is 1, the rest of the body is the same.
Now if I give this end a fixed displacement equal to a1 and the displacement of this end
of this node that is node number 1 that is Q1. The deformation in the spring delta will
be how much. Q1 minus a1, because a1 and Q1 will be in the same direction. The actual
deformation of this spring will be Q1 minus a1 and the potential energy stored in the
spring that will be equal to let's say the u of the spring will be equal to half kx squared
which is equal to half c into Q1 minus a1 whole squared. If now I look at this complete
system, the potential energy in this complete system would be this potential energy plus
the potential energy of the spring.
So now my complete potential energy, I will write that as pi which will be half c into
Q1 minus a1 whole squared plus half Q transpose KQ minus Q transpose F. And essentially what
we will say is that if I take this spring and I make the spring very stiff, as C approaches
infinity this system would approach this system. As this spring becomes very stiff, this system
will be almost the same as this system because then I will essentially get Q1 equal to a1
because this spring is so stiff and whatever deformation I giving to this end Q1 will also
go through the same deformation. So we will now solve for this system and then take the
limit as c tends to infinity. So I have this expression, again I will write this as half
of c into Q1 minus a1 whole squared plus I expand this term which is half of the same as what I have written yesterday Q1 K11 Q1
plus Q1 K12 Q2 plus so on up to Q1 K1n Qn plus Q2 K21 Q1 plus Q2 K22 Q2 Q2 K2n Qn minus
Q transpose time F which will give me Q1 F1 plus Q2 F2. So this is the expression we have
for the potential energy.
If I take this expression and again apply the constraint that del pi by del Qi is equal
0. So first I will differentiate with respect to Q1 and I will get del pi by del Q1 and
what will that be equal to. When I differentiate this with respective
Q1, I will get this term multiplied by 2, I will get this term so on up to this term
and this column I will get this term so one up to this term. And just as plus I will get
a differential with, when I differentiate with respect to Q1 I will get a term from
here. So the term from here will give me c times Q1 minus c times a1. The term from here
will give me this plus K11 Q1 plus K12 Q2 plus K13 Q3 so on up to K1n Qn and then this
term will give me minus F1 so this minus F1 will be equal to 0.
Similarly when I differentiate with respect Q2 I will get this term, I will get two times
this term, I will get this term and so on and on this side I will get this term so on
and this term, just like we differentiated last time. And this is going to give me expressions
like K12 Q1 plus K22 Q2 plus K23 Q3 so on till K2n Qn minus F2 is equal to 0 and so
on. This way we will get a system of n equations. Now here n equations, not n minus 1 equations
because I am also differentiating with respect to Q1. Last time I had put Q1 equal to a1
and I had eliminated the variable Q1, I am not doing that now. I will now have a system
of n equations and if I write them in a matrix form in this equation, in the first equation
the coefficient of Q1 is c plus K11 and all the other terms have coefficients which are
the same as the stiffness matrix K12 K13 K1n and so on. Similarly in this equation all
the coefficients are K12 K22 K23 and so on. So when I write it in the matrix form, it
will become the coefficient of Q1 is K11 plus c.
So the first term here will become K11 plus c. The coefficient of Q2 is K12. The other
terms will remain the same K12 K13 so on till K1n. If I take out the second equation, I
will get K12 K22 K23 and so on. Now this will become K12 sorry K21 K22 K23 and so on up
to K2n and this way we will continue we will get Kn1 Kn2 Kn3 up to Knn and this thing is
multiplied by Q1 Q2 up to Qn and this will be equal to, the force term here is this F1
and c times a1. So the first term here will become F1 plus c times a1. The force term
here is only F2, the term second term will be F2 and all the other terms will be the
same. So this is a system of equation that will I now get. This again we will write in
the similar manner K prime Q prime will be equal to F prime but now K prime is obtained
by taking the same n by n matrix but in the diagonal element I am just adding the stiffness
of the spring that I have taken c. And F prime I am getting by just adding c times a1 in
the force term in the first location. So this is the system of equations that we will have
and in this system of equation, as I take the value of c to be very large, if I take
c to be very large I will finally approach the case where Q1 will be equal to a1. We
can see that mathematically also.
So in this penalty approach, what we are doing is we are taking the stiffness matrix and
in the diagonal element for the constraint which is, for the node which is constrained
we will add a very high stiffness, very high value c and corresponding force term we will
modify it by adding c times a1. Mind you we are talking of a constraint which is Q1 is
equal to a1. So I will I am taking a very high stiffness c, adding that in the stiffness
matrix and adding c times a1 in the force matrix. And then I can solve the system of
n equations and get the different values of Q1 Q2 till Qn, Q prime is the same Q. So this
is the penalty approach. Next time I will take a small example of this penalty approach
and then see how it can be used.