Lecture 20: Beam, Plate, and Shell Elements II

Related Resources

Study Guide (PDF - 1.2MB)

Readings

Section 6.5

Examples

Problems 6.20, 6.21

References

Bathe, K. J., and A. Chaudhary. “On the Displacement Formulation of Torsion of Shafts with Rectangular Cross-Sections.” International Journal for Numerical Methods in Engineering 18 (October 1982): 1565-1568.

Dvorkin, E., and K. J. Bathe. “A Continuum Mechanics Based Four-Node Shell Element for General Nonlinear Analysis.” Engineering Computations 1 (1984): 77-88.

Bathe, K. J., and E. Dvorkin. “A Four-Node Plate Bending Element Based on Mindlin/Reissner Plate Theory and a Mixed Interpolation.” International Journal for Numerical Methods in Engineering 21 (February 1985): 367-383.

Bathe, K. J., and E. Dvorkin. “A Formulation of General Shell Elements: The Use of Mixed Interpolation of Tensorial Components.” International Journal for Numerical Methods in Engineering 22 (March 1986): 687-722.

Bathe, K. J., and P. M. Wiener. “On Elastic-Plastic Analysis of I-Beams in Bending and Torsion.” Computers & Structures 17 (1983): 711-718.

Bathe, K. J., C. A. Almeida, and L. W. Ho. “A Simple and Effective Pipe Elbow Element: Some Nonlinear Capabilities.” Computers & Structures 17 (1983): 659-667.

Lee, P.S., and K. J. Bathe. “Development of MITC Isotropic Triangular Shell Finite Elements.” Computers & Structures 82 (May 2004): 945-962.

Bathe, K. J., and P. S. Lee. “Measuring the Convergence Behavior of Shell Analysis Schemes.” Computers & Structures 89 (February 2011): 285-301.

• Download this transcript - PDF (English - US)

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 MIT.edu.

PROFESSOR: Ladies and gentlemen, welcome to this lecture on non-linear finite element analysis of solids and structures. In this lecture, I'd like to continue with our discussion of structural elements, and in particular, I'd like to concentrate our attention first on the formulation of beam elements. After that, I'd like to talk more about the formulation of shell elements, and then show you applications.

When we talk about beam elements, we might think about the usual Hermitian beam elements, in which the transverse placements cubically interpolated, and the longitudinal displacements are linearly interpolated. This is an element that is usually used in linear analysis, and very widely used-- you might have used it quite a bit yourself-- and it is also very effective in the analysis of linear response of structures. However, in my discussion now, I will talk about the isoparametric beam element, and how we formulate the isoparametric beam element, in other words, and how we apply it in the analysis of structures. So when I talk about the beam element, I really mean this isoparametric beam element that we are formulating in this lecture.

The isoparametric beam formulation can be very effective for the analysis of curve themes, for the analysis of geometrically nonlinear problems, and for the analysis of stiffened shell structures, where we couple the isoparametric beam with the isoparametrically formulated shell elements. The formulation is quite analogous to the formulation of the isoparametric degenerate shell elements that we discussed in the previous lecture, and so we can go a bit fast over the formulation of the beam element here.

I'm showing a typical beam element, rectangular cross section. Here we show three nodes. Actually, the beam element can be used with two nodes, three nodes, or four nodes, as we will discuss further later on. We notice that the depth off the element is ak, the thickness of the element is bk. And notice that at these nodes, we measure certain nodal point variables. and here I'm showing the nodal point variables that we do measure at the nodes. Three displacements, incremental displacements, and three incremental rotations. Notice that the description of the beam element follows much the same as the geometric description of the beam element follows, much the same way as we are dealing with the shell element.

We introduce director vectors. Now we have two director vectors. In the shell, of course, we had only one director vector per node. Now we have two. Here we show tvsk, the director vector at time t into the S-direction, corresponding to node k. Notice here tvtk. The director vector at time t corresponding to node k into the t direction. this t means direction t that t opting out means time t Notice that corresponding to this direction t here, we have the t isoparametric coordinate. And of course, this t isoparametric coordinate also runs into the thickness ak. Notice that corresponding to bk, we have the s-direction coordinate, and that s-direction coordinate, s-isoparametric coordinate, runs into the direction of the director vector, tvsk.

So these are the variables used to describe the geometry and the displacement of the beam element. We start with the general equation for the geometry at time t. Here we have txi. The coordinates of a material particle within the element at time t-- and remember, these are the coordinates measured in the stationary coordinate system. And these coordinates are obtained as shown on the right-hand side.

Let's have a close look at what we see on the right-hand side. We have the hk interpolation function. This is the interpolation function corresponding to the r-direction running through the nodes of the element, along the mid-surface of the element. And it's a one-dimensional interpolation function, the same one that we used to see for the truss element. These are the nodal point coordinates at time t. We are summing over all the nodal points, k from 1 to n, and n can be equal to 2, 3, or 4, as I mentioned earlier. Then we add to this term a term that comes in because of the sickness of the beam into the t-direction, ak, thickness off the beam into the t-direction-- here is a t-coordinate. The isoparametric coordinates going from minus 1 to plus 1. Here we have hk, same interpolation function as there, and here we have the direction, cosines of the director vector in the t-direction at time t.

A similar term, similar to this one, is added for the f-direction. s, the isoparametric coordinate, running from minus 1 to plus 1, bk being the thickness into the s-direction at node k, hk, same interpolation function that we have seen up here already, and these are the direction cosines of the director vector in the s-direction, corresponding to node k at time t.

This equation is very similar to the equation that we used in the shell element formation. Of course, in the shell element formulation, we interpolated here in a two-dimensional domain over the shell mid-surface. Now we only run along the neutral axis of the beam. And we had only one term in the shell corresponding to the director vector. There was only one direction, namely, the one normal to the mid-surface of the shell. Here we have two directions to deal with.

Since the displacements can directly be obtained from the geometry interpolation at time t minus the geometry interpolation at time 0, of course this geometry interpolation is obtained by simply substituting in the equation that we just looked at for time t, the time 0. Since we obtain tui, as shown here, this is the result. We simply apply our geometry interpolation at time t, and at time 0, subtract. This is what you get on the right-hand side.

Nodal point displacement from time 0 to time t, measured in the stationary global Cartesian coordinate system. Here change in the director vector cosines, or I should say, more precisely maybe, the change in the cosine of the angle director vector time 0, director vector time t. Of course, remember here we have direction cosines corresponding to time t. Here we have direction cosines corresponding to time 0. We subtract, and that goes in here.

