Lecture 10: Solution of Equilibrium Equations in Dynamic Analysis

Related Resources

Study Guide (PDF)

Readings

Sections 9.1, 9.2, 9.4

Examples

Problems 9.1-9.5, 9.12-9.14

• 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 Lecture Number 10. In this lecture, I would like to discuss a solution of dynamic equilibrium equations. And in particular, I would like to talk about the direct integration solution of the equilibrium equations in dynamic analysis.

These, of course, are the equilibrium equations that we derived already earlier. M as a mass matrix, C as a damping matrix, K as a stiffness matrix, U as a displacement vector, and here we have the velocities and the accelerations at the null points. R, of course, being the load vector and this load vector is now time dependent.

I also should draw your attention on the fact that I would like to talk about direct integration of these equations. We will talk about the modes of opposition procedures in the next lecture. By direct integration, we mean that we integrate these equations directly without-- that is, without-- a transformation of these equations into a different form prior to the integration.

In modes of opposition analysis, now you will see in the next lecture, we're actually transforming these equations first, and then we integrate the response. So the procedures that I'd like to talk about are explicit, implicit integration procedures that directly operate on these equations without a transformation and integrate, solve these equations.

I also would like to talk about some computational considerations with regard to these explicit and implicit integration procedures, the number of operations involved, the costs involved in performing these integrations and so on. And then there is an important consideration, namely the selection of the solution time step, delta-t.

We will see that when we integrate the dynamic equilibrium equations, we have to select the time step delta-t with which we perform the integration. And this time step has to be selected in order to preserve stability of the integration and to preserve of course or to obtain accuracy in the integration. Finally, I'd like to talk about some modeling considerations pertaining to the solution of the dynamic equilibrium equations.

Well, here on this view graph, I'm showing once again the dynamically equilibrium equations here, noticing once more that R, the load vector, is now a function of time. The same equations here have been rewritten once here, where FI denotes inertia forces, now time dependent, of course. The mass matrix is constant, but the accelerations of course depend on time.

Then this part here, which are the damping forces, has been expressed here as FD a function of time again. The C matrix is constant with time, but the velocities, of course, vary with time. So here we have damping forces. It's a nodes of the finite element system. And these are the elastic forces, KU, which are really due to the internal element stresses. And on the right hand side, we have the external forces.

So what we're saying is that we want to solve this force equilibrium over all times t. The inertia forces plus the damping forces plus the elastic forces, it's the nodes of the finite element system, of course, calculated in the virtual work sense, the way we've been talking about in the earlier lectures. The sum of these forces here must be equal to the externally applied nodal point logs.

Well, the procedure that we will be following is the following. We will consider time steps delta-t a part at which we want to calculate the response of the system. In order to do so, of course, we have to define our loads first. And here, I have schematic shown once the description of the load Ri, that is, the load at degree of freedom i as a function of time t.

Notice that we are having here this distribution of that load component as a function of time. The distribution is quite simple. In a computer program, of course, all we would have to put in is this point, this intensity here, that load intensity, at a specified time. That load intensity, at a specified time, and that load intensity at a specified time. In between, these intensities, the computer program would simply interpolate linearly.

This is the load at the component i, displacement component i. This is the load at displacement component j. Notice here, of course, we would have in general a different curve at degree of freedom j, then the curve is used at degree of freedom i.

I show here time steps, discrete time steps, delta-t1, delta-t2, delta-t3. These of course apply for the complete system, for the i component and the j component of displacements or loads, and for all other displacements and loads also.

Notice that they can vary with time. In many analysis, however, we simply pick one time step and keep it constant all along. The objective is the following. In the direct integration of the equations of equilibrium, the objective is knowing the initial conditions. In other words, knowing the displacement at time 0, the velocities at time 0. We want to calculate the displacements, velocities, and accelerations at these discrete time points.

We are matching here with the time interval delta-t1. And here we are changing the time interval to delta-t2. Knowing all the velocities, accelerations, and displacement at this time, we want to calculate the same quantities at the discrete time increments of delta-t2, and so on.

