Topics: Formulation of the displacement-based finite element method
- General effective formulation of the displacement-based finite element method
- Principle of virtual displacements
- Discussion of various interpolation and element matrices
- Physical explanation of derivations and equations
- Direct stiffness method
- Static and dynamic conditions
- Imposition of boundary conditions
- Example analysis of a nonuniform bar, detailed discussion of element matrices
Instructor: Klaus-Jürgen Bathe
Study Guide (PDF)
Sections 4.1, 4.2.1, 4.2.2
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 number 3. In the previous two lectures, we discussed some basic concepts related to finite element analysis. In this lecture, I would like to present to you a general formulation of the displacement-based finite element method.
This is a very general formulation. We use it to analyze 1D, 2D, three-dimensional problems, plate and shell structures. And it provides the basis of almost all finite element analysis performed at present in practice. We will see that the formation is really a modern an application of the Ritz/Golerkin procedures that we discussed in the last lecture. We will consider in this lecture static and dynamic conditions, but as I pointed out earlier, we will only be concerned with linear analysis.
On this first view graph, I've prepared schematically a sketch of a three-dimensional body. This three-dimensional body could represent, typically, a bridge, a shaft, a building-- whatever structure we want to analyze. And this three-dimensional body here is subject to the following forces-- it is subjected to concentrated forces, forces that have components, Fxi, Fyi, Fzi, at one point i, and there are many such points i. The body's also subjected to body force components, Fbx, Fby, Fbz. These are forces per unit volume. And we will see later on that we include in these forces the d'Alembert forces, when we consider dynamic analysis.
The body is also subjected to distributed surface forces, with components, Fsx, Fsy, and Fsz. These surface forces would be, for example, distributed water pressure in a dam, frictional forces, et cetera. So we have, basically, concentrated surface forces, we have distributed surface forces, and volume forces-- externally applied volume forces, body forces. The body is, of course, also properly supported. We have here, typically, a support that prevents displacements in any direction. And here, we have another such support. Here, we have a roller support, which prevents displacements only in this direction.
The body is defined in the coordinate system, XYZ, and notice that I'm using here capital XYZ's. And the displacements of the body are measured as U, V, and W into the capital X, Y, and Z directions. I'm using capital letters here to denote global displacements and global coordinates. We will later on in the finite element discretization also introduce small lowercase x, y and z's, and u, v, and w's to measure the displacement in the individual elements.
So the problem is, other words, that we have this body, this general structure, subjected to certain forces, properly constrained, and we want to calculate the displacements of the body, the strains in the body, and the stresses, of course, in the body.
Well, on this view graph, here, I listed the external forces once as vectors, here's our force FB, the body force per unit volume, with components into the x, y, and z directions. Here, we have the surface forces with components in the x, y, and z directions. And here, we have a typical concentrated force at that point i, with components FX, FY, and FZ again. The displacements of the body measured in the global coordinates are U, V and W, as shown here. Of course, notice that these U, V and W's are functions of the capital XYZ coordinates.
The strains corresponding to these displacements, which of course, are known, are listed here. And the three-dimensional analysis, we have six such known strains, from epsilon XX to gamma ZX. Of course, the last three being the shearing strains, and the first three being the normal strains. The corresponding stresses are listed here. Again, six components from tau XX to tau ZX. Of course, if we were to actually analyze a two-dimensional problem, such as the plane stress problem, we would only use the appropriate quantities from here and from there, as we'll discuss later on.
The starting point of our analysis, in which we want to calculate the stresses, the strains, and of course, the displacement also. The starting point is the principle of virtual displacements. Now, this is the principle which we already discussed very briefly in the last lecture. Remember, please that we can derive it by looking at the total potential of the system, which is given as the strain energy, minus the potential of the external loads, W, U being the strain energy.
If we invoke the stationarity of pi, and we use the essential boundary conditions, which are the displacement boundary conditions, then we can derive the governing differential equations of equilibrium, and the force boundary conditions, the natural boundary conditions, as I have shown in the last lecture. Well, we will not derive these boundary conditions and the governing differential equations in this approach, but rather what we do is we invoke this principle, we set del pi equal to 0, and that gives us the principle of virtual displacements. And it is this principle which is a starting point all our finite element analysis.
Let's recall once what does it mean. Well, here we have the body forces that I applied to the body, the surface forces that I applied to the body. These are externally applied loads and concentrated forces that are also applied to the body at the points, i. These forces are in equilibrium with the stresses, tau. Let's assume that we know the stresses, at this point. Then the principle states the following-- if we subject the body to any arbitrary virtual displacements, listed in here-- and I'm saying any arbitrary virtual displacements-- excuse me, that however satisfy the essential boundary conditions, and that means just the displacement boundary conditions. Then the work done by the loads, and that total work is given here. This is a virtual work because we are taking virtual displacements and subject the forces to these virtual displacements. Then the external virtual work done is equal to the internal virtual work done, which is obtained by multiplying the real stresses, which are in equilibrium with these extraordinary applied forces. Multiplying the real stresses by the virtual strains, which correspond to the virtual displacements.
So let me use here a different color. These virtual strains correspond to these virtual displacements, and of course, these virtual displacements over the body-- these are, of course, a function of x, y and z. These virtual displacements over the body give us also virtual displacements on the surface of the body, which are listed in here. So let us put another arrow in there. And these virtual displacement also give us concentrated virtual displacements at those points where we have concentrated load supply.
So once again, if we take the body and subject that body, who is in equilibrium under Fb, Fs, and Fi, with tau-- tau being the real stresses. If you take that body and subject it to any arbitrary virtual displacement that satisfy the displacement boundary conditions, then the external virtual work is equal to the internal virtual work. The internal virtual work being obtained by taking the real stresses, times the virtual strains, which correspond to the virtual displacements here, and integrating that product over the volume of the body.
And the external pressure work is obtained by taking the real forces, multiplying these by the virtual displacements, and integrating these contributions over the complete body. Physically, what does this mean? Here, we have again, our general body. Let's see once, pictorially, what we're doing.
Well, let's take a certain virtual displacement, which I depict here. Now, here, we have a boundary condition, so this point, P can only move over to there. It could not move this way because we have to satisfy, in the virtual displacements, the actual displacement boundary conditions. Here, the point cannot move at all, and here, this point can also not move. So a typical set of virtual displacement might look like that. Just sketched in here. This roller out has moved over there. So there's our new roller right there.
What I satisfy are the order displacement conditions. Only horizontal movement was possible here. No movement here, no movement here. The virtual displacement vector here is that one. u bar, for a particular point. And that is the point that I'm looking at. Then what the principle says, once again, is that if I take these virtual displacements, multiply them by the real forces, integrate that product over the total body-- that is my external virtual work, and that external virtual work shall be equal to the internal virtual work, which is obtained by taking the real stresses, which are in equilibrium with these externally applied loads. And multiplying these real stresses by the virtual strains that correspond to these virtual displacements, and integrating that product, as the internal virtual work over the whole body.
This is an extremely powerful principle and an extremely important principle, and provides a basis of our finite element formulation. In our finite element analysis, we are proceeding in the following way-- we say, well, let us idealize this complete body as an assemblage of elements, and what I've done here is to draw one typical element. This is an 8-node element, a brick element, a distorted brick element, to make it a little bit more general.
It is an 8-node element because we have four nodes on the top surface, and four nodes at the bottom surface. There's another node here. This element here undergoes certain displacements, of course. And what I will be doing is I will express the displacements in that element as a function of the letter coordinates system, x, y, and z. The displacement in the element being lower u, v, and w. If we idealize the total body as an assemblage of such elements-- in other words, there's another element coming in from the top, and another element coming in from the sides, from the four sides, and another element coming in from the bottom.
So if we idealize the total body as an assemblage of such brick elements that lie next to each other, et cetera. And we express the displacement in each of these brick elements, as a function of the nodal point displacements, of the displacements of the corners of the bricks, then we can, of course, express the total displacement in the body as a function of the nodal point displacement. And that is the important step in the finite element analysis. That the displacements in each of these sub-domains or elements are expressed in terms of nodal point displacements. These corner nodes, as shown here, for the brick element. And then since the total body is made up of an assemblage of such brick elements, we can express the total displacement in the body as a functional of these nodal point displacements. And invokes the principle of virtual displacements. Now let's go into the actual specifics.
Well, for element m, this might be element 10, m in that case would be equal to 10. Then we have the following relationship-- and this is the important assumption of the finite element discretization. We say that the displacements-- there are three displacements, U, V, and W, of course, now. For element m, U, V and W are listed in this vector u, are equal to a displacement interpolation matrix, Hm, which is a function of x, y, and z, times the nodal point displacements. And what I'm listing in this u hat vector are all the nodal point displacements that I've called in the finite element discretization.
For this brick element here, we have eight nodes, and 24 nodal point displacements. At each node, we have U, V and W. But notice, once again, there are other brick elements on top of it, on the sides, and below of it. And each of these brick elements, of course, has a set of such nodal point displacements. We notice, however, that the element below it here-- if I take my pen here, and draw in another element, we notice that that element has the same node as the top element. In other words, this node here is common to this top element and the bottom element. And that's where we have the coupling between elements. We will see that more distinctly later.
So what I'm doing here is I express the displacements of element m as a function of all the nodal point displacements, and I'm listing here in u hat these displacements for N, capital N nodal points. We have this vector. Now, in general, later on, we would simply call all of these components, Ui's, and so our u hat here, would be written in this way. Notice I use the transpose, the capital T here, to denote the transpose of a vector. So this UN is equal to that W, capital N. That's just for ease of notation.
This is our major assumption. This is the major assumption in the finite element analysis. We will have to define for each element this displacement interpolation matrix. We notice that when we define it, there will be many columns that are simply 0's because only certain displacements in this vector listed here, in this vector, really affect the displacement in an element. In other words, typically, for this element here, if we look at this node, then the displacement at this node do not affect the displacement in this element because this node does not belong to the element.
That is recognized by the fact that in this Hm matrix there will be many columns that are simply 0's. In fact, the only non-zero columns in this Hm matrix are those that correspond to nodal points in this vector, or nodal point displacements in this vector that belong to element m. Well, having laid down this assumption, and we will define, once again, later on the Hm matrix for specific elements-- 1D, 2D, 3D elements and so on. Having laid down this assumption, we can derive the strains and the strains are simply given by this relationship, with the Bm matrix is the strain interpolation matrix, of element m. Notice that the rows in this matrix are obtained by using the rows in the Hm. And differentiating these rows and combining these rows in the appropriate way. I show examples later on. There is no more assumption in this step. This Bm matrix is simply obtained from the Hm matrix. By recognizing what strains we are talking about, and by recognizing that we can simply use the rows here, differentiate them, linearally combine them, if necessary, to obtain the Bm matrix.
Of course, we also have to use our stress strain law to obtain stresses from the strains, and the stresses in element m are given as shown here. These are the strains that I talked about here already. This is our stress strain law, which can vary from element to element. I also introduce here an initial stress, which might already be in the body. This might be due to overburden pressure in an underground structure, as an example, and so on. So, this is our stress strain law, which we have to satisfy for the body, of course. Our compatibility conditions in the analysis will also be satisfied. The strain compatibility conditions are satisfied because we are deriving the strains from continuous displacements, within the element. We'll impose on to the different elements that they remain compatible under deformations. By that, I mean if we have an element coming in here and another element going out there that the elements underloading, when the top element here, and the bottom element, both of them are loaded. No gap is opening up here so that displacement compatibility between the elements is satisfied.
So if we look at the three conditions that we have to satisfy in an analysis-- the first one being the stress strain law. That is satisfied because we are using this equation. The second one being compatibility. That is satisfied because we using this relationship here to calculate our strains from the Hm, via the Bm matrix. And we are satisfying, of course, that the elements remain together, so no gaps opening up. We will be talking about it later on when we talk about the convergence requirements also of finite element analysis. And finally, our equilibrium condition has to be satisfied. That is the third condition where that equilibrium condition is embodied in the principle of virtual work.
Here, I have it once again written down. The equilibrium conditions are in this principle of virtual work. And as I stated earlier, that if this equation is satisfied for any and all the arbitrary virtual displacements that satisfy the displacement boundary conditions, the real displacement boundary conditions, then tau is in equilibrium with Fb, Fs, and Fi. Well, what we will be doing is we will be applying this principle of virtual displacements for our finite element discretization, which means that in an integral sense, we satisfy equilibrium. However, if we look into an element, then within the element, we will only satisfy the differential equations of equilibrium in an approximate way. We will not satisfy them exactly.
However, if we have a proper finite element discretization, and by that, I mean if we satisfy all the convergence requirements that have to be satisfied in order to obtain a valid solution, or a reliable solution in a finite element analysis, then we know that as our elements become smaller and smaller and smaller, we will be finally-- and always, of course, applying the principle of virtual displacements-- we will be finally satisfying, also, the differential equations of equilibrium locally within each element. So stress strain law is satisfied, compatibility is satisfied, both of them exactly. The equilibrium requirements are only satisfied in an integral sense, if we have a coarse finite element mesh, but as the finite elements become more and more, as we refine our finite element mesh, we will be satisfying the equilibrium requirements. Also locally, within an element, always closer and closer, and we will be approximating, or we will be getting closer to the satisfaction of the differential equation of equilibrium.
The first step now is to rewrite this principle of virtual displacement, in this form, namely as a sum of integrations over the elements. There's no assumption yet. All I've done is since our total body is idealized as a sum of volumes, namely the volumes over the elements, I can rewrite the total integral as a the sum over the element integrals. And that's what I have done here. Notice we have now here, m, denoting element m and we are summing over all of the elements. There's no assumption here yet. Now, however, I can substitute our assumption. Namely, that Um is given in this way, and epsilon m is given that way. Once again, this is actually not an assumption.
This is the major assumption. That epsilon m follows from this assumption. This equation follows from that equation entirely. Substituting now from here, these equations into the principle of virtual displacements, we directly obtain the following equations. Here, I have on the left hand side-- let's go through this equation in detail. On the left hand side, I have the following part. I have an integral over all of the elements, that is the integral here, summing over element m and integrating over each of the element. This part here, Bm transposed, times U hat T, is equal to the epsilon bar, mT. So notice that our epsilon bar m, is a virtual strain, is given in this way. Before, I talked only about the real strains. In other words, the bar was not there. I put that now into brackets here. This bar was not there.
Well, what we're doing in the finite element analysis is to use the same assumption for the virtual strains as we use for real strains. In other words, we use this equation without the bar and with the bar, on top of the epsilon and on top of the u hat. So here, we have our epsilon bar transposed. This part here is the stress. Tau m equals Cm epsilon m-- that is our cm. The Bm here comes in from the epsilon m. So Bm times U hat is equal to epsilon m, once written down again here. So this total part here is nothing else than-- let's go back once more to the previous view graph-- nothing else than this part here, than that part there. But with the initial stress, tau i, m being on the right hand side. Since we do know the initial stress, we put that one, of course, on the right hand side. It is a load contribution. And here, you see it as a minus sign because we put on the right hand side, and this part here. This part times this u hat, notice there's a big bracket here. That u hat bar here operating on that Bm transposed, there's a transposed here, too, gives us again the epsilon bar m transposed.
So let us look now at what we have on the right hand side. On the right hand side, I want to have discretized this part here. Well, what we do is we substitute here from our displacement interpolation, here from our displacement interpolation, and each of these integrals can directly be expressed, in terms of the nodal point displacements. And that's what we have done here. The Hm times the u hat bar is the u bar m transposed. Here, we have the u bar s. That is the u hat bar times the Hsm. Notice that this part, you should not have been written that far. In other words, there is the end. The u bar s m transposed only goes up to there and it embodies this Hs m transposed and the u hat bar transposed.
This part we talked about already. So what we have done then is to rewrite-- this is the important part-- is to rewrite this principle of virtual displacements, in which we had no assumption yet. We rewrote this in terms of the nodal point displacements and element interpolation matrices that we use for our finite element discretization. Now, we of course, have the assumption that the displacements within each element are given by the Hm matrix, the strains are given by the Bm matrix. So this is the result that we have obtained. And at this point, we now invoke the principle of virtual displacements. Since this principle here shall hold for any arbitrary virtual displacements that satisfy the displacement boundary conditions, we can now invoke this principle n times. And by that, I mean, once by imposing a unit displacement at the first displacement degree of freedom, and leaving all the others 0. Second time around, imposing a unit displacement at the second degree of freedom, all the others being 0. Third time around, imposing a unit displacement at the third degree of freedom, all the other displacements being 0, and so on.
And that really amounts to then saying that this vector here becomes an identity matrix. And similarly, this one here becomes an identity matrix. And therefore, we can take those two out, and our resulting equation then is simply what we have left. Taking these two identity matrices out, and that is our finite element equilibrium equation, or I should rather say that these are n finite element equilibrium equations, namely corresponding to the n nodal point displacements that we are considering. In general, what one does most effectively is to really derive these corresponding to all displacements. Having removed the boundary conditions, and then later on, one imposes the known bounded displacements. And that's what I want to discuss a little bit later in more detail.
So this is the result then that we obtained. And we obtained really in shorthand, Ku equals r. Where K is this matrix. It's the structural stiffness matrix. Notice that I'm summing here over the elements. This being here, the element stiffness matrix. Notice that this element stiffness matrix here is an n by n matrix, has the same order as this K matrix. However, we will recognize that a large number of columns and rows in this matrix are simply 0. In fact, all those columns and rows are filled with 0's that do not correspond to a nodal point displacement degree of freedom of element m. I show you later on some examples.
However, by using this Bm here, and making the element stiffness matrix of the same order as this total structural stiffness matrix, we can directly sum over all of the elements stiffness matrices. And that, of course, is our direct stiffness procedure, which already I pointed out to you earlier. The direct stiffness procedure means that we are adding the element stiffness matrices into the total stiffness matrix via this summation here. So the Km here, being what I have here in the blue brackets, must be, of course, of the same order as K in order to be able to do that method theoretically, at least. Later on, we will see that we indeed only work with the non-zero rows and columns in the K matrix, and then use connectivity arrays to assemble Km effectively into the actual K matrix.
For the RB vector, we have this part. Once again, we're using the direct stiffness procedure to add the contributions of all the elements in order to obtain the total RB vector. Once again, the rows now, or rather, the elements because this of course, is a vector here of n long now. Those elements that do not correspond to nodal point degrees of freedom will all be 0's. The RS vector, similarly, is obtained as shown here. Now, we of course, sum the element contributions as they arise from the surface forces and the RI vector is obtained as shown here, and the concentrated load vector is simply a vector listing all the concentrated forces in F. Notice that this HSM matrix here is directly obtained from this Hm matrix.
Sometimes one has difficulties visualizing what this matrix really is. Well, we will see later on if this is an element here, and our coordinate system, say, lies in that element this way, then the Hm matrix, of course, gives us a displacement within the element. Whereas the HSM matrix gives us, say, the displacement on this surface element, if it is that surface that we want to consider.
Now, to get the displacement on the surface of the element when we know the displacement within the total volume of the element, well, what we simply have to do is we have to substitute the coordinates of the surface in the Hm here to obtain the HSM. I will show you later on some examples.
Now, in dynamic analysis, of course, the loads are time dependent and if we are considering a truly dynamic analysis, then we have to include inertia forces. And the inertia forces can directly be taken care of, or can directly be included in analysis if we use the d'Alembert principle. Here, we have the body loads, which are the externally applied forces per unit volume. And if we split these up into those forces that are externally applied, and those that are arising due to the d'Alembert forces, as shown here. Then we directly have the inertia effect in the analysis. We now can, of course, express our accelerations in the element in terms of nodal point accelerations again, and we are using here the same Hm matrix that we use already for the displacement interpolations. If we substitute from here and here into the RB which I had written down here. If we substitute into this RB here, this equation for fB, then we directly can write down this equation here, M U double dot plus Ku equals R, where the M matrix now is obtained as shown here. Notice that this R vector now only contains this RB part. In other words, not anymore the fB here, but rather an fB curl.
Notice also that in this analysis now, or in this view graph, I've dropped the hat on the u. These are the nodal point displacements, these are the nodal point accelerations. I've dropped the hat just for convenience. I had already dropped it actually here also. We earlier had the hat there. Here, we had the hat still because I wanted to distinguish the actual nodal point displacements from the continuous displacements in the structure, or in the body. So here, the hat still being there. And here, I dropped the hat already, just for convenience of writing, and from now on, when we have this vector, U here, then that means that we're talking about the concentrated nodal point displacements, or the actual note point displacements of the finite element mesh.
I mentioned earlier that it is most convenient to include in the formulation all of the nodal point displacements, including those that actually might be 0. In other words, for our three-dimensional body, to make a quick sketch here. If we have here our support in an actual analysis, it is most effective to say well, let us remove the support and assign a node there with three unknown displacements. Once we have derived these equations of equilibrium, of course, we now will have to impose the fact that the displacements are 0's there.
And that is then done effectively, for example, as shown here. We have here the general equations, M U double dot plus Ku equals R. And what we are doing is we are listing the displacements and accelerations into vectors, U double dot A, U A, and U double dot b, and Ub, where the b components of the displacements and accelerations are known. And now they might be 0, as I showed here in this particular example, or they might be actual values that we want to impose. If they are known, well, we can look at the first equation, as shown here, and put all the known quantities on the right hand side, substitute for Ub and U double dot b, and we know then the right hand side load vector. Thus, we can calculate Ua, and U double dot a. Having now calculated the velocities, the accelerations, and displacements, we can go back and get the reactions. The reactions, of course, being Rb.
Here, I assumed that the displacements which we are talking about in the vector here are actually the ones that we also might want to impose. Well, in some cases, of course, we might have defined in our finite element formulations the U and V displacement, as shown here. But a displacement that we want to impose is actually this one here, namely, that one might have to be restrained, and this one here might have to be free. In that case, if our finite element formulation has used the U and V displacement, we have to make a transformation as shown here. A well known transformation from the U to the U bar displacements, and this, in a more general sense, is written down here once again. When we have many more degrees of freedom, our T matrix would look as shown here. It's an identity matrix with the little transformation matrix that I've shown here. The cosine minus sine sine cosine matrix. Now, put into the appropriate rows and column. The i'th column, j'th column, i'th row, and j'th row would carry these 2x2 matrix. And otherwise, we just have 1's on the diagonal. So this is the more general transformation that we are using when we have many more degrees of freedom than just the two that we want to modify.
Substituting from here into our equations of equilibrium M U double dot plus Ku equals R. We directly obtained this equation, where M bar now is shown here, K bar is shown here, and R bar is shown here. Let me mention here that this looks like a former matrix multiplication. In fact, two former matrix multiplications. M times T, And then the product should be taken times T transposed, pre-multiplied by T transpose. Well, in actuality, of course, all we need to do is combine rows and columns. The i'th and j'th rows and columns to obtain directly our M bar matrix, similarly for the K bar and the R bar matrices.
Another procedure that is also used in practice-- can be very effective-- is an application of the penalty method. In this procedure, we impose, basically, physically a spring of very large stiffness, where K is much larger than K bar i i. And then we supplement our basic equations that are shown here by this equation here. So if this K is much larger than K bar i i, and if we supplement this equation or add this equation into this equation here, then we notice that the spring stiffness will wipe out basically the other stiffnesses that come into this degree of freedom, and our solution will simply be that U i is equal to b, which is the one that we want. So this penalty method can be used to impose displacement degrees of freedom, and it really physically amounts to adding a spring into the degree of freedom where we want to impose a certain displacement. And it's important, however, that if we use that method that we are always dealing with single degree of freedoms being imposed.
And by that, I mean the following-- if we had a system, like this one here, and our original displacement degrees of freedom are these, then we first have to make the transformation onto that finite element element system to these degrees of freedom here. And now we add our spring in. In other words, we want to add our spring into this system of equations because now there's no coupling from this degree of freedom that we want to be impose into other degrees of freedom, through that spring. This spring only enters on the diagonal, and now it is a numerically stable process. However, if we were not to perform this transformation-- in other words, if we were still to deal with these degrees of freedom, and then add our spring in-- of course, that spring now would introduce coupling between these two degrees of freedom, and numerical difficulties may arise in the solution.
So basically, then in summary, if we do have degrees of freedom to be imposed, we first go through this transformation to obtain the M bar, U double dot bar, K bar, U bar equals R bar system of equations where the displacement that we're talking about are containing those displacements that we actually want to impose. We then can impose these displacements using the penalty method, which is this one. Or we can impose these displacements using the more conventional procedure, using, in other words, simply this procedural of imposing Ub and rewriting the equations into two equations. The first equation we solved for the Ua now, of course, in this particular case, we would now have all bars on there. And if we have done a transformation, and in the second equation then, we obtain the reactions.
Let me now go through a simple example to show you the application of what I have discussed. This is a very simple example, but a very illustrative example. In particular, it is also the example that we talked about already earlier in lecture 2, when we did a Ritz Analysis on this problem. In fact, our finite element analysis that we are now pursuing, using the general equation that I presented to you is really nothing else than a Ritz Analysis. And in fact, if you look at the earlier solutions that we obtained-- solution 2, in the Ritz analysis, corresponds to the finite element solution that I will be discussing with you now.
So here is a problem once again. We have a bar of unit area, from here to there, and then of changing area, from here to there. The length here is 100, the length here is 80. The bar in actuality, is supported here, but as I mentioned earlier, we remove that support in our finite element formulation, and introduce, in fact, a displacement degree of freedom there.
So here, I want to put down the first node. Now, the first step in any finite element analysis must, of course, be the step of idealizing the total structure as an assemblage of elements. And there are generally choices-- how many elements to take, what type of element to take, and so on. In this particular case, I know that there's a discontinuity in area here and for that reason, intuitively, I will put one element from here to there with a constant area. Also, let us consider for the moment, this also as being one element, and this then will correspond to the Ritz Analysis that we performed earlier.
Notice that we have here a bar of unit area, a bar of changing area. This total bar assemblage is subjected to a load of 100, a concentrated load of 100, as shown here. The only strains that this bar can develop are normal strains. In other words, if a section originally is here, that section we move over a certain amount and by that amount. That is U. In the coordinate system that we are using, the y-coordinate being in this direction, this is the y-coordinate here. This would be our U of Y. However, since we are dealing with two elements to analyze this bar, what I will do is I will introduce a little coordinate system here.
And I used little y in this particular case. So we put a little y here. There's also, for this element, a little y. And the area, in this particular element, is given as 1 plus y divided by 40 squared. So this is the changing area in this domain. However, also, remember please, if in our analysis, if an original section in this area was vertical like that, after deformations, due to the load here, it will still be vertical, and now our displacements in this element will also be given by a u. And that u is a function of-- if we look at this little y, if you use that little y of that y as this u here is a function of this y. This y corresponds to the y in this element, that y corresponds to the y in this element because we use different coordinate systems for each element.
Now, this is actually an important point that we can use for each element, a different coordinate system. We could use Cartesian coordinate systems for each element, different ones. In fact, that is most generally done. If we have specific geometries, we might use cylindrical coordinates systems for certain elements, Cartesian coordinate systems for other elements, and so on. This is an extremely important point that we can have different coordinate systems for different elements because that eases the calculation of the element stiffness matrices.
So here, our element 1, here our element 2. And the displacements that we are talking about are U 1 at this node, U 2 at this node, U 3 at that node. These three displacements shall give us the displacement distributions. Of course, only an approximate displacement distribution in the complete element mesh. Two elements make up our element-- complete element idealization or complete element mesh. Well, the equation that, of course I will now be operating on is this one, KU equals R. In this particular case, we recognize that we want to calculate our stiffness matrix, K. Here, we have two elements, so M, in this particular case, will be equal to 1 and 2. We don't have a body force vector, we don't have a surface force vector, we don't have an initial stress vector. However, we have a concentrated load vector. So what I will want to do then is calculate our K matrix, and establish our concentrated load vector.
The K matrix embodies the strain displacement interpolations, which are obtained from the element displacement interpolations. Now, as I already pointed out, we have now here, two elements. The first element shown here, second element shown here. Notice that the U1, U2 correspond to these two displacements, U1, U2. The U2, U3 correspond to these displacements, U2, U3.
Notice that there's a coupling between the elements because U2 is here a displacement of that element, and is here the displacement of element 2. The length of the element 1 is 100, length of element 2 is 80. Well, if we have two displacements to describe the displacement in an element, then we recognize immediately that all that we can have is a linear variation in displacement between the two end points, between these two nodes. And so the element displacement interpolations must involve these functions. For unit displacement at this end of an element, Y over L is the interpolation of the displacement. For unit displacement at this end of the element, this is the interpolation. Notice that the actual displacement, of course, goes into this direction, but I'm plotting it upwards to show you the magnitude of that displacement.
If we have established these interpolation functions, recognizing, of course, that for element 1, L is 100, for element 2, L is equal to 80. Then we can directly write down our H1 and H2 matrices. They are given as shown here. Notice that Um is equal to Hmu, where U is equal to-- I write down the transpose-- U1, U2, U3. Now, we notice that element 1-- let's go back once more for element 1. Only U1 and U2 influence the displacements in that element. And that is shown here. U3 has 0 and does not influence the displacement in the element.
Similarly for H2, U1 does not influence the displacement in that element. So we have these displacements. Notice also that I've written here, Um, of course, but that Um here, for our specific case is simply this displacement, Vm. We have only one displacement component. Taking the derivative of these relations here, we get directly the strains. Notice that here, we should have probably put an m there. This is the normal strain and we obtain these matrices by simply taking the derivatives.
Well, now we have the components that we need to evaluate the K matrix. And as I mentioned earlier, that is the total K matrix obtained by summing the contributions over the elements. This is coming from element 1, this is coming from element 2. Notice that the area is 1, Young's modulus of stress strains law. We are integrating from 0 to 100. This is here, B1 transpose, that is B1. This is here, B2 transpose, that is B2. This is the area that I pointed out to you earlier.
Evaluating these two matrices, we directly obtain these matrices. And notice the following-- that there's no coupling from the third degree of freedom into element 1. Similarly, there's no coupling from the first degree of freedom into element 2. In fact, what we will do later on is simply calculate the non-zero parts. We call these the compacted element stiffness matrices. And knowing these non-zero parts and knowing into which degrees of freedom they have to put in the assemblage phase to obtain the total stiffness matrix, we can directly assemble the stiffness matrix. In other words, if I know this part here and I know that the first column corresponds to the first column of the global stiffness matrix, the second column corresponds to the second column in the global stiffness matrix, then I can just add this contribution into this part here.
Similarly, I can simply add this contribution here into that part there, without carrying always these 0's along. And that is, of course, a very important computational aspect. However, in theory, we are really still performing these additions as shown here, we really still perform the additions as we pointed them out in the direct stiffness procedure. In other words, we still perform this summation, as I pointed out to you earlier. The important point, however, is that we now have established the K matrix, corresponding to the system.
Our R vector is simply, in this particular case, 0, 0, 100. Because we only have 100 applied at the third degree of freedom. We now have to impose that U1 is 0, we simply set U1 equal to 0 in the equilibrium equations, as I pointed out earlier. We solve for U2 and U3, and having obtained U2 and U3, we know the displacement in each of these parts, and we know, therefore, the strains and the stresses in each of these parts. The solution is plotted in the example that I discussed with you in lecture 2.
This example, really, showed some off the basic points of finite element analysis. Of course, we have to discuss much more how we actually obtain the Hm matrices for more complex, more complicated elements that I use in actual practical analysis. However, this is all I wanted to mention in this lecture. Thank you 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