Lecture 14: Solution of Nonlinear Dynamic Response II

Related Resources

Study Guide (PDF - 1.1MB)

Readings

Sections 9.1, 9.5

References

Bathe, K. J. “Finite Element Formulation, Modeling, and Solution of Nonlinear Dynamic Problems.” In Parter, S. V. Numerical Methods for Partial Differential Equations: Proceedings of an Advanced Seminar. Mathematics Research Center, University of Wisconsin-Madison. New York, NY: Academic Press, 1979. ISBN: 978-0125460507.

Bathe, K. J., and S. Gracewski. “On Nonlinear Dynamic Analysis Using Substructuring and Mode Superposition.” Computers & Structures 13 (October-December 1981): 699-707.

Ishizaki, T., and K. J. Bathe. “On Finite Element Large Displacement and Elastic-Plastic Dynamic Analysis of Shell Structures.” Computers & Structures 12 (September 1980): 309-318.

• 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 ocw.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 the methods of solution for the finite element equations that we are encountering in non-linear dynamic analysis. In the previous lecture, I pointed out that the techniques used for the solution of these finite element equations in non-linear dynamic analysis can be classified as data integration techniques, methods of multiple position, and substructuring. We talked in the previous lecture about these types of methods, and in particular, we discussed in detail an explicit integration method, namely the central difference method, and an implicit integration method, namely the trapezoidal rule. We talked about how these methods are applied, how time steps might be selected, what the restrictions of these methods are, and so on. I'd like to now in this lecture discuss with you the method of multiple position, and also substructuring technique.

After the discussion off these two techniques, I'd like to talk with you about some examples. As a first example, we want to consider the wave propagation in a rod. As a second example, I'd like to share with you experiences regarding the response of a 3 degree of freedom system. As a third example, then, we want to consider the analysis of a 10-story tapered tower. Example 3 and example 1 are really solutions of a linear system, but we learn quite a bit, even regarding non-linear analysis, from these examples.

As the fourth example, then, I'd like to share with you experiences obtained in the analysis of a pendulum that goes through very large displacement. As a fifth example, then, we consider the pipe whip response solution. This is of much practical interest. And after these examples, which I've prepared view graphs for, I'd like to consider with you some slides. On these slides, we will have the analysis of a control rod drive housing with the response of a spherical cap, a highly non-linear problem, large displacement, elastoplastic response, and then the analysis of a fluid structure interaction problem, a pipe test in which we are comparing our computer's results with experimental results.

I should point out right now that the details regarding these example solutions are given in some papers that I had also earlier. And please look at the study guide, where the references to these papers are given.

Let me now walk over here to the view graph that we have prepared regarding this lecture. And the first method of solution that I wanted to discuss with you is the method of multiple position. Quite a number of you probably have used this technique already in linear analysis. In fact, the word "superposition" somehow implies that you are dealing with a linear analysis. However, we can extend that method directly to non-linear analysis if we recognize that the mode shape vectors of the system at a particular time can be used as generalized displacements, or basis vectors, to express the response to any other time, as well. The method is effective when, in non-linear analysis, the response lies in only a few vibration modes, just the same way as we are applying multiple position also in linear analysis, namely, it's only effective really when a few vibration modes capture the total dynamic response. The same holds here, of course, in non-linear analysis. And when the system has possibly only local non-linearities.

Let's look at is equations that we're using in this multiple position technique. The basic starting point are the equations that we have seen before regarding dynamic analysis. But notice in this particular equation, I have assumed now that the damping matrix is not present, so we have no damping in the system coming from the damping matrix. There could, of course, be damping due to the material non-linearities, which are captured in the K and the F matrix.

Well, this is the set of equations that we want to solve. And usually, we just iterate on these equations until the right-hand side is 0, where of course we also would bring this one here onto the right-hand side. So think of this one here as a minus on the right-hand side, and if, then, this right-hand side, with M on that right-hand side is 0, meaning the data UK is 0, no more updating to the incremental displacements, then we have established equilibrium configuration corresponding to time t, of course, including inertia effects.

In multiple position, we don't want to directly to deal with these equations in this form. Instead, we want to transform them into a new form, and into a form such that we obtain less unknowns to deal with. Let's assume for the moment that tau is equal to 0 so that we are dealing with the solution of the equations using the initial stress method. We talked about initial stress method in an earlier lecture. Then using this transformation here, where phi R are the eigenvectors of this problem down here, And the xi are the generalized displacements corresponding to these eigenvectors.

Notice this equation is exactly the same that we're used to seeing in linear analysis. Of course, in linear analysis, we are usually not putting this t plus delta t up there. It's understood that the time dependency is in xi, but the space dependency is in phi i. And notice we are summing here from r to s. r would have to be selected by the analyst, and s the same way.

Having introduced this transformation with this eigenvalue problem here, and substituting for the U displacements in the general equilibrium equation, and using the same approach as we are used to in linear analysis, we directly obtain this set of equations. Now, notice here that in this vector x double dot, and in this vector x, delta meaning increment, we have only the components corresponding to the r mode up to the s mode. These are the generalized displacements corresponding to phi r up to phi s. Here we have phi r to phi s listed in this matrix, capital phi. Lower-case phi for these here. These are the columns in this matrix that is denoted as a capital phi.

Notice here we define an omega squared matrix, 0 elements in the off-diagonals. And on the diagonal, we have the frequencies of the system, of the actually linear system, because we use 0K in the eigenvalue solution. These are the entries for this equation up here. And notice that we have now here much less equations than the total number of finite element degrees of freedom. Typically, r might be 1, s might be 20, and then we would have 20 equations here, 20 equations only. Of course, r could also be 25, and s 30, in which case we have here a multiple position corresponding to the 25th up to the 30th mode. In any case, in general, we will have much less equations then finite element degrees of freedom that are listed, for example, in this vector here, R, I should say. The externally applied load corresponding to these degrees of freedom are listed in here.