So basically what we are saying is that knowing the initial conditions, we want to calculate in a step-by-step procedure, the displacements, velocities, and accelerations at these discrete time points. Of course, the real problem then reduces to the following. Given the initial conditions at the particular time t, we might call this time t, and by initial conditions, I mean we know the displacements, velocities, accelerations, we know all the quantities at time t.

Our objective is to calculate the unknown quantities at time t plus delta-t1 in this particular case. So this point is t plus delta-t1. In general, or in many cases, we keep a constant time step, and so we will simply talk about t plus delta-t from now onwards.

We should remember, however, the delta-t might change doing the response solution. I repeat once more, this is an important point, given the initial conditions or given the conditions at time t, our objective is to calculate the conditions, displacements, velocities, accelerations, at time t plus delta-t.

In the explicit integration, we proceed in the solution as follows. We use the equations of dynamic equilibrium at time t to obtain the solution at time t plus delta-t. This is important. In explicit integration, we use the condition, the equilibrium equations, the equilibrium equation, the dynamic equilibrium equation, and by that I mean, just to refresh your memory, this equation-- we use that equation at time t in explicit integration to obtain the solution at time t plus delta-t. Once we have the solution at time t plus delta-t, of course we can march ahead in the same way again.

In implicit integration, however, we are looking at the equilibrium equations at time t plus delta-t t to obtain the solution at time t plus delta-t. In other words, if you look at once more at the equations here, at the equilibrium equations, these equations would be applied at time t plus delta-t to obtain the solution at time t plus delta-t. This is indeed the fundamental difference between explicit and implicit integration.

In explicit integration, we are looking at equilibrium equations at time t to obtain the solution at time t plus delta-t. And in implicit integration, we are looking at the equilibrium equations at time t plus delta-t to attain the solution at time t plus delta-t.

In the following now, I would like to show to you an explicit and an implicit integration scheme. The explicit integration scheme that is abundantly used in practice is the central difference method. I will later on then discuss how this method is effectively used.

In the central difference method, we are applying, as I said earlier, the equilibrium equations at time t, because it's an explicit integration scheme. The two additional questions that we're using to solve for the unknown displacements at time t plus delta-t, is firstly an equation for the accelerations at time t. And notice that this equation involves as the unknown, as the only unknown, the displacement at time t plus delta-t.

In addition, of course, we are using an equation for the velocities at time t and again this equation only involves as the only unknown the displacement at time t plus delta-t. Substituting these two equations into the equilibrium equation at time t gives us one equation in the unknown t, displacement at time t plus delta-t. We can simply solve for it, as shown, now on the next view graph.

You can see here the equation that is obtained by substituting the relations for the acceleration and for the velocities in the equilibrium equation at time t. The unknown appears here. This is the displacement at time t plus delta-t, which is unknown.

And there is a coefficient matrix here, that involves the mass matrix and the damping matrix. And on the right hand side, of course, we are calculating a load vector. And that load vector, notice once again, corresponds to the conditions at time t. That's why it's an explicit integration scheme.

The difference operator, the essential difference operator that I've been used here, simple operators which you have seen on the previous slide, that involve only the displacement at time t minus delta-t and the displacement at time t. And notice that here are these two displacement vectors appearing.

Notice that this equation is cheaply-- or can be used in a very cheap way to calculate the displacement, the unknown displacement at time t plus delta-t, if the mass matrix is a diagonal matrix. And, say, the damping is neglected. In other words, this term not being there, this being a diagonal mass matrix, then the solution is very trivial. In fact, all we need to do is calculate the right hand side load vector and divide each component of that load vector by this coefficient here which, if m once again is a diagonal mass matrix, is simply Mii divided by delta-t squared.

And in fact, this is really how was the procedure is usually applied, namely, with a diagonal mass matrix and no damping. The damping usually is not taken into account, or at most, we have diagonal dampers also, and the mass matrix being a diagonal matrix, then the solution is very cheaply performed.

There's one further important point. This load vector of course is constructed in the usual way. This multiplication here is quite cheap if the mass matrix is a diagonal mass matrix. K times tU. Well, K here is of course a banded matrix. It's the stiffness matrix of the system, constant.

