Lecture 2: Basic Considerations in Nonlinear Analysis

Related Resources

Study Guide (PDF)

Readings

Section 6.1

Examples

Problems 6.1-6.4

References

Sussman, T., and K. J. Bathe. “Studies of Finite Element Procedures: On Mesh Selection.” Computers & Structures 21 (1985): 257-264.

Sussman, T., and K. J. Bathe. “Studies of Finite Element Procedures: Stress Band Plots and the Evaluation of Finite Element Meshes.” Engineering Computations 3 (1986): 178-191.

• Download English-US transcript (PDF)

The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high-quality educational resources for free. To make a donation or view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

PROFESSOR: Ladies and gentlemen, welcome to this lecture on non-linear finite element analysis. The title of this lecture is "Basic Considerations in Non-Linear Analysis." In this lecture I would like to discuss with you first the application of the principle of virtual work. We will spend quite a bit of time on the basic principle of virtual work, because it's such an important ingredient for non-linear finite element analysis. In fact, the displacement-based finite element method, which is widely used for non-linear analysis, is based on an application of the principle of virtual work.

I also would like to discuss with you and particularly emphasize some basic requirements of mechanics. There are a few very important requirements of mechanics that it is important to always keep in mind, and I'd like to share some experiences with you regarding these requirements.

And finally, to emphasize all the points that I have tried to make in the lecture, I like to present and share with you some experiences regarding some actual analyses. And I've chosen two analyses. First we will look at the analysis of a plate with a hole, and then we will look at the analysis of a plate with a crack.

Let me now go over to my view graph and start with a discussion of the principle of virtual work. The principle of virtual work, in essence, says that this integral here on the left-hand side is equal to what is on the right-hand side. And we should, of course, discuss now in detail what we have on the left- and on the right-hand side. On the right-hand side, we have this tR-- this is a script R, an unusual symbol, but some time ago, some long time ago, I decided to select this symbol for this quantity.

This script t, superscript t R, contains two integrals. The first integral is an integral over the volume tv. Notice it's the current volume of the body of tfiB, the body forces at time t multiplied by the virtual displacement, del ui. This product here is integrated over the volume at time t. We add to this integral another integral, which is an integral over the surface area of the body. And here we have surface forces tfiS per unit area, and these are virtual displacements, delta uiS. We integrate this product over the total area of the body. So this defines the right-hand side.

On the left-hand side, we have one quantity here, which is t tau ij, t once again referring, of course, to time. tau is a stress, and ij being the components of the stress. This stress is called the Cauchy stress. It's a force per unit area at time t. And on the left-hand side, we have also a quantity which is a virtual strain. Delta here meaning virtual, teij meaning the strain. And it's defined as shown here.

Now let us look briefly at what these quantities are. We have a partial of del ui with respect to the current coordinate txj, and we have a partial of del uj with respect to the current coordinate txi. Here I have prepared a view graph which separates out, once more, these important quantities. And on the top here in red, we see the virtual displacement and the virtual strain. Notice that the virtual strain here carries a subscript t, a left subscript t. This left subscript t means that it's a strain referred to the current configuration. del always meaning virtual, of course. And it's important to recognize the following. It's important to recognize that these are the virtual strains corresponding to these virtual displacements. That's the most important point to recognize.

In green here, we have some other quantities-- the current volume, the current surface area of the body. Current meaning always at time t. And here in blue, we see the externally applied forces. tfiB being the body forces, the volume forces, and tfiS being the surface forces per unit volume, per unit surface. These are externally applied forces.

So we have here quite different quantities-- namely, virtual quantities and real quantities. So let me go back to this view graph to discuss with you in more detail what the principle really tells us. It tells that the internal virtual work consisting of the product of the actual stresses-- t tau ij, which in general analysis are, of course, unknown-- times the virtual strains, which we impose, this internal virtual work-- there's an integration over the whole volume-- must be equal to the external virtual work. And the external virtual work, once again, is given here.

This equation must hold-- and this is an important point-- for arbitrary variations in displacement, or arbitrary virtual displacements that satisfy the geometric boundary conditions, that satisfy the essential boundary conditions, or you might say, the displacement boundary conditions. In other words, this is an equation that must hold for any virtual displacements and corresponding virtual strains. Remember, I just said that for virtual displacement, we always have corresponding virtual strains. This equation must always hold, provided the displacements satisfy the displacement boundary conditions.

Here we have, once again, the unknown forces, or stresses, I should say, forces per unit area. Notice these are Cauchy stresses. Here we have the word Cauchy. In an infinitesimal displacement analysis, we really only talk about one kind of stress-- the engineering stress. In large displacement analysis, large strain analysis, we have different kind of stress measures.