Here, same kind of term, but for direction s. Also the incremental displacements are obtained from this relationship here, and if you apply the geometry interpolation at time t plus delta t and at time t, you directly obtain this right-hand side. These are the increment in the direction cosine from time t to time t plus delta t of the director vector in the t-direction. Similar here for the s-direction.

Of course, we are used very well to deal with increment and nodal point displacement, but we cannot very well deal with these quantities here, vti and vsi. We want to deal with nodal point rotation. And so what we do is, we express this quantity and that quantity in terms of nodal point rotations, and that those expressions are given down here. These are the two expressions. You can simply multiply out, and you would immediately see that they indeed hold.

With these geometry and displacement interpolations established, we can now directly calculate or establish the strain displacement matrices corresponding to the Cartesian strain components, and then a simple transformation, standard transformation gives us the strain displacement matrices corresponding to the eta zeta psi directions. Of course, these are the directions with the components of strains and stresses that we need to deal with for the beam formulation. Notice that zeta here is a physical coordinate running into the direction, so to say, of the a thickness of the beam, and zeta corresponds to the t-isoparametric coordinate that we just talked about. Psi here corresponds to the s-isoparametric coordinate, and runs into the b width, so to say, of the beam. Eta here corresponds to the r-isoparametric coordinate.

Of course, these zeta, eta, and psi are actual physical coordinates that we are dealing with. We have to, of course, calculate the tau psi eta, the tau eta eta, and the tau vita psi stress components and corresponding strain components for the beam formulation. And the stress-strain law, corresponding to these components, is given here. Notice tau eta eta is obtained by taking the Young's modulus and multiplying the Young's modulus by epsilon eta eta, and so on. Notice here we have one share component and another share component corresponding to this column and corresponding to that column, and here we have k factor, the share correction factor, that we already used from the shell analysis formulation that we discussed earlier.

Notice that via the same approach, we can also directly obtain a beam element for elastoplastic and creep analysis. In this particular case, of course, the only stress component that would be non-zero again would be tau eta eta, tau eta zeta, and tau eta psi. The approach here is very similar to what we have discussed earlier regarding the shell element formation.

There is one important point in the beam element formulation. Namely that with the kinematic assumptions used so far for the beam element, we do not allow cross-sectional out-of-plane displacement. We do not allow warping. Our interpolation does not contain the effect of warping. In torsion or loading, however, we know that warping is very important, and for that reason, we need to amend our displacement assumptions by two or more displacement patterns, so to say. We allow these displacements to take place in the element. u eta is equal to alpha times psi beta. This is the exact warping displacement for infinitely narrow sections. And this term here, with beta a constant, is the exact warping displacement for a square cross section.

Now the actual section might not be infinitely narrow, might not be square, but we can easily test, and I will show you just now an example of that in fact, a superposition of these functions and the evaluation of alpha and beta by static condensation actually yields good results, even for cross-sections that are not infinitely narrow, and that are not exactly square.

Let's look at an example, linear elasticity. If we take a section and simply put it into torsion with b/a equal to 1, the analytical value given by Timoshenko for this k constant here in that formula is 0.141. And we obtain with this formulation exactly that same value. We need to get it, because we have embedded the proper interpolation for u eta into our formulation, and then u static condensation, of course, to condense out these two constants, alpha and beta, that I just talked about. For b/a, a very large, very narrow section, we also get the exact result in our finite element solution, the exact analytical result. In between, we have a small arrow, an arrow that, however, we believe is quite acceptable.

Let's look at one example that I'd like to show you the results of. You might be interested in this example. Here we have a narrow ring. Here you have a plan view of that ring with the dimensions given, and Young's modulus and Poisson's ratio given for that ring, subject to these bending moments. Of course, the twisting capability or capacity of the beam is of utmost importance now, has to be properly modelled under this loading condition.

The model that we use to analyze the problem is shown here. And you can see, at one side of the ring, we have fixed it, and at the other side, we prescribed a rotation, see eta x. We use here four node elements as shown. This analysis was actually inspired by us looking around in our office and seeing a magnetic tape lying there. And if you look around, you might have yourself a magnetic tape there. And you pick up this little ring from the magnetic tape. And here I'd like to just demonstrate to you briefly what we are doing. You see here, we have this ring, much alike of what I've shown you on the view graph. I'm fixing it on the one side, and on my right hand now, these two fingers, I'm going to put a rotation onto that ring, and I'm going to turn my fingers 180 degrees.

So please watch closely. This is what I'm doing. I'm just turning my fingers 180 degrees. And by turning my fingers 180 degrees, if you watch very closely, you can see that the top of the ring does not touch the bottom of the ring. There is actually a little gap in between, of course, because of the elasticity in the ring.

So this is what I've tried to do now on the computer. And what I'm particularly interested to see is the amount of moment required in my right hand, so to say, to be transmitted to the ring in order to twist it. Well, here we have the results for the computer solution. Notice here, we have a force deflection, curve, rotation here, up to 180 degrees. Moment plotted here in pounds, inch. And you can see that there's fairly linear behavior first, but then a highly non-linear behavior as the rotation becomes very large.

Here we use the TL formulation for the beam element. And once again, we rotated it 180 degrees, as I have shown you with my little experiment here. If you look at the beam, you see a top view like that, and the side view looks this way. I've pointed out that top and bottom do not touch each other. There's a little air gap in between.

I'd like to show you now another solution regarding the torsion in an analysis. And let me walk over to the slides here, where I actually have two slides regarding that solution. The interesting part of this problem is that we now, in this particular analysis, are using the same displacement interpolation functions for which I've shown you already-- two analyses in which the material remained elastic. The same interpolation functions are now also used to analyze an elastoplastic problem.

And the problem is described here, where we have now a square section subjected to a torque. The material data that Greenberg used are given here. And in our analysis, we used this material data, where the yield stress and the strain-hardening modulus are matched to this stress-strain law here.

Now notice what is going to happen is, we are increasing this torque incrementally, and we want to bring, of course, this section into the elastoplastic regime to see how the section will take the torque. On this slide now, we have plotted the torque vertically and the rotation horizontally, both quantities with some constants attached, as you can see. And the Greenberg result and the Adina result are indeed very close. Please look at the reference that is given in the study guide if you're interested in seeing and reading more about this problem.

Let me now walk back to the view graph and continue with the discussion of the material relating to isoparametric beams and shell elements. The next topic that I like to talk to you about is the actual use of these elements. Of course, the shell element that I'm here talking about, or that I will share some more experiences with you upon, is the shell element that I have been already discussing in the earlier lecture. So please refer back to that lecture with respect to this discussion.

