Topics: The two-noded truss element - total Lagrangian formulation
Instructor: Klaus-Jürgen Bathe
Study Guide (PDF)
Problems 6.16, 6.17
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 nonlinear finite element analysis of solids and structures. In this lecture I would like to discuss with you the total Lagrangian formulation of the two-noded truss element.
In the previous lecture, we considered the updated Lagrangian formulation of the same element. We started from the general continuum mechanics equation, and obtained the governing finite element matrices. I would like to do the same now for the total Lagrangian formulation. And we will see that to find element matrices are indeed the same that we obtained in the updated Lagrangian formulation.
Here we have the governing equation of the principle of virtual work that we derived in an earlier lecture. This, of course, is a governing equation that we apply to obtain the governing equation, finite element equations, for two or three dimensional analysis, and also for the truss element. It's a very general equation that we discussed earlier in great detail.
We use this equation now for the truss element. And the first point that we have to look at is what stress measure and corresponding strain measure do we need to deal with? Now for the truss, it is important to realize that the only non-zero stress is this stress here, t0 S11. In other words, only the 11 component is non-zero. Let us look at this fact in detail.
There's first of all, what we might call a mathematical explanation. Let us assume, of course, that the cross-section here is constant, the same way as we did in the discussion of the updated Lagrangian formulation. And here we have the element in its original configuration, and the element has undergone large displacements and large rotation to get into this configuration at time t. Notice the length is almost constant. We have an epsilon there. That, of course, stands for the strain that the truss has undergone, but that epsilon is assumed to be very, very small.
We will later on also look at a similar picture, but this red truss element at time t will have to be moved into this configuration here, starting, in other words, starting on the left-hand side here, we simply move this truss element here straight back so that this node lies on top of that node. In other words, if I just show it to you like that, it will have rotated like that.
Notice that as far as the derivation of the stiffness matrix of the truss element in this configuration, and as far as the element force vector derivation are concerned, it does not matter whether we look at this particular truss element in this configuration, or whether we move it back to make this point coincide with that point.
Well if we drive using what we have discussed earlier the deformation gradient for this particular element, we recognize that this is the general expression. Notice that X, of course, is equal to R times U, R being the orthogonal matrix, and U being the stretch matrix. Notice that in this X matrix now, this here really is the R component, and the U matrix is nothing else but the identity matrix with an epsilon added to the 11 component. And the product of R times U gives you then, of course, this right-hand side.
Since the truss carries only the Cauchy stress, tau11, corresponding to the tensile force in the truss, we can directly write down the Cauchy stress in the stationary coordinate frame as given here. This, of course, is nothing else but a Mohr's circle transformation. Notice that our truss element moves through the stationary coordinate frame. I pointed out earlier very strongly that we always deal with stationary coordinate frames.
So this is the expression of tau, the stress in the truss, but written now in the stationary coordinate frame. So we have sheer components and normal components. Of course, if theta's is equal to 0, then we have only here a 1, and all of these are the components drop to 0.
If we now use our general formula, that once again we discussed earlier, to evaluate from Cauchy stresses, second Piola-Kirchhoff stresses, we directly obtain this result here. Since epsilon is small, we set this equal to 1. And we have, indeed, that's the only non-zero stress component in this stress tensor is the S11 component.
Well there's also a physical explanation. And that physical explanation goes as follows. We have initially, say, this situation here in the configuration at time zero where 0S refer to the configuration at time 0 is equal to tau0, and that is equal to nothing else than the initial force in the truss divided by the area, 0 components, of course, here.
If we now think of first extending that truss element by a very small amount as shown here, notice that blue truss element here corresponds to the extended truss element. It lies really on top of the black, but in separating it slightly out so that you can see the different coloration. Then this blue truss element can be sort of corresponding to, and I'm showing it now up here, equal to a time T star, which is really a conceptual time. We are pulling out the truss first, and this means of course, corresponding to that T star conceptual time, we can write down this equation here.
Surely there are 0 components here. Now, we start rotating the truss about the left node, and the result is shown on this U graph. Now we have the RAD configuration, theta is here, and the stress tensor that we are talking about has not changed. Why has it not changed? Well because of this statement down here, the components of the second Piola-Kirchhoff stress tensor do not change doing a rigid body motion, in particular doing a rigid body rotation.
We discussed this statement earlier in an earlier lecture. It's a very important property of the second Piola-Kirchhoff stress tensor. And this means, of course, that the 0s which we had here in blue remains 0's here in red. Very important point.
Well having identified then that S11 is the only non-zero stress component that we really need to deal with in the general principle of virtual work, we can directly write down this relationship here. It's nothing else than the application of the general equation, which we looked at earlier, reduced to the special form of the truss element taking only the S11 component into account. See here, S11, of course, corresponding strains-- eta11, and here we have only the 11 components as well.
The next step is to look at what are the individual terms in that equation. Recognize that S11 is nothing else than tP over A. C1111 is equal to E. The volume, of course, is given as shown here. And of course, we also recognize that the stress and strain states are constant along the truss. Once again, because we're dealing with two nodes only, and we can only have a lean interpolation of displacement between these two nodes, meaning that the strains must be constant, and the corresponding stresses must be constant as well.
We substitute these things into, these relationships up here, into the equation that I showed you earlier and we directly obtain this equation. Notice that in this equation now, we want to, of course, express the strain terms e11, eta11 in terms of displacements, nodal point displacements. And that is achieved as usual via a strain displacement matrix. This is the relationship here corresponding to the linear strain increment, and that is the relationship corresponding to the variation on the nonlinear strain increment.
Notice again, and I've pointed it out a number of times already earlier, that we write this right-hand side in such way that the right-hand side is equal to the left-hand side. Therefore, these BNLs are defined such that the right-hand side is equal to the left-hand side.
Also remember, u hat here, this hat up there, the nodes, the [? discrete ?] nodal point displacements are listed in this vector. And here they are. Notice upper superscript refers to the nodal point. Lower subscript refers to the directions, the coordinate directions that we are looking at. Similarly, of course, for node two, 2 up there, and 1 and 2 down here.
We go back to our general continuum mechanics relations that we discussed earlier, and identify that the total increment in the Green-Lagrangian strain is given via this relationship here. Notice that this relationship is obtained from the general equation of this epsilon 0 IJ. Of course, I and J are now equal to 11 for the truss because S11 is the only non-zero stress that we just have been considering earlier.
So we'll go into our earlier view graphs, the information that we discussed earlier. To look up the equation on 0 epsilon IJ, we set I and J equal to 1. There is a k also in that equation if you look back. That k goes from 1 to 2 now because we're looking at the two dimensional motion of the truss. If you have a three dimensional motion, then, of course, we would have to let k go from 1 to 2 to 3. Here we have k only going from 1 to 2, that's why you see the 2 at the highest subscript here. And automatically, you would get this equation.
Notice that the linear part in here gives us the linear strain increment. And that the nonlinear part in here gives us the nonlinear strain increment. And if we take the variation on that part, we directly obtain this relationship here. We like to write this now in matrix form, and that is achieved as shown down here. Notice this part, of course, is nothing else than what we have up here written in matrix form.
There is, of course, in the total Lagrangian formulation the initial displacement effect. And that initial displacement fact lies in these derivatives here. This one and that one. And we need to evaluate them. We can evaluate them from kinematics, as shown in this picture here. Notice that we have moved the truss now, as I pointed out earlier. We would do to the nodal point to coincide such that, I should say such that the nodal points or nodal point one has not moved. Once again, that is quite OK. If we want to derive the stiffness matrix and force vector of the element, because those quantities are not affected by this rigid body translation.
Here we have theta, and notice these are displacements that we call delta tU2, delta tU1 with a negative sign there because this displacement would be measured positive in the opposite direction. And we can directly see that this differentiation here is nothing else than this value right here. From the kinematics. Similar, this differentiation is nothing else than sin theta, once again, from the kinematics.
So we have the initial displacement effect in these derivatives, and can evaluate those as shown by these functions here. We also want to, of course, evaluate these derivatives here, 0U11, and in fact, 0U2 comma 1. And we recognize that, as I pointed out earlier, the strains are constant. Therefore, this is a constant value, and in this differentiation it's really nothing else than this relationship here.
We take the difference in the nodal point displacement, and divide by L. Quite the same way as we did it for the updated Lagrangian formulation. You can, therefore, write down this differentiation, these two derivatives here, as shown on the right-hand side.
Now, we are ready to take these expressions and substitute into the right-hand side of the linear strain term. Notice that this gives us one term here, and another term there. Let's look at these terms. This term here comes from that part there. And this term here, down here, comes from that part there. Let's look at this second term a little closer. This here is, of course, the initial displacement effect. All of that I should say is really the initial displacement effect. But if these terms are 0 here, then, of course, initial displacement effects drop out.
So therefore, these two terms here are the important ones, and we can see that they are written in this vector here. We just derived them earlier. Notice that this matrix, this vector here, is given by this matrix there multiplied by the nodal point displacement vector. And of course, this total here gives us the initial displacement effect, which is captured in t0BL1.
The sum of t0BL1 and t0BL0, of course, gives us a linear strain displacement matrix. And that summation is performed on this view graph here. Notice we have written these terms once again. Of course, if you had to multiply out in order to get t0BL1, the sum of these two gives us this strain matrix, the linear strain displacement matrix. And we now notice that in the linear strain displacement matrix, cosines and sines of thetas of the angle theta appear. That is something to keep in mind.
For the nonlinear strain displacement matrix, we look at our nonlinear variation on the nonlinear strain term. Here it is on the left-hand side. And we substitute for the differentations the way we have been doing it for the updated Lagrangian formulation of the truss. And we recognize, of course, that this part here is given by what is underlined as blue. And this part here is given by, again what's underlined in blue. And then this part here is the nonlinear strain displacement matrix, the BNL. That part right there.
Well if we now use all this information, we directly obtain the stiffness matrix and stiffness matrices, and force vector off the truss element. We substitute into here from what we derived earlier, and directly obtained this expression.
And what's underlined here in blue is, of course, the linear strain stiffness matrix of the truss element. And notice that this linear strain stiffness matrix is identically equal to the linear strain stiffness matrix of the updated Lagrangian formulation truss element.
We use this term to obtain the nonlinear strain stiffness matrix. And we notice that this nonlinear strain stiffness matrix is identically equal to what we derived earlier in the updated Lagrangian formulation as well.
The force vector, the F vector for the truss element is given on this view graph and is obtained from this relationship here. Notice it's, once again, identically equal to the force vector of the updated Lagrangian formulated truss element as well. And so we notice that to find all the matrices that we obtain in this total Lagrangian formulation are the same matrices that we already had derived in the updated Lagrangian formulation.
Important point to recognize is that the coordinate transformation used in the U.L. formulation is contained in the initial displacement effect of the T.L. formation. Let's look at that a little bit closer. What this means is that in the updated Lagrangian formulation, we did the evaluation corresponding to the configuration at time t with a local coordinate system. And then after that derivation was performed, we transformed the nodal point variables to the global coordinate system. In the total Lagrangian formulation, we obtain the matrices directly corresponding to the global coordinate system through the initial displacement effect.
Now, the same can also be shown, of course, for other elements, and we have looked here, particularly at a two-noded beam element. And we have asked ourself the question then which method or which formulation is always more effective? If you use numerical integration, as we might want to use for the beam element, particular we have to use if we want to capture elastoplasticity effects.
Then clearly, the updated Lagrangian formulation must be more effective because what you are doing in the total Lagrangian formulation is to carry these cosines and sines theta terms along for every integration that you're considering. In the numerical integration, the product of B transpose CB, with the B containing the initial displacement effect now, is much more expensive than the product for the updated Lagrangian formulation where the B does not contain this initial displacement effect.
So in the updated Lagrangian formulation, the B matrix is much simpler. We do the numerical integration with a simpler B matrix along the beam element, and the details are given in this paper. And once we have done the numerical integration, we simply transform the nodal point displacements from the local coordinate system. The [? coord ?] system, we call it a [? coord ?] system when we discuss the updated Lagrangian formulation of the truss element, and we transform these nodal point displacements to the global system only once by this t transpose kt product.
So this means then that certainly for a two-noded beam element, the total Lagrangian formulation is not as effective as the updated Lagrangian formulation, and it is very natural to use the updated Lagrangian formulation because you get the same result in the end. For the two-noded truss, of course, the numerical operations involved are so close, so close, that one can really say the updated and total Lagrangian formulations are basically the same.
Let's now look at an example. Here we have a truss, a simple truss structure. And we perform the analysis of this truss structure using the U.L. formulation, the T.L., total Lagrangian formulation would give the same result, as I just discussed. And we want to particularly calculate the collapse load of this truss. We also want to test the MNO, the Materially Nonlinear Only formulation.
The reason for testing this formulation is to identify what differences do occur, and to understand more clearly the U.L. Formulation, as well as the MNO formulation. Notice that for this truss, the [? geo ?] matrix, the geometric data are given here. The material data are given here in this table. Notice that, of course, these nodes dawn here are fixed. Notice that we are putting a load that pulls on the truss, P, as shown. And that we have two degrees of freedom, U and V. So it's a fairly simple problem, yet there is quite a bit of meat in it when one looks closely at the solution. And that's what we like to do now.
We calculate first analytically, the limit load, the elastic limit load. In other words, the load for which the yield point is reached in one of the members. And in fact, for this particular load, the side trusses just become plastic, and this is the value of that load. The ultimate limit load, when the side trusses are plastic and the center truss is plastic as well, is reached at this level.
If we use an automatic load stepping incrementation, we will discuss that in one of the next lectures, we obtain this response. Notice P plotted vertically here, and V, the displacement corresponding to P, plotted horizontally. Notice that we come up here in the step-by-step solution and march along here.
There are two solutions given, actually. One with v1, very small, and one with v1 larger. This is here the parameters that we use to start a solution algorithm. We have to start the solution algorithm by imposing a small value of displacement to a node, and then the solution algorithm can perform automatically the total step-by-step solution. We will discuss that in one of the next lectures. And as you see, with the larger value here, we get these discrete points. No difficulty solving for this response. Notice also that there's no difficulty going beyond this analytic elastic limit load, shown here by the blue line.
If we were to look very closely at this slope here, we would see that it is slightly positive. You can hardly see it from this graph. This is the solution obtained with the U.L. formulation, with the U.L. formulation. And once again, the T.L. formulation would give identically the same result.
We now consider an MNO analysis. MNO means Materially Nonlinear Only. No large displacement effect included. No geometric stiffening effect included. Geometric stiffening, meaning the effect of the nonlinear strain stiffness matrix being included.
So we want to use this formulation, and once again, apply the automatic load stepping incrementation. If the stiffness matrix in that automatic load step incrementation is not reformed, in other words, we keep the original stiffness matrix, then almost identical results are obtained to the U.L. or T.L. formulation solution results.
However, if the stiffness matrix is reformed in the MNO formulation, and for a load level larger than the elastic limit load-- remember, the elastic limit load is quite a bit lower than the ultimate limit load-- then the structure is found to be unstable, as 0 pivot is obtained in the stiffness matrix.
I like to ask the question why is that so to you now? And ask you to just think a little bit about it. A few minutes. And then when I come back, I give you the explanation by looking at the problem once again.
I'd like to now address the question that I asked you before the break. Namely, why, if we use the MNO formulation do we get virtually identically the same results as in the U.L. formulation if we only set up a stiffness matrix that corresponds to a configuration below the elastic limit load.
If on the other hand, we setting up a difference matrix above the elastic limit load in the MNO formulation, then we cannot, so it's a problem. Well the answer is that as soon as the two trusses on the side have become plastic, the MNO formulation does not have any more stiffness corresponding to this degree of freedom, meaning the matrix is singular, and we cannot solve the system of equations.
Whereas if we look at the U.L. formulation, then we see that the geometric stiffening effect, or the KNL effect, the nonlinear strain stiffness effect, introduces the forces that are acting in these two elements into the structure into the structure assemblage as the stiffening part. And of course, is a KNL matrix corresponding to this element as well. So in other words, these three KNL matrices corresponding to these three elements, these three truss elements, introduce a stiffness into the total structure assemblage, meaning that we can still solve the system of equations.
We also notice that these KNL effect that we are carrying along in the U.L. formulation, when measured on the KL effect are rather small. And this means that the response of the U.L. formulation calculated with the U.L. formulation is basically the same as the response calculated with the MNO formulation. As long, of course, we use in the MNO formulation a stiffness matrix that corresponds to a configuration below the elastic limit load.
Well this completes this example. Let us look at another example. And in this example, I like to consider with you a pre-stress cable, as shown up here, that has initially a span of S. And we then prescribe on the right-hand side a displacement to this node to bring this node towards the left, meaning that the cable will sag through. We want to calculate the deformations of the cable, and in particular, we want to have the final deformations when S is equal to 30 meters.
Notice that the cable carries an initial tension of 500 Newton. And here we see the other data corresponding to the cable. This, of course, is a geometrically nonlinear problem. We want to use here the U.L. formulation. And the interesting part is that to solve this problem we have to use relatively small time steps, or load steps, because they're considering it to be, of course, a static solution.
We need to only allow small perturbations in the nodal point coordinates, because if we do will allow too large perturbations in the nodal point coordinates, we find that the out of balance loads at the nodes become very large and we have difficulties solving the equation. So we use many load steps with equilibrium iterations so that the configuration of the cable is never far from an equilibrium position.
In this particular problem, we used full Newton iterations without line searches. We will talk about this scheme in detail in the next lecture. We use also convergence criteria because we're iterating, we have to break up that iteration or stop the iteration at some point. And we will use these two convergence criteria.
We consider these convergence criteria in quite some detail in the next lecture. Notice that basically we are saying is that the out of balance loads shall be small enough, and here we have an energy criterion that also has to be-- this measure has to be small enough for convergence. Once again, we will talk about these in detail when we consider the solution of the equation.
Here we have a table that summarizes the results. In the first time step, puts the gravity loading onto the cable with a mass density that was given, of course, on the earlier view graph. And we notice that with this iterative scheme, we need to 14 iterations.
We then from time or load step 2 to 1,001 prescribed the displacements in 1,000 equal steps on the right-hand side of the cable. And we notice that we need less than 5 or equal to 5 iterations per load step. This might appear to be a fairly large number of load steps. Of course, one could now experiment and bring this down. If you try to bring this down, this number, of course, would go up.
And one could experiment around. But it shows what kind of considerations are necessary when you solve this problem, namely that you don't want to be too far away from the equilibrium configuration force of the cable in any one of the load steps, and if you are, then, of course, you will need a large number of iterations here, and your iterative scheme might not even converge.
Pictorally, the results are shown on this view graph here. Notice we have the undeformed cable here. And then with S equal to 55 meters, the cable takes on this configuration with S equal to 30 meters. We see this configuration. And if you were to try to bring this point in further, then to model the curvature down here, of course, properly, you would have to have a finite discretization. So this is probably the largest amount that you want to bring this point over with this particular discretization.
If you were to use a finite element discretization, this curve, of course, would hardly change. And this one here will change a bit down here. But these results are quite acceptable up to S equal to 30 meters with this coarse discretization. And we have learned a bit from this problem.
Thank you very much for your attention.