And the stress that I now talk about here is the Cauchy stress, which is the force per unit area at time t, which is, in other words, an actual physical stress. Later on, we will introduce another kind of stress, a very important stress measure, which, however, is not quite physical. This one here is a physical stress, the force per unit area. And that is the stress that you, for example, in a computer program, would like to get printed out. This is the force per unit area that you would design your structure with.

This here is, once again, the virtual strain. And notice that this strain is referred to the current geometry, the coordinate at time t. Therefore, stress at time t times virtual strain at time t-- that's what we're looking at here, and that gives us the total internal virtual work when we integrate over the original-- over the current volume, current, I should say, very important-- the current volume of the body. So if we have a body that moves through space and undergoes large displacements, large rotations, large strains, then this here is the quantity of the internal virtual work with the current stresses at time t in the configuration at time t.

Let us look at a schematic example, because there are so many important considerations. Here we have, in a three-dimensional Cartesian coordinate frame, denoted with x1, x2, x3 as coordinate axes, a body in which we are showing two material particles. Notice one of them in the body itself, in the volume of the body, rather, and one on the surface. The body is here supported, and this is the configuration at time 0.

Now this body moves through a certain amount. It undergoes large displacements, large rotations, large strains, and takes on the red configuration here. Notice that in that red configuration, we have certain body forces, which we already talked about earlier, and surface forces, that we also talked about earlier.

Notice that the particles have now moved to new configurations, of course, to new positions. The surface particle has moved there, and the particle within the volume has moved there. We denote this as tui, the displacement at time t, or the displacement from configuration 02, configuration t, for this particular particle. Similarly, here we denote this displacement as tuiS.

Now the principle of virtual work says that we will be imposing on to this red configuration the configuration at time t, a variation in displacement, which is the virtual displacement. And here we see now one such virtual displacement. Notice that it is shown by the blue line for the whole body, and it's a virtual displacement from the configuration at time t. Notice that the particle, which was originally here, or which was at time t here, I should more precisely say, has moved to this point in the virtual displacement. And similarly, the particle that was here at time t has moved to this configuration, to this point here.

Now, this is one variation. However, the principle of virtual work, of course, states that the left-hand side must be equal to the right-hand side-- I think you know now what I mean by left- and right-hand sides-- for any variation. So we have here another variation. All such variations, however, must satisfy the displacement boundary conditions, and those displacement boundary conditions are given right down here.

Now notice that the particle, which I looked at earlier already, has now moved to this configuration, or in this variation of displacement. Similarly, this particle has moved here. Now the principle of virtual work states once again that for any variation that is satisfying the displacement boundary conditions, the left-hand side, the internal virtual work, must be equal to the right-hand side, the external virtual work.

And these are two variations that we looked at. Of course, there are many, many more. In fact, we could think of millions, billions such variations. And we will, in finite element analysis, however, discretize the whole structure using finite elements, and then only allow certain sets of variations. Namely, those contained in the interpolation functions of the finite elements.

Well, let us now look at one important point. Namely that if we were to integrate the principle of virtual work by parts, we would obtain the governing differential equations of motion plus the force or natural boundary conditions, just like an infinitesimal displacement analysis. However, the difference, of course, is that we have referred all our static and kinematic variables now to time t, to the current configuration of the body.

Let's look now at an example. It's a simple example, in some sense, but it still illustrates many of the basic points that I just mentioned. The example is a truss stretching under its own weight. Here we have the truss in its original configuration shown in black. It has an original length L0. It's subjected, of course, to gravity, g. A typical cross-sectional area is given here, 0A. Notice that this bar, under its own weight, moves down to the configuration at time t here. tA is this cross-sectional area now. In other words, this area has now moved to that area. The length of the truss, of course, has changed, and in fact has changed to tL.

In order to analyze this problem, we have to make certain assumptions. The first assumption is, plane cross-sections remain plane. Second assumption-- constant uni-axial stress on each cross-section. And then if we make these two assumptions, then we have a one-dimensional analysis.

Using these assumptions, we then obtain directly, from the general principle of virtual work, the following. Here we have, on the left-hand side, the general left-hand side of the principle of virtual work, Cauchy stress, the virtual strains, integrated over the current volume at time t, which then becomes, on the right-hand side here, the one stress only, because we have only one stress in this bar, the vertical stress, times the virtual strain corresponding to that motion of the bar, namely, vertically downwards. And this product is integrated over the current area at time t, and of course, the longitudinal coordinate tdx.