One interesting point is that these elements can all be programmed for use with different numbers of nodes. For the beam, we can have two, three, or four nodes, as already mentioned earlier, and for the shell, we could have 4, 8, 9, up to 16 nodes, as I mentioned also in the earlier lecture. The elements with these nodal point configurations can be employed for analysis of moderately thick structures. Remember, shear deformations are approximately taken into account.

However, if the sink of the analysis of thin structures, then we have to be very careful, and only certain elements of the ones that I just referred to should be employed. For shells, we can only recommend to use a 16-node element with 4-by-4 Gauss integration over the mid-surface. This element is depicted here. 1, 2, 3, 4, 4 times gives us 16 nodes, and 4-by-4 integration indicated by the blue crosses.

For beams, we'd best use the 2-node beam element with 1-point integration along the r-direction, or the 3-node beam element with 2-point integration along the r-direction, all or with Gauss integration, or the 4-node beam element with 3-point Gauss integration along the r-direction. We have to ask, of course, now, why is that so? Well, the reason is that the other elements become overly and artificially stiff when we use them to model thin and curved structures.

Two phenomena occur-- shear locking and membrane locking. Most interesting phenomena. And I'd like to share with you now some information regarding these phenomena. The 2-, 3-, and 4-node beam elements with 1-, 2-, and 3-point Gauss integration along the beam axes do not display these phenomena, these shear- and membrane-locking phenomena. But please notice, I'm using 1-, 2-, and 3-point Gauss integration corresponding to the 2-, 3-, and 4-node beam element.

The 16-node element with 4-by-4 integration, Gauss integration, is also relatively immune to shear and membrane locking. However, we have to be careful here. The element should not be distorted for best predictive capability.

Let us study what this shear locking phenomenon is. And to do so, let us look at a 2-noded beam, a very simple beam, in just 2-dimensional action. So we have one node here and one node here, and the defomations of the beam are described by w1 and theta 1 at this node, and w2 and theta 2 at this node. Notice the deformations of the beam, of course, are given by w, the transverse displacement, and beta, the section rotation. The transverse displacement and section rotation, w and beta, are interpolated linearly, as shown here, via the nodal point displacements and rotations. Notice the linear interpolation, of course, is given here via this expression and that expression, and similarly here via this expression and that expression.

Now there's one important point, namely that w and beta are independent quantities. If you are familiar, as you probably are, with Euler-Bernoulli beam theory, you recognize that in Euler-Bernoulli beam theory, they are not independent quantities. In fact, beta is dw dr, or dw dx, x being the physical coordinate.

Here in our formation, isoparametric degenerate element, we have that beta and w, these two quantities, once again, are independent and linearly interpolated in this particular beam element. Of course, these two quantities are tied together by the physics, and that is by the shear deformation. So let's look at that next.

Here we have that gamma shear strain is equal to a del w del x minus beta. And this is shown here on the picture. Notice this is a neutral axis of the beam, which originally of course was horizontal. Here we have a vertical line shown in green, and the section rotation is shown at the angle beta here. These are the material particles shown in red that were originally, of course, at right angles to the neutral axes, and that were originally lying vertically up. In other words, this line here, this red line was a vertical line originally, and was originally at right angles to the neutral axis. That right angle is not preserved because gamma is there, and gamma is, in general, a non-zero.

This is here of course dw dx. So this picture here shows what we are calculating up here. And you can see that gamma ties together w and beta, which are, however, by themselves, interpolated independently.

Consider now the simple case of a cantilever subjected to a tip-bending moment. And we use one 2-node beam element to model the cantilever. Clearly at this end, the section rotation is 0, and w1, the transverse displacement is also 0. So these are the bounded conditions at this end of the cantilever. And if we use these bounded conditions and put them into our general equation, we get directly that beta is given as shown here, and gamma is given as shown here.

Now look at this equation carefully. We have a constant term and a linear varying term. The exact solution to this problem, analytical solution to this problem, if we use elementary beam theory, we would immediately recognize that the shear strain must be 0 for this element, or for this beam situation. The shear strain should be 0.

Now, if we look at this carefully, we can see that gamma can only be 0 if w2 is 0, and theta 2 is 0. But if theta 2 and w2 are 0, rotation and transverse displacement right here at this end, then the beam would have not deflected. Of course, that is not borne out in physics. We know that if we take a beam and put a bending moment on it, then we get a deflection here and a rotation at this end.

So you can see already from this equation here that there will be some difficulty. There will be some difficulty because gamma should be 0 physically, and this equation tells us that gamma can't really be exactly 0. Well, that is once more summarized in this sentence here-- clearly, gamma cannot be 0 at all the points along the beam unless theta 2 and w2 are both 0. But then also beta would be 0, and there would be no bending of the beam.

Let's look now at this point. Since for the beam, the bending strain energy is proportional to h cubed, elementary beam theory again, the shear strain energy is proportional to h-- elementary beam theory-- we see immediately that any arrow in the shear strains due to the finite element interpolation functions becomes increasingly more detrimental is h becomes small. Why is that the case? If h becomes small, then this number becomes very rapidly small. This number here become small, but this one goes much faster to us, very small, so to say. And if there's an error in the shear strain energy, that error, because of this proportionality sign, is going to be magnified relative to any error that we would see here, and it's going to be detrimental in our solution.

Let's look at an example. For the cantilever example, this very simple example that we just studied already a bit, the shear strain energy should be 0, and as h decreases, the relative error in the shear strain increases rapidly, and in effect introduces an artificial stiffness which we identify as locking of the model. This is the terminology that is quite widely used now. We say that the model locks.

Let's look at some results. In this table, we show the results for h/L with L equal to 100 for these three values. The analytical value for the rotation at the end where the bending moment is applied is given here. Bernoulli beam theory. The finite element solution, using, of course, exact integration-- I always assume I should maybe point that very strongly still out. We assume that we use exact integration of the K matrix. Meaning in this particular case, we're using 2-point Gauss integration along the length of the beam.

With that scheme then, we obtain with a h/L being 0.5. Well, an error, but something still reasonable. After all, we're only using one element. Of course, we are pretty far off, but still something reasonable. But now notice, as you decrease h, h/L, of course, becomes much smaller. Theta analytical goes up. The beam basically bends, is very flexible as it gets thinner. And the finite element solution also becomes more flexible, gives us a more flexible beam, but is much too stiff when compared to the analytical solution. Notice for orders of magnitudes difference. If you plot it, this is basically 0 compared to that value, and that means the beam locks. It doesn't really at all bend anymore as h becomes very thin.

