Home » Supplemental Resources » Finite Element Procedures for Solids and Structures » Linear Analysis » Lecture 6
Flash and JavaScript are required for this feature.
Download the video from iTunes U or the Internet Archive.
Topics: Formulation and calculation of isoparametric models
Instructor: Klaus-Jürgen Bathe
Lecture 6: Formulation and ...
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 dot MIT dot edu.
PROFESSOR: Ladies and gentlemen, welcome to lecture number six. In this lecture I would like to discuss with you the formulation and calculation of isoparametric finite elements. We considered earlier already, in lecture four, the calculation all finite element matrices. But in that lecture we considered the generalized coordinate finite element models. The generalized coordinate finite element models were the first finite elements derived.
However, I now want to discuss with you a more general approach of deriving the required interpolation matrices and element matrices. And this more general approach is the isoparametric finite element derivation. The isoparametric finite elements that I will be discussing in this and the next lecture are, in my opinion, the most effective elements currently available for plane stress, plane strain, axisymmetric analysis, three dimensional analysis, thick and thin shell analysis. These elements are being used, for example, in the computer program [? Aldena ?]. They are also used in other computer programs and represent a modern approach to the solution off structure problems.
In this lecture, I would like to talk about the derivation of continuum elements. In the next lecture we will talk about the derivation of structural elements, [INAUDIBLE] and central elements, beam elements. The basic concept of isoparametric finite element analysis, is that we interpolate the geometry of an element, and the displacements of an element in exactly the same way.
Let us look at the geometry interpolation. Here you see for three dimensional analysis, the interpolation of the x-coordinate within the element, where we use interpolation functions h i. Of course these interpolation functions are unknown and I will have to show you how we derive them for certain elements. The x i as a nodal point coordinate. The x i value for i equals 1, for example, is nothing else than the x-coordinate of nodal point one. And so on.
So the x i, the nodal point coordinates for all nodes, they are given, they are input in the analysis the way we have been discussing it in the previous lecture. And if we know the h i, we have a direct relationship for the x-coordinates within the element, as a function of the nodal point coordinates. Again, I have to show you how we obtain the h i functions for the various elements that we are using. Similarly, for the y interpolation and similarly for the z interpolation.
Having derived the h i functions we are using in the isoparametric finite element, the same functions also to interpolate the displacements. Notice that we're having here, the nodal point displacements u i, v i, and w i, and here we have the same interpolation functions that we already used in the coordinate interpolations. And if the number of nodes that are used in the description of the element-- and, in fact, we will see that N can be a variable, it can be equal to three, four, five, up to a large number of nodes. And in practice we generally don't go much further than 60.
I mentioned that I want to discuss in this lecture continuum elements. Well, the continuum elements that we addressing ourselves to are the truss element, the two d elements-- plane stress, plane strain, and axisymmetric elements-- and then, the 3 d elements for three dimensional and thick shell analysis. These are continuum elements. We call them continuum elements, because we only use u, here u and v, and here, u v and w-- the displacements-- to describe the internal element displacements. We're only using the nodal point displacements u, u v, and u v w, for all of the nodal points to describe the internal element displacements.
Structure elements, on the other hand, also use the rotations at the nodal point theta x, theta y, and theta z, in order to describe the deformations within the element. And I will be discussing structural elements in the next lecture, so, for the moment, we do not have any rotations at the nodes. Or what we will allow, u v and w displacements at the nodes. Typical examples are given on this Viewgraph. These are elements that are available in the computer program [? Aldena ?].
Here we have a truss element, a two-noded truss element. The only displacement of interest in the truss is the actual displacement here. Here we have a three-noded truss, or cable element. Here we have a set of two dimensional elements. Notice we have triangular elements here, a triangular element here. We have here an eight node element that is curved. We can construct curved elements in the isoparametric approach. And here we have a rectangular, eight node element, straight sides, in other words.
These elements are used in plane stress, plane strain, and axisymmetric analysis. In three dimensional analysis we might be talking about-- we might be using-- such elements such as the eight node brick. Each node now has three displacements, u, v, and w. This is here a higher order element, where we have, in addition to the corner nodes, we also have mid-side nodes. You don't need to have a mid-side nodes along all sides. For example here, I did not put a mid-side node as an example. In fact we could simply have mid-side nodes here, and no mid-side nodes anywhere else. I will show you how we deconstruct these interpolation matrices in the next few minutes.
Let us consider, as a very simple example, a special case. And the special case that I would like to look at is a truss element, which is two units long. In other words, the length from here to there, is equal to 1, and similarly the length from here to there, is also equal to one. I'm describing the element with a coordinate r that I set equal to zero at the midpoint of the element, plus one at the right end and minus one at the left end. This is special geometry and we will later on have to generalize our approach to more general geometries where the length is not equal to two and the element might even be curved. But for instructive purposes, let's look at this geometry first.
Similarly, for two dimensional analysis, I want to look at an element that has length two this direction and length two that direction. In other words, a generalization of the concept that we just looked at for the truss. Now we have two dimensions and in each dimension we have a length of two. I embed into this element an rs-coordinate system. And this coordinate system now has its origin here at r equals zero and s equals zero. The r-axis bisects this side, and the s-axis bisects this side. Notice that this is a coordinate system that is embedded into the element. We also call that the natural coordinate system, or the Isoparametric coordinate system. This element, two by two, would lie in space in an x xy-coordinate system as indicated here.
Now for three dimensional analysis we would proceed in the same way, then we would consider an element that is two by two by two units long into each coordinate direction. And then we would have an rs-axis and a t-axis coming out of the transparency in this particular case. Well, I want to look at the truss element, the special truss element, and the special two d element. And I want to show you how we can construct the interpolation matrices, the displacement interpolation matrices, h i. These are also the coordinate interpolation matrices. I want to show you how we construct them for the special elements, how we then can calculate the strain displacement interpolation matrices for the special elements. And then I want to go on and show you how we generalize to concepts to curved elements. And once we have discussed the two dimensional case, I think you can see yourself how the concepts are generalized to the three dimensional case.
Let's look then at our two-noded truss first. Here we have once more, the truss, and we have node one on the right hand side, node two on the left hand side. Our r-coordinate system starts in the middle of the truss. Notice what we want to obtain is that we interpolate our u displacement via the interpolations. U being equal to h i times u i, the h i are the unknowns. I want to show you how we obtain the h i. The u i are the nodal point displacements. In this particular case it would be here, u two and here we would have u one. I in other words, goes from one to two for this particular case. The h i has to be a function of r. For a given r, however, we can evaluate the u displacement if we have u i is given. In other words, when u one and u two are given, then for a given r, we can evaluate directly from this relation the displacement at that point r.
Well, with two nodes, we can only have a bi-- only a linear representation in the displacements. The h one, from this relation, must be this function. Because if we look at this, and we say let u two be equal to zero, we would simply have that u is equal to h one u one. Well, therefore, our h one must look this one, because u one has its full strength here. If we put u one equal to one, we would have u is equal to h one, and this would be the variation. Similarly, for h two, in this case, we would have u is equal to h two u two. H one is now equal to zero. If we put you two equal to one, as I have done in this particular case, then our linear variation is indicated as shown here , and the function h two is given right here. Notice that h two is equal to zero when r is equal to one. And h two is equal to one when r is equal to minus 1. So, this gives us h two. The two-noded truss, therefore, simply has this description here for the displacement and the coordinates. Remember that we're using the same h i for displacement and coordinates. It simply has this description where h one and h two are defined as shown.
Let us now say that we want to add another node, that we want to put another node right there. In other words, we want to go from a two-noded description to a three-noded description. On one of the earlier Viewgraphs, you could see a three-noded cable element. Well, the way we proceed in the construction of the h i functions is as follows.
The two-noded element simply had this description here. H one is given as shown here on the left side of the blue line. H two was simply this. And we knew that we could do no better than a linear description in displacements between two nodes. However, now we have a third node right here. And a third node means that we can use a parabolic description in displacements. Now h three must be equal to one at the third node and zero at both sides. Why? Well, because, remember we have u now equal to h i u I, where i equals one to three. And if we put u one equal to zero and u to equal to zero, we simply have u equals h three u three, and, therefore, h three is this function. I've written it down here when r is equal to zero its one when r is equal to minus 1 or plus h three is zero. So this is our h three.
However, if we now look back to our earlier description of h one and h two for the two-noded element, we remember that we had a linear variation here. That we had a linear variation here and here. This one was a linear variation of h one-- let me take here the green color to show once more what we are talking about-- that was our linear variation. And here, this one here was our linear variation there.
Well, now with the third node there, we recognize that, at this third node our h two-- our actual h two for the three-noded truss-- must be zero here. H one must be zero here. Well, how can we make it zero here? We can make it zero by subtracting from our two-noded truss h one, one half of this parabolic description. And that's what I have done here. Similarly, we have to subtract it here, because then, what we are doing is we are taking one half off this parabola, and we are putting it right on here. This is one half of the bottom parabola. We are putting it on there. And that brings this point back to that point. Remember this is equal to one. I want one half as a correction here, and that's why I take one half of one minus r squared, to bring this point back to there.
This total description then, this total part-- this part here, all of that together-- is our h one for the three-noded element, for the three-noded element. And similarly, we would have four h two. This total part here is this function here for our three-noded element. Now the important point that I really would like you to understand is that, we have started off with a two-noded element description-- h one only the linear part, h two only the linear part. We have constructed the interpolation function for the third node, and then we have corrected the earlier two-noded interpolation functions via subtracting a certain part of the third interpolation function, in order to obtain the new h one and h two for the three-noded element. And this is indeed the actual procedure that we can use very effectively in constructing higher order elements. We are starting off with the lower order element descriptions-- the ones shown here dashed with a dashed line, the linear part only-- and we add in the higher order description, and subtract the correction from the lower order description. The correction has to be subtracted to bring this point back to zero, because h one total must now be equal to zero here, h two total must be equal to zero here, and this way we have constructed our new h one and h two for the three-noded truss element.
Let's look at the four-noded element in two dimensions now. Here we use very similar concepts. These are the four nodes for the element. The description along a side is just like we have been discussing now for the truss, linear here and linear here. Notice h one can directly be written down. It has to be equal to one here and zero everywhere else. This is a function which is bilinear and which satisfies these conditions to be equal to one here and. Zero at the other nodes. It is a linear variation along this side, along this side, and also across the surface. In other words, for a given value of r, we have a linear variation across here too. How do we obtain h two? Well, we simply have to change the signs over here-- the plus signs here-- in an appropriate way. H two shall be equal to one here zero at all the other nodes, while we see that this function satisfies these conditions, let's put in r equal to minus 1. That makes this two, divided by four, gives us one half. S has to be equal to plus 1 at this point-- you get another two in here-- so h two is equal to one right there, and its zero everywhere else. Similarly we construct h three and h four. Indeed we can immediately observe that these signs-- plus in both cases here, a minus here, plus there, minus, minus, plus, minus-- these signs correspond to nothing else in the signs of the r- and s-coordinates of the nodal point under consideration. R and s is plus here, you have two plus signs here. R is negative here, you put a negative sign here. But s is positive here, you put a plus sign there. r and s both are negative here, so we have two negative signs here, et cetera.
Let us now see how we construct from this basic four node element, which really corresponds to our basic two node element in the case of the truss, how we can construct higher order elements. Well, we proceed in much the same way as in the truss formulation. Here we have a four-noded basic element and we have added a fifth node to it. Now let's look at this in detail. Since we now have added a fifth node here, we know that we can allow a parabolic distribution. In fact, we should allow a parabolic distribution of displacements along the side. Well, if its a parabola along the side for s being plus 1, we immediately see that this function here-- this s equal to plus one, this is equal to two, knocks out that one half, we have a one minus r squared along the side-- just the same function that we already had in the truss.
Well, along this direction we can only vary things linearly because we have two nodes only in along these directions. And this is the reason why we have to put in a one plus s here. The linear variation is given by one plus s. And we also notice that when s is equal to minus 1, in other words, we are looking at this side, this function is zero. When s is equal to plus one, as we pointed out earlier, this function here-- this part and that part-- gives us a one together, and we have simply a one minus r squared along here. So, here you can see these triangles that show the linear variation along this side. These parabolas here run really across here, but with a different intensity. The intensity off the parabola goes down linearly from one at this end, to is zero at that end.
This is our h five. Well, if we have this as our h five, we remember that our original h one of the of the four-noded element had a linear variation along here, and also a linear variation along here. We will now have to take this h one and subtract some correction from it. In order to make this point here have zero displacement, four h one. Well, what we will do is we take this h five and subtract a multiple of h five from h one. In fact, you can see since the original h one is equal to one half here, we simply have to subtract one half of this function to obtain the new h one. And this is what I have done on this Viewgraph.
The h one now is the original h one that we had. And we are subtracting one half of h five, which is one half of the parabolic distribution, to bring this point here down to zero in displacements. And that's what we have done here. The resulting function then is shown here. And our h one for the five-noded element is shown right here. Similarly, for h two-- this is the h two function for the five-noded element-- interesting to note that h three and h four are for the five-noded element, the same as for the four-noded element. Because this fifth node lies between nodes one and two. And there is no effect along this side and along this side here-- along that side four h three and h four-- so, we have our original functions also for the five-noded element.
Let us look now at a generalization of this concept. Here we have a typical nine-noded element. A very effective element for many types of applications. I already show it here in its curved form, but think of it, please, as follows. This is the x-axis. This is the r-axis. S is equal to zero along this side. S is equal to plus one along this side. S is equal to minus 1 along this side. R is equal to plus one along this side. R is equal to minus one along the side. And r is zero along this axis. So, in the r s description, in the embedded coordinate system, this element is still a two by two squared element.
Well, if we look then at the interpolation functions, for this element, they look as follows. Now maybe you have difficulty seeing all this information, so, please then refer to the study guide where you'll find this Viewgraph. For the four-noded element, we had these interpolation functions. If we want to deal with the five-noded element, what we have to do is, we add this interpolation function. And-- as I pointed out earlier-- we have to correct our h one and h two. But that is all of the correction that is required. H three and h four not corrected-- they are blank spots here, they are blank spots here. So, our five-noded element would have these interpolation functions now shown in red.
For if we added a sixth node-- and I now should go back to our earlier picture-- if I wanted to add, in addition to the fifth node, also the sixth node. Well, then, I have to put another interpolation function down here, which now is parabolic in s and linear in r. And I have to correct h two and h three again, these are the corrections. So, now I have in green here shown to you the interpolation functions off a six-noded element. And like that we can proceed by adding interpolation functions and correcting the earlier constructed interpolation functions as shown. Like that we can proceed to directly obtain the interpolation functions for the five-noded, six-noded, seven-, eight-, and nine-noded element.
In fact, its also important to notice that we could have this node, that node, this node, and that node, and just h nine also added. We could have, in other words, a five-noded element which looks like this. It has this node, that one, that one, and that one, and that one in the middle. Another five-noded element would be this one, this one, this one, this one with that node in there. So there is no necessity in having all the nodes below a certain number. But we can simply use four nodes, and then add whichever nodes we want to have into the element.
This is, then, how we construct interpolation functions. And once we had the h i, we directly can obtain the h matrix, the matrix that gives the displacements in terms of the nodal point displacements. Notice that the h matrix really is constructed from the hi's. And the elements of the b matrix-- the strain displacement interpolation matrix-- are the derivatives of the h i or zero. And I will show you an example right now.
Because we are using four, we are still looking at the special case of a two by two by two element in a truss case, only this two. Plane stress, plane strain, axisymmetry, two dimensional analysis, we have these two two's. In other words, we're talking about a two by two element. And in three dimensional analysis we would talk about a two by two by two element. In these cases, we have x equal to r, y equal to s, z equal to t. So the strains-- which are derivatives with respect to x the actual physical coordinates-- can also directly be obtained by simply taking the derivative with respect to r, and then similarly for s and t.
Let us look at a four node, two dimensional element. This is really the element that we have used in our earlier example of the cantilever analysis. Here we would simply have that u r s, v r s, are described as shown. Notice that in the first row we are really saying nothing else than u is a summation h i u i. Where i equals one to four because we have a four-noded element. H i u i. In the second row we are really saying nothing else than v being the summation of i equals one to four h i v i. And all I have done is I have taken these hi's and assembled them into a matrix form to obtain our h matrix.
This h matrix here-- the entries in that h matrix-- are dependent on the ordering that you're using here for u one, v one, u two, et cetera. With this ordering, these are the entries. Notice there are zeroes here, because the v degrees of freedom at the nodes have no contribution to the u displacement in the element.
Well, if we look at the plane stress case and want to construct our b matrix, remember the b matrix gives the strains in terms of the nodal point displacements. And we remember also that our epsilon xx is equal to a dy dx, our epsilon yy is equal to dv dy and our gamma xy is equal to a du dy plus dv dx. Well, if we recognize these facts. And we also note again that for the two by two element, r is identical to x, s is identical to y. Then, we can obtain the epsilon rr or epsilon xx by simply taking the derivatives of the hi's with respect to r. Notice here, since u is equal to the summation of h i u i, dy dr, which is equal to du dx, is nothing else then the summation of partial h i dr u i. And this part, which runs from one to four-- i going from one to four-- I simply put right in there.
I proceed similarly for epsilon yy. Now I'm talking about the v displacements, which are stored after the u displacements. That's why I'm looking here at the second column, the last column. And here we have dh one ds, because we are talking about the derivative with respect to y or with respect to s, which is the same thing. So here we have the entries for the epsilon ss part. Now for the strain part, we are talking du dy or du ds, dv dx, dv dr. And all we have to do in now is take this term and put it right in there, and take this term and put it right in there. And then the last row here gives us the shearing strength.
So this is our b matrix for the special two by two element. It is constructed in a very simple way. The h i are known, and if we had five or six or seven nodes, we would proceed in exactly the same way. All we would have to do is include additional columns in the b matrix that would give us the appropriate entries for strains generated by the additional nodal point displacements.
Let us now look how we can generalize these concepts to the element that is not, anymore, in the physical x and y space, two by two element. In the physical x and y space, this element-- four-noded element now-- might look as shown here. However, what we do is we still deal with the r- and s-coordinate system embedded on the element. We are on s one here, minus one s plus one here, r and s both minus one here, plus one here s minus one. So in the natural coordinate space, the rs space, we still have a two by two element.
The interpolation of the displacement therefore-- the interpolation of the displacements-- even for the distorted element, is exactly still as shown here. Provided we are entering always with the appropriate r and s. H one is a function of r and s. So if we look at a point in the general element, if we look at a point in the general element, here, for example, at such a point. Because that point has a specific r and s value. And if we want to find the displacement at that point, well, we would have to put the r and s value of that point into h one, h two, h three, h four, and that gives us then the displacement at that point, in terms of the nodal point displacement. So, the displacement interpolation for this element can still be done in the r and s space.
However, difficulties arise-- or additional considerations I should say rather-- arise when we talk about strains, because the physical strains that we have to deal with are derivatives with respect to x and y, and not r and s anymore. Well, so what we have to do is, use a Jacobian Transformation. What we want are the derivatives with respect to x and y. What we can find easily are derivatives with respect to r and s on the displacements, because the displacements are given in terms of r and s values. So, remember u is equal to some of h i u i and the h i is a function of r and s, so, we can directly find derivatives with respect to r and s of u. What we cannot find easily are derivatives with respect to x.
This is the relationship that we use. It gives us a transformation from derivatives of x and y to derivatives r and s. A question must immediately be in your mind, why do we not write down directly this relationship, which is given by the [INAUDIBLE]. If we want d dx of displacements, why not just use d dx being equal to d dr, dr dx, and so on. Well, the difficultly is that we cannot find dr dx very easily. We have x being this function of h i x i. This is the interpolation which we have to use now.
I mentioned earlier on the first slide that we interpolate displacement and coordinates in the same way. So, here we have a linear interpolation of the coordinates, from this node to that node, and in between here too. So we are using this interpolation here on the coordinates. And we can easily dx dr, but we cannot easily find dr dx. You would have to invert this relationship somehow so that we have r in terms of x. Well, it is easier, therefore, to write this relationship down, which is really the chain rule. This is also the chain rule, its a chain rule the other way around, which gives d dr being equal to d dx, dx dr plus d dy dy dr, if we multiply this out. And its this relationship that we can use effectively.
Well, this relationship in three dimensional analysis would involve a third row and third column. In one dimensional analysis we only talk about one by one. In other words, in one dimension analysis for a truss we just have that entry. In general we can write it in this way, where j is the Jacobian transformation from the xyz-coordinate system to the rst-coordinate system. And since we want these derivative-- because these derivatives give us actual strains, we have to invert this relationship.
Having constructed then, these derivatives in terms of these derivatives, which we can find just as we have done before, we can now establish the b matrix in much the same way as earlier. And since we now have h and b matrices for an element-- these are a function r, s, and t-- we would perform the integration. And now I'm referring to the integration of the stiffness matrix. Remember k is equal to b, transposed cb over the volume. Now notice that since the b matrix that we are using in here is a function of r and s in a two dimensional analysis. And r runs from minus one to plus one, and s runs from minus one to plus one. We now have to use a transformation also on dv to integrate over the r s volume. And that integration-- and that dv element is expressed as shown here-- and that is given to us from mathematical analysis.
So basically what we are saying here is that we are replacing this integral over the physical volume by an integral minus one to plus one, minus one to plus one. And that signifies from minus one to plus one over r and over s, if we had a three dimensional analysis we would have another integral sign here. B transposed now r and s function of r and s, cb, function of r and s. And then our dv, now in terms of r and s. So this is how we really do things. And remember this dv here, this dv, is that one there. This integration is effectively performed using numerical integration and I will discuss later on.
Let's look once at the Jacobian transformation for some very simple examples-- the Jacobian transformation for some very simple examples-- just to make things a little clearer. In this case we really have taken our two by two element and we have stretched it into the x- and y-axes. That stretching is giving us a three here and a two there because our two by two element has a length of two here, and six divided by two gives us three. Similarly here we have stretched the element by a factor of two. This relationship here is, in general, calculated-- the j is, in general, calculated-- as shown here. And the result of that, by putting these interpolations into there are these values here. Physically what these mean is a stretching, in this particular case, into the x- and y-axis.
Somewhat-- the case where we cannot directly-- not easily directly write down to j matrix is this one. Here we would go through the actual evaluation the way I have indicated it. In other words, we would go through this actual evaluation, substituting from here and of course for y also. And this would be the result Notice that we now have a stretching here of three, compression from a two length to a one length-- therefore we have a one half here-- and there's also an angle change that gives us the off-diagonal element here.
Another interesting case here, as an example, we have the same lengths here-- two same lengths here-- but a distortion in the element because this node two has come down from there to its midpoint. And the resulting Jacobian is given here. Now notice that that Jacobian is a function r and s, so the inverse, which is used in the construction of the b matrix would also be a function of r and s. A particularly interesting case is the one where we shift nodes to advantage. See here we have our original-- our three-noded element that we talked about already-- in the r space now. Its a truss element. And let's say that in our actual physical space, we have this node there, the three node there, and the two node there. The element in the actual physical space has a length of l. We have taken this node and shifted it to the quarter point of the element, to the quarter point of the element.
What I will show you right now is that by having done so-- by having taken this node from its midpoint and shift it over in the actual physical space to the quarter point-- we will find that the strain has a singularity here of one over square root x. This is a very important point which can be used in the analysis of fracture problems, because we know that in the analysis of fracture problems, we have a one over square root x singularity at a crack tip. And, if we want to predict the actual stress there or the displacement around the crack tip, it can be of advantage to use this fact, shifting nodes to quarter points in order to capture the stress singularity more accurately.
Well let me show you then how this strain or stress singularity comes about. If we look at this element here and we use our interpolation on the coordinates, this would be the result. Now, notice that we have substituted the x i values here. X one is zero, x two is l, x three is l over four. We have substituted those values and directly come up with this result. Well, we can see that this indeed is true. Let's put r equal to plus one n. In other words, the right hand side node-- for the right hand side node-- and we would have a two here squared, gives us four, goes out with that four. So at i equals plus 1 we have x equal to l, which is correct of course. Let's put i equal to minus one n, we find x is equal to zero. Let's put i equal to zero n we find x is equal to l over four. In other words, this has been a simple check in that we have the right interpolation-- geometry interpolation-- for this element.
Our j now is simply, dx dr, X is given here. If you take the differentiation of that you get the two in front that gives l over two times one plus r, this value here, in other words. Then our b matrix is constructed by the inverse of the j-- that is this one-- two times the r derivative of the interpolation functions. Of course here we talk only about one strain. Remember, in the truss, the only displacement of concern is the u displacement and the strain is simply epsilon xx, a strain into this direction also.
Well, this is our b matrix then. And if we take our h one, two, and three, and we differentiate these-- as indicated here with respect to r-- we directly obtain these functions here. If we now recognize that since we have x related to r here, we can also invert this relationship. We can write r in terms of x. And if we have done so we can take that relationship and put it right in here. Then we would get b in terms of x.
And that's what I have done here. The first line shows simply r in terms of x now, and I have substituted that r value into the b matrix and this is there is the result. Notice that we have a strain singularity. This is the first element here. The next element-- this is here-- the next element in b matrix, and that is the third element in the b matrix. Notice that we have in this element, that one, and that one the one over square root x, which means that we have one over square root x singularity at x equal to zero.
Well, this fact is used very effectively in fracture mechanics analysis. Assume that we have crack here and that we want to analyze the stress conditions around that crack. What we can do is we use now two dimension elements, as its a plane stress situation. We would put a two dimensional triangular element there and we shift the midpoint nodes. This is now very small here, but I hope you can still follow. You shift these midpoint nodes to the quarter point, just the way we have been doing it here. You are putting the third node the quarter point and the result is that at this crack tip we have a one over square root x singularity, using this element layout. And we know that in fact, there is a one over square root x singularity in linear fracture mechanics analysis. And so this is an effective way of capturing this singularity, and has been used or is currently being used very abundantly in practice.
The important point that I wanted to make really is that we can shift nodes in the element to our advantage. But, we really do that in specific applications such as fracture mechanics. In general we will see later when I talk about modeling of finite element systems, in general, it is most effective to leave the mid-side nodes at their physical midpoints. In other words for an eight-noded element, two dimensional analysis, we would put this mid-node in the physical space also, actually at the midpoint. We would not shift it. Then the element has good convergence characteristics into all directions and this is really how the element is used most effectively for general applications. However, in specific applications, such as fracture mechanics analysis, it can be of advantage to shift these mid-side nodes to pick up certain strain or stress singularities that we know do exist.
Now on the last transparency that I wanted to show you, I wanted to indicate something to you that I will be talking about in later lectures more abundantly, namely the fact that we're using numerical integration. B, for the k matrix as an example, is once again now a function of r and s. This part here is also a function of r and s. So we have here function of r and s. and we have here also a function of r and s. So this f matrix here is a function of r and s. Notice that the b also includes the inversion of the j, the Jacobian matrix. It includes the inversion of the j, because we had to construct the x and y and z derivatives from the r, s, and t derivatives.
So, what we do in practical analysis is that we use numerical integration to evaluate the k matrix. I have indicated that here schematically, if you look at this element here. What we do is, we evaluate the f matrix-- this is a matrix. In two dimensional analysis we would only run over i and j. I is this direction. J is that direction. In three dimensional analysis, which is in general analysis, we run i, j, and k this way. We evaluate the f matrix here at specific points, r, s, and t. T now being this axis. And then multiply that f matrix by certain weight constants and sum these contributions over all i, j, and k, in order to obtain an accurate enough approximation to the actual stiffness matrix.
The order of approximation with which we obtain the actual stiffness matrix-- or rather how closely the numerically calculated stiffness matrix approximates the actual stiffness matrix-- depends on, number one, how many integration points we are using and what kind of integration scheme we are using. These points here correspond to the Gauss numerical integration. In this case for two dimensional analysis, we would use a two by two integration. In other words, i and j would both run from one to two. K is not applicable and we would have altogether four evaluations off the Fs here. Multiply each of them by weighting factors, which has been derived for us by Gauss some long time ago. And summing up these contributions gives us a close enough approximation to the k matrix.
Of course, the question of how many points we have to use, what integration scheme we should use, is a very important one. We must use enough integration points to get a close enough approximation to the actual stiffness matrix and I will be addressing those questions in a later lecture. This is all I wanted to say in this lecture. Thank you very much for your attention.
MIT OpenCourseWare makes the materials used in the teaching of almost all of MIT's subjects available on the Web, free of charge. With more than 2,200 courses available, OCW is delivering on the promise of open sharing of knowledge. Learn more »
© 2001–2015
Massachusetts Institute of Technology
Your use of the MIT OpenCourseWare site and materials is subject to our Creative Commons License and other terms of use.