Notice that here we are integrating over the current volume, because we have a tA and a tL here. The right-hand side becomes, in this particular case, an integration over the current length of the gravity force, which involves the mass density at time t, in other words, the current mass density. Notice the mass density, of course, changes as the bar extends. g is the gravity constant. And here we have the virtual displacement, area at time t, and of course the differential increment tdx. Once again, we here also integrate over the length Lt at time t. Hence, if we now equate these two quantities, we directly obtain what's shown here.

Now it's important to recognize that in this equation, these here are the virtual strains corresponding to these virtual displacements. And once again, we apply the principle of course, at time t. So stresses at time t, meaning forces per unit area at time t, area at time t, and so on.

This virtual strain, also used to describe this quantity, is given right here. Notice that here we have the partial of del u with respect to tx. tx, the coordinate at time t. And this del u is exactly the del u that we see right there.

Now this principle must hold for any arbitrary variation in displacement, any arbitrary virtual displacement, that satisfy the displacement boundary conditions. And the displacement boundary condition here, of course, is simply that the displacement is 0 at the top of the bar.

If we use, now, integration by parts on this equation that I just discussed with you, we obtain this equation here. And that equation contains really two terms, two parts to it. First, an integral that is taken over the length L of the bar, and then a quantity that corresponds to the length of the position 0 and L, tL. And all of the sum of these two must be equal to 0.

Since the variation del u is arbitrary, except, as I mentioned, right at the top of the bar, we directly obtain that this equation must hold, which is really extracting this term here from the integral. And this, of course, is the governing differential equation. And also, we must have that at the lower part of the bar, this force or natural boundary condition must be satisfied.

Therefore, this is an example which shows, really, how the principle of virtual work contains the basic differential equations of equilibrium. Here it is just one, because we have only one direction for the motion. And also the stress boundary conditions, the natural or force stress boundary conditions.

The finite element application of the principle of virtual work proceeds as follows. Here, once again, virtual strains, virtual displacements in the volume of the body and on the surface of the body. These are the virtual quantities. The actual applied forces, externally applied forces, are these blue quantities. We deliberately chose a different color to emphasize that these are real forces, blue, applied to the body, and these are virtual displacements and virtual strains. The stresses here at time t, of course, are unknown, and those we want to calculate.

In the finite element method, what's shown in this big box, we interpolate displacements, virtual displacements. We calculate from those interpolations, strains, real strains, virtual strains. Via the stress-strain law, we got stresses, and we will see in the later lectures that in fact, what we are arriving at is a force vector tF and an external load vector.

Now, this tF here, shown in green, is a result of this stress t tau ij. And it's really the nodal point force vector that corresponds to these internal element stresses, t tau ij. This R here is a nodal point force vector that corresponds to these externally applied forces. This del UT, on left-hand side and on the right-hand side, is a vector of virtual displacements, nodal point displacements, and which, in the finite element analysis, we then set equal to the identity matrix, invoking the principle of virtual work in turn for each degree of freedom, just like we do in linear elastic analysis. This becomes the identity matrix, that becomes the identity matrix, and the result is, of course, simply that tF must be equal to tR.

The procedure that we are following in this box, of course, is of much interest in non-linear analysis. And what's happening in this box here, we will be talking about in some considerable detail in the later lectures.

But for the moment, let's now assume that the solution at time t is known. Hence the stress at time t, the volume at time t, et cetera, are known, and we want to obtain the solution at time t plus delta t. In other words, we want to move ahead 1 increment in time delta t and establish the solution at that time. The principle of virtual work applied now at time t plus delta t in exactly the same way as we just applied it at time t results in this set of equations. A vector of nodal point forces that is equivalent to the current internal element stresses at time t plus delta t now. And these are the nodal point externally applied forces at time t plus delta t.

Well, of course, we don't know t plus delta tF, and for the solution, we now say that t plus delta tF is equal to tF, which we do know, plus a matrix tK-- we call this as tangent stiffness matrix-- times del U. If we set this left-hand side equal to t plus delta tR, the externally applied loads, we directly obtain this equation. We are putting the unknowns on the left-hand side, and the unknowns, of course, are the incremental displacements. The tangent stiffness matrix is known. It's the matrix corresponding to the configuration at time t. The load vector is known, and tF is assumed to be known.

Now, this is here an out-of-balance load vector, and this set of equations then result into the solution of an increment in nodal point displacements delta U. This increment in nodal point displacement is added to the displacement at time t to obtain the displacement at time t plus delta t.

But-- and here is the big but-- by writing this equation up here, we had to linearize the problem. You see, we are dealing with a large amount of non-linearity that can come from the kinematics, from the material relationships and so on. And we had to linearize. We linearized here about the configuration at time t to obtain this term, and this means that down here in this equation, we only can have an approximation sign. In other words, what we have calculated are not the exact increments in nodal point displacement, but an approximation to it.