Well, there are some operations involved here. And if we look at this operation here in detail, we see here how it can be effectively performed. We have a stiffness matrix times the displacement vector here. And of course, the stiffness matrix remember is obtained by summing the element contributions.

So you can write this K matrix in this way. Since tU is independent of the elements, we can also take it into the summation sign. And if we look at this product closely, we find that is nothing else than tFm. It's a nodal point vector of loads or elements forces that correspond-- element elastic forces-- that correspond to the current displacements. And this tF vector is much more cheaply calculated than the Km times tU. And in fact, this is the way how we actually perform the calculation of this product here.

Now notice the following important further point. If we were to leave this K times tU the way it stands here, we would have to assemble a global stiffness matrix of the system. However, if we write that product this way, we never need to assemble a global stiffness matrix of the total system. In fact, all we need to do is calculate element forces and add these elements forces via the direct stiffness procedure the way I've been discussing earlier.

So this is a very effective way of calculating K times tU. We never need to assemble a global stiffness matrix. In fact, this is one of the main reasons why the central difference method is effective. It's this reason that we can calculate the KtU part, this times that, very effectively via this procedure. That is the first reason.

And the second reason is that if M is diagonal, we really never solve equations. The equation solving here is so trivial, that we might just say that we do not solve equations. There is no factorization involved. There's no factorization of a stiffness matrix or another coefficient matrix involved, because if M is diagonal, and as I pointed out, we simply divide the right hand side by Mii divided by delta-t squared.

Notice that if we start the integration process, we know of course the initial conditions. We know tU. And we do not know-- usually directly the t minus delta-tU. That is the displacement prior to the solution.

And to obtain those, we can use this formula, which is really the central difference method applied at time t, where we however now involve the accelerations corresponding to time 0. We involve the accelerations corresponding to time zero right here, and the velocities corresponding to times 0, and the displacements corresponding to time zero. And of course, these quantities must all be given as the initial conditions to the solution.

I pointed out already that we're using this technique mostly with the lump mass matrix. In addition I should also point out that we're using it primarily with lower order elements, with lower order elements. The reason for that is the following, namely, one of stability and accuracy of the central difference method. The method can only be applied if delta-t, the time step that I'm talking about here, is smaller than a critical time step. And that critical time step is given by the smallest period in the system divided by pi, 3.14 et cetera.

The smallest period in the system can be of course very small. And therefore the critical time step that we're talking about here can be very small. The smallest period in the system, just to give you an example, would be extremely small if you have a bar element, and you have some very short truss element in there. The shorter the truss element, the stiffer of course it is. AU over L is the stiffness of a truss element, and therefore, a very small period would be in that assemblage of elements.

Therefore, what is required really is to have a mesh that is more or less uniform, so that there is no artificially small time step required in the integration of the equilibrium equations. Because there is this condition that delta-t must be smaller-- we should really say must be smaller or equal-- to delta-t critical, but to be on the conservative side, we are selecting delta-t smaller than delta-t critical in practical analysis of complex meshes. Because this condition has to be satisfied, we're talking about a conditionally stable scheme. It is only stable when delta-t is smaller than that value.

If we were to select delta-t larger than delta-t critical, we would see that after just a few time steps, the solution-- it gets out of bounds. And after 10 or 20 or 30 time step, in a computer program, you would all of a sudden notice overflow of the numbers. So the solution wouldn't make any sense whatsoever anymore.

And this is an absolute condition that does have to be satisfied. There's no way you can use a time step larger than delta-t critical, and still get an acceptable, reasonable solution to your equations of motion.

In practice, of course, we have to estimate, in a conservative way, a time step. And that means we have to have an estimate also of delta-t critical. For continuum elements-- by continuum elements, I mean plain stress, plain strain, axis a matrix, three-dimensional element, the way we have been discussing it earlier.

A good way of obtaining time step delta-t is to use this formula here, delta-t being delta-L over C. C is the wave velocity through the system. [UNINTELLIGIBLE] is divided by rho, square root out of it gives us C. And delta-L is an element length that has to be defined. In fact, it's the smallest distance between nodes. The smallest distance between nodes if we talk about low order elements.

