Lecture - 2 advanced finite elements analysis

We saw various aspects of non-linear finite element analysis, what are the types of problems that can be solved using non-linear FEM. In that process we recognized certain very important things. One is that when we do a problem, of course, we have to decide whether we are going to do stress analysis problem or whether we are going to do heat transfer problem or fluid flow and so on; that is very simple aspect. But, once we decide that we are going to do a stress analysis problem, the first thing we have to do is to recognize what is the type of analysis? Whether you are going to do a static analysis or whether a dynamic analysis? That depends, of course, on the problem. That is the first thing that we have to decide. Now, it is not an easy decision. Many times, certain dynamic analysis you may have to simplify using static analysis and many static analysis problems can be very easily done using dynamic analysis. So, it is not that it is a very straightforward decision many times. Of course, you can do it. You can say, for certain class of problems, you can definitely say that they are dynamic and they are static. But, certain class of problems, though it may appear to be a static problem, for example many of the metal forming problems are today done or carried out using dynamic analysis; dynamic explicit analysis, basically because you also have to consider certain problems that may arise during static analysis. So, you have to know theory of non-linear finite element quite thoroughly in order to make that judgment, how you are going to do that. So, once you decide what is the type of analysis you have to do or you are interested in, you also have to look at whether the problem is going to be a small deformation problem or a large or finite deformation problem. That is the next thing that you will decide. Then, afterwards you will decide what the material property is that you will use? One of them you have to choose. Of course, along with it you will choose the type of element that you will be employing in order to solve that problem; whether you are going to do a shell or whether you are going to do a membrane or whether you are going to do or whether you are going to use a solid element and so on. So, these are some of the other things that you will also consider apart from the type of analysis and so on. Of course when there is contact you will do also contact analysis and that will complete to certain extent, to certain extent, there are other issues that you have to look into before you can actually do a finite element analysis, but this to a great extent completes the picture. Now, let us go into the theory of finite element or non-linear finite element analysis. The theory of non-linear finite element analysis also moves in a very, very similar fashion. What is that movement or what is that we are or how is that we solve a non-linear finite element problem? Basically if you, if you look at it, in very loose terms I can say that non-linear finite element analysis consists of a number of algorithms; a number of algorithms, algorithms to solve a non-linear equation. This can be, the most popular among them being Newton-Raphson scheme. You have other schemes as well, BFGS and so on, but 90% of the problems are solved using Newton-Raphson Scheme. So, this is an algorithm for solving the non-linear equation. So, this can be employed, whatever be the problem. Now you have for solving, N-R or Newton-Raphson Scheme for solving the non-linear equation. Now, you have algorithms or I would say formulations, which would take into account incompressibility, formulations for incompressibility. Again, whenever you face incompressibility constraints, whatever be the type of material property you can use this kind of formulation for incompressibility. It can be selective reduced integration or it can be a mixed formulation or whatever it is, you can employ that technique right across. Similarly we will have certain considerations including kinematics for taking into account large deformation. The way we are going to or in other words we may look what is called an updated Lagrangian; we are going to look at total Lagrangian and so on, which would, which becomes important when we are looking at finite deformation. Again this is common to whatever be the type of the material behavior or whatever be the type of element that you are going to do; they become or they are common to it. So, you pick up one algorithm from here. Of course, there has to be certain some small compatibility issues, to pick up one from here and one from here and you can say whether total Lagrangian or updated Lagrangian and so on that also can be used to solve. This can be looked at as a very general algorithm procedure to solve problems independent of type of material behavior and of course you have contact and contact is looked again from a very independent angle and contact algorithms can be used for solving certain problems. Of course, it is not that, just it is a plug and play; not that you just pluck one from here, one from here, one from here, from here, whatever be the consideration combine it and get an algorithm; it is not like that. There are considerations. For example, in contact not all algorithms can be used for large deformations. There are algorithms which can be used only for a class of problems. When sliding is large, then you cannot use it and so on, but still once that class is clearly put, what is the class of problems that can, that this algorithm can solve, then that is what you should look for, in order to pick that up to get to your solution. With that in mind, let us now quickly go into the solution procedures for the non-linear problems. In fact, if you had noticed, I had not really, rigorously defined non-linear problems. I know I have, I have appealed to your feel to say what non-linearity is. I have not yet defined mathematically what non-linear problem is. All of you know looking at it, yes, if it is non-linear elastic it has to be a non-linear problem; if it is going to be elastoplastic, you know that it is going to be a non-linear problem. But, I have not rigorously defined what non-linearity is. Where does that term come into, from where that term comes into picture, especially in finite element analysis? The definition is not very difficult. All of you know the very fundamental definition or fundamental equation for finite element analysis. What is the equation for general finite element analysis, most fundamental equation? No; ultimately what you solve? Yeah; K d or K u is equal to F. This is the fundamental equation we solve and when we solve that, all of us know from our earlier fundamental finite element that determination of K is one of the major issues in the element formulation and that we see that this K is independent of the, what is the assumption that we make ultimately when we get this formulation? K is independent of the d that it goes in. So, K is not a function of d. You solve the problem only once; you apply the load, you solve this K, use the K to solve this equation for d by whatever be the procedure. Right from conjugate gradient or use Gauss elimination, whatever it is, you solve this problem. The difference in this case is that this K is not a constant. In a non-linear problem, this K is no more a constant. So, this K depends upon the deformation; this K depends upon the deformation, so, I cannot, in other words, I cannot determine. Once I say K say, a function of d, say displacement and then put d is equal to F without even doing finite element analysis, without even knowing finite element, you are going to ask me how are you going to determine K? It is not. So, you have an equation which becomes non-linear and you cannot solve it in the same fashion as you did here. Yeah, the question is will K depend only on d and how it depends on d? This depends upon the problem. In fact the major issue here is, like in the previous case calculation for K, calculation of K for different elements is an issue in the linear finite element. In this case, calculation of K for different material behavior is an issue. The whole issue here is the calculation of this K. But, we are not going to write it in the same form; we are going to change this form, nevertheless. But, the issue here again is that formulation of stiffness. We are going to call that in a slightly different fashion; we are going to give a new name to it. We are going to call this as tangent stiffness. The reasons for it will be clear as we go along, but this K is a function of deformation and how it is, what function is it and how does it depend upon it and how will it vary, all these things we will do. Is it clear? So, right now I am not going to define them. The question that you may ask is, is it only function of deformation alone? It may also be a function of damage. Suppose you are looking at damage mechanism, say for example you are analyzing concrete or composites and then there is damage in the material. Then, K definitely would vary. You can, of course, say that damage is due to deformation and so, in one sense K is a function of deformation. That is why I said that K can be function of d or another parameter which is again a function of d and so on. Is that clear? We are going to, anyway we are going to see quickly, the whole idea is to quickly get in and see how we are going to solve these problems and get to know the whole issues of non-linear analysis by using a one dimensional case. So the first and foremost lesson, yes? In the linear analysis we have K independent of d. Yes, K may depend upon time; K may depend upon time, but ultimately time; no, no; it would, because if you say that K depends on time because of deformations, again it becomes the non-linear problem. As long as K is not one factor, as long as K is not one factor and K depends upon deformation by some means, we are moving into non-linear problems. We will, we will see that as we go along how K would depend upon and so on; just may be it will take a couple of classes or more. Let us, let us first look at the procedure. The first difference, the major difference between a linear problem and a non-linear problem is that we are going to resort to an incremental iterative approach, iterative approach for solving this non-linear problems; incremental iterative approach. Now incremental iterative approach has come about because of the type of algorithm we are using; of course, Newton-Raphson scheme and so on. But, let us first understand what we mean by this term. Then we will look at the equations which govern incremental iterative approach. What it simply says is that, since K is a function of d, I cannot apply the force, all the force straight away. What it simply says is that I have to apply the force in an incremental fashion. That means that if I am given 1000 Newtons, I do not apply all the 1000 Newton at one go. I will apply 10 Newton; I will apply 10 Newton or 100 Newtons initially, make the body or component to come to equilibrium with that load, then increase that load may be to 200 or 300, then make the body again to come to equilibrium and so on. This is what we mean by an incremental approach, where I keep incrementing the load. What is iteration? As and when I increment the load, I iterate or in other words get better and better displacement quantities or displacement measures, so that the internal forces are equal to the external forces or in other words the body is under equilibrium. So, I keep iterating, so that the displacements are good enough for me to develop stresses which can hold the external forces. This is what we mean by an incremental iterative approach. In other words, we have two big loops. What are these two big loops? One is that I will keep increasing the force or external forces; force incrementation, till we reach the final figure and within this force incrementation, I have an iterative loop; an iterative loop, so that for the force given, I get a good value or good in the sense that correct value of the displacements so that equilibrium is maintained. That is what is meant by an incremental iterative approach. Is that clear? Yeah, so now, see whatever it is, ultimately in reality, the question is in reality we do not apply loads like this. May be if I want to apply 100, then I will apply 100. Whatever it is please note this carefully and please note what the difference between a static problem and a dynamic problem is. Right now, we are talking about static problem. When I say static problem, then time is not a factor. We do not consider inertia forces. So, I apply this 100 Newton or 1000 Newton, whatever it is and I do not invoke that inertia force; ma does not exist. It is so small that you neglect it and keep applying the load and your final interest is at the point when, whatever load you had applied is under equilibrium with the stress. Is that clear? So, it will not be any different as long as the process is static. But, when will that have a difference? Suppose it is dynamic. I go and hit, then it is a dynamic process. Then, I cannot do it obviously as a static problem. In that case, I will automatically make a mistake if I do it as a static problem. Is that clear? With that in mind, let us now develop a simple equation using Newton-Raphson procedure and first do a one-dimensional case to understand this problem. What is that we have? We have two things. One is an external force and another is an internal force. These are the two things that we have. Let us say that the external force is given by f and internal force is given by P and we know from our earlier classes or earlier study that P, which is the internal force can be written as, do you remember? No; from a finite element perspective you can write this down as, yeah; no, this is B transpose sigma d omega or d v or a d a, whatever it is d omega, B transpose sigma d omega. In fact, if you remember, it is this sigma which when replaced by means of stress strain relationship and then strain displacement relationship, gave rise to our K, because sigma can be written in terms of d epsilon; epsilon can be written in terms B u, so B transpose d B u is what we got and that is what we obtained k. So, B transpose sigma d omega is the original or the internal forces that would be developed due to external load f. At any point of time, obviously, f minus P should be equal to zero. Note that I am not correct in writing just P. Actually it has to be vector quantity. So, in order to distinguish that from a scalar quantity, I will put a small squiggle underneath to just say that these quantities are vector quantities. Obviously I cannot write it in bold letters, so I am just putting that and you know what the size of these vectors are or if it is a second order tensor you know what these things are and so on. So, f minus P is equal to zero. This is the equation. As I already told you, it is not possible for me to get the internal force straight away, because of the fact that sigma now is not just d epsilon; that d itself is a function of epsilon and so on. So, I cannot write it like K, like what I had done at that point of time. I have to or I will definitely end up with an error or residue. This is what we call as residue. The whole idea here is to make this residue go to zero; actually not zero, but to a very small quantity. The residue, what is residue? Residue is actually the error that is the difference between what has been applied and what has been realized. So, that residue has to be a very small quantity. Now before we go further, a small introduction to terminology, because especially when you use finite element software packages, the software packages do not or may not, I should not say do not but, may not use just increment of load but they would talk in terms of time. Most of the softwares talk in terms of time, application of or incremental load or application of loads in terms of time. Now, note that this would, this is a point of confusion, basically because as our friend here asked now time, is it a factor, it is a static problem. Many users of finite element get confused when they see time as one of the factors, because in static problem where do I get time, so how do I handle? Actually the time that you have in a static problem, it is just a pseudo time. There is no, you know the interpretation is not the same as that you see in a dynamic problem. It is just a pseudo time. In other words, time is a carrier of force. What do I mean? I can say that let me apply all the load in one second, say 1000 Newtons in one second. Instead of saying that at every step you apply 100 Newton, which becomes cumbersome when you program it, what the softwares do is to say that increment it in terms of 0.1 second. So, when you say that you increment by 0.1 second, 10% of the load is applied in the first step. So, 100; instead of saying increment it by 100, second step increment it, increment the load by another 100 and so on, you say that let there be a time step with delta t is equal to 0.1. It is much easier to handle. So whenever we talk about time in this static analysis, understand that they are used as a carrier for the force. That is all, nothing more to it. Is it clear? It is usual to call that as say t and that t is incremented and so it is usual to put a subscript here as n. So, t would start at t1, like what you do in a dynamic analysis our t can start at t0; then, t0 can increase to t1, where t1 is equal to t0 plus delta t and so on. So, you can write down that equation as tn+1 is equal to tn plus delta t. That is the time factor. So, n plus 1 is an incremental step. Is that clear? It is an incremental step. Say for example, if t is equal to 1 second, as I told you it does not matter; you can have t is equal to 2 seconds and say let my time step be 0.1. It is the same as saying that t is equal to 1 second and time step is equal to 05; it is the same. Here for example if I have, if I start with t is equal to zero, which you usually do, then tn is equal to zero. So, t1 becomes zero plus delta t which is say 0.1. So, t1 becomes 0.1. t2, which is an increment over t1, now is 0.2 and so on. Note that this delta t need not be a constant or no software will have this delta t to be a constant. See, actually given a choice, you are not interested to do this. You already realized that this is going to take enormous time, computing time. What is that your are actually doing? You had a nice equation in the linear case, one nice equation like this. You are going to develop similar equation for every increment and you are going to solve this 100, 1000 times; may be 10,000 times you are going to solve the same thing. You can immediately see the enormity of the problem and you can immediately see the type of time you are going to spend before the computer before you can get the result. So, the smaller the delta t, the worse you are, because you have to solve that problem so many more times. So, most softwares will not keep this delta t a constant. What would they do? They would keep galloping this time step; they would keep galloping. That means that they would start with 0.1. If things go well, then they would watch for may be one more step. Then they would raise the delta t to be 0.2. If things go well that means your convergence is good, then they gallop that again may be 1.5 times or 2 times or whatever the software decides. Then they will keep galloping it. So, delta t is not a constant. In the same fashion, they will do the other thing also. Suppose you put 0.1 and in the iterative process, I said already that in the iterative process I have to get convergence. What is convergence? That means your internal loads or internal loads developed by your stresses which is P and external loads they match. That means the error residue is very small. If this residue is not small even after very many number of iterations, usually you put a limit. Usually it so happens that, if it does not, if it does not converge within say, 16 iterations or 20 iterations, you are not going to get convergence that easily. That means your time steps are large. So what the software does usually is to cut down delta t. If it had used 0.1, then it would realize quickly that look at 0.1, I am not getting convergence; I will cut it down and go back to 05. Why this convergence does not happen, we will see that as we go along. So, the delta t is not a constant, it can gallop; when things go well it can increase or it can be made to reduce or cut down. We call this as cut downs in time steps and so that would happen. It is not correct to write just delta t alone, but put one delta tn as well. Is that clear? Now, I will bring that to this equation and write this as psin+1, this will be fn+1 Pn+1 should be equal to zero. This is the equation I have to solve and this equation is solved by, so this is not going to be zero first; I will remove that. This equation is solved by an iterative procedure, very standard technique. This iterative process is a very standard technique for solving any non-linear problem and Newton-Raphson, probably many of you would have studied this in your numerical analysis course; that is exactly what we are going to apply, nothing very different from it. So, what do I do in order to apply Newton-Raphson technique? I start with psin and get a better estimate of psin by what is called as linearization procedure or looking at the first term, linear term of the Taylor series. That is all we are going to do. So, this is called classically as linearization, so that I can write the psi say, k plus 1; look at that k, I am introducing it here. I will just come back to it in a minute is equal to psi k n+1 plus dow psi by dow u. Why u? It is obvious that this and this are functions of u or d however you want to call it; though this, we had called this as d, d is the matrix, u is the displacement and d is the matrix which consists of displacement. We will follow that u. So, obviously this is a function of, this may or may not be a function; for the time being let us say that, that is a constant. It may or may not be a function of u. P is obviously a function of u, because P is a function of sigma, sigma is a function of u and so on; epsilon epsilon is a function of u and so on. So, this is function of u and psi is a function of u; so, dow psi by dow u and so that k remains, k remains into du, this is n plus 1, du k, du k at say n plus 1 is equal to zero or you can in fact, it is better to write that as n because that will be easier to follow later. This is the crucial or crux of the whole of non-linear finite element, nothing else. If you understand this equation quite clearly, then you can understand the whole of non-linear finite element analysis. Whole game is to solve this equation. Of course there are, this will not be actually equal to zero. Why, because there are other higher order terms; we have left them out and so this is a truncation that we have introduced and we do not know the pattern of psi and so on. So, we cannot write that equal to zero. It is approximately equal to zero. In other words, these k's indicate, what? The iterations. We march in two regions; one is first in n that is increment, another is in k, iteration. Is that clear? So, solution of this is very straight forward, it is not very difficult. How do I do that? Look at that carefully. Dow psi k by dow un+1. Look at this. So, what is dow psi? Look at this and look at this. What is this? f, I had already said that it is independent of, yeah minus dow P by dow u. So this term is minus dow P by dow u. In other words, in fact this itself can be used to solve, so, this can be written as dow psi into du which is increment. What is this? This is the increment in u with respect to the, this particular iteration is equal to minus of this quantity here, or sorry, the residue in the previous case. This combined with this equation, I get or this can be looked at as minus dow P by dow u at n plus 1 and call this as KT which is the tangent stiffness, what is called as the tangent stiffness matrix; tangent stiffness matrix and that when substituted back into that equation gives rise to what? This equation, when you substitute it here gives rise to an equation of the form KT into du, let me remove all these things for clarity, you can put it back again; KT into du is equal to psi that is P minus f; P minus f. Of course, you can write this, please note that instead of this, this minus can be removed, here this minus can be removed, KT can be just defined as dow P by dow u. Either way you will get the same result or in other words, you can start with P minus f, in which case you will, actually if I define KT to be just dow P by dow u, if I define, if KT is dow P by dow u, then what would be the equation I will get? KT du is equal to f minus P. Anyway that minus sign and plus sign is quite clear. It is very simple; where you start with will matter or will decide what the type of signs we have. So, KT, the tangent stiffness matrix is the core issue in non-linear finite element analysis. Like you had K which is the core issue, in non-linear finite element analysis here, KT is the core issue. How am I going to solve or how am I going to get KT? Now, once I get KT, my worries are solved, because I can solve this equation and get du K and build up on du K. In other words, I will just remove this and say that once I solve this equation, I can get down to writing u; u at n plus 1 at say K to be, that is at the end of the Kth iteration once I solve this, to be how do I write this? un plus sigma i is equal to 1 to K du K n. Please again note that this n plus 1 and n, there can be confusion, because people may start, put here n, here n minus 1 and so on. I think that should not be confusion, because when we write down completely the procedure, things will become very clear. This is what we would also call that as delta u; delta u k that we are going to calculate. That is that sum of, the sum of du's, sum of du's as we proceed in iterations. It is quite, I hope it is clear. Any question, is it clear? So, once I know this u k n plus 1, and I know un, of course. What is un by the way? un for the nth time, the converged value of course when I start, it is going to be zero. If I am in the 5th increment, un is the converged value at the end of the 4th increment. So, this is the end of the 4th increment. Here suppose, I am, I am in the 3rd iteration of the 5th increment; suppose I am, I am, in the 3rd iteration of the 5th increment, then, what are the factors that will be involved? u, no, it will not be u2, it will be u4. u4 plus, I had already calculated du for this. So, this will be du 1, du 1 at the 4th that is 4th to 5th; I am going from 4th to 5th, so that one or if you want, you can write that as 5. How you do that is in your thing and then plus du 2 plus du 3 that will give you 5th; the previous converged plus what you have calculated in this current increment, this iteration. It is customary to put n. If you want you can also change it to be to n plus 1 and say that this is 5. That will not change any way the result as long as you understand things. Basically why because you may start with 1; so, t is equal to 1 or n is equal to 1 is where you may start, in which case u 2 will become u 1 plus the results of the first increment. So, you will get 1 here. You can put either 5 or 4. That will not cause any confusion as long as you understand that in the current increment as I go from 4th to 5th, what are the changes that I see in du's? Now, once I come to this step, it is very simply. Then, I calculate epsilon based on these results, strains. Then, I use constitutive equation to get stress and then what do I do? I would check or I would calculate p. From here, I would ultimately get sigma actually; sigma, which would now be tested against f. From this sigma, I would go and get P say, let me call this as P5 at the end of the third iteration. Now, I would check how good is this P 3 5 when tested against external force? What is that fn+1? That is f5; f5. How good or what is this difference or what is the residue? Exactly; f5 is going to be a constant. f varies only with respect to n, f5 would be constant; so, what really happens is that I will keep increasing this P, so that I will approach f. Let us give a very graphical procedure for this whole thing in a minute, so that we understand what we mean by this whole procedure. Once I see that this difference here is small, is small, I say that small is again a word which we have to be careful; we will define it. So, once we understand this whole thing, we will put down this procedure very neatly then, you will understand what we mean. Small means it is not zero, but it is a very small quantity; what is small is also important. Suppose I do a problem in Newton meter, the residue may be one value. Suppose I do the same problem in mm, the residue may be another value and so on. It has to be independent of what you are working on. So, you have to be very careful in defining what this small is. We will define that again and so, when once it is small, we say that convergence has been achieved. Then, we go back and go to the next iteration; 6th one, 6th one, increase the load exactly. If not, if this fellow is not small, what do I do? I get back to here, next iteration; that means I recalculate. What do I recalculate? No; look at that carefully. That is what is that called as? Correct; it is the tangent stiffness that you recalculate. Then, you again calculate du solving this equation and then again do it. Now, again look at the problem. Sigma is the stress which you can calculate; yes, because sigma is related to, we calculate, from sigma we calculate P because we had defined, no, no; it is not minus. This will lead to, sorry; this will lead to, this difference or residue calculation. Is that clear? Note this again. In fact, the problem dimension is tremendously increased. Now, every time I do this in a Newton-Raphson scheme, I form tangent stiffness; forming tangent stiffness takes time. Then, I solve this equation. In every iteration, in every iteration, I solve an equation or in other words, every iteration is just equal to that of linear problem. So, it is something like a linear problem. So if it takes, if linear problem takes a second, this you can understand how long it would take? Is that clear? So, that is why non-linear problems take tremendous time, not only incremental approach, but also an iterative approach. So, tangent stiffness every time when you form it, that takes time. Yeah; you can, you can modify it a bit. We will do that graphically in the next class. You can modify this procedure a bit; call this as modified Newton-Raphson. You will understand that, once I, once I put this whole procedure for one dimensional case in a graphical form. So, essentially what is that we did? Just to summarize what we did, we looked at residue. Residue was defined as the difference between the external and the internal or the internal and the external forces, difference between them. Then we linearized this residue, arrived at an equation by truncating to the first term, Taylor series approximation. From this we developed or we came back to our very familiar, simple simultaneous equation, defined tangent stiffness, solved it. For every iteration, you keep a watch on this residue. If the residue happens to be very small, we said that we are happy with it, go to the next step and so on. We keep repeating the steps until we finish the problem or when t becomes t total say, once second what we have defined the Yes; so, what is changing tangent stiffness, what do you understand by tangent stiffness? This is exactly what we are going to see in the next class. We will see it starting with a graphical procedure. We will see what exactly do we mean by this term tangent stiffness?