More generally, then, we repeat that process. We apply the same process in an iterative way to obtain a very accurate solution. And this is done in the following way. On the right-hand side now, we have a load vector, the one that we discussed already before, minus the nodal point forces corresponding to time t plus delta t in iteration i minus 1. I should say at the beginning of iteration i, or end of iteration i minus 1. This gives us an out-of-balance load vector. We're calculating an increment in displacements, delta Ui. And that increment in displacement is added to the displacements which were known already. At the end of iteration i minus 1 corresponding to time t plus delta t, we knew these displacements. We add the increment that we calculated, and we obtain a better guess, so to say, for the current displacement at time t plus delta t.

Now of course in this iteration here, we need some initial conditions, and the initial conditions are given down here. The t plus delta t F 0, in other words, when i is equal to 1, then we need t plus delta t F 0 is given by tF. The displacement t plus delta t U0, when i is equal to 1, is given by tU. These initial conditions are used in the iteration here, and we notice that this set of equations in the first iteration, when i is equal to 1, reduces to the equations that I already discussed with you with the previous view graph.

Well, if we look at this set of equations, we recognize the following. Nodal point equilibrium is satisfied when the equation t plus delta tR minus t plus delta tF, i minus 1, is equal to 0. When the externally applied loads are equal to the nodal point forces corresponding to the stresses at time t plus delta t, when these two vectors are equal, or rather, when R minus F is equal to 0, then of course we have nodal point equilibrium. Compatibility is satisfied provided a compatible element mesh is used. The stress-strain law enters in the calculation of the tangent stiffness matrix and the nodal point force vector F, which corresponds to the current stresses at time t plus delta t, iteration i minus 1.

Well, we will get back to these considerations just now. But most important is to recognize that we have to calculate F t plus delta ti minus 1 very accurately from the displacements t plus delta t Ui minus 1. The general procedure is depicted on this view graph. Here we have the displacement at time t plus delta t ui minus 1, which give us strains. We can directly, via differentiations of the displacements, get strains, which gives us stresses if we enter here with the appropriate constitutive relation. Of course, we have to use the constitutive relation corresponding to the problem that we want to solve. We get stresses, and the stresses, then, give us this nodal point force vector that I talked about.

Notice that in this relationship here, we have to perform an integration, and that integration is expressed here, where we obtain the stresses at time t plus delta t in iteration i minus 1 by taking the stresses at time t and integrating the stress-strain law times the incremental strains from time t to time t plus delta t, where we use, of course, the strains at time t, and the strains at time t plus delta t corresponding to that iteration. This step is most important, because you have to calculate the F vector correctly so that when you say that this equation is satisfied, that indeed you have obtained the right configuration for the finite element mesh that you're considering. So it is most important to calculate the F correctly. Of course, the tangent stiffness matrix should also be calculated appropriately, and we will talk about that as well later on.

In all this discussion, I should point out we assumed that the nodal point loads are independent of the structural deformations. In other words, they vary only as a function of time. Here's a small problem that schematically shows what we mean. The load R is always acting vertically, never mind how much this beam has displaced. And that load variation is shown here on the graph as a function of time.

We satisfy the basic requirements of mechanics if we perform this finite element analysis appropriately. We satisfy the stress-strain law because we are evaluating the stresses correctly corresponding from the given strains. We also satisfy the compatibility requirements provided we use compatible elements, and of course we satisfy the displacement boundary conditions.

Equilibrium we also satisfy corresponding to the finite element nodal point degrees of freedom. But there is one important point. This means that we only satisfy global equilibrium in general. We satisfy equilibrium locally only if we use a fine enough finite element mesh.

In other words, I'd like to distinguish now between satisfying global equilibrium and local equilibrium. Local equilibrium means that we want to satisfy the differential equations of equilibrium at every point in the mesh. And that is only achieved if we have a large enough finite elements for the given problem-- in other words, if our discretization is fine enough.

A check on whether we are satisfying local equilibrium well enough is to look at the stresses along the boundary and see whether the stress boundary conditions are satisfied. That is one check. Another check is to test whether along boundaries of the elements, there are small stress jumps only. In other words, whether the stresses from one element to another element do not jump too much.

We call that the stress jumping because as we will just talk further, when we calculate the stress at a node for an element, we get, of course, one value. And if at that same node, the stresses are calculated using another element, you would get, in general, another value. This stress difference should not be very large, and of course, if we have a homogeneous material, if we have other conditions satisfied as well, in other words, if typically the exact solution would have a continuous stress variation, then the stress jump should be very small for an accurate solution. And this is a good check, an indication of whether we indeed have obtained a valid solution.