What I mean by that is the following. If we have a mesh here that looks like that. And I deliberately show elements that are not equal in size, then delta-E in this particular case would be this distance here. That is our delta-L, sorry, our delta-L is that distance there.

This, of course, as I have here in the heading, assumes that we have low-order elements. That means only corner nodes in the elements. For higher-order elements-- let me put down here one higher-order element, where we have these nodes now, say an 8 node element. In this particular case, we take the smallest distance between nodes, again, just as in this case. And in this particular case, it would be, say, this distance here. That is our smallest distance between nodes.

However, we have to also divide this distance by a relative stiffness factor. And why does that relative stiffness factor come into the picture? Well, it does come into the picture, because remember that the stiffness at an interior node is larger than the stiffness at a corner for the element. The amount it is larger depends on the particular element used. But since there's a lot of stiffness here, we have to put here a relative stiffness factor into the calculation of delta-L, the effective lengths that we're using in this particular formula.

Now, for this element here, for example, the relative stiffness factor would be conservatively just equal to 4. But for other elements, and particular also when there's node shifting involved and so on, the relative stiffness factor has to be calculated, has to be estimated. And we should always remember that it should be estimated here in a conservative way, because a smaller delta-L gave us a smaller delta-t. And then we are then conservatively applying the fact that delta-t shall be smaller than the delta-t critical, value that I have been referring hereto.

The method is used mainly for wave propagation analysis. The reason being that we have to use a fairly large number of time steps to predict the solution. Large, because a time step delta-t is usually very small.

And in wave propagation analysis, many modes of the system are being excited. So we really want to pick up a large number of modes in the mesh. And the small time step integrating accurately many modes in the mesh goes together with the fact that we want to predict many modes in the mesh or a large number of modes in the mesh, and therefore the method is effectively used for wave propagation analysis.

The number of operations are proportional to the number of elements that we are using and the number of time steps. Of course, this is an approximate proportionality here that I'm talking about. This fact derives directly from the consideration that we talked about here. Notice that as we include more and more elements in the solution process, of course, we have to sum over more elements here.

And as we go more and more time steps, we have to apply this equation as many times also more and more. The important point is that we don't assemble a K matrix, so the bandwidth off the system, the bandwidth of the system does not appear in an operation count. It's the number of elements that do appear in the operation count.

The next procedure that I like to talk about is the Newmark method. And in the Newmark method, this being an implicit integration scheme now, we proceed as follows. We apply-- because it's an implicit integration scheme-- the equilibrium equation, the dynamic equilibrium equation at time t plus delta-t. Here it is given the additional formula that we're using is this one here for the velocity at time t plus delta-t, and this one here for the displacement of time t plus delta-t.

Now notice that with this set of equations, with these three equations, knowing the displacements and velocities at time t and accelerations at time t, we have with these three equations, three unknowns. The velocity at time t plus delta-t, the acceleration at time t plus delta-t, and the displacement at time t plus delta-t. These are the only three unknowns that appear. There are three linearly independent equations. And we can simply use them to calculate these three unknowns. Of course, we want to do so in a very effective way. And I will show you how we proceed.

The important point is that with this scheme, we are using an implicit approach, as I defined it earlier. And when we do substitute from here into there, we will find that on the left hand side, we have a coefficient matrix that involves a bandwidth even if the mass matrix is diagonal. And this coefficient matrix here is written here down as K hat. This coefficient matrix invoice this K matrix.

I should say maybe a few words to these formulae. These formulae can be obtained by Taylor Series expansion, around the conditions at times t plus delta-t. They have been proposed originally by Newmark. That's why it's called now the Newmark method. Newmark in particular used the values of delta equals 1/2 and alpha equal to 1/4. And then it's a constant average acceleration scheme.

The constant average acceleration schemes use these values as I pointed out already earlier. The method is unconditionally stable. This means that the time step delta-t does not to have to be smaller than a certain value, as it was the case in the central difference method. The method is primarily was for analysis of structural dynamic problems. I will talk about why this is so just now.