This behavior is, of course, also observed when you take more than 1 element to model the beam. Here we show more elements than just one, and notice that each of these elements should still carry a constant bending moment. If it cannot do so, it develops, in other words, these furious shear strains. We are going to have a very stiff model, even when we use that many elements.

Let's look at this particular example close off. And here we have an example beam, in other words, cantilever beam, subjected to the bending moment. l equals 10 meter. Square cross section, height 0.1 meter. We use the 2-noded beam elements to model this cantilever, and we plot the tip deflection as a function of the number of elements. A most interesting graph.

Look, here is the beam theory solution. And here we plot the number of elements-- here we talk about 400 elements, because there's times 10 to the 2 down here. 400 elements at this point. 600 elements, 800 elements, 200 elements. At this point we have 100 elements. The height of an element is equal to the length of the element. And at this point, we still don't even get a good approximation to the beam theory solution. We need basically something like 400 to 500, 600 elements to get close to the beam theory solution. However, if we take a very large number of elements, we actually do converge, but we need a lot of elements, as you can see here.

Of course, it is quite impractical to use this kind of 2-noded beam element to model thin structures, then beams. Once again, I must point out that we use 2-point Gauss integration into the r-direction, and therefore, we're dealing with the exact K matrix corresponding to the interpolation that we have embedded into the element.

In fact, if we look in this regime down here, we find that as shown on this view graph, the solution really becomes beta very slowly as the number of elements is increased. Notice that the solution really doesn't increase in accuracy much at all as we go from 10, 20, to 30 elements.

A remedy for the 2-node beam element is to only use 1-point integration along the beam axis. This then corresponds to assuming a constant transverse shear strain. However, by using 1-point integration, we still integrate the bending strain exactly. Because remember, beta is linear, del beta del x is constant, and 1-point integration will pick up that term exactly. If we look at the results, now once again, of this problem, cantilever problem, the cantilever beam is subjected to an end moment, tip end moment. We find now which h/l, as shown here, theta analytical, shown here, and our finite elements solution with just one element, using one element, of course now with 1-point integration in r-direction, gives us excellent results. In fact, the exact results, as you can see here.

The 3- and 4-note beam elements, evaluated using 2- and 3-point integration, are similarly affected. In fact, they give us also the exact results for this problem that we've considered. We should note now that these theme elements, based on what one might reduced integration, are actually reliable elements, because they do not possess any spurious zero energy modes. They have, in other words, only six zero energy modes corresponding to the actual physical rigid body mode that the elements should contain, should be able to represent, in its three-dimensional space. So they only have those six zero eigenvalues in this 3D analysis that they should have, and no spurious zero energy modes.

The formulation can actually be interpreted as a mixed interpolation of displacements and transferred shear strains. And this is really the key for developing all the low-order effective shell elements. I will get back to that in just a little while.

So far, we have talked about shear locking. I mentioned earlier also that there is a phenomenon of membrane-locking. Membrane-locking means that in addition to not exhibiting erroneous shear strains, the beam model must also not contain erroneous mid-surface membrane strains, particularly when we analyze curved structures. The beam elements that I just mentioned, 2-noded beam with 1-point r-integration, 3-noded beam with 2-point r-integration, 4-noded beam with 3-point r-integration, do also not membrane-lock.

Let's consider here, once again, a simple example-- a simple curve beam subjected to a bending moment, tip-bending moment, fixed at the other end. Here is an angle alpha that the beam spans out. The exactly integrated 3-noded beam element, when curved, does contain erroneous shear strains and erroneous mid-surface membrane strains. In other words, shear locks and membrane locks. And that is shown in this table here, at least to some extent. What I can't show is how much there is shear-locking and how much there is membrane-locking. You would see that both effects are very significant.

If you were to run, for example, this beam element with a computer program, using, for this beam element, 3-point integration into the r direction and printing out the stresses, print out all the stress components and look at the shear component and the membrane components, in other words, the normal stress components, along the neutral axes. If you do so, you would see that they are very significant shear stresses and very significant membrane stresses, which of course should be 0 for this problem. And these are, of course, causing the locking effect.

Here we have the results in terms of just rotations, h/R now, with R equal to 100, given here, analytical value given here, finite element solutions, 3-node element with 3-point integration, you can see we have a locking phenomenon right here. Once again, how much there is membrane, and how much shear locking you would see by actually running as a problem and looking at the shear stress and the membrane stress. And the finite element solution with 2-point integration gives excellent results. So this is, of course, the element that we would be using in engineering practice.

Similarly, we can study the use of the 4-node cubic beam element. And here we find that for h/R, these values, analytical value once more listed here, we find that the 4-node element with 4-point integration gives also very good results, and the 3-point integration gives also excellent results.

So we note that the cubic beam element performs well, even when using full integration. It is not susceptible to membrane and shear locking. But notice, we have our third-point node exactly at the third point of the elements. Once we start shifting these nodes, then we would also not get as good results as shown here. We would not get these good results once you shift the node at the third point away from their physical third point locations.

Now, considering the analysis of shells, the phenomenon of shear and membrane locking are also present, but the difficulty lies in that this simple reduced integration approach that we're using for the beam elements cannot be directly applied, because the resulting elements contain spurious zero energy modes. For example, if you integrate the 4-node shell element just with 1-point integration, contains 6 spurious zero energy modes. In other words, it has altogether 12 zero eigenvalues where it should only have 6. The 6 of them are spurious modes. And such spurious energy modes can lead to very large errors in the solution that, unless we have a comparison with accurate results, are not known, and therefore, it can be very dangerous to use such elements. We want to have elements that are reliable in general applications, and these elements should only contain the actual physical rigid body modes, no spurious zero energy modes. And of course, should also have high predictive capability, not have any shear or membrane locking.

For this reason, we can only recommend the 16-node shell element with 4-by-4 Gauss integration. The cubic beam element result that I just showed you regarding the analysis of the curve beam already indicated that the cubic element does quite well. In fact, that is also borne out by the 16-node shell element. Using 4-by-4 Gauss integration on the shell mid-surface, we have an element that is relatively immune to membrane and shear locking, but I mentioned already that the element should not be distorted. And if it's not distorted, it performs best.

I'd just like to share now with you some thoughts regarding development that we have recently pursued. Namely, recently we have developed elements based on the mixed interpolation of tensorial components. This sounds very complicated. Actually, if you look closely at what we're doing, it is not that complicated. The elements, however, do not lock in shear or membrane action, and also do not contain any spurious zero energy mode. And of course, that is the key of success of an element.