Of course, the F vector is the vector of nodal point forces corresponding to the element stresses at time t plus delta t, and at the end of iteration, k minus 1. This matrix times that vector gives us a vector which has the same length as this vector has and that vector has. So we are not dealing anymore with all the degrees of freedom that are listed for which we have the forces here in the solution of the equation, but with much less degrees of freedom, the degrees of freedom being equal to the number of mode shapes that are being employed.

An interesting and most important point is that in non-linear analysis, we have a coupling, due to this vector appearing in here, between the different mode shape equations, or modal equations. We cannot solve, as we usually do in linear analysis, one equation of the modes of a position after the other. In linear analysis, you would typically solve the dynamic response equation in mode 1, then the one in mode 2, then the one in mode 3, and like that, you'd sweep over all the generalized degrees of freedom, one after the other, for the total response from time 0 to the end of the response time that you're interested in. In non-linear analysis, we have to go step by step still forward in time, solving all equations, all modal equations at time t, then all modal equations at time t plus delta t, then at time t plus 2 delta t, and so on, because there is this coupling, and because we don't know f unless we have all the information coming from all the mode shapes.

Well, a typical problem that you might want to solve with this approach is the pipe whip problem. Here we have a very long, slender pipe. There is here a restraint that has a gap. And suddenly, this end of the pipe is subjected to a very large force. The pipe undergoes dynamic motion, a wave travels along here, goes plastic, hits this stop, et cetera. It's a very complicated problem, and this approach of multiple position is it cheap way of obtaining a solution to this problem.

Of course, the accuracy of the solution depends on how many mode shapes you include. If you only include a few mode shapes, the solution is quite inexpensive, but the accuracy is not very high. If you include as many mode shapes as there are finite element degrees of freedom, then of course you get the exact response. All you have done is transformed from one basis, the finite element degrees of freedom basis, to the mode shape basis. But then, of course, the solution is very expensive, and you don't want to ever do is that. The point is that you want to perform such analysis if you only need to include a few modes in the response solution. We would actually look at this problem, or such a problem, a little later, when we discuss example analyses.

As a second technique that can also be very useful in non-linear dynamic analysis, I'd like to discuss with you the substructuring technique. This procedure is used is implicit time integration, and the important aspect is that prior to the analysis, we are setting u this effective stiffness matrix corresponding to the linear degrees of freedom only, and we are statically condensing out all the degrees of freedom that will remain linear throughout the response solution.

It is quite effective for local non-linearities. Here, for example, we have a tower structure, and the only non-linearity that this structure sees, say, is the slip down here. This typically would mean that we have a large number of linear degrees of freedom that we can deal with very effectively, using the substructuring approach. The non-linear degrees of freedom here, of course, require iteration and updating of matrices, and those we deal with in the usual way.

Let me show you in more detail what I mean here. Let's look schematically at an example. Here we have a 10-story building, and that 10-story building is discretized, idealized, as shown here. We also want to include the elasticity of the foundation, and that is modeled very simply here as an assemblage of springs.

We deal here with two substructures, shown in red for the substructure degrees of freedom and in black for the master degrees of freedom. This is also part of the master degrees of freedom. Notice that here we have one substructure-- in other words, this part here is nothing else than that part right here. That's what is one substructure. And of course, you see here another substructure. But we have a repetition here-- in other words, this substructure here is identically equal to that substructure. This, of course, is the fact because of the particular design of this building. The upper stories look the same way as the lower stories. That's the reason why we can say, this is one substructure model which can model this part of the structure and that part of the structure, as well.

If we look at the stiffness matrix, or the effective stiffness matrix of this structure that he would encounter in a direct integration solution using implicit time integration, schematically, we see the following. Notice here we have the entries coming from the substructures, and here we have entries from the master structure. This coupling here corresponds also to master degrees of freedom, and similarly, this coupling corresponds to master degrees of freedom.

Here we have the displacement vector, the incremental displacement vector. And let's see here that these master degrees of freedom of course correspond to these element entries here. These master degrees of freedom correspond to these entries, and these master degrees of freedom correspond to these entries.

Now what we can do is condense out all the linear degrees of freedom effects. And that is achieved as follows. We recognize that the total effective stiffness matrix corresponding to time t is made up of this part here, which is the linear stiffness matrix plus the total mass matrix. I should point out here that this K matrix contains all the linear contributions, where that M matrix here corresponds to all the contributions from the masses at the master degrees of freedom and the substructure degrees of freedom. Because we're assuming here that the total mass matrix is linear, we can put all the mass effects right in there.

And then since this is not the total effect of stiffness, we add the effect of the non-linear element stiffnesses, right here. And in this particular analysis, this might be, for example, the non-linear springs down there. The non-linear springs modelling the foundation on the structure.

Well, this can be written in this form. And if we perform static condensation on this combined matrix, we arrive at this matrix here. Notice here we have a condensed degree of freedom vector carrying only the master degrees of freedom, and here we have the entries corresponding to the master degrees of freedom in the tK hat matrix. We put a c on there, meaning condensed. After static condensation, this is a matrix that has been reached. These black squares here correspond to as the black squares that you saw earlier. Of course, they have now been changed a bit because of the static condensation, because the effect of the linear degrees of freedom have been carried into these entries here.

The major steps in the solution are as follows. Prior to the step-by-step solution, we establish the K hat, the linear effective stiffness matrix containing K plus 4 over delta t squared m. And this K hat matrix is then used to statically condense out all internal substructure degrees of freedom to obtain this K hat c, c standing for condensed matrix.