And the number of operations that we're talking about here are an initial factorization of the K hat matrix. This is the K hat matrix, And that initial factorization involves this K matrix, so the bandwidth of the K hat matrix is the same as the bandwidth of the K matrix. And therefore the operations here are one half NM squared.

And for each solution, in time, because this equation here has to be solved for each time step, you're starting with delta-t, and then we go to 2 delta-t, 3 delta-t, 4 delta-t, and so on. And since this solution is to be formed for each time step, however, remember the factorization is only done once because this K hat matrix is independent of the time.

Therefore, we have a forward reduction and back substitution in each time step. And the number of operations in each time step are 2nm. t is the number of time steps here. t is the number of time steps here.

Well, let us look then at the accuracy considerations that we have. I mentioned already the method is unconditionally stable, which means that the delta-t time step that we're using can be any amount, and we will never have a blow up of the solution, a blow up of the solution due to round-off errors that I involved.

In the central difference method, of course, the explicit time integration scheme, I mentioned that delta-t has to be smaller than a critical time step. So in the implicit Newmark integration scheme, what we can do is we select our of time step based solely on accuracy considerations. Let us consider then these equations that we want to solve, and I'm now leaving out the damping term in order to obtain some understanding of how we would select the time step, delta-t.

Of course, remember the larger the time step, the less normal operations I involved in the solution. So we want to select the large time step. However, if we make the time step extremely large, then surely we cannot expect any accuracy in the solution. In fact, if delta-t is very large, we simply predict the static solution. This is one property of the Newmark method. If you were to use a computer program with the Newmark method, if you're interested in it, why don't you try, say, for example, using ADINA, in which we have the Newmark method, and just put delta-t very large, and all you would be getting out is a static solution.

Well, since we're interested in the dynamic response, we have to therefore select delta-t small enough to integrate accurately the solution response. And these are the equations that I want to now consider, and I want to show you how can we select delta-t rationally and to obtain the required accuracy. The displacement vector here can be expressed in terms of the mode shape vectors times xi values.

I will be talking about how these phi i vectors are calculated in the next lectures. This involves the solution of the eigenvalue problem written down here. The free vibration eigenvalue a problem corresponding to this equation.

I like to defer the discussion at this point, because we will be talking about it later on in detail. All I like to say is that this transformation here, phi i being a vector time independent and xi being a scalar, dependent on time. This transformation can be used in this solution here, or can be used rather to transform this system of equations into a different form.

And in that transformation, we use the fact that when we calculate phi transpose K times 5-- notice these are capital phis. I've crossed bars here, capital phis. This is a lower phi, no crossed bars top and bottom. This capital phi here stores all of the lowercase phi vectors, and their n of such vectors.

If I use this capital phi and calculate this product here, I get a diagonal matrix. In fact, the diagonal matrix here stores on its diagonal the frequency squared, the free vibrational frequency squared. Also phi transposed M phi gives us an identity matrix. These, once again, the vectors in this matrix, capital phi, are these phi i vectors which are the eigenvectors of this problem. And I will talk later on about how we calculate them.

Well, if we apply this transformation-- whoops, let me put it right once more-- if we apply this transformation to this system of equations, and we use this property right here, we obtain n decoupled equations of this form.

X double ought i plus omega i squared, xi equals phi i transposed times R. And they are in such equations, in other words I run c from 1 to n. In fact, if we use multiple position analysis, we actually perform this transformation.

What I now want to do is simply refer to it theoretically. I do not want to transform it. You're talking about direct integration of the solution of equations, in which there is no transformation involved. So I do not want to go through this transformation.

However, I want to refer to it theoretically. And it gives me a tremendous insight into what happens. So basically, what I can say theoretically is that the direct step-by-step solution of this system of equations, using the explicit integration scheme, the central difference method, or the Newmark implicit integration scheme, corresponds completely to the direct step-by-step solution of n equations, of this form, where the displacements here, U, are given in this way. x being now a function of time, so if I differentiate here, your double dot would be obtained by differentiating xi also twice, with respect to time.