We will use later on in the example solutions a 4-node element that is based on this mixed interpolation of tensorial components. We call that element the MITC4, Mixed Interpolation Tensorial Component with 4 nodes element, and I want to discuss this element just briefly with you now. We also have developed an element that we call MITC8 element. It's an 8-node element based on mixed interpolation of tensorial components. A most interesting development. If you're interested in that development, please refer to the paper that is given, or that is being referred to, in the study guide.

I'd like to discuss with you now the MITC4 element, and then show you example solutions. However, we have, in order to be able to do that, to change now the reel, so let us just do so, and then continue with this discussion.

So let us look then at the MITC4 element briefly, the element that I mentioned before we had the break. The element is described using four nodes. We have the same isoparametric coordinate system that we are usually having. Notice here director vector, same description as for the shell element that we discussed in the earlier lecture. Each node has 5 degrees of freedom, 2 rotational degrees of freedom, and 3 translational degrees of freedom. Once again, the way we have been discussing it in the earlier lecture. We use the element for analysis of plates and for the analysis of moderately thick shells, and also thin shells. And that is the important point. This element is directly applicable in a very effective manner to thin shell analysis.

The key step in the formulation of the element is to interpolate the geometry and the displacement as earlier described, as we discussed in the earlier lecture. But to interpret the transverse shear strain tensor components separately, and these interpolations are selected judiciously, to tie, then, the intensities of these components to the values evaluated using the displacement interpolations.

Let's look at what I'm saying here now pictorially. The strain tensor interpolation that we are using is shown here in blue on the black element. These are the rt transfer shear strain tensor components. This is one interpolation function, and this is the other interpolation function. Notice that if we look at the r- and s-axes, the rt strain tensor component is constant in r, but varies linearly with s. Here, too, constant in 4, varying linearly with s.

Now it's these intensities, one here and one there, that we tie to the nodal point displacements and rotations. Here for this intensity, we use these nodal point descriptions, and for this intensity, we use these nodal point descriptions, or nodal point degrees of freedom.

So if you recognize that these are, it seems, unknown quantities that enter into the formulation since we're using these two interpolations-- this one is an unknown quantity that enters into the formulation, that one is an unknown quantity that enters into the formulation-- if you recognize that, then you have to ask yourself, what do we do with these unknown quantities? Well, we eliminate these unknown quantities by expressing them in terms of the nodal point degrees of freedom, these two nodes here and these two nodes. And thus we eliminate these tensor component intensities. And in the K matrix and the F vector, we end up having then only the nodal point displacement degrees of freedom. And that's the key part.

For the ft transfer shear strain tensor components, we use these two interpolations. Now we notice that we have a linear description in r and constant in s for both of these two components. We proceed, of course, to eliminate this intensity and that intensity in terms of the nodal point degrees of freedom. Same way as we are proceeding here.

This element, the MITC4 element, has only 6 zero eigenvalues, no spurious zero energy modes. It passes a patch test. And this is a very important point-- that the element process the patch test. What do we mean by the patch test? Well, the idea in the patch test is that any arbitrary patch of elements should be able to represent constant stress conditions.

Let's see how we perform that patch test. Here we take an arbitrary patch of elements. Some of the elements in that patch of elements would be geometrically distorted. And we subject this patch to the minimum displacement rotation boundary conditions to eliminate the physical rigid body molds, and then constant boundary tractions corresponding to the constant stress condition that is tested. These two items we apply to the patch of elements, and then we calculate all nodal point displacement and element stresses. The patch test is passed if the calculated element internal stresses and nodal point displacements are correct.

In other words, let's go through this thought once again. What we're doing, really, is we're looking at this structure. In essence, we are taking out of this big mesh of elements in the structure a certain set of elements. We call that set of elements a patch of elements. We apply to this patch of elements the minimum boundary conditions to eliminate the rigid body molds and traction along the boundary that should result in constant stress conditions within the elements. If in fact, our solution for that patch of elements gives us constant stress conditions and nodal point displacement that correspond to these constant stress conditions, then the patch test is passed. and this means that we would have an element that ultimately, as the mesh is made finer and finer with that element, the solution will converge to the correct solution. That's the patch test, and the most important test, particularly for shell elements, when we have these complications of various kinds of interpolations in the element formulation, and so on.

Well, here, schematically, I show such patch of elements. x1-axis, x2-axes, flat patch of elements. Notice the elements not distorted, not rectangular elements. We take certain material properties, assign a certain thickness. We make it a thin patch of elements. The plate that we're looking at here is, in essence, thin. Notice thickness 0.01, width 10 in both directions. And so if the elements would lock do to shear or membrane locking, the way we discussed earlier, we would certainly see it by subjecting this patch of elements to bending moment boundary conditions.

The membrane test would correspond to these stress conditions. The bending twisting test would correspond to these stress conditions. These are the externally applied forces, tractions that the patch of element is subjected for bending twisting tests, here for membrane tests. And once again, we adjust the rigid body modes from the patch of elements by setting the appropriate number of degrees of freedom, nodal point degrees of freedom, equal to zero, and then subject that patch of elements to the tractions corresponding to these stress conditions, the six ones that I'm showing here. And we measure the internal stresses, and find whether these internal stresses indeed correspond to the constant stress situation.

Also, since only those nodal point degrees of freedom have been removed on the boundary that take out the rigid body modes, we should also see that the other nodal point degrees of freedom, of course, have taken on the exact analytical value corresponding to the constant stress condition that is tested. This one, that one, that one, or the other three here. That's the patch test.

Well, let us now to look at some examples regarding the elements that we have discussed in the previous lecture and this lecture. And the first example that I like to look at with you is the analysis of a spherical shell with a hole. There's a hole here, and you see the sphere here. Here's the backside of the sphere, shown by a dashed line. And this sphere is subjected to concentrated forces, shown by these arrows. Radius given here, thickness given here, and the material property of the sphere given right there.

The first step is to select director vectors. We talked in the previous lecture very heavily about the fact that we use director vectors, and we have to define these director vectors. The initial director vectors, or the director vectors for the initial configuration. And we in fact, in this particular case, can quite easily generate the director vectors for each node. The director vector for each node, in this particular case, is chosen to be parallel to the radial vector. So if we look, this is the skin of the shell from the midpoint out via the radial vector, we would see the director vector right here in two dimensions, of course, x and y.

