Module 2 Lecture 1 Finite Element Method

In the previous lecture, we had talked about the Rayleigh--Ritz method, which was one of the ways of getting an approximate solution to the boundary value problem of interest. What we had done is we had taken a simple problem of a bar subjected to an end load P of size 10 and an uniformly distributed body force, f(x) of intensity x. At x equal to 0, we had fixed the bar; that is this was the fixed boundary condition and at x equal to 1, we had P applied; that is, the force was applied at x equal to 1. For this problem, we wanted to get an approximate solution using the Rayleigh--Ritz method. So how did we proceed? We took the so-called functional which was nothing but the total potential energy for the structure, which was given as half integral 0 to 1 u prime squared dx minus P, u evaluated at 1. Remember, here we are taking the material to be such that EA is equal to 1. Minimization of this pi u, that is the first variation of pi, gave us integral 0 to 1 u prime variation of u prime dx minus P variation of u at x equal to 1 was equal to 0. Then we went ahead and used a series representation for u. We said we will take the so-called two terms solution u2x was equal to u1x plus u2x squared; this is the two-term solution that we took. Why did we take it of this type? Because, we wanted this u2x also to satisfy the specified geometric conditions at the point x is equal to 0. And as we see that here the x and x squared vanishes at the point x equal to 0. In order to find the coefficients u1 and u2 of this series solution, we went back to this variation form that we have written here. In that we took delta u was equal to delta u1 into x plus delta u2 into x square. When we substitute that now the delta u is something that is under our control. It is the variation of our virtual displacements that we talking about. I can choose first since it is under our control. We will set delta u2 is equal to 0; delta u1 is equal to 1 and we will get one equation. Then in the second one, we are going to set delta u1 is equal to 0 delta u2 equal to 1 and we will get the second equation. This we have done in the previous lecture; so I am not going to go there again. After doing all this we will get the solution to the problem that we have defined. Let me again write it; we will get integral 0 to 1 u2 prime into x prime dx is equal to integral 0 to 1 x into x dx minus plus 10 into x evaluated at 1. This is first equation we get and second equation we have is u2 prime x squared prime dx is equal to integral 0 to 1 x into x dx plus 10 into x squared evaluated at 1. So out of these two equations substituting for the representation for u2 prime I will get the solution to this problem. And what we got in the last lecture was u2x was 10.5833x minus 0.25x squared. This solution that we obtained for u2x is very close to the exact solution to the problem which was equal to 10.5x minus one sixth x cube. When we plot this we will find that this u2 x is quite close and for engineering accuracies that we desire this is a good enough solution. If we went ahead and took three term solution instead of u2x; if you took u3x, this was equal to u1 into x plus u2 in to x squared plus u3 into x cube. Then we see that by going through the whole procedure here we will get the exact solution back. What the Rayleigh--Ritz method has done for us is that by taking these appropriate functions which we had called as our basis functions, which are polynomial functions, x, x square and x cube we have been able to recover the exact solution for this particular problem exactly using a three term solution. Is this picture as rosy as what we have said using this example? The answer is no. Let us look at the second problem that we have posed. So here again we said that P is 10 units; at this point I am going to apply a concentrated load of size 20. This is the point x equal to half, this is the point x equal to 1, this is x equal to 0. If I go ahead and write the total potential energy for this problem, it will be 1/2 integral of 0 to 1 u prime square dx minus 20 u evaluated at x equal to 1/2 minus 10u evaluated at x equal to y. If I take this first variation of this pi, this will be equal to integral 0 to 1 u prime variation of u prime dx minus 20 variation of u at x equal to 1/2 minus 10 variation of u evaluated at x equal to 1 this whole thing is equal to 0. Again if I take a two term solution as an approximation to the u of the problem, so u2x again as done for the previous case is equal to u1 into x plus u2x squared. So again we go through the same steps as we had followed for the previous problem, and we can obtain the coefficients u1 and u2 of the exact solution of the approximate solution that we have obtained. For this problem let me first write the exact solution. Exact solution of this problem is equal to if we look at it, it is 30x for this part and 10x plus 10 for this part. If we look at this exact solution, the exact solution is piecewise linear. That is it is linear in the part 0 to 1/2 and linear in the part 1/2 to 1 with the slopes different in the two parts. If I look at the two terms solution for this problem, let us see what it is 35x minus 15x squared. Look at this solution it is certainly not a piecewise linear, it is a quadratic polynomial. If I look at the value of the solution at the point x is equal to 1, u of 2 at the point 1, this is equal to 20. If I look at the value of the exact solution at the point x is equal to is equal to 1 it is also 20. So things look to be quite good; answer is no. Look at the derivative of this solution; u of 2 prime evaluated at x is equal to 1 will be equal to 35 minus 30x evaluated at 1 which is equal to 5. What is the derivative of the exact solution? The derivative of exact solution is 10. We have an error of 5 in the derivative of the approximate two terms solution that we have obtained and this error is 100 %. In most of our computations we are not interested in the value of the solution; we are mostly interested in the derivatives of the solution because we want to obtain the strain as well as the stress information out of this computation. So we see that this is a disaster as far as the numerical solution concerned. Now again if I go back to the point 1/2 and obtain the u2 prime at x is equal to 1/2. This will be equal to, as we have obtained it will be 35 minus 30x evaluated at 1/2. So this will be 35 minus 15 which is equal to 20; this is one number. What about the derivative of the exact solution? Exact solution if I look at it is equal to the derivative of it 30 for x equal to half minus epsilon that is I am very close to x but coming at it by the left hand side. And it is equal to 10 at the point x is equal to half plus epsilon. That is, I am coming to the point x is equal to half from the right hand side. This is my point x equal to half derivative in this region is 30; derivative in this region is 10. What has the approximate solution given me? It has given me the average value of these two derivatives. Again the numerical solution is not able to capture the jump in the derivatives which is inherent in the approximation that we have made. Are these problems of interest to us or are these simply artifacts that I created to show that the Rayleigh--Ritz method does not function? Unfortunately, most of the engineering problems that we are interested in do have this kind of a feature. For example, I may have a very general problem of a bar with multiple concentrated loads coming due to fixtures which may be attached to this bar. So I may have concentrated load F1 here, F2 here, F3 here. Similarly I may have another problem where I will have a bar I can gave it end load P no problem but the bar is now made of different materials that is our one material here, another material here, another here another here. I may have EA1 here, EA2 here, EA3 here, and EA4 here. This may be a problem that is of interest to us, in all this cases that is in both these cases if I look at the exact solution it does something like this. That is the solution is so-called piecewise smooth and continuous. This is what the engineering problems have as a feature; because continuity is required otherwise my specimen is going to break. So most such problems where I have dissimilar materials or points loads applied at certain locations the solution is going to behave like this. I could have also put some distributed loads no problem, on the structure. That is not going to change the nature of the solution that is at these points of transition of the material are the points where the consecrated loads are applied, we will have a change in slopes that is inherent. These are the solutions that we are interested in and we have seen that our so-called two terms solutions that we have obtain for the model problem of interest that the second model problem that we had put; the two terms solutions something like this. This is the exact solution and if I look at my approximation, it does something like this. The question is what if I took three terms solutions? Three terms solution gives exactly the two term solution back that is cubic part of the solution if you go ahead and solve it is 0. Taking higher and higher terms really does not solve the problem. How do we solve this problem? For this very simple thing that we can do is, we can change the definition of these basis functions that we are using to represent the solution. If we remember that this we had written as i equal to 1to N ui phii of x. Till now we had taken these as polynomials. What if I go ahead and do the following thing for the module problem two that we have taken. This is x is equal to 0, x is equal to 1/2, x is equal to 1 what should I do? I do the following that now I construct these functions phii in such a way that this is my so-called phi0 of x this is my so-called phi1 of x and then this is my phi2 of x. So what have I done? I have put these intermediate points at x is equal to 1/2 in my domain; this is the point where I have applied my concentrated load. This phii(s) I am going to define as so-called piecewise linear functions; these are piecewise linear, and continuous. If I look at this points, these points is these functions have a very special property that is phi0 at the point x is equal to 0 is equal to 1, phi1 at the point x is equal to 1/2 is equal to 1, phi2 at the point x is equal to 1 is equal to 1. And if we look at these pictures that these functions at all other points go to 0; that is, if I call this as point x0 this is x1 and x2 then phi0 at 1 at the point x0and 0 at the point x1 and x2, phi1 is 1 at the point x10 at the point x0 and x2 and phi2 is 1 at the point x20 at the point x0 and x1. This is my construction of these functions which are going to be used represent the series solutions. Let us see if this is going to help. If I have to use these functions, what are these functions if I go to the representation? If I do now, I will still call it u of 2 (x) you see the reason why this is equal to u0 phi0 of x plus u1 phi1 of x plus u2 phi2 of x. What do we want our functions to do, because we are still in the regime in the Rayleigh--Ritz method? We want this u of 2 at the point x is equal to 0 is equal to 0 that is my solution, my series solution has to satisfy the specified geometric boundary conditions. If this is equal to 0 at the point x is equal to 0 what do I know? u1 is equal to 0, phi1 is equal to 0, phi2 is equal to 0, phi0 is equal to 1. So this becomes equal to u0. I can knock out the u0 from here. I am left with two terms that is why I wrote the two terms solutions. I have to obtain the coefficients u1 and u2. How do I obtain the coefficients? Again I go back to definition of the total potential energy of the structure and we again look for the values of this u1 and u2 which minimize the total potential energy that is I am looking for the first variation of the pi which here is the function of u(2) is equal to 0. So what do we get? We will get integral, let us take the second problem only. 0 to 1 u2 prime phi1 prime dx is equal to F phi1 evaluated at the point x is equal to 1/2 plus P phi1 evaluated at the point x is equal to 1. Similarly, I will get for the second one, u2 prime phi2 prime dx this is equal to F phi2 evaluated at the point x is equal to 1/2 plus P phi2 evaluated at the point X is equal to 1. I have simply replaced instead of the polynomials a generic definition of these files. These are the two equations that we will get; now I will substitute u2 as u1 phi1 plus u2 phi2. If I substitute these things then what are the equations I am going to get? I am going to get u1 integral 0 to 1 phi1 prime squared dx plus u2 integral 0 to 1 phi1 prime phi2 prime dx this is equal to F phi1 evaluated at the point half plus P phi1 evaluated at the point 1. This is the first equation. Similarly, I will get for the second equations integral 0 to 1 u1 into phi1 prime, phi2 prime dx plus u2 integral 0 to 1 phi2 prime squared dx this is equal to F phi2 evaluated at x equal to 1/2 plus P phi2 evaluated at x is equal to 1. We see that these are essentially two algebraic equations in terms of simultaneous equation, in terms of the unknown coefficients u1 and u2. Let us now define these functions phi1, if you look at these functions phi1 what does it do? It is linear in this part from 0 to 1/2 linear in this part from 1/2 to 1 and it is continuous at the point x is equal to 1/2. It is very easy to define this function phi1 x, this is equal to if it is linear here in this region and it has the value 1 at the point x is equal to 1/2. Then if I make it 2x, then you see that it satisfies these conditions at x is equal to 1/2 at 1 x is equal to 0 it is 0. In this region it is equal to 1 here 0 at the point x is equal to 1. This one will be this region if I want to define it, it will be 2 into 1 minus x look at it. At x is equal to 1/2 this number is equal to 1/2; 2 in to 1/2 is 1 at x is equal to 1 which is 0. So these is the linear, this is also linear so the phi1 is piecewise linear and it is given by this one. Similarly if I want to define phi2 of x, phi2 of x if you remember it is like this. What is the value of phi2 of x in the region 0 to 1/2? If you see this line these is 0 in this region, it is this. In this region it is linear and it takes the value of 1 at the point x is equal to1 and 0 at the point x is equal to 1/2. If I make it as this, 2x minus-1 in the region 1/2 less then equal to x, less then equal to 1, you see that at the point x is equal to 1/2 this expression is 0 at the point x is equal to 1 this expression is 1.This is our representation of phi to x. These two representations we have to plug back in our equations that we have written. The two simultaneous equations - evaluate these integrals and then solve the corresponding problem. I will do the first equation. If you look at these you have the first equation u1 into if I write it integral of 0 to 1/2 phi1 prime which is dx plus integral 1/2 to 1 dx plus u2 into integral 0 to 1/2 phi2 prime dx plus integral 1/2 to 1 integral phi1 prime phi2 prime dx. This will be equal to F is equal to 20 for us into phi1 evaluated at the point x is equal to 1/2. So phi1 at the point x equal to 1/2 is equal to 1 plus P is equal to 10, phi1 evaluated at the point x is equal to 1which is 0. Let us look at these terms that we have obtained. So for us the phi1 is equal to what? In the region 0 to 1/2 phi1 is equal to 2. I will get phi1 prime is equal to 2. So 2 squared is 4, integral from 0 to 1/2 so it will be 4x. It is going to be 4 into 1/2 which is 2; plus in the region ½ to 1. What is phi1? We obtained that phi1 was equal to 2 into 1 minus x. So phi1 prime is equal to minus 2, minus 2 square is 4, 4 into x is equal to 4 into ½ which is again 2, plus I will have the u2 part 0 to 1/2 phi1 prime is equal to 2 what is phi2 prime?. Since phi2 is equal to 0 in the region 0 to ½ phi2 prime is also equal to 0 so this is equal to 0, plus in this region ½ to 1. In this region both phi1 prime and phi2 prime are non zero. What is the phi1 prime in the region ½ to 1? It is -2. phi2 prime what is it equal to? It is plus 2. It will be equal to what we have minus 2 into 2 into 1/2 this is equal to 20. I will get this equation will become highlight this equation. It will become 4u1 minus 2u2 is equal to 20. I will call this equation (a). Similarly if I went ahead and did the exercise for the second equation I will end up getting for the second equation minus 2 u1 plus 2u2 is equal to F into phi2 evaluated at the point x is equal to 1/2. But phi2 at the point x is equal to 1/2 is 0 plus P in to phi2 evaluated at the point x is equal to 1. It is phi2 at the point x is equal to 1 is 1. This is my equation (b). If I write it in matrix form what will I get? I will get 4, -2, -2, 2. This into u1 u2 is equal to 20 and 10. So I have to solve this matrix problem. It is very easy to invert this matrix. I will get after inversion: u1 u2 is equal to... this is equal to 1/2 into 20 is 10 plus, 1/2 into 10 is 5. So this will be equal to 15 here. 1/2 into 20 is 10 plus 1 into 10 is 10. I have been able to solve this matrix problem to obtain this two terms solution coefficients u1 and u2. My solution u2 of x is equal to 15 phi1 of x plus 20 phi2 of x. What was our exact solution if you remember? u exact of x was equal to 30x in the region 0 less than equal to x less than equal to 1/2 and 10 plus 10x in the region 1/2 less than equal to x is less than equal to 1. If I go to this point x is equal to 1/2; at the point x is equal to 1/2 u exact is equal to 15, at the point x is equal to 1/2 what is u2 equal to? At the point x is equal to 1/2 you see that phi1 is equal to 1 phi2 is equal to 0. So u2 is equal to 15 which is equal to u exact at the point ½. Similarly if I go to the point x is equal to 1, at the point x is equal to 1u2 at the point 1 is equal to 15 phi1, phi1 at the point x is equal to 1 is 0, phi2 at the point x is equal to 1 is 1. This will be equal to 20. And from this expression if you see it becomes u exact at the point x is equal to 1. At both this points x is equal to 1/2 and x is equal to 1 these 2 solutions are the same in between they are both piecewise linears. So what do we know that to define a linear we need two points; at both points these two solutions are the same. So my u2 of x is identically equal to ux x. That is using these functions phii that we have redefined this piecewise polynomial functions piecewise linear functions I have been able to capture the exact solution to the boundary value problem that we have been interested. So is this the solution? The answer is yes. This is one way we can obtain the solution to problems which have concentrated loads or material interfaces. Let us see what is so nice about this functions phii x that allows us to get exact solution to the second problem that we have posed. See there are certain properties that these functions have to satisfy; this is something that we have to keep in our head. Whenever we want to construct such piecewise polynomial functions we have to ensure that certain basic properties are satisfied by these functions. One is that well they are piecewise polynomial, and they are continuous; that is, in the whole domain these functions are continuous. Secondly, in this problem what have we done? We have taken these functions such that values continuous derivative is not continuous at the interface of the two pieces that we have taken. These kinds of functions, which do not satisfy continuity of derivatives at certain points, are called C zero functions. Fourthly, what these functions have to do is that they have to form something called a complete set. What does this completeness means? So if I look at what we have done earlier x0 is equal to 0; I will just take only one part of the domain. In this I had drawn this function phi0 and this function phi1. What completeness means is that the linear combination of these functions phi0 and phi1, should be able to represent any polynomial that we take such that this polynomial is linear. That is I should be able to find these unique constants a0 and a1. Such that I should be able to represent the polynomial b1 plus b2 x the linear polynomial in this region exactly using in the linear combination; this is called completeness. If we look at this particular example and let us say we knock off b2 that is a look at the constant. Then I should be able to represents the constant also exactly in terms of the linear combination of this phi0 and phi1. How will I do it? Very simple example is that if we look at phi0 and phi1; sum them up; sum of these two functions is equal to 1 at all points in this region x is equal to 0 to x is equal to 1/2. Very simple situation we have is, a0 phi0 a1 phi1 is equal to b1. Trivially we get that is a0 is equal to b1, a1 is equal to b1. I am able to represent even a constant in this region. Another very important property that these functions have to have which we have outlined earlier is linear independence. What does linear independence mean? That if I take any combination of these functions, so I will do a0 phi0 plus a1 phi1 a2 phi2 and I set it equal to 0; everywhere in the interval 0 to 1 because these functions are defined in the 0 to 1 so that is where I am going to concentrate. In the interval this I want to be equal to 0 for all x lying between the points 0 to 1. What does the linear independence requirement tell us? If this is so then all these coefficients a0, a1, a2 should trivially come out to be equal to 0. Is it so for our functions that we have taken? Since this has to value vanish at all points in the interval 0 to 1, let us take the three specific points x0 x1 and x2 that we have taken. They have to vanish there also; x0 is equal to 0, x1 is equal to 1/2, x2 is equal to 1. What happens at the point x0? If I look at this expressions phi1 is equal to 0, phi2 is equal to 0. At the point x0. At x equal to x0 this expression becomes a0 this is equal to 0. Similarly, at the point x is equal to x1 what happens? My phi0 is equal to 0, phi2 is equal to 0, phi1 is equal to 1. This expression a1 becomes a1 is equal to 0, at x is equal to x2 again by the same token phi2 is equal to 1, phi1 and phi0 are 0; so a2 becomes 0. This is exactly what we needed for linear independence. This functions phiis that we have constructed are indeed linearly independent piecewise linear, complete, continuous and we say that belongs to Czero. Now the question is what kind of a job do these functions do in approximating the exact solution of the first model problem that we have taken. Because, yes, we did a good job what about the first one? Let us go back to the model problem; here we have P is equal to 10. This fx is equal to x. For this module problem, how do we go about finding the solution? The same process that we had done for the problem with this concentrated loads; the only difference we will have is that in this case I will get integral of 0 to 1 u prime delta u prime dx is equal to integral 0 to 1 fx delta u dx plus P delta u evaluated at x is equal to 1. If we remember that our delta u is again equal to, if I take a two term solution, delta u1 phi1 plus delta u2 phi2. Plug everything in to our expression and taking first our delta u is equal to phi1 then taking delta u is equal to phi2. I get the two equations for the two term solutions. It should be obvious to you that the right hand side of the equation remains same as what we had obtained for the case with the concentrated loads; why? Because the left hand side of the equation, remains the same as far as the concentrated load; simply because the left hand side remains unchanged. For both the problems, only the right hand side which corresponds to the load effect of the load externally applied forces, that are going to change. What we can do is we can again pose the problem like this. I am simply going through the steps. It is integral 0 to 1 x phi1 dx plus 0 because phi1 at the point x is equal to 1 is equal to 0 and this one will be the second part will be x phi2 dx plus P . Why? Because phi2 at the point x is equal to 1 is equal to 1. The load vector that we obtain is different from the problem with the concentrated loads we have taken but the so-called the stiffness matrix that we have defined much earlier remains unchanged. If I go and evaluate the load vector, I will get again; solve for u1 and u2, u1 and u2 becomes equal to 251 by 48 and 31 by 3. If I compare with the exact solution of this problem you will see a very curious observation that u exact at the point x is equal to 1/2 is equal to 251/48 which is equal to u1 which is equal to two term solution evaluated at the point x is equal to 1/2. Similarly u exact at the point x is equal to 1 will be equal to 31 by 3 which is equal to u2 which is equal to two terms solution evaluated at x is equal to 1. We see that the exact solution matches the two term solution that we have taken using these special phi(s) that we have defined at the points x is equal to 1/2 and x is equal to 1 which we have called as points x1 and x2. If I plot the functions, here is x, here is u; you will see that u exact will do something like this and this is 1/2 this is 1 your two term solution will do something like this. What if I want to improve the accuracy of the solutions? If we see this solution certainly in between here it is not great that is the gap between the two term solution and exact solution is quite high. What if I want to reduce the gap? So for that what do I do? I will introduce more points in the domain. Earlier I had these three points x0, x1, x2. I am going to add another two points such that the divide I will call it x2 this point is x3 this point is x4. What they do? They divide each of these earlier pieces into 1/2. I am going to define these functions phiis over these pieces. This will be my phi0, this will be phi1, this will be phi2, this will phi3, this will be phi4, phi0, phi1, phi2, phi3 and phi4. If I go ahead and again represent the so-called four term solution in terms of u1 phi1 plus u2 phi2 plus u3 phi3 plus u4 phi4 and solve the problem what I will see is the following; I will get the new solution that will look like this. The four term solution is closer to the exact solution as compared to the two terms solution. In principle what I can do is I can keep on refining this sub division that is to keep on adding; these points will divide the previous sub interval into 1/2 and I can keep getting hopefully closer and closer to the exact solution. This is something that we call Convergence. In principle at least it looks like that this kind of a series solutions does converge to the exact solution. If I take more and more terms in this series that is I take this finer and finer distribution of points with respect to which I am defining these functions. Now these points we have given a name; they are called Nodes and the region between two consecutive nodes, the interval the sub interval, is called an Element. These nodes and these elements together connect the partition of the initial domain of interest into smaller pieces. What these nodes do? They define extremity of the elements. For example, if I am looking at element 2 here; this element 2 will have extremity at points x1 and x3. Element 2 will have these two points as extremity and so on. What we have done is we have laid the foundation of a method by which we can define better and better approximation to our boundary-value problem of interest using certain principles of partition of the domain into smaller and smaller elements defining these extremities of the elements which are called the nodes. This whole concept of breaking a domain into smaller pieces or sub domains is what we call the Finite Element Method. And we see these phi(s) that we have defined are only non zero in small neighborhood of the point with respect to which they are defined. That is the two neighboring elements which are joined at this point xi, this function phii is non zero; elsewhere it is 0 . These are also called basis functions which have local support. That is they are non zero only in a small part of the domain, 0 everywhere else. In terms of these basis functions we can construct a series solution to the problem of interest and hopefully the series solution if I take finer and finer partitioning of the domain will converge to the exact solution. The whole process of partitioning of the domain that is assigning these nodes is called Meshing. The process by which we are converging that is, by adding more and more nodes this kind of process is called h version of the finite element method. What we are going to concentrate is in the so-called h version of the finite element method and will see whether we can guarantee conversion of the approximate solution to exact solution. We will see that this kind of the basis is able to capture the solutions to problems with continuously distributed loads or concentrated loads or material interfaces, provided we put the so-called nodes at the interfaces. If I have a domain like this, let me take a very pathological example. Here I have some material; here I have a concentrated load. Then the partitioning that I do that is these nodes that I put will be here, it can be here, and it has to be here. One node has to be at the interface of the two materials. This is one material, this is another material, another node I can put here and another node has to be at the point where the concentrated load is applied and so on. This where I can define the partitioning of the domain and based on the portioning of the domain since this phii(s) defined with respect to these nodes I can define what are these functions phii and so on. Go ahead and construct my approximate solution to the problem of interest. What we have seen is that the Rayleigh--Ritz method does a good job with polynomial approximate functions for domains with one type of material with loading; which is nice and continuous, then the Rayleigh--Ritz method does a good job. If I have dissimilar material or concentrated loads the Rayleigh--Ritz method using polynomial approximation fails. That was the motivation for introduction of these piecewise polynomial functions which constitute the basis of the finite element method that we are going to develop in the next lecture using the simple one dimensional module problem. We will go ahead and give the detailed formulation of how to construct this matrix problem in order to determine the unknown coefficients. We will solve that problem and see how the solution looks like, what is the accuracy of the solution, and when the solution will do certain nice things and when will the solution fail. So that will be the starting point of the finite element method.