Well, I like to now demonstrate to you some of these considerations by showing you the results of two example analyses. But it is best if we show these examples on a different reel, on a second reel of this lecture. Thank you very much.

We just discussed some issues regarding how equilibrium is satisfied in a displacement-based finite element analysis. We discussed briefly that global equilibrium is always satisfied for any finite element mesh that you have selected. By global equilibrium, we mean equilibrium at the nodal points of the finite element mesh. However, locally, equilibrium is only satisfied if the mesh is fine enough. And an indication to see whether a mesh is fine enough is given, for example, by whether these stressed boundary conditions are well satisfied, and whether there are no stress jumps between elements.

I would like now to share with you experiences relating to two problems. These are actually linear elastic analysis problems that exemplify how we are looking at stress jumps in order to choose an appropriate finite element mesh. The first problem is the analysis of a plate with a whole. The plate has certain geometric dimensions, of course, and it has also material data. The material data are given right there. Notice that we also give material data for an elastoplastic analysis. We will not perform now the elastoplastic analysis. We will consider it later. For the moment, we just need these two material data. The| geometric data for the plate are given here.

Now, this plate with a hole has a doubly symmetric structure. In other words, there are two symmetry lines, as shown here. So we only need to look at 1/4 of the plate in the analysis. Notice that the plate is subjected to a pulling force up here and similarly down here.

The purpose of the analysis is to accurately determine the stresses in the plate, and of course, once again, assuming that the loads are small enough to warrant a linear elastic analysis. The considerations that are valid in linear elastic analysis regarding the choice of a mesh are of course also valid in a non-linear analysis, and we more easily can discuss them now, just using a linear elastic analysis as an example.

Using symmetry, as I pointed out already, we can just focus our attention on 1/4 of the plate. Here is that quarter of the plate. Here is the hole. Notice we have roller boundary conditions here, roller boundary conditions there, and the load is applied up there.

Accuracy considerations in displacement-based finite element analysis are that first of all, of course, compatibility must be satisfied. And that is achieved by using a compatible finite element mesh. The material law must also be satisfied, and that is achieved by just using the appropriate material data. Of course, in non-linear analysis, this item needs a lot more attention, but we are performing now a linear elastic analysis, so this is a very simple condition, and it's satisfied by simply selecting the proper Young's modulus, Poisson's ratio for the model. Equilibrium, once again, is locally only satisfied if we choose a fine enough mesh. Otherwise we only approximate equilibrium locally. And as I pointed out earlier, we can observe the equilibrium error by plotting stress discontinuities.

Here we have the results of a two-element mesh for this problem. Notice this quarter of the plate, idealized using just two 8-node elements. This is the line z equal to 0, and we will later on look at stresses along this line. This is a line y equal to z. We will also look at stresses along this line.

But we will plot the stresses along this line here as a function of y. You should keep that in mind. In other words, the first stress quality that we will see will actually start somewhere here, because we are plotting this stress as a function of y.

The displacements, somewhat amplified after load application, are shown here. Notice the maximum displacement here, uz, is 0.0285 millimeters. The maximum stress is at this point here, 281 MPa, megapascal.

Well, if we plot the stresses along the line z equal to 0, and we're evaluating the stresses, I should also say, at the nodal points of the elements, then we obtain these values here. One value there, one value there, one value there. And once again, we are plotting the stresses along the y-axis, measured along here, distance, and the stress magnitude is measured along here. Tau zz. These are nodal point stresses. And we simply have drawn a smooth curve through these nodal point stresses. The applied stress, far away from the hole, is 100 MPa.

If we plot the stresses along the line y equals z, along this line here, we notice that we have a choice of choosing either the nodal point stresses from the element below it or from the element above it. Well, the nodal point stresses from the element below it are given by the crosses here. The nodal point stresses from the element above it-- and notice, these are the nodal point stresses at these nodes here are given by this point.

I should say, actually, that we are plotting here the principal stress, the maximum principal stress. We only look at one stress. Of course, you could plot such curves for all three stress components, but here we're looking just at sigma 1 being the maximum principal stress.

Now there we can see the stress discontinuity. We see a stress discontinuity right at that point. Let's see once, what does this stress continuity result from? Notice that there is a node here-- and I just like to put it in one color-- there's a node here coming from this element. And there's another node coming in from that element. Of course, this red point corresponds to this mark here. So this node here gives us a stress in the stress graph given by that symbol. This blue point there is a node corresponding to the element with the cross, and there we get also a stress value, and that stress value corresponds to the cross mark on the graph.