Once again, this is an important, very important, fact. I look at these equations, and I say that's a solution of these set of coupled differential equations. It's completely equivalent to the solution of n equations of this form with this transformation used. No dots here or dots there. In other words, this equation holds for displacements without the red dots and for the accelerations with the red dots. It is equivalent. This solution is equivalent to that solution.

However, in direct integration, I never perform actually the transformation. All I know is that when I operate on this system of equations, using the Newmark method or the central difference method, I, in fact, operate on n such decoupled equations. That's important.

And this means that for an analysis on our direct integration scheme, we can in fact analyze our direct integration scheme by looking at these equations, instead of these equations. You see if I were to analyze my direct integration scheme using these equations, I would have to deal with this huge matrix K, the huge matrix M. All I now need to deal with, in the analysis of the direct integration scheme, is with n such simple one degree of freedom equations.

And in fact, I don't even need to deal with n such equations, for the analysis purpose, I only need to deal with one equation, where I make now omega a variable and r a variable. So my analysis on this equation, using this particular time integration scheme that I want to analyze, is completely giving me all the information that I need to know when I use my direct integration scheme in the solution of this equation. For analysis purposes of my direct integration scheme, I can simply analyze my direct integration scheme using this equation and get all the relevant information from that analysis.

This is an extremely important fact. Once again, I do not actually do perform the transformation. We don't do the transformation. Only in theory for analysis purposes. I'm looking at the solution of this equation using the Newmark integration scheme. And I get all the relevant information on the accuracy, stability, et cetera of the integration scheme by looking at the solution of this equation with that direct integration scheme.

Now, even in the analysis of this equation, of course, remember-- or in the use of this equation, remember that r is available now. And I have some choice what to put in there and so on. But the important point is that what I want to identify is the characteristics of the integration scheme.

And the characteristics of the integration scheme are really quite well seen already if we look at a very simple case, namely this case. I set r deliberately equal to zero. And I look at the case where I have initial condition x, the displacements being one, the velocities being zero. this means, of course, that the accelerations are minus omega squared. If I use the Newmark method, the one that I just talked to you about, on that system of equations, then I can plot errors that reside in that integration. And what I've plotted here the percentage period elongation, pe divided by t times 100, which is this amount here that we see when we integrate the differential equation that I showed you via as the Newmark method with different delta-t over t values.

On the bottom here, we are plotting delta-t divided by t, the period t, of course, being given for the particular oscillator that we're looking at. We are starting off with one. The exact solution, I should take better now another color, the exact solution would look something like that, finishing off here. That would be the exact solution.

And what we're seeing now is an amplitude decay. It's already there. It would come in right there. That is the exact solution. We'll stop here. The period is of course this t here, as shown. This would be the exact solution for one period.

What we are seeing as an amplitude decay, a drop in the numerical solution, and period elongation in the numerical solution. And the percentage period elongation is plotted here. We notice that the Newmark method corresponds to this curve here. So with delta-t over t being 0.10. We have about-- let us just look it up here-- we have about 3%, 3% period elongation for the Newmark method. I'm also showing other techniques here, which are also implicit techniques, the Wilson Theta method with Theta 1.4 has about 5% period elongation for delta-t over t being 0.10.

Another point of interest is of course also is the amplitude decay. And that amplitude decay is plotted on this graph here, percentage amplitude decay as a function of delta-t over t again. Delta-t over t.

And notice that the Newmark method does not appear. In fact, the Newmark method has no amplitude decay. The solution is right down here.

No amplitude decay for any delta-t over t, whereas the Wilson methods and the Houbolt method is also shown here, have the amplitude decay shown by these curves. Notice that the Wilson method has less amplitude decay then the Houbolt method. The period elongation for the Houbolt method was also shown on this view graph here.

So the Newmark method does have period elongation, but it does not have any amplitude decay in this solution. I should also mention this, of course, as one important point, I forgot to mention that I'm looking here at the Newmark method with delta 1/2 alpha equal to 1/4. This is the constant average acceleration method.