So we generate these director vectors. In ADINA, that can be done automatically. And then we would select the boundary, displacement boundary conditions. If we look at what happens, for example, in the zx phase, using symmetry conditions, of course we would only look at a particular part of the shell. The part of the shell, in fact, maybe I should show you that one first, that we want to look at is shown here. Notice here we have the part of the shell that we analyze because of symmetry conditions. We have here a symmetry boundary and there a symmetry boundary. Notice that we have x, y, and z, the coordinate system used here, and that this part of the shell really lies in the part of the coordinate system that we're looking at. In other words, zy, y here, z up there, xz, x here, z comes out there. So here we have y equal to 0 on this face, and on this face here, we have x equal to 0.

And what I'd like to now talk briefly about is about the selection of the boundary conditions along these two faces. Well, these boundary conditions, if we look at the zx face, in other words, this face here, we see on this view graph here the x face-- here is the shell. Here we would have our director vector at time 0, and at time t, that director vector would have moved as indicated by the red here.

Notice that with this movement, the only rotation that we can have is the theta y rotation. And that is now shown on the next view graph. Notice theta y here is free, a free rotation. Theta x must be 0. Theta z must be 0. And of course uy must be 0, because these are material particles, and in particular that one there is a material particle on the zx phase. So uy is 0. uz and ux are have free. The free degrees of freedom, once again, uz, ux, and theta y.

Well, to impose, then, the boundary condition, we would proceed as shown for this material particle along all the nodes on this face, and similarly, of course, on the other face-- the face on the zy. To prevent the rigid body translations, we would also have to take 1 degree of freedom out, namely one z degree of freedom out, one displacement degree of freedom corresponding to the z-direction out. And that then would basically give us all the boundary conditions that we have to apply for this problem.

So let's look at this view graph here to show the mesh that we are Using. And we show here that at this node, we set uz equal to 0 to prevent the rigid body motion into this direction. The linear elastic analysis results for this problem are given here. Notice displacement at the point of load application is 0.0936. The analytical value is shown here. Pictorially, the original mesh is shown in black, and the deformed mesh is shown in red.

The important point of this problem, and I'd like to go back to it once more, is really that when we look at this mesh, that we define the director vectors at each of these nodes-- and that can be done automatically, they can be generated-- that at each all of these internal nodes, we use 5 degrees of freedom, we have two locally aligned alpha and beta degrees of freedom, the one that I talked about in the earlier lecture, at each of these internal nodes. I like to distinguish the internal nodes here from the boundary nodes. On the boundary nodes, we of course have also only 5 natural degrees of freedom, but we assign now 6 degrees of freedom there. Meaning that the alpha and beta degrees of freedom are rotated into the three Cartesian coordinate axes so as to have three rotations at each of these nodes on the symmetry faces. And then we can very easily impose the symmetry boundary condition that I just talked about. That's the important point of this problem, really-- how we model the boundary conditions, and what degrees of freedom are assigned to each of the nodes.

Well as a second example, I'd like to consider with you the analysis of an open box, a 5-sided open box. Here we show the box upside down, or open site down. The open side, in other words, is down here. And the box lies on a rigid, frictionless surface. On top here, we have pressure applied to the face. The box is modelled using shell elements.

The point of this problem is to discuss with you how we are modelling this box using the concept of director vectors, 5 degrees of freedom, 6 degrees of freedom, that we discussed in the earlier lecture. So what we need to do here is choose initial director vectors, choose 5 or 6 degrees of freedom for each node, and choose the appropriate boundary conditions. The way we proceed there is as follows. We recognize that we could deal with director vectors at every one of the nodes, but if we do so, at a corner, you would have no unique director vectors. You could, of course, assign there a mean normal. In other words, if you look closely here, my hands show the corner, and the pointer shows the mean normal at that corner. And that mean normal, mean in quotes, would be the direction of the director vector. Let me show it from up like that. Here you have the right angle between my hands, and you have the pointer giving you the director vector.

That would be one possibility. However, if you use ADINA, then, in this particular case, it is more effective to not input director vectors. Because if you don't do so with ADINA, then for each node, ADINA will generate automatically a mid-surface normal vector. In other words, if no director vector has input for a node, then ADINA generates for each element connected to the node a nodal point mid-surface this normal vector at that node from the element geometry.

Now, what does this mean? It means that if you use this option, there will then be as many different nodal point mid-surface normal vectors at the node as there are elements connected to the node. Now we will have to look at this statement much more closely for this particular example.

In this particular case, using the option of automatic generation of element nodal point mid-surface normal vectors, we find that at a node away from an edge, we will end up with just one mid-surface normal vector. Because see, this is a flat surface, and the program will calculate for this element and that element and that element and this element, corresponding to that node, the same normal vector. And that is the normal vector that the program will then use.

However, if we look at an edge, along an edge, and pick a node there, then at this node, there will be 1 normal vector calculated for this element and that element at that node, shown by this red arrow, and one normal vector at this node calculated for this element and that element, shown by this red arrow. Of course here, you would have three normal vectors. One on that phase, that phase, and this phase.

Well, with these normal vectors, then, given, you recognize that there will be no stiffness corresponding to the rotation about this vector, because these elements don't carry stiffness in this rotational degree of freedom. We have to be careful with that, and we have got to make sure that we solve the problem with this rotation, the rotation corresponding to this vector here or about this vector, being deleted in the model. Whereas at this node here, we assign 6 degrees of freedom to the model, because the rotation corresponding around this vector will have stiffness, around this vector will have stiffness, and if you think of a rotation in the direction of that vector but at this node, about this direction, you will have also stiffness. So we need 6 degrees of freedom being assigned at a node that corresponds to an edge or corner. We have 5 degrees of freedom to be assigned to all other nodes.

This is once more summarized here. 5 degrees of freedom at all the nodes inside a face. 6 degrees of freedom along an edge, and of course also at a corner.

Let us look now overall at the boundary conditions that we have to assign to the box in order to be able to solve the problem. Here we show now the box open side up. So this is really the side on which the box is placed, but we just switch the box around so that we can more easily look at what happens on this boundary here.

Well, we have nodes that we might call internal nodes along here, nodes that are not corner nodes. And these nodes would have the following boundary condition. First of all, all the rotational degrees of freedom shown by this arrow, this double arrow, must have been deleted. This I just discussed. It is due to the fact that there is no stiffness corresponding to this rotation at that node, because all the elements coming into that node don't have stiffness corresponding to this rotational degree of freedom.