We note that the actual effective stiffness matrix including the non-linear effects is then obtained by taking the linear effective stiffness matrix and adding the non-linear stiffness effects to it. Notice, once again, this matrix is obtained from this matrix here, as I mentioned earlier. We use this matrix to statically condense out all the linear degrees of freedom.

The solution is then obtained in each time step as follows. We update the condensed matrix K hat c-- there should be no t here-- K hat c for non-linearities. That means we're adding the tK matrix to it, corresponding to the non-linearities that have occurred in the system. We establish the complete load vector for all degrees of freedom, and condense out the substructure internal degrees of freedom.

This is important to recognize. We have to calculate the complete load vector because of the intertia effects on the internal substructure degrees of freedom. The intertia forces on the internal substructure degrees of freedom carry forces to the master degrees of freedom that we have to take into account in each time step and in each equilibrium iteration. We then solve for the master degrees of freedom, using the small system size, the tK hat c, the velocities, accelerations, and calculate all the substructural degrees of freedom-- velocities, displacements, and accelerations. These are required here to calculate the load vector corresponding to all degrees of freedom.

Schematically, the solution procedure, therefore, proceeds as follows. For each time step and each iteration, we start with the displacement velocities and accelerations at time t. We calculate the effective load vector, corresponding to time t plus delta t, the way I've been talking about it earlier. And then we condense that load vector, as shown here, to the condensed load vector, c meaning condensation. We then, from this load vector, calculate with the effective tension stiffness matrix, the t plus delta t Uc vector, and this vector is used to calculate the displacement, velocity, and accelerations of all the degrees of freedom at time t plus delta t. Of course here I'm showing only the process for a typical time step. If you iterate, as we usually do in non-linear analysis, of course, then you would have to also introduce the iteration counter here in these vectors.

Notice that using this solution procedure, we get the same results as if we had not performed the substructuring. In other words, there is no additional assumption that we're introducing into the solution. You get the same results. However, the solution cost is reduced, because we are dealing with one substructure, a typical substructure, only once regarding the factorization-- the LDL transpose factorization and static condensation that has to be performed.

Let us now look at some example solutions. And the first example that I'd like to consider with you is the solution of the wave propagation in a rod. Here we have a uniform freely floating rod, and we apply to the left end of that rod a force that varies as shown here as step load of 1000 newtons. The material data for the rod are given here. We want to perform a linear elastic analysis, but although this is a linear analysis, we still will learn quite a bit from it regarding non-linear analysis.

We want to consider the compressive force at a point at the center of the rod, right there. And the exact solution for this force at that point is plotted here. You have, once again, apply 1000 newtons, and here we have the time scale. Notice up to 1/2 t star, there's no force at point a, because the applied force in fact has to propagate through the rod. Then we have the full force up to 3/2 t star, and then no force again, where t star is the time for the stress waves to travel through the complete rod. This is the analytical solution, well-known, well-established, and we want to compare with this solution when we perform our finite element analysis.

For the finite element analysis, we use a mesh off 2-noded truss elements, as shown here, and we measure the force at the midpoint of the total structure by the compressive force in this element. We use in this analysis first essential difference method, and recognize that the time step that we use has to be small or equal to the critical time step, as we discussed in the previous lecture. And that critical time step is given by Le/c. Please refer to as the previous lecture notes, and you would see this formula there, which is nothing else than t star divided by the number of elements. We recall that if delta t is greater than delta t critical, we would produce an unstable solution.

We also need to deal with the initial conditions for this problem. And we recognize that if the initial displacements are 0 and R is given, then we can calculate the initial acceleration, as shown here. Using now a time step equal to the critical time step, we obtain the correct result, and I pointed that out in the earlier lecture already, that we in fact obtain the exact solution for any number of finite elements to model the rod. And here is a finite element solution, compared with the exact solution. A wonderful solution, right on the exact solution.

Using, however, now a smaller time step-- and this is frequently of surprise-- using a smaller time step, we get a much worse result. Here we use a time step of 1/2 delta t critical. Half the time step that we used before. And just solution looks actually very bad right here.

Frequently one sort of intuitively thinks, well, if I use a smaller time step, I should get better results. My costs go up, and surely I should get better results. Well, it turns out that using the explicit time integration scheme for this type of problem, this wave propagation problem, with a smaller time step, more cost, higher cost, you get worse results. That is really something to keep in mind.

Now let us consider the trapezoidal rule. The trapezoidal rule-- we discussed this method in the earlier lecture. And I'd like to refer you once again to the notes corresponding to that lecture, where we have identified that a stable solution would be obtained with any choice of delta t. Remember, this is a linear analysis, and the trapezoidal rule is unconditionally stable.

We could use a consistent or lumped mass matrix. And in this particular case, we use a lumped mass matrix, because we also like to compare our results with those that we obtained using the central difference method.

If we use the trapezoidal rule with the time step that we employ in the central difference method-- in fact, the critical time step that we used there-- and we also use the initial conditions computed as shown here, the way we did it for the central difference method solution, we find that the solution is very inaccurate when compared to the exact solution. Here, of course, we also use the 10-element mesh.

If we use 0 initial conditions for the trapezoidal rule, we get this solution. Not much different, really, from the solution that I just showed you. So whether we use 0 initial conditions or the initial conditions corresponding to the equation MU double dot is equal to R, which we have seen on the previous view graph, makes very little difference when using the trapezoidal rule. It makes a lot of difference using the explicit central difference method. And for that reason, in fact, we have to use the proper initial conditions obtained from MU double dot equals R at time 0 when we use the central difference method.

Using now twice the time step of what we just used, namely, delta t equal to 2 times delta t, critical of the central difference method, this solution is stable, but very inaccurate. It is stable because we're using an unconditionally stable operator. We discussed it in the previous lectures. Note that this is the solution we're obtaining when compared to the exact solution.

