Lecture - 10 advanced finite elements analysis

In the last class we were looking at how to deal with incompressibility. In that process, we saw that one of the things that we can do is to do a selective reduced integration. One of the techniques what we saw was to modify D, but the other technique which we saw was to modify B. One of the comments at the end of the class, which is very important, is the way we write B. How do we write B? Because, this also depends upon two things; one is the order in which you write epsilon. Whether you write it as epsilon 11 22 33 12 23 31 or in some other order, 11 22 33 31 or whatever it is. So, depending upon that, B will change. That, that is very fundamental, may be all of you know it; please take care of that. The next one what also is important is the order in which you write u. You can write u in different orders, say for example, you can write this as u1 v1 u2 v2 and so on and call them as u1 u2 u3 u4 u5 u6 u7 u8, for a say, a 4 noded element. Of course, already you know that the size of this depends upon 2D or 3D. These are issues which I just want to bring it to your notice or you can even change this order. You can say u1 u2 u3 u4, then v1 v2 v3 v4; you should be consistent in writing it, when you implement it. So, that will settle; that is why there you might have felt that there is a difference in the way I wrote B. So, that is very obvious; as long as you know that B is that S into N into u and as long as you know this is what goes in, you can adjust B; that is not an issue. The second thing is that we said that this KT can be decomposed into two, into two parts. One is what we called as B dilatational plus deviatoric transpose say DT B dilatational plus B deviatoric. There was a question that it can be, there can be cross terms. Of course, there are going to be cross terms. There are various ways in which people have implemented this. In order to look at it very cleanly as a selective reduced integration, we do not take into account the cross terms. The question was am I correct in doing it? Yeah, you are not correct. If you remove the cross terms, you are not correct in doing it. But, what people have done over the years is to adopt their own procedure, so that their results are met, in the sense that the incompressibility condition is satisfied. In other words, what we call as B-bar method has taken different shapes, no doubt about that. In one of the implementation, people feel that they can write this P as B transpose say, sigma into d omega and say that sigma can be split into deviatoric and the dilatational part and use only B dilatational separately and B deviatoric separately. In fact, many of the flow formulations in plasticity which is used for, which is used for doing metal forming problems, they use this in a very simple fashion. They will not look at the cross terms, they will take B dilatational BT B dilatational deviatoric separately and then they add this together. They get the results, it works. But, strictly speaking you are not very correct in doing. That answers your question. So, you cannot just like that neglect. There are other ways in which people have done it. What are the other ways in which people have done? In implementations, there are so many implementations. The other way is to calculate B at a reduced Gauss point and call this as B bar. So, instead of evaluating it at say, 4 Gauss points which you do for this, you evaluate B dilatational at say one Gauss point. This is the standard procedure of use. So, this and this part, both of them are evaluated at say single Gauss point. Many people together they call this as B bar method. So, B bar transpose DT B bar and so on. The essence of the story is that the B term is split, is split into two parts and that the culprit, dilatational part is the one where maximum emphasis or maximum attention is given. This is what I explained in the, initially when I said that there are two parts to K matrix itself; when I look at it from D, aspects of D and that one fellow who is the culprit, who is going to make that infinity is the guy who is going to make, when nu is equal to say 0.5, make the, one of the terms to be infinity. But, having said this, then there may be a doubt. Is it that is it ad-hoc procedure? Is it that ad-hoc procedure came out of some sort of a justification from what I gave in the last class or does it have any variational basis to the whole thing? After all finite element may be looked at as an ad-hoc procedure, especially when you look at finite element from the perspective of the development, early development in the 50's, but it is actually not ad-hoc. It has a very sound mathematical basis. So, same question you can ask here. Does it have any variational basis, which forms the foundation stones of finite element analysis? The variational basis for this comes from what is called as mixed formulation. We are going to see how, when I develop mixed formulation, now, how closely we are into this kind of formulation? Now, what is mixed formulation? Mixed formulation ultimately leads to a variational form based on what is called as Hu-Washizu, Chinese and Japanese scholars who came together, who simultaneously proposed Hu-Washizu variational forms. There is a very interesting book by Washizu, if you have access to it please look at that - Variational methods in Elasticity and Plasticity. But, let, let us now concentrate on this. All over, formulations so far are based on displacements; displacements. Now, the current formulations what we are going to use are not only based on displacements, but also are based on stress, either stress or pressure; sigma or p as well as the volumetric strains, which is denoted by epsilonv or theta in different case. It can be either p or sigma, this or this and epsilonv or people can use also what is called as J. We will come to what is, what this J is in a finite deformation case later. Now, what is this mixed formulation? When we go to u, we go to what is called as an irreducible form. You cannot go below u; simply talking, you cannot go below u. u is the variable to which you can go down. When you are in sigma, it is not, we are not talking about an irreducible form. It is a reducible form. From sigma, I can go to strain; from strain, I can go to u. So, it is not an irreducible form. So, u is one where it is irreducible and straight displacement formulation. In mixed formulation, we have forms which are both reducible and irreducible, but interestingly we will remove the dependence of these later. We will see that we will remove those dependence, but simply means that you have three things. You will or in other words, you are going to treat them as if they are also, these two guys here, as if they are also variables like u. Is that clear? Having said that, let us develop this theory and see on mixed formulation. It is a very, it is not very difficult; it is quite simple. Implementation may be slightly difficult, but the theory is not very difficult. We follow Taylor and Zinkevichs as we develop this mixed formulation. Yeah, one of the, I think one comment is in order. There are other types of formulations as well and probably one of the comments which was passed in the, in my earlier course is that the strains have to be enhanced in order to meet certain demands, say for example, while bending. An ordinary element, a quadrilateral element cannot take bending and hence you have what are called as enhanced strains to take care of that, so that it can bend. You can have the, though you have 4 degrees of freedom, it can still take the shape as it bends. In other words, an element which is like this, in bending or rather, sorry, it takes like that shape on bending. It is not a nice shape which you would like to get, but you would like to get a shape which is like that and so on. So, that kind of shape is possible by enhancing the formulations or enhancing the strain calculation. Epsilon is not only a function of u of my original finite element, but also is enhanced such that I can handle such situation. So, epsilon consists of two terms; the regular epsilon terms, say let us say, R plus enhanced terms. So, this formulation is different from mixed formulation. You have also combination of these two as enhanced mixed formulation or mixed enhanced formulation. So, you have different combinations as well. Is that clear? So, now let us start our procedure. Let us start with the definition of what is called as m. So, let me define m to be 1 1 1 0 0 0. This is standard way in which Zinkevichs and Taylor defined ultimately the equations in terms of the matrix forms. In fact, just to illustrate that it is a matrix form, it is important that we write this term, so that Id matrix is defined as I minus one-third of m n transpose. This is very helpful to us to split the deviatoric part; the deviatoric part and the dilatational part. That is I is split into deviatoric and dilatational part, so it is very useful. So, we define an Id matrix as I minus one-third d. Remember that we already said that epsilon is equal to S u and sigma consists of two things. We will see what it is in a minute, but we can as well write that in terms of two parts here before we go to sigma dilatational deviatoric part and write epsilon to consist of Id into S u, first part of it, plus one-third m epsilonv. That is the dilatational and the deviatoric part of epsilon. You can tell me; look at this. What would be the dilatational part? Yeah, 1 by 3 m epsilonv; yeah, that is the two parts, we are splitting that into two parts and writing it. m is, what is m? m is this. It is I mean it is just as I told you, do not put a physical meaning to m. It is just a way of writing, so that you can write the whole formulation in a matrix notation. There is no physical meaning to m. It is just like a unit matrix you have; 1 1 1 0 0 0 1 0 so on. This is a way of writing the indicial notations or the tensorial notations into a, writing them into a form which is more amenable for implementation. That is all or else it becomes difficult to implement. Now, let us look at sigma term. Sigma can be split again as Id sigma; I am using this sigma with a v on top, because I am going to, please note all, all of them are matrices. I am not going to, you know the situation and you can write it. I am not going to every time put that as a matrix; may be we can look at this as an underscore. On the board when you write you want to put all the matrix form, it becomes bit complex. So, sigma is again split into two categories; that plus m into p. Epsilonv, of course you know, is i i plus j j plus k k, sorry, 11 22 plus 33 or epsilonii. So, again sigma is split into two terms the deviatoric term and the dilatational term. This sigma here is a guy who participates in the constitutive equation, in the constitutive equation and depending upon the way you write the constitutive equation, it can be written as say delta sigma and some function of delta epsilon. Remember, this we derived in the last class, the relationship between delta sigma and delta epsilon as DT. But, why we have written it like that is because, this formulation what we are going to develop is not only applicable for our small strain plasticity which we are dealing now, dealing with now, but are also applicable for a variety of other situations where incompressibility is an issue. Any other constitutive equation, incompressibility is an issue, you can use this. But, we are still in the small strain regime, so, you cannot use this formulation for large deformation. You have to make some amends if you want to use this for large deformation. Is that clear? Having said that let us look at the variational or rather the Galerkin formulation of this. I am not going to write the Hu-Washizu principle; we will, may be we will, go into details later. Now, let me look at the second step and I will explain what this means, this step means. Just to, you know, evoke your memory we said that you can write down, in the last course we said that you can write down, a variational form and then you can write down the minimisation of the variational principle which we have wrote it as delta I and said that the variational of this or the minimisation of this variational principle leads to virtual work principle. So, you can either go to virtual work or to variational form and we also made a comment that many times variational form may not be available. So, you stop with the virtual work principle, because you need a constitutive equation, if you have to go to the next step and especially in a situation like plasticity, you may not be correct to go to a variational form straight away and write it; may not be very mathematically very correct, because you do not have a strain energy function. Most of the variational forms are written in terms of strain energies and you do not have strain energies, relationship like what you have for elasticity, so, you may not be very correct. But, the next step where we write down the Galerkin form, virtual work form, is fine for many of these things. So, let me write down the virtual work form for this. Now, what I am going to write down also involves the inertia part of the force. So, if you want to include it, you can include it or you can neglect it as well. The first term is that inertial force rho u double dot d omega. That is the first term plus delta of S u transpose. Note what is, delta S u transpose is actually delta epsilon transpose, nothing else, sigma d omega. These are the two terms which happens to be available for us as the internal forces is equal to the external virtual work which happens to be delta u transpose b d omega plus delta u transpose t theta or d gamma. No; these are, these are see, the internal virtual work is equal to the external virtual work, our underlying principle. Last time we did that, so these two left hand side and the right hand side same, but the main difference is that we do not stop with this. We write down two more equations for p as well as for the volumetric strains. So, let us write down the second equation for volumetric strain. So, that is written as delta epsilonv, integral delta epsilonv, note this term carefully and let me have your comments, into one-third m transpose minus p d omega is equal to zero and next term, let me write down the next term; watch this carefully and tell me your comment, delta p into m transpose del sorry S u minus epsilonv d omega is equal to zero. Look at those two terms. Can you pass some comments? Delta epsilonv, epsilonv is the volumetric strain and delta p is the variation in pressure. We said that, we are going to have three things. Like we had delta u in the first case, here we had delta u in the first case; we should also have variation in p as well as in epsilonv. So, those two also enter into the picture. Delta u is the virtual displacement. Delta epsilonv, delta epsilonv is the variation in the volumetric strain. Delta p is the variation in the pressure. Please note one thing carefully. I mean, again I am commenting; we did that in the last class, because it is good to remember these things. In delta u you call this as virtual displacement. Now, you may wonder what this virtual displacement is. See, it is just that mathematically you may understand this equation, but as engineers many times you are more comfortable to look at terms like displacement, force, pressure, energy and so on. So, just to give you a physical flavour to the whole mathematics, we call this as virtual displacement, we call this left hand side as internal virtual work and call that right hand side as external virtual work and so on. In actuality, if you look at it, there is no need for you to call delta u to be a virtual displacement. You can just leave it as, as it is. If you want, even you can call it as virtual temperature; it does not matter, as long as you understand that it is a variation of u. Variational calculus does not depend upon you calling this as virtual displacement. Let us now concentrate on these two, these two equations. Can you comment on this, how we got these two equations? Very, very simple, straight forward as we did; if you want, you go back and look at the formulations which we did for our variational formulations, which we did for this. How we started from equilibrium equation and came here; you see it. Can you comment on it? Weighted residual of the equation, but, but how did I put it as ? You are absolutely right, it is a weighted residual. But, how did I put it as zero? What did I do? You look at that slightly more carefully. Absolutely; quantity inside the bracket should be satisfied at all points and this has to be equal to zero. Actually, what is that we did? p is equal to nothing but one-third of sigma. So, this is nothing but zero multiplied by some quantity. This is exactly what we did, if you remember on our variational form itself. When we got the variational form, this is what we said; we said we can, we can play; either from equilibrium equation you can go to the, to the variational form, from variational form or Lagrange equation you can come back and so on. In one side, one route when I took what I did was exactly the same; anything multiplied by zero is equal to zero. It has to be satisfied at all points and in other words, when, when I discretize it, as our friend here said that it is a weighted residual of this quantity. Look at that. Here again, this is epsilonv, of course, I hope you know. Here again it is exactly the same. One-third here and m transpose epsilon which happens be epsilonv. So, this quantity has to be zero that multiplied by variational form of delta p. So, these are the three equations which we now use in order to develop the finite element solution. Is there any question? You can combine them together and write also Hu-Washizu variational principle, but I will defer that till we develop finite strain plasticity, sorry, finite strain elasticity and plasticity. At that time we will write this term, because now it is very easy to understand at this stage. But, may be I will leave it as an exercise. In the same fashion as we did in the first course, when we converted this variational form into a functional, you can try to convert this into a functional form. It is very simple. Hu-Washizu, though it is a brilliant idea, has come from very, very simple concepts; nothing very difficult to understand. Now, watch the next step what I am going to do. Now, I am going to do exactly what I did for N or for u, for p, as well as for epsilonv. In other words, I am going to write u is equal to, my good old friend, shape function into say, u hat. What is u hat? You remember, yesterday we said that, correct, they are the displacement at the nodes. I am going to do the same thing for pressure. I am going to say that p is equal to N p hat and we are going to do that for epsilonv. So, epsilonv is equal to N, sorry, let me do that as manipulative, epsilonv hat. But, so, these are shape functions. But these shape functions, of course all of you know, depends on how many variables that you use. The shape functions if they have to be the same, then the number of variables for u, p and epsilonv have to be the same. But, if they are going to be different, then the shape functions have to be different. So, they are going to be different due to certain very important conditions that we require. So, we write them separately; we write them as Nu, Np and Nv separately or in other words, it is not necessary that all of them are the same. What does it mean? It simply means that if I have an element or triangular element, whatever it is and at the nodes I am going to keep the degrees of freedom for displacement, it is not necessary that at the same node I will keep for p and at the same node I will keep for epsilonv as well. They, they can be at different places, may be at the centroid only. No, no; you would not have studied these three things separately. It is, you use please note the difference. You will use the same shape function for all u's. We are not touching them. Whether you are going to interpolate u1 u2, it does not matter. There you will use exactly the same shape functions. We are not talking about that. We are talking about a situation where I have pressure as well as epsilonv to be also interpolated. They are going to be different, they can be different. Again as I told you, this is due to certain conditions which I require and we, we will call that later as LBB condition and so on. Let us not get into lot of, lot more mathematics at this stage; let us just look at the end result. Is that clear? So we will now make Np is equal to just Nv alone and leave N as it is. Any confusion? No, no; number of shape function, Nu is what I am keeping for displacements. No; N's depend upon, of course, u's, u's; number of, number of degrees of freedom. So, whether, the whole question, I am repeating, is whether I have the same degrees of freedom for u, p and epsilonv. I can have p to be constant in one element. In other words, N can be 1, epsilonv can be 11, which would, which means that I am keeping this to be a constant. I am going to evaluate it, such that Np is equal to, have to have a single value. Is that clear? That is what we mean. So, this is the most confusing step. We are also putting down a condition that Nv is equal to Np. 11 plus 22; 11 plus 22 plus 33 is always there; depending upon whether you have, you know, 33 is equal to zero or 33 is there and so on. See, epsilonv does not change; 11 plus 22 plus 33. Of course, it will have. Depending upon whether the incompressibility is equal to zero or not, so epsilon, once you have for an elastic case, you have epsilon11 and epsilon22, you will have a volumetric strain; of course, you will have. Now what do we do? Go back to and look at, remember what we did? I have to write down all the variations in terms of N's. That is in other words, u and epsilonv and epsilon and delta epsilonv. So in other words, delta u is equal to how do I write that? Nu into delta u hat; please write that down. Delta p Np delta p hat and delta epsilonv is equal to Nv into delta epsilonv hat. Now, what is my next step? Remember my previous steps, I calculated what is called as B. Now, I have not done that yet. So, what is the B that I will get or what is the relationship between strain and displacement? You remember, I had already written down epsilon. We had epsilon, we had two terms. We just now wrote that epsilon is equal to, what is it that we wrote? Id into S u plus one-third m epsilonv. So, in terms of the discretized quantities how will this appear? How will this appear? First term will be Id into, instead S u this will become, B u hat, same fashion as you write down the B matrix, plus the second term; second term we had as one-third m epsilonv. So, what will that become now? Nv epsilonv hat. So, one-third m Nv epsilon v hat. See the procedure; please note this procedure is the same. Follow this step by step. Remind, I mean keep in mind what you did in the last class. What is delta epsilon? So, what am I doing? What I am doing is very simple. I am looking at where all I have to substitute things there, I am going to keep writing them in that fashion; then go back and substitute. Delta epsilon is, of course, I can write that as say, same thing, just I will change it; delta epsilon here and delta u hat and delta epsilonv hat. What is now sigma? Remember what we did for sigma? What do you do? Just substitute now in terms of my discretized quantity. So, sigma has again two things. What are the two things or what are the two terms? Id into plus m into p. Now, what will happen to p? Np p hat, m Np p hat. Having known all this, substitute these things into my first equation. I am going to get, obviously I am going to get three equations; I am going to get three equations. Substitute this; I will give you two minutes. Please look at this and this, substitute it and let us see how, what, what you get? Please look at this equation. Of course, you should have written these things in your note book. So, these equations, these equations, substitute it here. Let us see what you get? Of course, delta u transpose you have to delta u hat N transpose. Note that, this is the inertia term. If you are doing dynamics, this term will be there; u double dot term. If you are not doing dynamics, you can, term is out. So, please look at these two terms and write that term. Let me write this ultimate result. I am not going to say what this is; of course, all of them are matrices is equal to f. Let me see, you come up with an equation for say P and M. It is very simple. It is, there is nothing there; P is nothing but, look at that equation, there is, there is no change. P is nothing but B transpose sigma; no more than that and M is so, P is equal to from where do I get? No; this I get from P, I get from here; delta epsilon transpose. So, delta epsilon transpose is delta u hat transpose; delta u hat transpose will go out. So, this term I can write this down as delta u hat transpose integral B transpose sigma d omega. So, P remains the same, B transpose sigma d omega. My M, M will remain the same as you have it in a dynamic case. What is that? Rho, substitute it here; rho is scalar quantity, so, N transpose this, this will become, this will go out; delta u hat transpose N transpose, u double dot you have to substitute; so, this will become N u hat double dot. So, M is now, that is the quantity, that is the quantity. M is rho N transpose N d omega; your regular M, nothing has happened to it, is equal to f. f will consist of these two terms. So, this will become N transpose B d omega N transpose t d gamma; will be the two terms for the external forces f. So, the first equation is the same. Now, there is no change in it. Now, I am going to leave that second thing for a minute and let us see what you get out of the second equation? Now, this is clear. Yeah, delta u hat transpose will go off, because obviously you know that from your earlier class, then because it is valid for every delta u hat transpose, so that will go off; so, it will not be there. Let us now look at the second equation. Write this down for the second equation. Let me write this in Zinkevichs terminology as Pp, the, the final equation minus C p hat is equal to zero. See, we are writing these things at the element level; understand that. So, now what is Pp and write that down as well as the third equation; that is the first equation we know, second equation and third equation is written as minus C transpose epsilonv hat minus E u hat is equal to zero. Now, I will give you two minutes; let me see how many of you come up with what Pp is. What is the procedure? Very simple; in order to get Pp I substitute, in the same fashion as I did in the first equation I substitute, all these discretized forms into this equation and rewrite this first equation. Is that clear? Let us see what Pp is now. So, one-third term will be there; very simple. One-third, one-third term will be there. So, this will be delta epsilonv transpose. So, here you will get Nv transpose from here. If you want you keep this here. So, you will get, in delta v, delta epsilonv you will have, Nv rather. Then, I think we should have written that as N delta epsilonv transpose, so that you can take that out. So, N transpose one-third will be there, of course, you can write that out, N transpose. Then, m transpose that will be there. For sigma hat, what do you do? You can keep that as it is. Yeah, substitute it. Nv, of course, d omega. Yes, this may be a worrying term. We will see how we can eliminate that. So, that is the Pp term. Obviously I have, just substituting it straight away. In a proper matrix notation this should have been transposed, so that is why we kept that. Then only you can multiply both. Then, how does the second term come up? How does the second term come up? What are the terms that will be there? That is C term. Is this clear? That gives me Pp. Yes, so, Nv transpose; yeah, this will be there; this or this fellow will go out. So, Nv transpose. For P, I am going to substitute for Np, so, integral Nv transpose Np d omega. That is all. So, this will now become Np p hat. So, C now will become integral Nv transpose Np d omega. p hat; p hat, how did p hat come? Because, p is equal to Np p hat, so, that should come there. Now, third equation; having done all these two, look at the third equation and substitute. Third equation is slightly more involved and let me write, I, I think I have written down that here. So, your whole job is to find out what C transpose and E is; slightly more involved, especially E; you have to define to me what that term E is. It is not that mathematically very difficult; one or two step of algebraic manipulation will give you what E is. But, the procedure is direct. Substitute this back into that equation and determine. So, that is going to take few minutes. We will stop here, please try this out. We will carry on with the third equation in the next class. Any question? We will stop here and continue in the next class.