This degree of freedom is free. That one is free. This rotation is free. That rotation is free. This translation is also deleted. That is due to the fact that the box can't move vertically. It's placed on the rigid surface.

If we look at a node that is a corner node-- and let's pick a more typical one. Here we have one. Then at this node, we only take out, as you can see, the z displacement degree of freedom. Because, once again, the box lies on the rigid flat surface. But notice that these rotations, all three rotations, are left free. They all carry stiffness the way I already discussed it. Also these translations are free. Typical corner node, another typical corner node, the same action as this corner node.

Now here we see something in addition. We see that this degree of freedom is taken out, the same way as over there. But we also see that this degree of freedom is taken out. This one is taken out in order to prevent a rigid body rotation about my pointer here. That's why we take this one out. And we take this one out and that one out to also prevent the rigid body rotation that I just mentioned, and to prevent the rigid body translation. See, this translation has to be taken out as well, and that translation has to be taken out, and this rotation has to be taken out. I repeat, this translation, that translation, and that rotation, all of those have to be taken out. And that is achieved by knocking out this degree of freedom, that one, and that one.

So these are the boundary conditions that we apply to the nodes on the edges that are lying on the rigid surface. If we perform now a linear elastic static analysis of the box with a uniform pressure applied to the top off the box, using our 4-node MITC4 element, here is a typical element shown, this is [? the measure ?] that we use, we find that the deformations are magnified, quite highly magnified, looking like that. Very reasonable deformation. We could not compare here with any other solution. But it's still a very interesting demonstrative example, because of the way we are dealing with the director vectors, which actually we don't define. We let the program generate mid-surface nodal point vectors. And because we could discuss the use of the shell element degrees of freedom, freedom is the imposition of boundary conditions, and other aspects.

This, then, brings me to the end of what I wanted to discuss this you related to the view graph, and I'd like to now share some further experiences with you that I documented on the slides. So let me walk over here and let us discuss the information that is even on the slides.

The first slide here relates to the analysis of a simply supported plate under uniform pressure. This plate is modelled as shown here. We only model 1/4 of the plate because of symmetry conditions. And as you can see here, n is equal to 2 for 4 elements. We analyze this plate and we calculate the center displacement, which means the displacement at that point, and compare that center displacement to the analytical solution.

Notice that in this particular case, l/H is equal to 1000, so it's a very thin plate. And certainly we would get very bad results if we had a shear locking phenomenon in the element that we're using. As I pointed out earlier, the MITC4 element does not have any shear locking. And as you can see now here, as we increase N, as shown in this axes, we get excellent convergence behavior to the analytical solution.

Now the next slide shows the analysis of a very similar plate. Again, a square plate, l/H still being 1,000. Simply supported plate. But now the plate is subjected to a concentrated load. And in this particular case, we see the solution results obtained with the MITC4 element as a function of N-- same number N that I referred to in the earlier slide. We see that we get also very good convergence behavior in this particular case. Not quite as good as for the uniformly distributed pressure case, but still very good as well here.

The next slide shows the plate now clamped at its edges. Not any more simply supported. l/H still 1,000, so still a very thin plate. And here we once again calculate the center displacement as a function of N, the number of elements used. Remember, N equal 1 means 4 elements. Remember N equals 4 means 16 elements. Once again, we see excellent convergence behavior to the exact analytical solution.

And finally, this was a case for uniform pressure. Now we look at the case for a concentrated load applied at the center. And here we see this convergence behavior. Of course, this is quite a complicated, difficult problem to analyze, because of the concentrated load and the clamp boundary conditions.

Here in this analysis, we used nice square elements. A very important consideration in a practical analysis is to identify how well an element performs when it becomes distorted. And the next slides now relates to that question. Here we now model 1/4 of the plate using mesh 1, shown here, and using mesh 2, shown here. Notice the elements that we're now using in these analyses is still the MITC4 element. Each of these elements is an MITC4 element, but it's a distorted element.

And our objective is really to identify what displacement predictions we obtain at the center, at point c of the plate, when the plate is subject to uniform pressure, and the plate is also simply supported. And what displacement predictions do we obtain for these two measures at these two points? How much do they deteriorate when compared to the usage of nice square elements?

Well, the displacement predictions are given here. Finite element over Kirchhoff theory results. Mesh 1, 0.93, mesh 2, 1.01. Notice mesh 2 has just 1% error. Mesh 1 has 7% error. But with mesh 1, we have a very coarse mesh, and highly distorted elements.

It's of course also of interest to see what moment predictions one obtains, and here we compare the finite element moment versus the Kirchhoff moment at point c. And you can see 15% error for mesh 1, 2% error for mesh 2. Actually also very good results for moment prediction, as well. I should also emphasize, please, that l/H is equal to 1,000, so we are still talking about this very thin plate.

The next slide now shows how the element behaves in the analysis of a skew plate. It's also now a distorted element mesh, or rather, the elements used in the meshing of this plate are distorted elements, but they're a little bit more regular distorted. In fact, they are parallelogram elements. Notice, please, that this is a top view, so we are seeing a plan view of the plate. And--

SPEAKER: Pick up lecture 20, tape B, slide 20-7, time code 20:33:18. Take one.

PROFESSOR: The next analysis that I like to consider now is the analysis of the skew plate. Here we show a top view of the skew plate with an angle of skew of 30 degrees. So it's a highly skewed plate. And this is the kind of mesh we want to use. Here we show the 4-by-4 mesh of MITC4 elements. Notice the elements are distorted, but sort of regularly distorted, because they are having a parallelogram shape.

The boundary conditions on the plate are simply supported boundary conditions, which I should point out are modelled using the isoparametric degenerate plate elements, shell elements, that we talked about. These boundary conditions are modelled by simply setting only the transverse displacement at the nodes, on the boundary equal to 0. In other words, the displacements corresponding to the z-axis, pointing out of the slide, are set equal to 0, but the rotations on the boundary nodes are left free. This is a very interesting point, most interesting point. And if you like to read up on it, please refer to the reference that is given in the study guide related to this slide.

Well, we use these kinds of meshes to analyze the plate. Notice the material data are given here. The thickness of the plate is given here. B is equal to 1, so we have a thickness with ratio of 1/100. It's a fairly thin plate. And we are considering uniform pressure applied to the plate.

The next slide now shows analysis results. 4-by-4 mesh, 8-by-8, 16-by-16, 32-by-32 meshes were used. We compare our finite element solution at the point C, the midpoint of the plate, versus the Morley solution at that point. And you can see that we really have obtained very nice conversions to the Morley solution. Notice here, we're looking at the principal moments, actually the maximum and the minimum moments, at that midpoint of the plate again, point C, and we're comparing finite element solutions with Morley solutions. And here once again, quite good convergence for the moment, as well.