Now we look at the solution where we have half the time step off the critical size, corresponding to the central difference method. And you see here, the solution, again, very inaccurate when compared to the exact solution, but of course stable, because the method is always stable for any time step delta t. The question now could, of course, be, and should be, what happens if we take more elements to discretize the whole structure? Here we use the 10-element model. Why not use more elements?

And here we now look at the solution obtained with the essential difference method when we use 100 elements to model the total bar. Here the exact solution in red, and that exact solution is reached if you integrate the equilibrium equations with this time step. We discussed that earlier. That exact solution is always obtained for any number of elements that you use to discretize a structure, provided you use delta t equals delta t critical. If you then use a smaller time step, in fact, delta t equals 1/2 delta t critical was used here, we obtain a solution that is not the exact solution, but as you can see here, it oscillates for this time step this amount about the exact solution.

So we have much the same observations that we have seen already earlier, namely, when the time step is decreased in the central difference method solution-- from the critical time step, here we see a decrease-- the solution that one might expect to be better actually deteriorates. Of course, it couldn't be better than the exact solution. But we see here a deterioration, and a quite significant one.

If you use the trapezoidal rule, this is the solution that you obtain with the time step equal to the critical time step of the central difference method. Again, here, we're using the 100-element mesh. The solution is certainly better than what we have seen for the 10-element mesh, but it does not give us the exact solution-- in other words, a [UNINTELLIGIBLE] wavefront as we would like to obtain.

Of course, this is a very difficult problem to solve, because we have infinitely steep wave fronts here and there, basically, to model. And an overshoot would obviously obtain it, but we would like to have said overshoot and undershoot, with this oscillation here as small as possible, of course, with the solution scheme that we employ.

There is one more important consideration that I'd like to bring to your attention. Here we use a 2-dimensional element mesh to model the bar. Notice these are 4-node elements, plane stress elements here to model the bar. And in this particular case, we use the central difference methods. The thickness is given here, the material data are given here. For this solution, notice that this distance here is half the distance that we have here. And that is something that leads to a phenomenon that one has to be aware of.

Notice our the critical time step that we identified for the one-dimensional wave propagation here would, of course, be calculated depending on this length. However, we have to allow for the fact that there is also a possibility of wave propagation this way, and the actual critical time step that we have to use in this problem is given by this length, by the shortest element distance. Please refer to the previous lecture, in which I discussed this already. You have to use the shortest distance between nodal points within an element to determine the critical time step.

So this relationship here does not hold anymore, because this distance is less than that. And if you were to use delta t equal to t star over 10 elements, which was a critical time step that we used when we solved the truss element model, then we would find that the solution diverges. In fact, in element 5, which was the element that we always used to measure our stress, we find that tau zz, the absolute value of tau zz, is greater than 1000 newtons divided by the area at this time. Notice that it was 1000 newton that we applied. So a stress that should be 0, really, because tau zz should really be 0, it's the stress in the transverse direction, the wave travels this way, and tau zz is acting that way. This stress should be 0. We find, actually, that it grows, and it grows because we are using a time step that is larger than the critical time step that should be used for this problem.

Let us now look at another example, a second example from which we can also learn quite a bit. Here we have a very simple 3 degree of freedom system. These are the degrees of freedom, x1, x2, x3. Notice we have masses m here, and we have a linear spring connecting these masses, a non-linear spring connecting these masses, and a linear spring connecting this mass to the rigid support.

For this problem, the linear spring has is distance here, and the non-linear spring is characterized by this force displacement relationship. Notice there's a kink going from 1 to 100. So the spring gets suddenly, at a particular distance, or rather deformations, I should say-- displacement x2 minus x3-- at a particular deformation, the spring becomes 100 times as stiff as prior to that point.

Well, we can calculate a critical time step based on this spring stiffness here, the linear part, or rather, this spring stiffness for small displacements. And that critical time step for this problem is given as 1.11 seconds. As long as the response of this system is such that we would, in this non-linear spring, never go around this kink, we would deal with a linear system, and this would be the proper critical time step to use.

However, as soon as we are going around this kink here, due to the fact that the forces that I applied, or the initial conditions that we are applying to these degrees of freedom, are large enough, as soon as we go around this kink, our new critical time step is based on this spring stiffness, and if this given right down here. This is a new critical time step. So if we now perform essential difference method solutions of this problem, and if the solution brings us into this regime, it is very important to use, as a critical time step, 0.14 seconds. It is very important, even if we're only in this regime for just a few steps.

You see, what might happen in the response is that most of the time, you're in this regime. But then all of a sudden, for just a few steps, you're getting into this regime, and because of this high stiffness, you're thrown back into this regime. If your time step that you use in the time integration is then larger than that value, you are having an instability for just these few time steps, and that introduces an error in the solution which, if you don't have a comparison with an exact analytical solution, will of course be an error that you not necessarily will see, directly. And that's one of the important points.

You see, in linear analysis, when you use the essential difference method in linear analysis with a time step that is larger than the critical time step, you will, after a few time step integrations, see that the response becomes very large, and the instability is obvious. The instability is really obvious when you apply the central difference method with a time step that is larger than the critical time step to the solution of the system.

However, in non-linear analysis, as I mentioned in the earlier lectures, and you are now well aware of it, in non-linear analysis, the frequencies change. The critical frequency being the largest frequency in the system changes. The smallest period of the system changes, as we are seeing here in this view graph, also.

And it is important that your time step is always smaller than the critical time step based on the highest frequency that you can ever encounter in the solution. If that is not the case, then you might have only an instability for a few time steps, and the solution may not have time to blow up, or give you very large displacements. So the instability will not be obvious. Instead, the instability over these few time steps only introduces an error in the solution, an error of which you really don't know the magnitude. And after these few steps have been completed, you're using now a time step that is again stable, giving a stable solution.