Well, let us look then at how can we actually perform the solution? How can we select an appropriate time step delta-t? Here we notice that the dynamic load factor that we are knowing that we can derive for the solution of this differential equation here-- here I included damping-- the dynamic load factor plotted along here as a functional of p over omega. p being the exciting free frequency, omega being the frequency of the system. The dynamic load factor looks as shown here. If we have no damping, you're talking about this curve.

Now notice the important point. The dynamic load factor, of course, is 1 for static response, and it really gives us a magnification of the static response. In a static response, we would have this solution. And if the p over omega value is say 2, we would have that maximum solution in the dynamic response analysis. If p over omega is equal to 0.5. we would have, for zero damping, this maximum solution in the dynamic response history. Or rather this would be the maximum dynamic response that you would see in the dynamic solution. This is how the dynamic load factor is defined.

Of course, at p over omega equal to 1, we have resonance. In practice, we have always some damping, and if you look at this amount of damping, this would be the maximum with p over omega equal to 1, we would have this as the maximum response that is ever measured in the dynamic response solution.

Well, if we look at this graph, we recognize of course that as p over omega goes against 0, in other words as the exciting frequency, as over the frequency of the system goes against 0, we really talk about static conditions. As p over omega is very large, we are talking about no response whatsoever. And the important dynamic response really lies somewhat in this region here with that being the tail end of it. But the important dynamic response really lies in here, that we have to include in the solution.

Out here we have basically a static response. Let us look at two simple cases. For p over omega 0.05, we can see that the static response here shown in a dashed line, the dynamic response showing as a solid line. They're basically the same. In other words, a dynamic response being equal to the static response, practically equal to the static response. We are down at this range here.

As another example, if we have p over omega equal to 3, we are down over here. Then the response, the static response is shown here in the dashed line, and the dynamic response is shown as the solid line. Very little response here at all. Because as p over omega goes to infinity, there is no response for the system at all anymore. Basically, what we're saying then is that the frequency of the applied loading is so fast, that the structure can;t just move. It just cannot move at all. And we don't have any response in the structure.

Well, the important point is that if we know the frequency content of the loading, then what we have to do in the use of the Newmark method is to calculate a delta-t small enough to pick up the loading, of course, because that is given in terms of this great steps or this great input points to the computer program, and the time step has to be small enough to integrate the dynamic response of importance accurately. The static response, the static response is contained in the solution anyway, because the K matrix is on the left hand side of the solution.

So the complete solution process therefore can be summarized as follows. Identify the frequencies contained in the loading, and using a Fourier analysis, if necessary. Choose a finite element mesh that accurately represents all frequencies up to about four times the highest frequency omega u contained in the loading. The highest frequency contained in the loading is really the cutoff frequency that will give us a time step. And that then, the time step is in fact obtained via these consideration.

Perform a direct integration analysis. The time step delta-t for this solutions should equal about 1 over 20 Tu, where Tu is 2 pi over omega u, which is the highest frequency contained in the loading. What we are saying basically is that we want to select a finite element mesh to represent the highest frequency in the loading correctly, accurately, and then we want to select a time step in the data integration solution using the Newmark method that integrates the highest frequency contained in the loading also accurately, or the response in that highest frequency contained in the loading also accurately.

And this then gives us basically a step-by-step procedure for the selection of the finite element mesh and the times step. Notice that the mesh selection and the time step selection both of course hinge on the structure and the particular frequency content of the loading. That's important.

This is the modeling procedure that we would use for a structural vibration problem. We would use this time step here in the implicit integration, because we get all the accuracy we want. In an explicit integration, we would have to use a smaller time step. The sentence down here "or be smaller for stability reasons" really refers to the explicit time integration.

In the Newmark method, we would not use that condition. We would simply use Tu, being 2 pi omega u and delta-t being 1 over 20 times Tu. So in the Newmark method, I blank this out, and this is here the delta-t that we would be using.

And just to summarize once the modeling of wave propagation problem here, here what we would do is we assume that the wavelength is Lw. The total time for the wave to travel past a point is then tw equal Lw over c, where c is the wave speed. Assuming that n time steps are necessary to represent the wave, we get our delta-t and the effective length, which I talked about earlier, of the finite element would be then given by this formula.

