Lecture - 5 advanced finite elements analysis

In the last class, we were looking at a very simple non-linear problem and we had worked out for the first increment the complete iterations that are required so that the convergence results are with us. Then, we started with the next increment and we had completed one iteration of it. In other words, what we were trying to do is that we have a big loop which we called as the iteration, sorry incremental loop. There is another loop inside which is called as the iteration loop. So, we had finished one here and we had finished this and saw that convergence was attained. We are in the second loop. Without wasting much time, let us look at the final results and make certain very important observations which would go a long way in seeing or in understanding the non-linear analysis. Now, let us come to the iteration two for the same problem. That is this is increment two, iteration two. I have calculated already the tangent stiffness matrix. Then du 2, note that du 2 is this increment in u for this particular iteration. Sigma of du's was the ones which gave us delta u. Then, I calculate the u, epsilon, sigma and calculated the error. T obviously . Then, we go ahead from iteration 2 to iteration 3 using the error to be the load or the right hand side; remember that is what we did. Now, what I do again is to calculate this stiffness parameter, actually stiffness matrix. In this case it happens to be just one number,because of the type of problem we are doing and then we calculate du 2 for the third. Look at that. It is the third iteration for the second increment. Already I know that you would have been exasperated looking at this problem on doing this. So, you can imagine that if there is a simple problem like this with only one element, just one elementand not a very complicated stress strain relationship and if this is going to take the kind of calculations that are involved here, these are the calculation are involved for this simple problem,imagine the amount of time that is required to solve a much more complex problem using the non-linear analysis. We calculate u, calculate epsilon, calculate sigma and calculate the error, go to the iteration; hopefully this is the last iteration. Calculate KT and repeat all these steps so that ultimately, thankfully the error is quite small. There are some very interesting observations you can make. Number 1, look at KT's. In fact, compare KT's with what we have got in this say, iteration even in this increment can compare them with the KT's which you had got in the previous increment. Have a look at that. What is it that you see? What was our first KT? If you remember, what was our first KT? Two point, 2.4; o,sorry! 4; yeah;it was 4.0, if you remember right, yeah. So, 4.0; we started it from 4.0. Look at the way KT's dropped; 4.0, 1.552 and now we have got 0.816. It is not a very good trend and this is what happens many times in finite element analysis, when you look at these kinds of problems. At one time this is going to be so small;it going to be so small that convergence becomes difficult. In fact, if this happens to be a multi dimensional case, K would become condition. So, you cannot apply whatever step you want. In other words, what it physically means for this problem is that the slope of this curve is very nicely sitting here with the value of 4 is now dropping; it is dropping. Once it drops like that, obviously KT which depends upon this ET, of course, you know the sigma versus epsilon; KT, which depends upon the slope of this curve, it is going to become smaller and smaller and note that when I calculate the second step, this KT happens to be in the denominator here. In other words, in actual multidimensional case that is two dimensions or three dimensions, I have to calculate KT inverse. Though you may not explicitly do it, you can say that you, you have to calculate that in order to solve for du and this calculation, number 1, poses lot of problems when K becomes singular; of course, you cannot do it and if it is, even if it is not singular and even if i condition, this kind of thing pose problem so convergence becomes difficult. So, these are the difficulties that we are going to have. But, of course, I am not saying that every problem will end up in a situation where you cannot solve. It would so happen that thankfully when we have lot of elements and the deformations are e, somehow KT would still keep going. In other words, when deformations are not the same at every point, the stiffness matrix that are formed due to these elements would be such that it will be better than what we have seen till now. Here there is only one element. You have to understand that many times one element problems are more difficult than, you know when we have more than one element, because of the inhomogeneous type of deformation. What is inhomogeneous type of deformation? That means the deformation is not the same in every element. This is what we call as inhomogeneous behavior, if the behavior itself is different from one point to the other. They are not the same; two points are not the same. Yes, in other words, the error and hence the number of iteration are required depends upon the shape of the curve. In other words, the question is, is it that we have taken four iterations and if I have the same step and then, in other words, same step in sense that same size of loading in the next step, I go to step 3 say, I put another 2.4, will I take more number of iterations? Of course, I will take more number of iterations, because as I move towards or these positions where the ET happens to be smaller and smaller, my time or my time step or my load step, both of them being the same, have to reduce. So, for example if you are going to do linear problem, yes, I can look at a linear problem like a non-linear problem as well. In other words, I can say that I would go ahead and do convergence test in a linear problem. Obviously you will converge in the next step, because when you move from say this step here to here, you would require only two iterations. Why the second iteration, because I have to calculate the error and compare it with the first one. So, I would get a beautiful answer in the first iteration itself. Second iteration will have error equal to zero and I will say that I have converged. Many packages till very recently look at even a linear problem like a non-linear problem. Is that clear? Yes, what is the question? Yes, epsilon is 3.998. There is very small increase in the, sorry, 0.3998; not 3.998, 0.3998; yes. Now, having understood how non-linear problems work, let us now go back to elastoplastic case. There are very interesting questions and see how these questions can be answered from an algorithm perspective when we do finite element analysis. Any question on what we did? Now, if you notice in this case what we did there was one step to calculate stress. The problem is not going to be very simple when I do elastoplastic case. Just for simplicity we assume a bilinear case, where we say that we have only two slopes that operate and that a well defined yield point exists. Now, please note that this is a simple assumption. We can definitely modify it. We are going to do that later, but in order to understand the issues I am going to just look at this as a bilinear problem. I had already made lot of comments about this slope. It has to be, I mean this is exaggerated; it has to be very close to the y-axis and so on. Now, what are the issues? 5 The issue is that when I have this incremental loop and iterative loop, at one point I have to calculate the stress and this stress calculation which is called stress update in the literature or stress integration, certain difficulties. What are the difficulties? 230 Assume that I am at this point, at this point here, when I am in n and I have to move ahead in stress. Assume that my deformations during that step is such that the epsilon that I get make me move from the elastic region into an elastoplastic region. Means that let us say that that is what I get as delta epsilon and that when I move my stress, I have to move from this point to the other point. So, what is the difficulty? The difficulty is that I have one slope or one type of behavior till I reach this yield point and another type of behavior as I move out of this yield point. So, this problem though that it looks very straight forward in a one dimensional case, because of the fact I have drawn very clearly and shown you this,is very and involves lot of issues when you look at it in a multiaxial case. In fact, this had prompted quite a few researchers to work on stress update algorithms and a series of papers have been published from the, from about may be late 70's to the late 80's. In fact, though the first paper on the stress update appeared way back in 1964, there have been a series of papers which had looked at the stress update algorithm very closely. Now, what are the other issues? There are other issues as well. I have to calculate the stress such that it lies on the stress strain curve. In a multiaxial case, we would see that this means the stress would lie at the yield surface. I cannot get a stress at this point for that strain, so I have to see that I lie on the curve. This is what we would call as consistency condition and consistency conditions also have to be satisfied you will understand the problem exactly now as we go along. So, I am going to take off from this point. 516 In other words, I have already calculated delta epsilon till that previous step and we will look at one dimensional problem and deviate from here and see how we calculate the stress. Clear? 537 The first step is that we invest all the strains, what all we have that is delta epsilon, which we have calculated here, as if it is elastic. This step is called as the elastic predictor step. Is it clear, elastic predictor step. The first thing is that we will go to elastic predictor. Of course, before that we had the other steps. If you want, we can write that down also so that you are clear about it. So, let me write down those steps: a, I calculated psiand then b, I calculated delta u and I calculated delta epsilon. That is the step which I did. Then c, I do an elastic prediction. I made some comments on delta epsilon, if you remember, if I remember right I said that I would all the time start again and again from the previous converged value of u and then calculate delta epsilon not based on du and up, but I would base my calculations on delta u and then put that as delta epsilon. This is what we said. Keep that in mind, we will make more comments on it as we go along. So, the first step is that I do an elastic prediction. What is elastic prediction? That means that I get an estimate. We will call this by different names. People say that it is a trial stress or an elastic estimate; by saying that this delta sigma, we will say sigma E, what sigma E is in a minute,delta sigma E that is elastic estimate of delta sigma to be E into delta epsilon. This is one dimensional case, so, justI am putting that as E. In fact, there is no harm even removing this, because it is for illustration purposes; we are going only look at the stress strain curve. What I have done is to say thatI am going to invest all these things in the elastic case. What do I do? Now, I calculate sigma E to be sigma at the previous converged step plus delta sigma E. Is that clear? N,. We are there, because I have extended the slope directly. I mean, though it is very, very special case there are other cases that exists. So, we have to exhaust all the cases as we develop this algorithm. So, I am going to, I am not going to write down all the steps here, but it will be clear as we go along. As we go along it will be clear, but I would like you to follow these things and then may be when we do a problem later, we can put down all these steps in a more form. Now, there are number of cases which may happen at this point. What are the cases? Though I have just drawn a figure like that, this may not be the case. This may be one of the cases that may exist; there are many other cases that may exist. There are two things that are possible. One is that, the previous step that is sigman can be in the elastic region like what we have said or it can itself be in the plastic region, which means that when I have a curve like this,I may start my whole sigman sitting at that point. This may be my sigman. Note that, I am going to call this as initial yield, because that yield point is going to change. Is it clear? So, there are two possibilities. I may be in the elastic region or I may be in the elastoplastic region; that is my first thing. My second possibility or second issue which results in two possibilities are that this elastically predicted case, elastically predicted case may be such that I may still be in the elastic region, may be still in the elastic region or I would have gone to a plastic region, in this case. In this case, I may march ahead in the elastoplastic region or I would have come back and I would have been in an unloading curve. So, there are two cases initially and each one has two other cases. So, we have to look at all the cases in order to develop the algorithm. Is it clear, any question? No, correct; doesthis case arise only when we are looking at situation where we are going from an elastic to plastic? That is exactly what I said; going from elastic to plastic is only a special case. General case, the general case is such that I can go from elastic to elastic, I can go from elastic to elastoplastic, I can go from elastoplastic to elastoplastic or I can go from elastoplastic to elastic. How do I do that? Unloading; so, I have to consider all the four cases. Yes; we have not yet come to KT definitions. That is one of the things that we are going to find. You will see that there are lot of small tricks that we are going to do. In other words, the question is what will happen if I sit exactly there? In most of the problems, it is almost impossible to sit exactly at that point. To give you a small into program development, what we usually do is we take that point and give a small band. S,if I say that the yield strength of the material is 300.2 newton per mm squareor 30.02 kg f per mm squared, sorry mpa, in terms of mpa, it is say 300.2 or 30.02 kg per mm square, then it is impossible for me to exactly hit that point, because we are doing a numerical process and it would so happen thatyou may be either slightly below this mark or above that mark, so we will define a small band there for ..... You can choose this as a,so it is not, that is exactly what we do. Is it possible to have a smooth variation? That is exactly what we do. We say that yielding takes place usually when it is in a small band. We will discuss this especially when we come to a multiaxial case. You know, this is going to be issue there, small issue, implementation issue; we will discuss that at that point. So, let us understand first the broader aspects here. So, these are the four cases. In order, in order to handle these four cases, first what we do is to have an elastic prediction like this. That will settle us, that will settle our position; that will settle where we stand in one of these four cases. Let us see what we have to do when we get into one of these four cases. Of course, the problem is very simple if we go from elastic to elastic. It is not a problem at all. But, what happens when we move from elastic to elastoplastic and plastic to plastic? We will right now see that the elastic to elastoplastic is a more general case ofelastoplastic to elastoplastic. The other case elastoplastic to elastic, which is unloading, is also very similar to the first case; it is not a problem. In other words, what do we do to find out these cases? Simple; check sigman and see whether sigman is less than, we will include his case, sowe will say that less than equal to sigmay. So, we will say that it yields when you just cross sigmayand so if this is the case, then we say that it is elastic. In order to do a stress update, we need to check this as well. Now, let us say that let us consider one case, most difficult of these cases. Let us say that this results in elastic, I will say, yes in other words, that is the first thing. Then, what do we do? We check sigma E and see whether sigma E results in a, what, in the plastic state or elastoplastic state. Why amI calling this as elastoplastic? Is there any, why is that I am not just saying plastic state? I can say that as well. No. The point is that it is not just coming back to elastic case. That is one thing, but more important is that I am not, I am not assuming that the elastic strains are negligible, absolutely. I am not assuming like that and I am saying that there is an elastic strain which witake into account as we go along. So, what do we do? So, we find out with sigma E whether sigma E has now crossed the yield point. So, the sigma E is greater than sigmay at n. I specifically put sigmay at n basically because, I can have a loading, unloading and reloading. Note that carefully. For example, coming back to this figure, coming back to this figure here, you see that I can go from this point to this point. That is fine that is very easy to understand. But, you can have a situation where you would like to look at a loading pattern where there is a loading, there is an unloading and then reloading. That means that I will go like that, I will come like this and then again go like this. Again go, in the sense that from this point I would go to this point. Then my sigmay is not this initial sigmay, but what has been updated? In other words, I have to update the yield strength of the material as well. It is very important point that yield strength of the material is also updated. When I update the stress, this goes along with it. Is that clear? That is why I said sigman. S,from here to here the procedure is exactly the same. This is the reason why we say that the behavior of elastoplastic material depends upon the history,the history of loading. What do I mean by history of loading? I have loaded here like that, unloaded, reloaded; that is the history. There is no one unique relationship between stress and strain. I cannot say if this is the stress and then that is the strain, no. This can be one, I unload, reload; that can be another and so on. There can be a number of situations or number of strains which would satisfy the stress. Is that clear? That is why I said that we will compare sigma E and see whether that is greater than sigmay. Let me say that it is, yes. Yes; absolutely. In other words, the points will always lie in a curve. This is what is called consistency condition. So, it will lie only in a curve in this case. In the multiaxial case we are going to define, already you know what is yielding. So, we are going to define an yield surface and we are going to say that a point will lie either inside the yield surface or on the yield surface. A point cannot lie outside the yield surface; it has no meaning in what we call as rate independent plasticity. I am introducing lot of terms, as I explain. Note these terms. I said history dependency or in other words, the stress or the strain depends upon the history and I now said, I used a term called rate independent. In other words, the behavior that we are right now seeing, in what you are looking at, is a behavior where strained rates have no role to play, the strained rates have no role to play. It is only the plastic strains which have a role to play. Where are the situations when strain rates have an effect? Yes, correct, but in solids, you haveanswered that fluids have an effect, viscoplasticity. So, whenever we have a situation where the temperatures are high, then we go over to a constitutive equation or relationship which is rate dependent. We call this rate independent. We call that as rate dependent or viscoplasticity where strain rates have an effect on the yield strength of the material and further behavior after . In viscoplastic materials, the situation is slightly different, let us not worry about it right now. But, for the case we are considering, for rate independent case, the material point or the stress, state of stress at a point should lie, should lie on the yield surface. We will come back to this may be after a couple of classes; may be next class or after that we will come back to this and let us see what we have to do in this case. Let us come back to this figure. What do we do with this? Can you have any clue? That is I am going from here to here. Yes; the first point is that I have to change a function at a point which lies between the two. In other words, my calculations are okay till this point, but after that my calculations are not going to be correct. So, I have a region. That is anyway that is epsilon, delta epsilonat a particular iteration. Note that carefully, we are in one iteration; one iteration. We have departed from that one place. So, if this is delta epsilon, then and this is delta sigma E, a part of it has to be corrected and another part need not be corrected. Let us say that this part, this part here which has to be corrected is called R in to delta sigma E, sorry, that is this part. That is the whole here, let me redraw it, so that remove all this andredraw it.That is the total delta sigma and that is the extra delta sigma that is there and that is what I have tocorrect. I? No; let, wone minute; let me draw a bit. That is the total. What is this?This is delta sigma E and that is the part which I have to correct; I have to correct. Let me call this as R into delta sigma E. So, 1 minus R which is this part, delta sigma E I need not correct, I can leave it as it is.I? So, my first job is to calculate R. In this case, R is less than 1, so that sum of this and this will give me delta sigma E. Now, this is my sigma E that is the total elastic prediction and that is my sigman and that is my sigmay.Those are the things I know. How do I now calculate R? Very good; R is equal to sigma E minus sigmay; sigma E minus sigmay divided by, yes, sigma E minus sigman. So that I can calculate. So, sigma E minus sigmay, I will get; I will get. Sigmay is already known. It is the yield strength of the material which is the previous step. So, this is already known. From this only I am going ahead. So, this is known to me and this is what I calculated. So, all these things here are known to me. So, I have to now bring this point down. I do not know this point. So, this is what I have to correct. Now, how do I correct that? The first step is, the first step is, I have calculate delta epsilon or split this delta epsilon as well. That is the first step. So, in order to do that, let me put some positions here say, let me call this A B C D E. Now, that is this total is my delta epsilon and this is my delta epsilon say, ep, elastoplasticand this is my delta epsilon E. Yeah;this is a straight line, A B D is a straight line. My whole idea is that if I calculate delta epsilonep, then I can calculate this point very easily by calculating ET at this point,in which case, this case it is only one value and multiplying by delta epsilonep. Is this clear? Now, how do I do that? How do I calculate delta epsilonep? What is that I know? I know R, I know this, I know this height, I know this height; that is the clue. I know this height AC, I know AE. So, how do we calculate? Yes; what is, what is the principle? That means that ABC and ADE are similar triangles. So, from this I can calculate delta epsilonep. So, BC by AC is equal to DE by AE. So, in other words AC by BC is equal to AE by DE. Is that clear? So, I can now calculate delta epsilonep. What is this AE by, sorry, AC by AE, can bring the other side? AC by AE that is R; R by 1 minus R you can write that term or you can sorry, AC by no, AE is delta sigma E and AC is equal to R into delta sigma E. So, AC by AE is equal to R, you can calculate that. Now, with that you can calculate delta epsilonep, which happens to be now equal to what, R into delta epsilon. Now, once I calculate delta epsilonep, my job is very simple. I will just remove it here and write that down. I multiply this delta epsilonep by ET and write delta sigmaep is equal to ET into delta epsilonep, so that sigma at n plus 1 for this iteration is equal to sigma at n plus delta sigma into1 minus, delta sigma E into 1 minus R plus delta sigmaep. Now, let us, let us recapitulate. There are lot of things we did. Let us recapitulate what all the things we did. First we are looking at one case, elastic to elastoplastic; that is the most complicated case. Now, what we said was because of my elastic predictor, I went from D to A which is not the correct position. I have to now correct my A. I realized that I have crossed the point, but I need not correct the whole of this jump; the whole of delta sigma E I correct. I need to correct only a part of it. Let me say that this part what I have to correct is R, so, R, so that this correct to be corrected part AC is equal to R delta sigma E and the part which I need not correct, which is CE, is 1 minus R into delta sigma E. So, I see that my calculation of R is very straight forward and that is equal to AE or AC by AE. I know AC because, I know this, I know sigmay and I also know AE because, that is nothing but sigma E minus sigman. So, I calculate R. Once I calculate R, I know this, no problem. But, I have to do some correction here. How do I do that correction? I find out of the total delta epsilon which I have completely invested in the elastic case, what is the part which actually is the elastic case and what is the part which is the elastoplastic case? The procedure is very straight forward when you look at the geometry of this figure and you see that by similar triangles what I have to do is to multiply this R and delta epsilon to get delta epsilonep. So, once I get delta epsilonep, I then multiply by ET to get delta sigmaep. So, once I get delta sigmaep, then I can calculate sigma. That is this plus this plus this. Note that, this that is sigman plus this, this height which is 1 minus R into delta sigma E plus, ,let me call this position as say F, CF. That CF is what I will get by ET, which is the slope of the curve into delta epsilonep, which is this. So, this multiply the slope of this curve which is that that is this curve will give me that position. Now, this is actually very easy to understand after I explain it, but has lot of problems. What if it is not bilinear, where the slope of the curve changes in the elastoplastic region, slope of the curve changes? Then, how do I calculate? I have to rely on my calculation, an explicit calculation ,I have to rely on explicit calculation, which means that I have to look at the slope at B and so on. In other words, this is, this comes under a category of what is called as an explicit technique.Though I am not very happy right now because of the issues I have raised, nevertheless for many problems this seem to work and when does it work? Especially when the load steps and hence the corresponding delta epsilons are not very large, this technique seems to be okay. We have to make some adjustments. We will see that what adjustments we have to make? But, many times when we are looking at problems, the time steps are small, this explicit time integration schemes are, sorry, explicit stress integration schemes are perfectly alright. Now, what kind of trouble can give? It is very simple to understand what kind of trouble it gives and it is very easy to understand it especially when I have a 1D curve. Now,look at this curve. This is no more bilinear curve. It is no more a bilinear curve. Now, I take a slope at this point. Let us say that, I mean let me it more flatly, so that the differences are clear. It is just an exaggerated case, let me take like this. Now, it is not very easy to find out the slope there, let me still find out the slope at this point. Let us say I am moving from this point and a elastically predicted and I have gone to this point. So, when I correct it, where will I correct? Yes; line which I have drawn, this point. I will only correct to that point. Let me zoom that out, so, I have two curves; that curve and that curve and I will now instead of falling back to this point, I will fall at that point. That means my consistency condition will not be satisfied. Now, what is happening here? Why is this happening? These are issues which wifocus in the next class. Yes; because we cannot; can we solve directly for and get this point? Yes, then equation becomes implicit. We will explain that in the next class and other issues which are involved with it, we will see in the next class. Yes, exactly. So, ETis changing. See, the whole idea here is first to explain to a method, then look at all the issues that come with the method. So,we will stop here and we continue in next class.