That error which we have accumulated, of course, is carried forward, but you have not seen a blow-up of the solution. And that's one important consideration when you apply the central difference method in non-linear analysis. Of course, the important point, the sort of ground rule that you have to observe, is that delta t must always be smaller than the critical time step of the mesh where that critical time step changes throughout the solution.

Let's look at this example here. What happens if we violate that condition? I show you here the solution with delta t equals 0.15 seconds, compared to the solution with delta t 0.05 seconds. Notice that the response of the right math here for both of these time steps is slightly different, but we can say that for this time step, we don't get displacements that are very much different than we obtain with this time step.

If we look at the response of the center mass, there are more differences in displacements. And if we look at the response of the left mass, in fact, we have even larger differences. But the important point is that if you only have this solution here corresponding to the time step 0.15 seconds, which is larger than the critical time step, you would certainly not expect that this solution was unstable. It in fact was only unstable for a very few time steps. There is an error accumulation, fairly large here. In fact, very large. You can see that these displacements are only about 50% of those displacements. But you cannot see directly an instability here if you don't know this solution.

Let's look at the force in the center truss. Here we list the time, here we list the solution for delta t, 0.05 seconds. And here we list the solution for delta t 0.15 seconds. And you can see here that at the beginning, the forces compare quite well, and then the error becomes larger and larger. In fact, here we have an opposite sign in the force. This being here, say, the accurate solution with a small enough time step, and this being, of course, a grossly inaccurate solution. But it is not evident that this is a grossly inaccurate solution if you don't know this solution. And that's the point I wanted to make with this example.

The next example that I like to discuss with you Is the analysis of a tower. The tower is shown here, subjected to a pressure that is induced by a blast. Here we have some geometric and material data regarding the tower. The applied load on the tower is shown here. Notice it's a load that we assume to be instantaneously there, and then to decay, as shown here, over 200 milliseconds. The purpose of the analysis is to determine the displacement and velocities at the top off the tower, to determine the moments at the base of the tower. And in this particular analysis, we use the trapezoidal rule and the lumped mass matrix to solve the problem.

We must make a number of decisions here, and two decisions are most important. We must choose a mesh that is, in other words, we have to establish an appropriate finite element discretization for the tower, and we have to choose a time step for the time integration. The choice of mesh and the choice of the time step are two very closely related items, and that's what I would like to discuss with you in this example in more depth. Of course, the time step and the mesh must both depend on the loading that is supplied.

Let's look at some observations. The point is that the choice of mesh will determines the highest natural frequency, and of course, the corresponding mode shape, that is properly, accurately represented in the finite element analysis. The more elements you use, just as a rule of thumb, the higher the frequency that can be properly represented, that can be accurately picked up, with the choice of mesh we're using. The time step that you're using will of course determine how high a frequency you're integrating accurately. In other words, you want to use really enough elements to have a mesh that is accurate, and you want to use is small enough time step to integrate the frequency of the mesh that you need to integrate in order to have an accurate response prediction.

I used the word "accurate," and we have to now of course look closely at what I mean by that. It is most effective, really, to choose the mesh and the time step such that the highest frequency, accurately integrated, is equal to the highest frequency accurately represented by the mesh. To choose the mesh and the time step, we look at the applied loading, in other words, the physics of the problem, and we represent this applied loading as a Fourier series. If we do so, that Fourier series will display the important frequencies that are to be accurately represented by the mesh.

Here for this particular loading that we're looking at, we did a Fourier analysis, and you can see that with this Fourier analysis, we include terms now in our first analysis, in case one, up to fn 17 hertz, and in case two, in our second analysis, up fn equal to 30 hertz. If we truncate our Fourier representation of the loading at 17 hertz, we obtain this approximation here to the actual applied noting. If we truncate our Fourier representation at 30 hertz, we obtain this Fourier approximation to the actual loading. As you can see, surely, as expected, if you include more terms in the Fourier approximation, your approximation to the actual loading becomes better.

Now having identified that we want to perform the analysis with two meshes, first there's a mesh that gives us frequencies accurate up to 17 hertz, and then with a mesh that is accurate up to 30 hertz, we have to look for those meshes. And we do so by establishing three meshes that are shown here. Here we have a 30-element mesh. Notice that one element spans from the bottom of a floor to the top off the floor, and across from one column to the other. A 60-element mesh, in which we simply have subdivided each element here in this more coarse mesh into two elements, and a 120-element mesh, in which we again have taken one element here and represent it as two elements. We use these meshes here to calculate frequencies. And on this view graph here, you see the frequencies calculated for the 30-element mesh and the 60-element mesh.

Notice mode numbers in this column. 30-element mesh results here, 60-element mesh results here, in each case, hertz. We notice that up to this red line, the difference between the frequencies is less than 1 hertz. In other words, the difference between the frequencies obtained from the 30-element mesh and the 60-element mesh is less than 1 hertz, and that is also right at the 17 hertz level. So we can say that the 30-element mesh really represents the frequencies up to 17 hertz quite accurately.

You should also remember, please, that since we are using a lumped mass matrix in the analysis, there's no bounding theorem saying that the frequencies of the 60-element mesh must always be below the frequencies of the 30-element mesh. For a larger mesh, we can actually have frequencies that are smaller or larger than for the more coarse mesh. That of course also applies for when we compare the 60-element mesh result with the 120-element mesh result later on.