So let's go back to the graph once more. Once again, same location, because this blue node lies really on top of the red node. They are put together. And yet although we're looking at the same material particle, we have this stress discontinuity. And this stress discontinuity tells us that the mesh is really not fine enough.

So we go to a finer mesh. And here now we have selected a 64-element mesh. Once again, 8-node elements, as shown here, for that quarter of the plate with the hole. Here is the typical 8-node element. This is the undeformed mesh, and when we apply the loads, we obtain this deformed shape. The maximum stress is now 345 MPa, and the displacement right at that point here is 0.0296 millimeter.

Well, we can now once again plot stresses. We can plot stresses along this line-- or I should rather show it here, along this line, and along that line, same way as we have done before. And let's look at the stresses along the line z equal to 0. Of course, y is the coordinate along which we are plotting the stresses, and you can see now here a stress jump, a stress discontinuity. This is where two elements adjoin each other at one node. Otherwise, we see very little stress discontinuity. In fact, too small to be visible by eye.

If we plot stresses along the line y equal to z-- and here I should be more specific. In fact, we only plot, once again, the maximum principal stress. We find that we get these curves here. Now notice that at that location, we get four distinct stress values. The reason being that at this node, four elements couple into that node, and that gives us four stress values at that node.

If you go further away, we find that there are very few or very small stress discontinuities. So however, since we have still one stress discontinuity there at that one location that I pointed out, we went ahead and analyzed an even finer mesh, 288-element mesh. This mesh is shown now here, 8-node elements, once again, in the undeformed shape, and here in the deformed shape. Notice that the maximum stress is now 337 megapascals and the maximum displacement here is 0.0296 millimeters.

We can again plot stresses along this line and along that diagonal line as a function of y. And if we do so, we find that along the line z equal to 0, we find for tau zz virtually no more stress discontinuity. A nice, smooth stress curve. And along the line y equal to z, for the maximum principal stress, also almost no discontinuity. You can see here a slight discontinuity, but that is really a very small one that in engineering analysis, we generally can live with.

Of course, these stress jump plots are very useful and very helpful in identifying whether a mesh is actually a good mesh. However, we should realize that in order to identify whether a mesh is really good, we would have to plot stress jumps along many lines in the mesh, and that can, of course, provide quite some cost. It also means that we have to interpret a lot of data.

And one perhaps more convenient way to look at whether a mesh is a good mesh is to plot pressure bands. We have been lately using this procedure to identify whether meshes are good, and I'd like to just share that latest experience with you. You plot bands of constant pressure, where the pressure is defined by taking the mean of the stresses tau xx, tau yy, and tau zz. Of course, pressure meaning a negative sign in front here.

The two-element mesh gives us this pressure band plot. Notice that these are bands of 5 MPa difference. In other words, here you have a jump from this point to that point of 5 MPa, and once again, from here to there again. 5 MPa. These are the pressure band plots for the two-element mesh. Here we have it as the same kind of plot for the 64-element mesh.

Let's look at this pressure band plot a little bit more carefully. Here we have, as you can see, again, differences of 5 MPa. This is a non-dark area, unshaded area, provides a band of 5 MPa. That shaded area provides a band of 5 MPa.

Now, if we have a good mesh, then the pressure band should be continuous. In other words, this band here should be continuous, and there should be no pressure band break. Let's look at this one here. It's nice and continuous, however, there's a break. And so there is no continuity in stresses here, which is quite simply displayed by the pressure band picture shown here.

For a continuous stress field, we would have no breaks in the pressure bands. We have selected here 5 MPa. Of course, you could also select 1 MPa, and that would be then a much tighter criteria. You could also select 10 MPa. Of course, then that would be a very loose criterion, and 5 MPa for this particular problem, since you can see, we're getting quite a number of bands along here, is a reasonable way to proceed.

If you go to the 288-element mesh, the pressure band plots look like this. Beautifully smooth. You can see here the smoothness in the pressure bands, no discontinuity at all. And this shows really that this mesh certainly satisfies the criterion of smooth pressure bands, smooth pressure variations, over the whole mesh. And that is a strong indication that the stresses are smooth, and therefore that the mesh is really an acceptable mesh.

We see, therefore, that stress discontinuities are represented by breaks in the pressure bands. As the mesh is refined, of course, these breaks in the pressure bands disappear, or certainly become smoother until finally they disappear. The stress state everywhere in the mesh is represented by one picture, and that is, of course, the beauty of this whole procedure. If you plot stress jumps, stress component jumps, then you have to plot for each component along all the lines where elements join each other, these components of stresses. And there is a voluminous amount of information that is being generated.