So the procedure here is somewhat different in that if we have a wave that looks like this say, we now select a certain number of times steps, delta-t, to pick up the wave accurately. And this being the length Lw, the time for the wave to travel past the point is tw. And delta-t then is the number of time points, or time steps, that are used to represent the wave accurately.

And once we have the time step, we select the effective length of an element. And then, of course it depends on whether we have a low-order element or a high-order element, that will give us then the procedure of how to apply the effective lengths criteria.

I like to finally now summarize the complete solution process. This is a table that I put together for you to show how the complete solution process for explicit and implicit time integration can be very effectively implemented in a computer program. We form a linear stiffness matrix K, the mass matrix M, the damping matrix C, whichever applicable. We calculate the following constants.

If we select the Newmark method, these are the constants that you would use. Let's not go into details. These constants are given in the textbook. They are given in the ADINA manual, if you are using the ADINA computer program. They are widely published. We simply cast them in this form.

In the central difference method, we would use these constants. So there is some initialization involved. Then in the central difference method, we also want to calculate our first displacement, delta-t U, as shown here. And in the Newmark method, we don't need to do that, because we simply march ahead from time t to time t plus delta-t. You only need the condition at the conditions at time t to obtain the solution at times t plus delta-t.

In the central difference method, remember we need two previous values. We then form an effective linear coefficient matrix. In implicit time integration, this is the coefficient matrix. The stiffness matrix being there, as I showed to you earlier. The mass matrix and damping matrix being there.

If we look at this a0, we notice that a0 is 1 over delta-t squared and there's a constant alpha here. The Newmark constant. But it's 1 over delta-t squared. That is the important point. So as delta-t gets larger and larger, larger and larger, this effect goes against 0. The same holds for the a1 constant. As delta-t gets larger and larger, this effect goes against 0.

And this is the reason why the Newmark method really goes down or collapses to the static analysis, if we have very large time steps. In explicit integration, we use this coefficient matrix. And now for each integration step, or for each time step, we perform the following calculations. We calculate the t plus delta-t R hat. And implicit time integration, the Newmark method, in explicit time integration, we calculate this vector.

Notice at tF here, no k times tU. k times 2u is equal to tF. In implicit time integration, having calculated for each time step this load step, we now use the equations K hat times t plus delta-t u, equals t plus delta-t R hat to calculate t plus delta-t u. Of course, we would triangulize the K hat matrix prior to going into this loop. We do it once and forward.

In explicit time integration, we use this matrix or this vector here, I should say, this vector now, and we can calculate directly for the t plus delta-t U also. The two equations that I used are summarized here. In implicit time integration, this is the equation that we use.

Notice we only need to factorize this K hat matrix once. We do it outside the step-by-step integration, as I pointed out earlier. And then this gives us the increment in displacements. In explicit time integration, this is the equation we use. If M and C are diagonal, it's a trivial operation as I pointed out earlier.

With this then, we have the new displacements. And on the last view graph, I show you the formulae that we are then using to calculate the accelerations in the Newmark method from the incremental displacements and variables that we knew already. Having now calculated the accelerations, we can get the velocities at time t plus delta-t. And of course, we have already our displacement at time t plus delta-t. In the central difference method, we use our basic equations that we already applied earlier in the development of the method to calculate the velocities and the accelerations at time t.

This completes the procedure that is effectively used in a computer program for the integration of the dynamic equilibrium equations. I like to just now mention that I've placed heavy emphasis on the Newmark method in an implicit integration. Of course, there are other implicit integration formulae available, the Wilson method, the Houbolt method, and so on.

The Newmark method has shown, has proven to be quite effective. And it's certainly a good example to demonstrate to you how the actual procedure is implemented and how an implicit integration scheme is being used. Other implicit integration methods would be used in much the same way.

This is all I wanted to mention to you, discuss with you in this lecture. Thank you very much for your attention.

Free Downloads

Video

• iTunes U (MP4 - 121MB)
• Internet Archive (MP4 - 121MB)

Caption

• English-US (SRT)