So we say that for the 30-element mesh, we have accurately represented now the frequencies up to 17 hertz, and therefore, there is the mesh to use. We calculate the time step via is this relationship here as being 0.033 seconds, being the cutoff frequency, and delta t then is given here as 0.0017 seconds. This mesh now should be integrated with this time step, and then we can say is that the mesh up to 17 hertz represents the frequencies quite accurately, and we have integrated the frequencies up to 17 hertz quite accurately, and therefore, we have an accurate solution of the problem up to 17 hertz-- and what I mean by this "up to 17 hertz" is that also in the Fourier representation, we in effect only have represented the load by the Fourier representation up to 17 hertz, in effect. Of course, in the actual step-by-step solution, we also will carry terms above 17 hertz. But these terms are assumed not to contribute to very much to the response, because we are not picking them up very well in the mesh and in the time step selection.

If we proceed now with the same ideas to also look at the 60-element mesh, then we perform there also the frequency analysis for the 60-element mesh and the 120-element mesh, and we find that up to 30 hertz in the 60-element mesh, we gain represent all the frequencies accurately. Notice that above the red line, below these frequencies here, the difference between these frequencies is less than 1 hertz. Above this line, the difference between the frequencies is larger than 1 hertz, and that is our cutoff line.

Well, if we want to perform now the direct time integration of the 60-element mesh, we proceed as for the 30-element mesh. We calculate our cut-off frequency, our delta t, the proper time step, 20 steps for the highest frequency or the smallest period that we want to integrate accurately. That's the rule of thumb we're using. And with this time step, we perform the analysis.

Notice a smaller time step would be accurately integrate frequencies which are not accurately represented by the mesh. Therefore, there is no point choosing a smaller time step. A larger time step would not accurately integrate all frequencies which are accurately represented by the mesh. It is not good to choose a larger time step, because your mesh is chosen to represent the frequencies up to 30 hertz quite accurately, and to get all of what you can get out of the mesh, so to say, you want to take a time step that is in accordance selected with the mesh. And that times step selection gives you this result here, 0.003 seconds.

Let's look now at the actual response that is calculated. Considering the horizontal displacement at the top of the tower, denoted as u, here, we find that in fact the 30-element mesh and the 60-element mesh give virtually the same response. If we look at the horizontal velocity at the top off the tower, we find that there are some small differences between the 30-element and 60-element mesh. If we were to look at the accelerations, these differences would of course be even larger. But the purpose of this analysis was to calculate the displacement and velocities at the top of the tower.

Let's look at the moment, which was also one requirement for our analysis, the moment right down here at the bottom off the tower. And here we see that the 30-element mesh gives the smooth result. The 60-element mesh gives the zig-zagged result. But overall, really, basically a good response prediction already with the 30-element mesh.

If we look at the displacements of the tower, we have amplified these displacements here. The 30-element mesh at 200 millisecond looks like that. The 60-element mesh. looks like that. At 400 milliseconds, the 30-element mesh looks as shown here. The 60-element mesh looks as shown there.

Regarding these analyses, we now can summarize some comments. We saw some high-frequency oscillations when we looked at the moment reaction response of the 60-element mesh. However, this high-frequency response, which actually can be identified to be 110 hertz, is probably inaccurate, because our mesh did not represent frequencies accurately up to that level. The obtained solutions for the horizontal displacement and the velocity at the top off the tower are really virtually identical for the 30- and 60- element mesh.

Let us now look at another example. This is the example of a pendulum. A large displacement analysis of this pendulum is performed, and basically what we're looking at here is, if I can just have your attention here, is at the pendulum that originally is horizontal, it's pinned right here. And we let this pendulum go through very large motions, and pendle around as shown here. If always should come back to the horizontal position, unless we have damping in the system. Since we don't have damping in the system, the pendulum should just keep on going forever.

We have a frictionless hinge right here, and this is pictorially represented here on the view graph, a frictionless hinge right there. The pendulum is modelled by a truss element of stiffness EA. The length is given here. The mass is concentrated right there. So there's no distributed mass here. All of the mass of the pendulum is concentrated here, half of it here, and half of it there.

The initial conditions of the pendulum are listed here. The initial angle is 90 degrees theta, so the pendulum is up here, as I explained already. And initially, it has 0 velocity. And now we let the pendulum go, and as I said, it should spring up and down here, of course, under the gravity loading.

The dynamic response that we want to calculate is obtained using the trapezoidal rule. And we are using full Newton iterations to establish the equilibrium during every time step. The convergence tolerance that we use, and that I have to refer you now to our earlier view graph regarding convergence tolerances. The tolerance that we're using here is ETOL equal to 10 to the minus seventh. We talked about what ETOL is in an earlier lecture.

In our first analysis, we chose the time step to be equal to 0.1 seconds. And in that case, the following response is obtained here. We expect that solution is shown as a solid line, and these circles here show the calculated finite element solution. Notice 90 degrees theta here minus 90 degrees here, and then we return to 90 degrees.

The solution looks quite fine, looks stable. And in fact, if you look at the textbook, you would find that there I'm giving the solution up to here. Everything looks good. However, if we plot the strain in the truss, we find that at larger times, the strain oscillates very much, and in fact, this oscillation grows. This clearly signifies an instability in the solution process. The instability is unchanged when we tighten our convergence tolerance or when we use the BFGS method instead of the full Newton method.

We should recall here that the trapezoidal rule is only unconditionally stable in a linear analysis. Here we are talking about a non-linear analysis, and you cannot select any time step magnitude delta t and expect a stable solution. So what we have to do now is to use a smaller time step, and we choose delta t equal to 0.025 seconds, and repeat the analysis. Theta degrees here as a function of time. And once again, a very nice, smooth response solution, starting with 90 degrees, going to minus 90 degrees, coming back to 90 degrees, and so on.

Of course, even with a larger time step, you got a smooth solution for theta. What really told us about the instability in the solution was our look at the strain in the truss. So we look now also at this strain, and in this particular case, with this time step, we see a smooth solution for the strain. You should please observe that the scale on the strain is here different from the strain scale that we had in the earlier view graph. If you refer back to the earlier view graph, you see that we had a different scale there. But we have a smooth solution here, in fact, a solution the way we would expect it.