Here we have just one picture that basically gives us a lot of information about the smoothness of the stress situation. Of course, the pressure band is plotted, maybe plotted by a computer program. Actual magnitudes of pressure are not given, but we are not really interested in seeing those from these pictures. What we would like to see are simply whether the pressure bands are smooth, and whether, therefore, the mesh is an acceptable mesh.

Let's summarize the results. Here we now performed three analyses of this plate. One two-element mesh was used, a 64-element mesh was used, and a 288-element mesh was used. The number of degrees of freedom are given here. Notice that the number of degrees of freedom goes up very rapidly as the number of elements goes up. The relative cost is given in this column, and you can see here that we have normalized the cost to the cost of the analysis using the 64-element mesh. So the two-element mesh is very cheap, 8% of the cost, and the 288-element mesh solution is quite expensive.

The displacement at the top of the plate is given in this column, and we notice that with a 64-element mesh, we really obtain a very good displacement prediction. In fact, there is no change going to a larger number of elements. And the stress concentration factor, interestingly enough, increases first and then slightly decreases. Notice this is not an infinite plate-- it's a finite plate, and therefore you do not get the value of three, for example, that you would see in an infinite plate, for a hole in an infinite plate.

So the two-element mesh cannot be used for stress predictions is certainly a conclusion. It's not accurate enough. The 64-element mesh gives reasonable results for stresses, and certainly very accurate results for the displacements. The 288-element mesh is probably overrefined for linear elastic stress analysis. However, for other types of analysis, for example, an elastoplastic analysis, if you want to have very accurate results, we might want to use, actually, this fineness of mesh.

One question that one, of course might have is, why did he not use 9-node elements? We used 8-node elements. Here we show now a mesh of 9-node elements. Notice each of these elements now is a 9-node element. Would the solution significantly improve if we had used 9-node elements? In other words, one 9-node element for each 8-node element that was used before?

Well, the answer is no. The solution does not improve very much. If you look at this column here, let's look at the 64 8-node elements solution, in which case we had 416 degrees of freedom. With a 9-node element, we will have 544 degrees of freedom. The displacements at the top are unchanged between the 8-node and 9-node element results. The stress concentration factor is very little effected here. Only in the fourth digit, as a matter of fact. And the stress jump and pressure band plot do not change significantly, either.

Well, the second example that I'd then like to discuss with you is the analysis of a plate with a crack. We will once again look at different meshes for this problem, and we will also discuss one additional aspect-- namely that for certain types of answers that you're looking for, you actually can use sometimes very coarse meshes. Whereas if you're looking for the answer of a much more difficult question, of course you will have to use a much finer mesh.

Here we have a plate with a crack. Notice the plate is subjected to a pull of 100 MPa right, and at the left-hand side, same pull. The crack is lying right there. Notice it's not a a symmetric crack. This is a crack tip A, that's the crack tip B. To model this plate, we have to really model half of it. We have to model half of it, because the crack is not a symmetric crack.

The material data for this problem are given down here. E, Young's modulus, nu, Poisson's ratio, Kc is the fracture toughness. And here we have the fact that the thickness of the plate shall be just 1 centimeter. It's a plane stress condition that we're looking at.

And the question that we want to address here is, will the crack propagate? Well, let us look at some background information. Since we want to perform a linear elastic analysis and we want to apply linear elastic fracture mechanic theory when considering whether the crack will propagate in the plate, we calculate a stress intensity factor, and we have a mode 1 stress condition, so we calculate K1. K1 determines the strength of the 1 over square root r stress singularity at the crack tip, r of course being the radius from the crack tip. And as very well established in linear elastic fracture mechanics, when K1 is greater than Kc, the crack will propagate. Kc is a property of the material.

The computation of K1 is performed by this here. It's equal to the square root out of E times G. G is obtained as the partial differentiation of pi with respect to a, where pi is a total potential energy of the structure, and a is the area of the crack surface. G is known as the energy release rate for the crack.

This is quite well established, and what we want to do, of course, in the finite element analysis, is to calculate this G value. The question is, how do we do that? Well, in this finite element analysis, each crack is represented by a node. And if you look at a cut through the plate where this is here the crack, we have this here as the node at the crack tip. This node in the finite element analysis is, so to say, moved forward a differential amount, so that the change in crack area is given via this little bit here. Now that the original crack tip location is here, the new crack tip location is right there.

Then we can, in this process, evaluate del pi over del a as del pi del l times 1/t, where t is a thickness of that plate. Of course, this here is only conceptually being done, that we're actually moving the crack. What we're really doing is the differentiation of del pi with respect to l, where l is the length of the crack and del l is the differentiation of l into this direction here.