Now in this analysis, we use the kind of mesh, this 4-by-4 mesh, that we had on the previous slide. And one can ask, of course, is that a good mesh to use? Well, we deliberately wanted to use this mesh in this analysis to test the capability of the element. Of course, one could use a different mesh. And on the next slide now, I'm showing the results obtained by using just four elements, a 2-by-2 mesh, to model 1/4 of the plate. And you can see, with this 2-by-2 mesh, we get already excellent predictions of the displacement at the point C. This is where point C is, once again. Excellent prediction of the displacement at point C when compared with Morley solution.

The moment is quite poorly predicted, still. The two moments that I defined earlier. If we go to a finer mesh, 4-by-4 mesh, the moments, of course, are much better, and so is the displacement prediction. But the displacement prediction was already goods, so there is not that much of an improvement. But in general, one can say, this here, certainly for displacement predictions, is a much better mesh to use than the earlier meshes that we used for the earlier table results.

Now I'd like to look at another problem, the analysis of a cantilever subjected to an end-bending moment, and undergoing very large displacement. Here are the data given corresponding to the cantilever, an elastic structure that we're looking at. And we will be measuring v, phi, and u, and u, v, and phi would take on large values. In fact, we will actually normalize u and v to the length of the cantilever.

The next slide now shows the results obtained. v/l, u/l, and phi over 2 pi as a function of the moment applied. We're using here moment parameter, but of course, the moment increases along this axis, and here we are plotting u/l, v/l, and phi over 2 pi. Notice that even for very large displacement, v/l being 0.7, with 1 cubic element, cubic in the direction of the cantilever, excellent results. This is, in other words, the use of the cubic element, in general the 16-node shell element, for these types of analyses. And it shows that the 16-node shell element is really quite effective.

The next slide now shows how we would apply the MITC4 element. Here we use just two elements to model the cantilever, bending moment applied here, clamped boundary conditions here. And you can see that here we are plotting u/l, w/l, phi over 2 pi, as before, moment parameter along here. And you can see that we're getting really excellent results just with 2 elements for w/l 0.6, or certainly up 0.5, very close to the analytical solution. The analytical solution is given by the solid line.

The next slide now shows the result obtained using three MITC4 elements to model the same cantilever. Notice these are now the nodes that we also have used when using the cubic element in our first result. And note this excellent comparison with the analytical solution for u/l, w/l, and phi over 2 pi. Very large displacement under very large displacement conditions.

The next slide now shows the analysis of an I-section clamped at this end and subjected to a torque at this end. This is actually quite a complicated problem if you think deeper about the problem and recognize that you have torsional conditions, restrained warping here, [UNINTELLIGIBLE] torsion here, then of course you have bending in the beam, et cetera, et cetera. It's not a very easy problem if you look deeply into the problem. And we wanted to obtain some relatively good prediction for this structural response when the structure goes into the elastoplastic regime. Here we have the material data used in the analysis. You can see perfectly plastic conditions by using tangent material modulus equal to 0.

The slide here now shows one model that we use, namely the use of altogether 9 isoparametric cubic elements. These are rectangular elements. Notice this is 1 isoparametric element. This is the mid-surface or the neutral axis of that one isoparametric cubic rectangular section element. We are modeling the whole I-beam by assembling the I-beam using 9 such elements. Of course we discussed earlier that these elements contain the proper warping conditions. We have supplemented our isoparametric interpolation by warping functions, and those now of course will be activated in picking up the torsional response, or modelling the torsional response. This is a free end, that's the fixed end. This was one model that we used, and in addition, we used another model, which is shown here-- a model of shell elements, 9-node shell elements, 9 such elements. Of course, a better response would have been obtained by using the 16-node element, but that of course would have been a more expensive analysis.

So in this particular case, we wanted to test this model versus the isoparametric beam model that I just showed you. Notice one interesting point, how we are modelling this connection here. You see a blown-up figure of it. At this point node c, we have just 5 degrees of freedom. No rotational stiffness about this vector here. So that is not a degree of freedom. Here we have 6 degrees of freedom, 3 translational degrees of freedom and 3 rotational degrees of freedom. Here we have 5 degrees of freedom. No degree of freedom, no stiffness corresponding to the rotation about this vector here. That's one of the important points that we discussed earlier already.

Well, to proceed with the analysis, we actually use the longitudinal element dimensions shown here. Notice that we use a shorter element at the fixed end than here at the free end, because we anticipate heavy plasticity here that we wanted to model quite accurately. Once again, the torque is applied here, and theta x is measured here. And on the next slide now, we show the variation of the torque as a function of the theta x [? rate ?] rotation. Notice here you have the isobeam model solution, and here you have the shell model solution. We compare our results with a sand heap solution and a Merchant's upper bound solution. These are simplified analytical solutions to this problem.

Our results here are somewhat close, but not very close. Of course, one could make a much deeper study of this problem, a very interesting problem. But the way we now already have concluded our study and the way I show you the results, see, I think it's already quite valuable, because we have learned a bit regarding the modeling of this type of problem.

Well, the next slide now shows a final problem for which I wanted to show you some solution results. Here we have a cylindrical shell supported on rigid diaphragms at the end, and the shell is subjected to uniform pressure. We measure here wB, the displacement at that particular point B. The data corresponding to the shell are given here. Notice that we model the shell as an elastic, perfectly plastic material, and notice that in the shell, we modeled only this part of the shell, this quarter part of the shell, using a 9-by-9 uniform mesh of the 4-node MITC4 element. The results are now given on the slide shown here. Pressure applied plotted along this axis, wB plotted along this axis, and this is a non-linear solution response we have obtained using the MITC4 element, and we compare our solution with the solution given by Krakeland.

I might point out that in this analysis, we used our automatic load-stepping scheme that we discussed when we talked about the solution of the non-linear finite element equations in an earlier lecture. Very good comparison up to this point with Krakeland's solution. We did not go much further than that solution given, because we could not compare with any other solution anyway.

Well, this brings me now to the end of what I wanted to discuss with you regarding shell elements. Of course, please remember there is a lot more information that we could have discussed and that we could have studied together, but we only have two lectures, and this is the information that I wanted to share with you in these two lectures. I hope you enjoyed it. Please also refer to the papers that are given as references for these two lectures, particularly for the slides, where you will find much more information, many more details regarding this material that we discussed. Thank you very much for your attention.