If we don't iterate, then we obtain a solution that looks as shown here. In other words, instead of getting this oscillatory response that we would like to obtain, we obtain really a highly damped solution response. This, of course, shows the importance of iteration. We have to iterate in the solution of a dynamic problem. It's very important to integrate in order to establish equilibrium for every time step, for every discrete time that we're considering. In this particular case, we see that we obtained a highly damped solution. In other cases, the solution might blow up. So it is important to iterate in the nonlinear dynamic solutions when you use an implicit time integration scheme.

Although the solution obtained with the equilibrium iteration is highly inaccurate, the solution is, of course, stable. And that is also seen when one looks at the strain plot. Notice here we're plotting strains in the truss element as a function of time. No blow-up of that strain, no violent oscillation. It is, in fact, going up here, and then becoming this rather small oscillating, as shown here. So we don't have an unstable solution, but a highly inaccurate solution.

Let's look at the next example now, the analysis of a pipe whip problem. Here we have a long pipe, dimension given here and here, that is subjected suddenly to a concentrated load at its end. There is also, below the pipe, a restraint. That restraint is only active once this gap is closed. Notice there is a 3-inch gap here between the pipe and the restraint.

Of course, this restraint, in an engineering design, is installed to prevent violent emotions of the pipe at this point. What our objective here is, to determine the transient response of this pipe. Notice that under this very large load, the pipe very rapidly becomes plastic. And the difficulty of this problem really lies in that the pipe going plastic becomes very flexible here, but the restraint suddenly hitting the pipe, or the pipe hitting the restraint, rather, means that a very large difference is suddenly introduced at this end.

The finite element model we used is shown here. 6 summation beam elements that of course can go plastic. There's more of the elastoplasticity in the system properly. And a truss element from this node to that node. This truss element contains a 3-inch gap.

The material properties for the pipes are listed here. Notice the material property corresponding to the elastoplasticity. We assume no strain hardening and for the restraint. Again, the restraint also is modeled as an elastoplastic material. The analysis performed, in this particular case, using multiple position with two modes and direct time integration. We want to compare the response that we predict using the multiple position procedure, which I discussed earlier in this lecture, with the response that we obtain using a time integration analysis.

For the time integration analysis, we use the trapezoidal rule. And in fact, even in the multiple position, we also use the trapezoidal rule. And for the mass representation, we have chosen a consistent mass matrix. In this particular case of convergence tolerance, ETOL was used with this value down here.

If we perform the eigenvalue solution of the initial system, in other words, of the elastic pipe in the initial state, we obtain these natural frequencies, 8.5 hertz and 53 hertz. These are also the mode shapes. And of course, these frequencies will change as time progresses in the elastoplastic response, and also the mode shapes will change, particularly once their restraint is hit. However, what we are saying is that we believe we can capture the elastoplastic response of the pipe, approximately, by using this mode shape and that mode shape as generalized displacement patterns to represent the total pipe response. Of course, there are the heavy approximations that we reduce all of the pipe degrees of freedom to adjust 2 generalized degrees of freedom.

To integrate a response, we have to choose a time step. And here we have delta t 1/20 of the cutoff period that we want to use here. Notice that we have only two modes in this problem that we want to consider, with the frequencies that I showed you just now. And this equation gives us, as a time step, 1/1000 of a second. Of course, once again, this estimate, or this choice of time step, is based solely on looking at the linear response of the pipe. In other words, assuming that the pipe remains elastic and that the restraint is not there either.

If we look at the results obtained from these two analyses, we find that the multiple position response for the displacements looks as shown here, the direct time integration response looks as shown here. Notice that this is here the distance of the gap, so at this particular time, the gap is closed and the restraint comes into action. We note from this graph that the multiple position solution is really very close to the direct integration solution for the displacement at the tip off the pipe.

Now let us look at the moment at the built-in end of the pipe. And here we have the moment response plotted as a function of time for the direct integration solution, as shown here, and for the multiple position solution, as shown here. We notice that the maximum moment occurring at the built-in end in the pipe is almost the same when predicted by the multiple position solution and the direct integration solution. There's some small difference, but the difference is certainly not very large.

However, if we look at what happened at earlier times, we see that there is a large difference between the results obtained using the multiple position solution and the direct integration solution. The reason is that the response here at earlier times is very complex, and the multiple position solution cannot pick that complex response up very well.

However, it's one important point you should keep in mind. And that is, if we had used more modes in the multiple position solution, we would come closer to the response obtained, predicted with the direct integration solution. In fact, if the number of modes used in the multiple position solution is equal to the number of degrees of freedom in the finite element mesh, then the multiple position solution gives you the same response as you calculate with the direct integration solution. I mentioned that earlier in the lecture already, and the reason, of course, is that in that case, all you would have done is a change of bases from the finite element displacements to the mode shape bases.

So if you use, even with this initial stiffness matrix that you have used to calculate the mode shapes, even with that initial stiffness matrix, if use all the modes, as many modes as there are degrees of freedom, in this multiple position solution, you would get the same response prediction with the multiple position solution as we're seeing here, with the direct integration solution. And that means we can be quite assured that as we increase the number of modes in the multiple position solution, we will get closer and closer to this complex response that is predicted using the direct integration solution.

In general, one, of course, does not know the error when you use the multiple position solution, and there must always be a question of how large is the error, for example, that we're seeing here. And that can only be assessed by a more accurate solution. However, I believe that's the multiple position solution is of much value in the initial design analysis stages when you want to obtain fairly inexpensive solutions for the response, and you want to just understand how the system acts that you are analyzing and designing. The multiple position solution really yields relatively inexpensive response solutions of the system which may be very helpful in the initial design stages of that structural system.

