Topics: Generalized coordinate finite element models
Instructor: Klaus-Jürgen Bathe
Study Guide (PDF)
Sections 4.2.3-4.2.5, 4.3, 4.4.1
The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high-quality educational resources for free. To make a donation, or view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.
PROFESSOR: Ladies and gentlemen, welcome to lecture four. In the last lecture, I presented to you a general formulation for finite element analysis. In this lecture, I would like to talk to you about the derivation of specific finite element matrices.
The formulation that we will be using, or the derivation I should say, that we will be using, leads us to generalized finite element models.
I will be talking later on about what we mean by generalized. Let's remember, let's recall what we arrived at in the last lecture. We had that the element stiffness matrix for element m is obtained by this integration here where we integrate the product of a strain-displacement matrix times a stress-strain matrix times the strain-displacement matrix over the volume of the element.
We also had derived the expression for a body force vector, B. m, again, denoting element m.
And here, we integrated over the volume of the element the product of the displacement interpolation matrix times the body forces in the element.
We also had derived a surface force vector of the element in which we integrate the product of the surface interpolation matrix times the surface forces applied to the element.
So in essence, we see then if we look at these expressions, we see that we need really a B matrix, a strain-displacement matrix for the element. Of course, we need the stress-strain law for the element. And that will really be given to us from the serial strength of materials.
We also want the displacement interpolation matrix for the element. Once we have this displacement interpolation matrix, we can directly obtain the interpolation matrix for the surface displacement of the element by simply substituting into this matrix here, the surface coordinates.
In summary then, we need the H matrix, the B matrix, and the stress-strain law for every element that we have in the assemblage.
Notice that the stress-strain law, of course, must depend on what kind of element we are looking at. How we have mechanically idealized our system. In a truss structure, we will have simply the Young's modulus in here. In a plane stress idealization, we will have a plane stress law. In a plate bending idealization, we will have the plate bending stress-strain law in there, and so on.
Well, I want to talk today about these matrices. How do we obtain these matrices?
Once we have talked about that subject, I finally want to also spend a little bit of time talking about the convergence of the analysis results. Of course, in finite element analysis, we must be concerned about the accuracy of the analysis results that we are obtaining, and about the convergence of the results to the actual analytically, theoretically accurate result that we want to obtain.
Well, let us then first talk about the derivation of the H, B, and C matrices. And the H, B, and C matrix, of course, depend on the particular problem that you're looking at.
Let us therefore categorize the different types of problems that we encounter in structural analysis. Here we have a truss structure, and a truss structure is really an assemblage of truss elements. And in each truss element that lies between the nodes, we have one-dimensional stress conditions. In other words, the only stress that we see on the section AA would be the tau xx stress. And that is the only stress, no other stress is involved.
Of course, just going a little bit ahead, we immediately would also say that the only strain we would be interested in is epsilon xx and the stress-strain law is simply Young's modulus.
The second kind of analysis that we might be performing is a plane stress analysis. In plane stress conditions, the only stresses that are non-zero are tau xx, tau yy, and tau xy.
Remember in lectures three, I had listed all six stress and strain components. And what I'm doing now is I take out just those stress and strain components that are really applicable for the kind of problem that I want to look at.
Plane stress situations we will encounter in the analysis of a beam, for example. Here's a simple cantilever beam subjected to surface loading p. And if we look at an element here, it is redrawn here. We would find that tau xx would be a stress into this direction. Tau xy is the shear stress, and tau yy is that normal stress. The cantilever, of course, can be subjected to a surface loading, a bending moment, tip loading, whatever loading. The stress components, the only stress components that we assume of course, act in the structure would be these three stress components. In fact, of course, in the actual physical situation, we would find that around here, depending on the specific support conditions that we are using, we would have a much more complicated stress situation than the plane stress situation assumed.
However, if we are only interested in the tip displacement here, or the displacements away from the support condition from the support. Similarly, if we are only interested in the stresses and strains away from the support, certainly it is accurate enough to-- it yields an accurate enough model to assume a plane stress situation.
Another plane stress situation occurs here, for example, where we have a sheet with a hole. If the sheet is thin, the only stresses that we would see on an element, such as shown here, would be z stresses, the plane stress stresses.
A next category of analysis that we would encounter in practice is a plane strain condition. Plane strain condition means that in the space xyz shown here with the displacements uvw, the displacement w is 0. It's identically 0 if our x, y-coordinate system lies in the plane of a typical section of the total problem. And let me show you one such typical-- or describe to you one such difficult problem.
Here we have a dam that is very long. In fact, we might assume it infinitely long. And it is supported on this face here, so that w, the displacement, normal to that face is 0. And the same support conditions are at the other face. If the dam is subject to uniform water loading, water pressure loading all along its side, then by symmetry it follows that w, the displacement, at any section of the dam, the displacement normal to that section must be 0.
Therefore, we have w 0 meaning directly epsilon zz 0. The strain into the z-direction is also 0.
Similarly, the shear stresses, tau zx and tau zy are 0. And the only stresses that we are left with which are non-zero are those shown here. These are the stresses acting on a typical slice of the dam. A slice say, of unit length would be subjected to just these stress conditions. Therefore, in that array of six stresses and strains that I talked to you about in lecture three, we really only are concerned-- you only need to calculate these four components of those six that I talked about there.
Another stress-strain condition that we are interested in frequently in practice is an axisymmetric analysis condition.
This case we have an axis of revolution for the structure, and here we look at a cylinder, which is axisymmetric in geometry about that center line here, the axis of revolution. And the loading is also assumed in this particular case to be also axisymmetric. So typically, we are talking about a cylinder, circular, of course, when looking from a plane view subjected to internal pressure, which is the same all around. And of course, the geometry is also axisymmetric when we look at the axis of revolution, the center line of the cylinder.
In this particular case, the only stress conditions that we need to calculate are shown here for an element in that cylinder. We have tau xx, tau yy, tau zz. This, of course, being the hoop stress. And the shear stress, tau xy. These are the stress conditions that we need to calculate for this particular case.
Similarity, because we are talking about an axisymmetric condition, an axisymmetric condition where really, every one of these elements behaves in the same way, we can look at one radian of the total circumference in the analysis. So this is an axisymmetric analysis, which we also are concerned with frequently in practice.
Finally, the very large number of analyses-- in a very large number of analyses, we consider plates and shells.
If we look at a plate, an element in that plate, we would find that in general we would have the tau xx, tau yy. These are the normal stresses. And tau xy shearing stresses in the plane of the plate. We would also have transfer shearing stresses, tau xz and tau yz. Sorry, tau yz, shown here as the transfer shearing stresses.
Of course, these transfers shearing stresses give us shear forces. And these stresses here give us membrane forces plus bending moments and twisting moments.
Elementary theory of plates shows us, tells us, that these are the significant stresses in a plate. The important point is that in a plate, as well as in a shell, the stress normal to the plate is 0.
The same holds for the shell. The stress normal to the shell is 0, and we have the stress components. However, we have to now remember that the xx and yy axes are aligned with the curvature of the shelf. Are aligned with the curvature of the shell. And then these are the stress components that are of interest in the analysis of shells.
In many cases, when we look at the analysis of plates, the membrane forces-- the total membrane forces are 0. In which case, of course, we assume that tau yy has a linear variation through the thickness of the plate. As stipulated, for example, in the Kirchhoff plate theory.
This much then for the static components, or rather the stress components of plates and shells when we look at the kinematics of plates and shells. Just to remind you, we assume that there's a normal to the mid-surface. And that normal in Kirchhoff plate theory remains during deformation normal to the plate surface. And all material particles on that normal AB remain on a straight line. These are the two things. I repeat, all material particles on the line AB remain during deformation on the line AB, and the line AB remains normal to the mid-surface in Kirchhoff plate theory, which does not exclude shear deformations.
However, if we use a theory that does include shear deformations, we would still have that the material particles on line AB will still lie after deformations on the line A prime B prime say. But the line A prime B prime-- B is still a straight line-- is not anymore normal to the mid-surface. In this case, we have included shear deformations.
We will see in the next lecture that, in fact, it is simpler to use a theory including shear deformations to derive effective finite elements.
Well, this then is a review for the kinds of stress components and strain components that we are interested in. Of course, that information which I discussed so far is available from the strengths of material theory and continue mechanic theory. And as an engineer, however, we have to be able to decide what kind of stress and strain situations you want to analyze. And there, of course, is already one very important assumption. A very important assumption in laying out or starting with the right model for the actual physical structure. Whether we want to model the actual physical structure as a truss assemblage, as an assemblage of plates, beams, et cetera, all of those situations can be covered in the finite element theory. But we have to make the right assumptions right from the start in starting with the right model.
So here I listed ones. In one table, the particular problems that we might encounter. We might have a bar. The displacement components then being simply u.
In a beam, which I have not given earlier, but here we have a transverse displacement, w.
In plane stress situations, we have two displacements, u and v.
In a plane strain situation, u and v.
In an axisymmetric situation, u and v.
Three-dimensional situation, we have u, v, and w.
In plate bending, simply the transverse displacements.
In analysis of shells, in the analysis of shells we would, of course, have the transverse placements and also the membrane displacements included. We will talk about that in the next lecture.
These are the displacement components for the particular problem that we are looking at. And correspondingly, we have particular strains. The strains that we are talking about are in the bar, simply the normals strains.
In a beam, the second derivative of the transverse displacements. Here we have kappa xx defined. This is kappa xx defined.
In plane stress, we have these three components.
Plane strain, these three components.
Axisymmetric these three components, and so on.
And it is here, of course, where the engineer has to decide what kind of problem he has in hand when he actually-- when he looks at an actual physical situation.
Corresponding to these strains, we have also stresses. And here we are listing the stresses.
Notice that for the beam, we are talking about a bending moment, which we might consider to be a generalized stress.
In plate bending, we talk similarly about bending moments, which you also might consider as being generalized stresses.
Otherwise, we have the actual stresses that correspond to the strains.
Having decided what stresses and strains we want to look at-- in other words, what kind of problem we have at hand, we also have to select the appropriate material matrix. Again, this one is given form strengths of material.
Here, for a bar, a truss structure, we simply have Young's modulus.
For a beam, we have the flexural rigidity.
In plane stress analysis, this would be the material law that you would find in many textbooks on strengths of materials.
Well, having therefore resolved the problem of how to establish the stress-strain law, or the stress-strain matrix c, we now go on to discuss the construction of the displacement interpolation matrix.
For a one-dimensional bar element, we would have simply-- we can assume simply a polynomial with unknown coefficients, alpha 1, alpha 2, alpha 3. And this is simply a polynomial nx. x being, of course, the coordinate that is being used to describe the displacement for the element.
For two-dimensional elements, here we're talking about plane stress, plane strain.
Axisymmetric analysis, we would have, similarly, these assumptions here. alphas for the u displacements and betas for the v displacements.
For the plate bending element in which we only discretize the w displacement, the transverse displacement, we would have this assumption here.
For three-dimensional solid elements, similarly, we have these assumptions.
Notice that these gammas, alphas, betas are unknowns for the element. And we will actually relate these coefficients to the nodal point displacements of an element.
If we look, therefore, in general, at the description of an element, we really have on the left-hand side here, the actual displacements, the continuous displacements in the element being given by matrix phi, which contains the constant and x, x squared, x cubed terms, and so on. And this alpha vector here contains the generalized coordinates, the alphas, betas, gammas.
Typically, therefore, just to focus our attention on a particular case, we might have if we say that u of x for a truss structure is equal to alpha 0 plus alpha 1x plus alpha 2x squared. In this particular case, we would have that we can write this, of course, as 1 x x squared in matrix form times alpha 0 alpha 1, alpha 2.
And here, I'm talking about this being equal to the phi, and this being equal to the vector alpha. So just to give you an example for this particular case.
If we now want to evaluate the alpha values. Of course, in two-dimension analysis as I shown earlier to you, we have also betas here. In three-dimensional analysis, we also would have gammas in this vector. If you want to evaluate these alpha values, then we can simply use this relation here, or that relation here, and apply it to the element nodal points. To the displacements at the nodal points of the elements.
If we do that, in other words, we substitute into here, into the phi matrix, the actual coordinate of the nodal points, and on the left-hand side, of course, the displacement at that nodal point, we directly generate this relation.
We can then invert this relation to obtain alpha. And here, we have now the generalized coordinates in terms of the nodal point displacements. This then would really complete the evaluation of the continuous displacements u in terms of the nodal points displacement u hat.
If we have laid down this assumption, we can directly also obtain the relation for the strains of the elements. We have decided what kinds of strains we want to look at that we have to include. And by proper differentiation of the rows here, we directly generate this relation. Of course, this is a [UNINTELLIGIBLE] relation then for the material of the element. And we then directly have tau in terms of alpha. But alpha was already given in terms of u hat here.
In summary, therefore, we have our displacement interpolation matrix here, given as phi times a inverse. This is obtained by simply substituting from here, substitute from here, directly into there. And you get that u is equal to H u hat. And our H being phi a inverse.
Once again, the u is the continuous displacement in the elements. u hat has nodal point displacements in the elements.
Similarly, if we substitute from here into there and recognize that our epsilon strain is equal to B times u hat, we directly obtain that B must be equal to E times A inverse.
Well, let us now look at a particular case, a particular example to demonstrate how we proceed for a particular case.
Here I'm looking at the analysis of a cantilever subjected to a load. And I've idealized that cantilever as an assemblage of four elements. Each element is a four nodal element. So for element 1, we have nodes 1, 2, 4, and 5 of the global assemblage. Notice that all together we have 9 nodes. And once again, 4 elements.
The cantilever structure is SiN and we believe that it is appropriately modeled using a plane stress analysis. In which case, the only stresses that are of concern are the stresses in this plane, tau yy, tau xx, and tau xy. So these other stress components here are 0.
If we look now at a particular element to focus our attention on one element, because all these elements are really the same. Basically, the same. Let's look at element 2.
And we notice that element 2 being a four nodal element is shown once more here with the displacements U1, V1 at the local nodal point of the element 1. V2, U2 being the displacements of the local nodal point 2. And here, V3, U3 similarly for nodal point 3. And V4, U4 for nodal point 4.
We have taken this element of the assemblage and what we want to really do is derive the displacement and strain interpolation matrices for this element corresponding to U1, V1; U2, V2; U3, V3; and U4, V4.
If we have done that, we recognize, of course, that U1 corresponds to a global displacement. Similarly, V1 corresponds to a global placement. Then we can directly construct the displacement and strain interpolations corresponding to the global displacements. But the finite element procedure is really to look locally at an element. Locally at an element. And work in terms of local displacements. Later on then, we relate the local displacements to the global displacements.
The advantage of proceeding this way is that we really-- in looking at one element, we look at all four simultaneously.
The procedure then is, once more in summary, we want this matrix here, the H2 matrix corresponding to element 2. The vector U lists all of the displacements, all of the 18 displacements.
However, we will look at the element 2 with local nodal point displacements, forgetting about the global nodal point displacement. And our objective is to derive a relationship of u equals H times u hat, where u hat lists U1, V1; U2, V2; U3, V3; U4 and V4.
Once we have this H, which is, of course, a 2 by 8 matrix, because there are 2 displacements, u and v. The u and v displacement for the element. And there are 8 entries in u hat. Therefore, this H matrix must be a 2 by 8. Once we have obtained this H matrix, we can construct directly our H2 matrix that we talked about earlier.
Well, let's proceed then in the way that I discussed in general form earlier. Since we have a four nodal element, let me just sketch it once more here. Since we have a four nodal element, and we are talking about u and v displacements, we really can only have four constants, four generalized coordinates. Generalized coordinates corresponding to the u and corresponding to the v displacements.
Remember that the only unknown u displacements are these 4. And these 4 u displacements will give us alpha 1, alpha 2, alpha 3, and alpha 4. The same holds for the v displacements. This is a bilinear approximation of the u displacements in the element.
We write this in matrix form as shown here. Noting that this capital Phi, notice this is a capital Phi because there are these cross bars. It's easy to confuse the capital Phi with the lowercase phi, which doesn't have the cross bars. So this capital Phi was the cross bars here is defined in here. Where the lowercase phi is simply an array of 1, x, y, x, y.
This relationship here is nothing else than that relationship, but written in matrix form.
The alpha lists all the generalized coordinates. The generalize coordinates. And this is the reason why we call this the generalized coordinate finite element formulation for this particular element.
Let us now define the u hat vector, which I refer to earlier. This case, the u hat vector involves 4 u nodal point displacements and 4 v nodal point displacements.
Now if we apply-- and this is really one important point. If we apply this relationship here, which of course must hold all over the element. In particular, it must also hold for the nodal points of the element.
If we apply this relationship for the nodal points of the elements, we can apply it eight times because we have 8 nodal point displacements.
We directly generate this relationship here. u hat equals A times alpha, where A involve the coordinates of the nodal points. Hence, we have our H Phi times A inverse the way I have derived it earlier. So this is a very systematic approach of deriving the displacement interpolation matrix. It consists of evaluating the A matrix. And of course, we assume that that A matrix can be inverted. This is, for this element, certainly the case, as long as we are talking about the rectangular element or not very highly distorted elements. We can for the four nodal element, invert this matrix and directly get this relationship.
Let's look now, once more back at the analysis that we wanted to perform. Remember that we looked at this element here, and we have now obtained the H matrix, the displacement interpolation matrix for this element. Once more, here is the H matrix. In fact, this H matrix here, since we have no superscript attached to it yet, really holds for every one of the 4 elements that we are talking about. Provided we put into the u hat here the appropriate nodal point displacements of the elements. We can directly say that this matrix here holds for any one of the elements.
If the elements, if the 4 elements are not completely identical, of course, the A inverse will change from element to element.
Having gone through that step and having looked at element 2 once more, we directly recognize how to establish now the H2 matrix.
Remember our u1 will be related to the global U11. Our v1 will be related to the-- or will correspond I should rather say, to the global U12. Let's look at why this is the case.
Why did I write down U11 here and U12 there? Well, let's look at this point there. Which is of course, the local nodal point 1. The local nodal point 1, right upperhand corner of element 1, corresponds here to nodal point 6.
Well, the global nodal point displacements of nodal point 6 are U11. Let me put them in. And U12. Why is that the case?
Well, if we start numbering here-- and I'm going to include now, I'm going to include now the displacement of these nodal points, although we know that they are actually 0, that boundary condition we can impose later on.
So if I start numbering here, I would have U1, U2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12. So the global nodal point displacements are U11, U12 corresponding to nodal point 9. And those are the ones that appear right here corresponding to the local U1 and the local V1.
Well, if we proceed in the same way for all of the other nodal points for element 2, and these are the results. You can look these up in the textbook figure 4.6, or in the accompanying notes. We can directly construct the H2 matrix for the element. And now the H2 matrix, of course, involves as the vector the global nodal point displacements, U1 to U18.
The H matrix for element 2. And if the 4 elements are identical, we will have the same H matrix for all the elements. The H matrix here looks as shown here.
Notice that we have a one quarter there. It's a little arrow here. The bracket should've been there. Bracket should have been there. So we have a one quarter outside of the bracket. And we have as the first entry, this one here. This corresponds to U1. This is the entry corresponding to U1. I left the entries corresponding to U2, U3, U4 out. This is the entry corresponding to V1. And I left out the entries corresponding to V2, V3, V4.
This is our H matrix here. And now we basically want to blow it up to obtain the H2 matrix corresponding to element 2.
Well, if we recognize as we just did that U1, the local U1, the local V1, correspond to capital U11 and capital U12, and I just showed you that that is indeed the case. We would simply take the first element here, this element here, and put it right there. This element here would go right there. This element here would go right there. And this element here would go right there. In fact, H1 5 is 0 as you can see here. And H2 1 is 0. So H1 1 is this term here. And H2 5 is that term here. Those are two non-zero entries in the complete H matrix.
Well, this is then the way how we actually can proceed to construct the H2 matrix.
I pointed out earlier in lecture three already that we do not perform this step computationally. In fact, what we do computationally is we stay with this H matrix. We calculate a compacted stiffness matrix, a compacted load vector, which the stiffness matrix will be of order 8 by 8. And then we assemble these element matrices with identification arrays and connectivity arrays into the global structural stiffness matrix, into the global load vectors. So we really don't perform this step of blowing up the H matrix. And I will show you in the next lecture how we actually perform the computations of assembling element matrices.
However, for theoretical purposes, it is neat if you understand how we can construct this H2 matrix from this H matrix. And of course, similarly, we could construct H1, H3, and H4 from this basic H matrix.
Once we have done that, we can directly write that k being equal to the sum of the km. In other words, the global structure stiffness matrix is equal to the sum of the element matrices, stiffness matrices, where each of these km's has the same order as the structure stiffness matrix k.
Let us now look briefly at the development of the strain displacement matrix. In plane stress conditions, we have three entries. The strains are given as shown here. Elementary strengths materials tells us that these are the strains. And I showed earlier that the B matrix is constructed by multiplying the E times the A inverse. The A inverse was already used in the calculation of the displacement interpolation matrix. The E matrix is obtained by appropriately differentiating the entries in the displacement interpolations. And this differentiation is performed on the phi matrix right here. Or rather, on these entries here. By simply differentiating these entries here, we would obtain the strains and writing those in matrix form, we obtain the E matrix here.
Once we have the B matrix, which in this particular case is 3 by 8 matrix for a typical element, we can, of course, construct the B2 matrix for element 2 just in the same way as the H2 matrix. Of course, we would again, blow up the B matrix from a 3 to 8. 3 by 8 I should rather say, to a 3 by 18 now because we have actually 18 degrees of freedom in the complete structural model.
This then completes what I wanted to say about the construction of the displacement and strain-displacement interpolation matrix for the element, the H2 and B2 matrices. However, let me just spend a few moments on the construction of the surface displacement interpolation matrix. In other words, the interpolation of the displacements along the surface of an element.
Here, assume for example, that this is a surface of the elements that we are concerned with. And that there would be some loading applied to that surface of the element. So that we need, in other words, the HSM matrix the way I have been pointing it out earlier.
What we would do to obtain this matrix here, is really simply use our earlier results for the HS-- sorry, for the HM matrix. And here we have the H matrix and the H2 matrix for element 2.
What we would do is we would use this matrix, or this one of course, this matrix is contained in this one. And simply, substitute the particular coordinates along the surface.
Now, going back once more to the element under consideration. If we are looking at this surface, then y is constant and equal to this distance here. So knowing y, we would simply substitute that value of y into this function here, and thus obtain the HS matrix.
So, in general, therefore the HS matrix is simply obtained from the H matrix, or the HS2 matrix here would simply be obtained from the H2 matrix by substituting the coordinates along the line, or along the surface that we are considering.
Let us now to look at the overall process of the finite element solution. Here I have summarized what we all are concerned about. Of course we have the actual physical problem that we want to analyze. We have a geometric domain that we have to discretize. A material that has to be represented, a loading and the boundary conditions. All of those have to be represented.
In the mechanic idealization, the idealization step that is a very important step that I referred to earlier, we would idealize the kinematics as being a truss, plane stress, three-dimensional, Kirchhoff plate, kinematics. The material isotopic elastic, the loading, and the boundary conditions. And this mechanical idealization really gives us a differential equation or equilibrium for one-dimensional stress situation, a bar. We will have an equation such as this one.
In the finite element solution then, we operate on this equation basically. We want to solve the governing differential equation of the mechanical idealization.
The errors that we are therefore, concerned with are discretization errors by the use of the finite element interpolations. And then, a number of other errors. Numerical integration in space. If we use numerical integration, there would be other errors introduced. In the evaluation of the constitutive relations, again, errors would be introduced. The solution of the dynamic equilibrium equations, by [UNINTELLIGIBLE] time integration, mode superposition further errors would be introduced.
If you solved the governing equilibrium equations, ku equals r, by iteration there would be further errors introduced. And of course, finally there are round-off errors on the digital computer.
What I want to do now in discussing convergence, I want to say that all of these errors are negligible. In other words, we are basically evaluating the integrations in space exactly. There is no round-off because we have an infinite precision machine, and so on. So all the errors that we are looking at now are the discretization errors,
And if we look only at the discretization errors, we can say some-- make some very fundamental remarks about the convergence of the finite elements scheme.
I've summarized here on the blackboard some of these very important concepts.
When we talk about convergence, let us assume first now that we are having a compatible element layout. In other words, we're using a compatible element layout. And then, we know that we have monotonic convergence to the solution of the problem-governing differential equations provided the elements that we are using contain:
Firstly, all required rigid body modes. This means, for example, that for in a plane stress analysis, the element must be able to undergo three rigid body modes, two translations, horizontally and vertically, and a rigid body rotation.
And the element must be able to represent the required constant strain states. In a plane stress analysis, three constant strain states. Pulling this way, pulling that way, and the constant shear.
If these two conditions are satisfied, we have monotonic convergence in a compatible element layout. What do we mean by a compatible element layout?
Well, here I have show a small picture showing two elements, a white and an orange element. And these elements are compatible because we have three nodes along this line for each of the elements. So the white element displacements basically vary parabolically along this line, and so do the orange element displacements. No gap can develop between these two lines.
Here we have an incompatible element layout. Two nodes that describe only a linear variation for the white element. But three nodes for the orange element describing a parabolic variation in displacement. Therefore a gap can open up and this is an incompatible element layout.
Well, let us assume then that we have a compatible element layout fist of all.
Then, typically for an analysis such as this one, a simple cantilever subjected to a tip load, the convergence that we would observe would be the following. As we increase the number of elements that idealize the cantilever, the displacement measured here would monotonically approach from below the exact displacement. What do we mean by exact displacement?
Notice I put it into quotes. Well, the exact displacement that I'm talking about here is the tip displacements that we would calculate from the problem-governing differential equation.
Of course, that problem-governing differential equation also corresponds to a mechanical idealization of the actual physical situation. And depending on which problem-governing differential equation we choose, we would have a different dashed line.
But assuming now that we have made up our mind which one to choose, which mechanical idealization to choose, and we are choosing a corresponding finite element idealization, this would be the convergence that we would observe. And it would be monotonically from below. In other words, our finite element model, actually is too stiff. It is stiffer than the actual physical structure.
Well, why I think we can intuitively understand that quite well if we just recognize that we are constraining the displacements within each element to only be able to basically, take patterns on that are contained in the element. In other words, for the four nodal element here, we really have for the U displacements and for the V displacements, only 4 displacement patterns each that the element can undergo. So we are constraining the exact displacements to be really, approximated by these displacement patterns. And that makes the actual finite element model stiffer than the actual physical structure.
If we have an incompatible layout, and I showed earlier what we mean by that, then, in addition, every patch of elements-- by patch of elements we mean small assemblage of elements-- must be able to represent the constant strain states. It's not just that each element must be able to represent the constant strain state, but each patch of elements.
Also, then if this holds true, then we have convergence. But non-monotonic convergence. By non-monotonic convergence, I mean the following.
We might get this result for a certain set of elements. And as we increase the number of elements, we are oscillating about the correct result, the exact result in quotes again, and finally we converge. This is non-monotonic convergence because we can have solutions lying on either side of the exact result.
Well, let me now show to you some of the facts that I just talked about on this view graph. What do these mean?
Well, I said already that each element in every case must be able to represent the rigid body modes.
Here for a plane stress analysis of a cantilever beam shown here, we have 3 rigid body modes. And as I mentioned earlier briefly already, the element must be able to go over, go up as a rigid body and rotate as a rigid body.
Why is this the case? Well, if we look at an analysis of the cantilever where we have not a tip load, but a load say, at this end here, we know that the bending moment looks this way. Hence, these elements over here must be able to freely rotate and displace. In other words, we know that the exact condition, the exact physical condition, there should be no stress here. The stress should be 0 in this element. And this means that the tip element here must be able to translate freely as a rigid body and rotate as a rigid body.
If this were not the case, then we would get stresses here. And of course, our analysis or finite element model can never predict the actual physical conditions accurately or closely, even if we take an infinite number of elements. So this somehow explains to you intuitively why the rigid body more criteria must be satisfied.
How about the constant strain state criteria? Well, if we assume that we take more and more elements-- and let me draw in here such elements. For example, now we have gone from a certain number of elements to a larger number of elements.
These elements, as we keep on refining the mesh, would become very small, and would go down. In fact, each element would go down to a point. But at a point we only have one constant stress basically. So the element must be able in the limit, to represent the constant stress at a point. And therefore, we have the constant stress criteria also.
One way of finding out whether an element satisfies these criteria, the constant stress criteria and the rigid body mode criteria is to calculate its eigenvalues.
Here, we have taken ones. A four node plane stress element and calculate the eigenvalues. Notice the lowest three eigenvalues are 0, representing the rigid body modes. And then, we have a flexural mode here, the fourth eigenvalue. The fifth eigenvalue is again the flexural mode. And here is the constant sharing mode. And here are the two constant, one compression and tension modes. So these are the constant strain states that the element can represent therefore. And earlier I showed you the three rigid body modes. So surely, this element therefore satisfies both of the criteria.
Let's look now at what happens when we don't have a compatible mesh. We have an additional condition that every patch of elements must also be able to represent the constant strain states. Only then we will have convergence. But we will not necessarily have monotonic convergence. We can only talk about non-monotonic convergence.
Well, here I've taken a patch of elements. Each element by itself satisfies the constant strain criterion and the rigid body more criterion.
In this case here, we use a compatible element layout. Compatible because along each side, the adjacent elements shares the nodes. And we are subjecting this patch of elements to a total stress here of 100 Newton per square centimeters. Of course, the stress in each of these elements would be simply a stress this direction of 1,000. Excuse me, 1,000 not 100.
The other stress, the vertical stress and the shear stress is 0 for this patch of elements.
Let us now perform the following experiment. Let us say that we are assigning two nodes, two different nodes to element 4 and 5 at this location.
Here they share the same node. Let us do the same at this end. And let us say that node 17 only belongs to this element 4. Node 18 belongs only to element 6. Node 20 and 19 belong to element 5. In other words, node 17 can basically gives this displacement pattern. And 19 gives this displacement pattern. An incompatibility can reside here. And similarly, along this line.
If we subject this patch of elements to the stress state of the external forces that we used here too, we would find that the stresses at the points A, B, C, D, E shown along here, which here were, of course, exactly equal to 1,000 for this stress tau yy and 0 in this direction. And the sheer stress was 0. If you look at the tau yy at A, B, C, D, E now the incompatible element mesh layout, we would find the following stress result.
Notice that the stresses now are not 1,000. In fact, they vary drastically from 1,000. 359 instead of 1,000. This means that this patch of elements with the incompatibility here, does not represent exactly the constant stress state that it should represent because we are subjecting it to a constant stress state.
Therefore, this patch of element does not satisfy the requirement for convergence. Not monotonic convergence and not non-monotonic convergence. It just should not be used in an actual practical analysis. This completes what I wanted to say then in this lecture. Thank you very much for your attention.
This OCW supplemental resource provides material from outside the official MIT curriculum.
MIT OpenCourseWare is a free & open publication of material from thousands of MIT courses, covering the entire MIT curriculum.
No enrollment or registration. Freely browse and use OCW materials at your own pace. There's no signup, and no start or end dates.
Knowledge is your reward. Use OCW to guide your own life-long learning, or to teach others. We don't offer credit or certification for using OCW.
Made for sharing. Download files for later. Send to friends and colleagues. Modify, remix, and reuse (just remember to cite OCW as the source.)
Learn more at Get Started with MIT OpenCourseWare