The quantities del pi del l may be efficiently calculated by an algorithm that we have just recently been publishing. Of course, there are different algorithms available as well, but we would like to use this procedure to calculate data pi over l. And you can look up how it's being done in this particular paper.

Let's look at the results. Finite element analyses have been performed using a 17-element mesh, and this is the mesh here. We call this a coarse mesh. Notice the crack is right there. Notice that we have only two elements above the crack, to the boundary of the plate, and we have here once again also only two elements below it.

Notice that these nodes around the crack tip B have been moved to the quarter point. And similarly, they've been moved here to the quarter point, meaning that this node is not at the midpoint of this element side, but at the quarter point of that element side. And similarly for the other nodes around the crack tips.

If we look at the results, the stress results, maybe plus the stresses on the line of symmetry, we find that tau yy varies as shown here, then of course increases very rapidly near the crack. This is where the crack is. Increases very rapidly, then drops down to very close to 0. Of course, it should be 0 right in here. And here we have, coming from the other side to crack tip A, this very rapid variation in stresses.

These are the stresses tau yy. If you look at the pressure band plot-- I introduced you to those kinds of plots in the earlier example-- we find that the pressure jumps are really larger than five MPa. For example, you can see it right here, you can see it right there, you can see it right there.

So it's a coarse mesh, and certainly the stresses are not very well represented, not very well picked up with this mesh. However, if we look at the stress intensity factors that we wanted to estimate in this analysis, because they will tell us whether the crack will propagate or not, we find that for the stress intensity factors, for example, KA, we get an excellent solution. For KB, the solution is not as good, but perhaps still acceptable.

We now go to a finer mesh, a much finer mesh. 128 elements are used in this mesh. And we can see once again here, 1/2 the plate. Here is point A, crack tip A. Here is point B, crack tip B. Notice the crack lies in between. And we will be plotting stresses along this line. This is a line of symmetry here, and we will plot tau yy, which is the stress across here, as a function of z.

But let's first look a little bit closer at the mesh. If we look at one detail closer to the crack tip A, we see this is the mesh that we use, and we can go still further and show the elements just around the crack tip, right A. Here is A. Notice that we use triangular elements, these are isoparametric degenerate elements, and that we have moved, once again, the nodal points from the midpoint to the quarter point. Because we know that if we do so, we are generating in the finite element the 1 over square root r singularity that we would like to have in the analysis, because we know that the stresses vary as a function of 1 over square root r. So it's of benefit to shift the nodes to the quarter points, as was done in this example.

Now we look here at the stresses. Stresses tau yy, as I pointed out, is a function of z. And you can see here a stress that very rapidly grows, then on the other side comes down to practically 0 in the cracked area, and on the other side, coming again now from the right here, a stress that varies very rapidly up close to tip A.

We also can see here that the stress variations are nice and smooth. Of course, that indicates to us that the mesh is really a good mesh, but we can also look at our pressure bands again. We find that the pressure jumps are quite small if we go not too close to the crack tip. In other words, here we have some pressure band variation, some pressure band breaks, but we have still quite smooth bands that tell us that the mesh is quite fine.

If we go very close to the crack tip A, then we find that the pressure jumps are, of course, large right here. You can see here breaks in the pressure bend. Right here, for example, a typical example, if you look very closely here, there's a black band coming into a white band. So there's clearly a break in the pressure band.

So here, then, in this area, the stresses are not continuous, and would also not be very accurately predicted. However, we know that our scheme of calculating the stress intensity factor works for even a coarser mesh, and for this fine mesh, we've got indeed excellent results, as you can see for KA and for KB. Both are very well predicted.

The point of this analysis is really that the degree off refinement needed in a mesh depends on what kind of question you have, what kind of problem you're solving, and what kind of question you're asking. If you're only interested in predicting displacements accurately, a coarse mesh may do. If you want to calculate stress intensity factors, also a coarse mesh may do. If you want to calculate lowest natural frequencies, associate mode shapes in a dynamic analysis, once again, a coarse mesh might very well do.

However, if you want to calculate stresses very accurately, then you may need a fine mesh. And in general non-linear analysis, we need accurate stresses, because the stresses influence the material relationship that has to be used as we go through the incremental solution. And therefore, frequently we need, in non-linear analysis, a finer mesh than what would do well in a linear analysis. So we should keep that in mind, that in fact what might be a sufficiently fine mesh in linear elastic analysis may not be sufficiently fine in a non-linear elastoplastic geometrically non-linear analysis.

Of course, how fine the mesh has to be in a general non-linear analysis depends very much on many criteria, criteria that we will still discuss in the next lectures of this course. But I hope this gave you a bit of an overview of what are some important considerations in linear as well as non-linear analysis. Thank you for your attention.