This completes, then, the discussion of our examples for which I have prepared view graphs. I'd like to now share with you some further experiences that are summarized on this slide. And let me go over here, where I have the first slide projected. I should point out, before we discuss the information on these slides once more, that we will go fairly rapidly through these slides, and that not all details pertaining to the solutions that we are looking at here are given on these slides. Please refer to the papers-- the references are given in the study guide-- in order to find all of the details pertaining to the solutions that I'm now discussing with the slides.

Here we consider now the analysis of a control rod drive housing, which is modelled using four beam elements and point masses. Notice that the structure is subjected to motion up here. The non-linearity in the system lies in the fact that there's a gap between the tip of that beam here and the spring on the right-hand side and on the left-hand side.

On the next slide, we see now the solutions obtained. We obtain a solution using direct integration, and using a multiple position analysis, including two modes. These two solutions are compared with a solution that was reported at some earlier time by Peterson and myself.

Here you see the comparison of the results obtained. Notice there is a difference here between the solution that we reported earlier, that Peterson and myself reported earlier, and the solution that we generated now. But that difference is due to the fact that Peterson and myself used, at that time, a different solution scheme. The response solution that we obtained with the multiple position analysis is very close to the response solution that we calculated using direct integration.

Now the next analysis that I'd like to discuss this you involves the analysis of a spherical shell subjected to a uniform pressure. And that pressure is suddenly applied as shown down here. Notice this is a step pressure constant in time. The material data of the shell are given here. We want to calculate the elastoplastic response. And some geometric data are given here.

Notice the shell was idealized using 10 8-node axis-symmetric elements. It's a spherical shell, so we use an axis-symmetric model to model the shell. We use Newmark integration with delta and alpha as given here. We chose these values because Nagarajan and Popov had chosen these same values, and we wanted to compare with their solution. We used 2-by-2 Gauss integration and a consistent mass matrix for the problem. The time step was 10 microseconds, and we used the total Lagrangian formulation in the analysis.

The next slide now shows the response that we predicted. Here we see the ADINA solution, a dashed line, compared with the solution by Nagarajan and Popov given by the solid line. Notice we are plotting here the apex displacement of the shell as a function of time. Now in this particular case, we use the consistent mass matrix, 2-by-2 integration, and specific solution parameter, and we want to do now change these parameters in order to see the effect of these parameters on the solution that is predicted.

The next slide now shows again the apex displacement of the shell as a function of time, but using once a consistent mass matrix, and once using a lumped mass matrix. And in this particular case, we use the trapezoidal rule, which corresponds to the Newmark integration with delta equals 0.5 and alpha equals 1/4. Notice that there are only slide differences between the lumped mass solution results and the consistent mass solution results.

Another parameter that one might want to vary is the integration order. Here we have the response, using 2-by-2 integration, which we saw already earlier. Actually, we didn't quite see this response earlier, because earlier we used different delta and alpha parameters. Very close to these parameters, however.

Anyway, here we have the 2-by-2 integration results, and then we have the 3-by-3 integration results and the 4-by-4 integration results here. Very small differences between these solutions, but it's interesting to perform these solutions, and really assure oneself that there are small differences. We know that shells are rather sensitive structures, and we want to make sure about the modeling that we're using for the representation of these types of structures.

The next slide now shows the analysis of a problem for which we had experimental results to compare with. Here we have a pipe, a very rigid pipe, and a flexible pipe, a nickel pipe here. And these pipes, both being straight, are totally filled with water. At the end off this pipe here, a pulse gun sets off a pulse, and that pulse travels through the water, all the way along. And of course, there's an interaction between the pipe and the water. In other words, it's a fluid structure interaction problem.

This pipe here is very rigid, so that we model it as infinitely rigid. You will see that just now on the next slide. Notice here we have pressure gauges indicated. The pressures and the strains in the pipe, in this nickel pipe, the flexible pipe, were measured by the Stanford Research Institute. Here we have the nickel properties, we assume an elastoplastic material, and the water properties.

On the next slide now, we see the pressure pulse as a function of time. Of course, this was given to us for our analysis. And this slide now shows the finite element model that we use. We used an axis-symmetric model. The center line is right here. And notice that the water here in the rigid pipe was represented by 81 elements. Notice the rigid pipe was modelled as infinitely rigid, because these are rollers on an infinitely rigid surface, so to say. Notice here is a nickel pipe representation, and here we have water elements again, and a thin layer of water elements here.

On the next slide now, we see the response predicted by ADINA and the experimental results, the ADINA response shown as a solid line and the experimental result shown as a dashed line. Of course, the experimental results reported by the Stanford Research Institute. Notice we are plotting here the pressure as a function of time for one of the pressure stations. And here you'll see the same information for another pressure station, and for another pressure station, and one more pressure station-- the pressure stations that I pointed out to you earlier.

Notice that here finally we compare the hoop stress predicted by ADINA with the experimentally observed hoop stress as a function of time at one of the stations within the pipe. Notice that for the hoop stress here, we have four experimental results corresponding to the measurements that were taken at the north pole, the south pole, east pole, and west pole on the pipe. In the ADINA solution, of course, we use an axis-symmetric model, and that axis-symmetric model only gives one hoop strain. And that hoop strain, as you can see here, lies quite nicely within the experimental results obtained by the Stanford Research Institute.

I should also point out that the comparison of the ADINA results with the experimental results is really very good, except that we don't pick up these high-frequency oscillations that have been observed in the experiment. Well, the ADINA finite element mesh that we used in this analysis is not able, really, to represent these high-frequency oscillations, and of course, that is certainly one of the reasons why we could not hope to pick these frequency oscillations up.

This, then, brings me to the last slide that I wanted to share with you in this lecture, and also to the end of this lecture. Thank you very much for your attention.