Topics: Use of elastic constitutive relations in total Lagrangian formulation
Instructor: Klaus-Jürgen Bathe
Study Guide (PDF)
Sections 6.6.1, 6.6.2
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 the previous lectures, we discussed quite general kinematic formulations and numerical algorithms that we use for nonlinear finite element analysis. We also referred to and used a stress strain matrix C, but we did not discuss this matrix in any depths. This is the topic of this lecture and the following lectures.
The stress strain matrix C is used in the calculation of the K matrix, and is schematically shown on this graph. And of course, the stress strain matrix is also used in the evaluation of the force vector F, as shown here. The stress strain matrix, of course, will be different for different material behaviors. It will be different for different kinematic formulations that we're using. And all these aspects we now have to discuss in this and the following lectures.
We notice that when we do finite element analysis, we have the choice of different kinds of formulations. And we talked about the formulation that takes only into account infinitesimally small displacements. But of course, material nonlinear behavior, we actually call that formation a materially nonlinear only formation, an MNO formulation. We may have or use a formulation that takes into account large displacements, large rotations, but only small strains.
This, for example, would be a total Lagrangian formulation. And we have also formulations that take into account large displacements, large locations, and large strains. And here we use, once again, the total Lagrangian formulation or an updated Lagrangian formulation. We discussed these formulations, of course, earlier in the earlier lectures.
The applicability of material descriptions is also very well categorized in these three categories. In other words, there are certain material description constitutive relations that hold for infinitesimally small displacements and small strains, and that are not directly applicable to large strains.
If, however, we have a material description that is applicable to infinitesimally small displacements and small strains, then we will see that same material description can directly be used also in the large displacement, large rotation, but small strain analysis, provided we use the proper formulation. We will discuss this aspect just now. In fact, the formulation that we are using in this category of problems is a total Lagrangian formulation.
Let us very briefly once more recall that we have the materially-nonlinear-only formulation, as I just mentioned, the total Lagrangian formulation and the updated Lagrangian formulation, which we discussed in earlier lectures. And that really kinematically, these formulations here kinematically, for the analysis of two dimensional, three dimensional solids, using isoparametric finite element formations, we include all large displacement, large rotation, and large strain effects.
Of course, whether the formulation will actually be applicable to model a particular problem, will also depends on whether we use and have the right constitutive relation for that problem. And that is an aspect that we will concentrate on now in this and the following lectures. But notice that kinematically these formations include all nonlinearities.
If we look very briefly at various material descriptions, we have here typical models that I used and examples of applications, materials, practically all materials behave elastically as long as the stress and the strains are small enough. A hyperelastic material would be, for example, a rubber material. A hypoelastic material would be the concrete material model. And elasto-plastic material would be applicable to model metals, soils, rocks under high stresses.
A creep material model would be applicable to model metals at high temperatures. When the temperatures are high, then metals are subjected very easily to creep. A viscoplastic material model really takes into account the elasto-plasticity and the creeping effect of the material and is used frequently to idealize or to represent polymers and metals.
We will talk about a number of these material models in this and the following lecture. We will not go into depths of concrete material models and rock material models, but still whatever we are talking about is applicable to a wide range of material models.
Let us first look now, in this lecture, at the elastic material behavior. Elasticity means that for a particular strain, you are given a particular stress. The history does not enter in the solution. And a linear elastic stress strain relationship is shown here.
Notice for the strain, tE, we are given a stress t sigma. And it's a unique stress. It is a unique stress. And of course, the history of the formation does not enter at all.
Notice that in this particular case, we would have that typically t sigma here is equal to E times te. te being of course this value here, t sigma there. Notice that an increment in stress is given by d sigma equal to E times de. Linear elasticity, a linear elastic material description, means E is constant.
For nonlinear, as a nonlinear elastic material description, we might have schematically in a one-dimensional analysis, this relationship here, the red curve. Notice that, once again, for a particular strain, we pick up a unique value of stress. And this relationship here shows that t sigma is given tC times te. tC is the slope of this line here.
Now, an increment in stress is obtained by this relationship here. d sigma equals C times de. Notice that C is not the same as tC. This is here, relating the total stress to the total strain. Here we talk about an increment in stress related to an increment in strain. It's important to note that C is not the same as tC. Unless, of course, we have a linear elastic relationship.
In a computer program, we would represent such a curve typically by a series of straight lines, as shown here. And the input, of course, would be these points here, and the corresponding points on the stress curve. Notice that if you take a large number of such straight segments, you can really approximate a continuous curve such as this one here quite well.
We can generalize this elastic material behavior using the total Lagrangian formulation, as shown here. Notice that the total stress, second Piola-Kirchhoff stress-- we talked about this stress measure in an earlier lecture-- is given by the stress strain law times the total strain. Green-Lagrange strain-- once again we talked about this strain measure in an earlier lecture.
Incrementally, we have this relationship. The increment in the second Piola-Kirchhoff stress is related to the increment in the Green-Lagrange strain via this constitutive law. Notice that this tensor or the components of this tensor are not the same as those components in general.
This material description is frequently employed in two types of analysis, the usual constant material moduli for infinitesimal displacement analysis. Or for rubber-type analysis, rubber-type materials, I should say, where of course this matrix is not constant. And that matrix is not constant either.
I like to now spend time on discussing these type of analyses. The use of the constant material moduli means, for example, for an isotropic material, that we use these matrices given via this relationship here. Notice lambda and mu are the Lame constants, that you're probably familiar with, and data ij is the Kronecker data, defined as shown down here. These are the same constants that are used also in linear analysis.
So no new parameters are introduced here. As examples, here we have the stress strain law for the 2-D plane stress analysis. I've written it in matrix form, and this zero here means total Lagrangian formulation. And notice that there could also be a t up here. In other words, we use the same matrix for the increment in the stresses and for the total stresses.
So I could very easily put a t up there, and we notice then that the last row corresponds to this relationship here. It's important to recognize that on the right hand side here, we have the Green-Lagrange strain, t0 zero epsilon 12 plus t0 21. In other words, we take the total sharing strain into account here. And that gives us then this row and that column, of course.
For 2-D, access symmetric analysis, we would use this stress-strain law, once again, for the increments in the stresses and for the total stresses. For an orthotropic material, we would typically have this stress strain law in plane stress, plane strain analysis. And there you see again, the zero.
Notice Ea, Eb are the material and nu ab are the material moduli corresponding to the a and b, directions and similarly the sheer modulis Gab. The material being orthotropic means Ea is not equal to Eb. This once again is the material law that we are quite familiar with from linear analysis.
Now let's look at the sample problem, the sample analysis. And I would like to choose a very simple problem to just demonstrate some of the aspects of using a total Lagrangian formulation with a constant material law. Here we have a bar, with cross sectional area A bar. We use the material constants E and nu, the ones that we are used to you in liner analysis.
And we want to use this stress strain relationship. We assume that the area A bar remains constant when we pull the bar out, or when we push it in. This E curl here is given via this relationship here, E being the Young's modulis nu being the Poisson ratio.
In tension, the bar would be pulled out as shown here. Notice the area, the cross sectional area, remains constant. And the displacement is measured as t delta. In compression, the bar would be pushed into, as shown here.
Notice here now t delta and the bar has been pushed in as shown by the blue outline. If we now apply this stress strain law, to the problem and develop a relationship between the force required here to pull the bar out or push it in, by delta, we obtain the following.
Here we have the Green-Lagrange strain. And of course, the right hand side involves derivatives that we have been discussing in the previous lecture. Notice t0u1,1 is nothing else than this relationship here. It's given by the current length and the original length of the bar. Substituting into here, we get directly this right hand side.
We also remember that we can calculate or we have a very clear relationship between the Cauchy stress, which is the force per unit area, which we are interested in calculating here. And the second, Piola-Kirchhoff stress, this is the inverse of the deformation gradient.
We substitute for here, the Cauchy stress of course being simply the force divided by A bar, the inverse of the deformation gradient being given by this relationship here. And therefore we can evaluate the second Piola-Kirchhoff stress in of the externally applied force, the physical force, divided by the constant area, cross-sectional area, A bar, and the original length and the current length.
We therefore have now the Green-Lagrange strain and the second Piola-Kirchhoff stress, in terms of current lengths and original lengths and force applied-- and of course the cross section area goes in there. Substituting these relationship and that relationship into our stress strain relationship of the previous view graph and using this fact here, we directly obtain this red curve in the-- which gives us the force applied to the bar as a function of the displacements at the end of the bar.
Now notice that this is a nonlinear curve-- the expression is actually given here-- and that it is quite unrealistic, certainly in this region here. It's not a realistic material description for large strains.
You see here at this point, we have a 0 force to be applied when the displacement, the end of the displacement, is minus the original lengths off the bar. Quite unrealistic, of course. And therefore, we recognize that really this use of a constant E and nu, going in here, constant Young's modulus and Poisson's ratio, going in here. It's really not a proper or realistic way to proceed for large strains.
The usual isotopical and also orthotopic material relationships, meeting constant E, nu, Ea, et cetera, the way I just showed it on a few previous view graphs are really almost exclusively employed in large displacement, large rotation, but small strain analysis.
We should now recall one very important namely that the components of the second Piola-Kirchhoff stress tensor, and of the Green-Lagrange strain tensor are invariant very under a rigid body motion of the material. And this is most important because we use this fact to apply material descriptions that are applicable to small strain analysis directly also to large displacement, large rotation analysis. You see, E nu, Ea are material constant that applicable to small strain analysis. And we can directly them some also in large displacement, large quotation analysis, but small strain analysis using this fact.
And this means-- another way of saying it is that the actual straining increases only the components of the Green-Lagrange strain tensor through the material relationships. Let's look at this equation here. We have here the nodal point force vector being equal to the stress strain-- the strain displacement matrix times a stress vector, the second Piola-Kirchhoff stress vector.
Now this vector is invariant under rigid body rotation. And of course, this vector is calculated using the constitutive relation. This strain displacement matrix takes into account the rotation.
Let's look at this fact a little bit closer. Here we have pictorially an element, a four-node element say, that originally is in this configuration, the black configuration, and we are pulling it out and sharing it. If this length here is 1 originally, and that length here is 1 originally, then this distance here corresponds to the Green-Lagrange strain in the 1:1 direction, in the one direction. In other words, it's equal to 1 0 epsilon 1 1. Similar, this distance here is equal to 1 0 epsilon 2 2, and the total sharing is equal to the sum of these two components here.
This is the deformation that brings the element, say, from the original configuration to state 1. Now, going from state 1 to state 2, we rotate the element rigidly, as shown here. Then we know that this distance here is still 1 0 epsilon 1 1. Same as what we have had here. This distance here is the same as what we have seen here. And the total sharing is still the same as what the total sharing was here.
Notice the Green-Lagrange strain components corresponding to the stationary coordinate axes. 1 and 2 have not changed. Because here this 1 and that 2 refer still too the stationary coordinate frame.
On this view graph, we summarize this information once more for small strains, we have of course this relationship here, which means that these strain components are much smaller than 1. And we can directly say that the second Piola-Kirchhoff stress in configuration 1 is given by the Green-Lagrange strains in configuration 1. Notice the second Piola-Kirchhoff stress is approximately equal to the Cauchy stress in state 1.
Now since state 2 is reached by just a rigid body motion, we find that this holds, as I already pointed out on the earlier view graph. We also find that this holds. And of course, the Cauchy stresses in state 2 are obtained by a rotation of the Cauchy stresses in state 1.
But the second Piola-Kirchhoff stresses and Green-Lagrange strains, the components of these two tensors are invariant under the rigid body motion. And that means that we can directly apply material relationships that are applicable for small strains and infinitesimal displacements. Those that we are, say, already familiar with from an MNO analysis. We can directly apply these to large displacement, large rotation, but small strain analysis, by simply using Green-Lagrange strains instead of engineering strain components, and second Piola-Kirchhoff stresses instead of the engineering stress components.
Applications are really all large displacement, large rotation, small strain analyses of beams, plates, and shells. We can use for this analysis frequently simply the 2D and 3D elements that we discussed earlier. But we may want to use beam and shell elements, and those we will discuss in later lectures.
We also have a very large application area here, namely, the linearized buckling analysis of structures. We have already performed the analysis of shells, or shown how we analyze shelves and beams using continuum elements. You might recall it from the earlier lectures.
Here you have just schematically how we would model say a frame using plain stress elements. And here you should see schematically how we are modeling using axisymmetric continuum elements and axisymmetric shell. Here we use 3D elements to model a general shell. But as I pointed out, we will talk in later lectures about actual beam and shell elements that are frequently much more effective in the analysis of beam and shell structures.
This completes the first large item that I wanted to talk to you about, namely, the use of constant material moduli, E, nu, Ea, Eb, et cetera in large displacement, large rotation analysis, but small strain analysis. We found that the applicability of constant material moduli to large strain analysis is really very questionable.
For large strain analysis, we use frequently a hyperelastic material model. And that is very widely used to analyze rubber materials. Here we have the relationship that is basic to this description. We have the Piola-Kirchhoff stress, given in terms of the partial of the strain energy density function, W, this respect to the Green-Lagrange strains. So we get one component here, by taking this t0W and differentiating it with respect to the corresponding component of the Green-Lagrange strain. And of course this right hand side gives us this tensor t0Cijrs.
Incrementally, then, we have an increment in the second Piola-Kirchhoff stress is given by this stress strain law here, involving a tangent material description and an increment, a differential increment, in the Green-Lagrange strain. This 0Cijrs is nothing else but the second partial derivative of W with respect to these Green-Lagrange strain components.
W is a strain energy density function per unit original volume. And of course that has to be defined for the particular material that you're looking at and will involve certain material constants. For rubber, we generally use this kind of description here, where W is a function of I1, I2, I3, and the Is are the invariants of the Cauchy-Green deformation tensor. We discuss the Cauchy-Green deformation tensor in an earlier lecture. Of course, it's given in terms of the deformation gradient.
So I1 is equal to this. I2 is equal to this relationship here on the right hand side. And I3 is simply the determinant of the Cauchy-Green deformation tensor.
For example, one material law that is quite widely used is the Mooney-Rivlin material law, where two constants enter here. And of course, we have here I1 and 12. These are material constants that have to be determined depending on what rubber you're dealing with.
And the incompressibility constraint is imposed onto the formulation by having I3 equal to 1. Of course, rubber is an incompressible material, at least almost an incompressible material, and is modeled as an incompressible material. And so we want to have I3 equal to 1.
In general, if you deal with incompressible materials, it is effective to not use the displacement-based finite element formulations the way I've been discussing it in the previous lectures, but to amend that formulation, you should really extend that formulation to have an effective formulation to deal with incompressibility in the analysis. Except for the plane stress case, we can directly use what we discussed earlier.
The formation for incompressibility, or for analysis of incompressible materials, such formulations of course are very valuable but the discussion of such formulations really deserves another lecture.
Well, we like to now show how this material law is used in the plane stress case. There we can directly employ the formulation that we discussed in the earlier lectures. Time 0, for example, we would have a piece of material as shown here. At time t, this piece of material has deformed, as shown here. Notice plane stress of course means stress is 0 through the thickness.
For this two-dimensional problem, the Cauchy-Green deformation tensor is as shown here. 0 components here. And we imposed the incompressibility of the material by simply setting the determinant of this matrix equal to 1. And this gives us then a constraint on T0C33.
Basically, we are saying physically that we constrain the material to shrink into this 33 direction such that for the applied stresses in the plane of the material, the material acts as incompressible. And that's what's being done here mathematically. We can now evaluate with this constraint I1 and I2. Here you see the expressions. And if we have I1 and I2, we can go back to the material law, substitute and directly octane obtain our second Piola-Kirchhoff stresses, in terms of the Cauchy-Green deformation tensor components.
And of course, we also have a relationship between the Cauchy-Green deformation tensor components and the Green-Lagrange strain components. Those we use. That relationship we use, and that actually brings into it this relationship here. And if you performed then the differentiation, the final result is shown on this view graph. This is the stress strain relationship that gives the total stresses in terms of the Cij components.
You similarly proceed to also calculate the tangent relationship, the tangent stress strain law. As I pointed out earlier, you need to now take second derivatives of W with respect to the Green-Lagrange strain components. And if you also use the relationships between the Green-Lagrange strain components and the components of the Cauchy-Green deformation tensor, you directly arrive at this expression.
And of course now, you would have to substitute as we did earlier for calculating the total stress strain law, in other words, the stress strain giving us a total stress as a function of the total deformations. You proceed in the same way, and obtain the tangent stress strain law.
Let us now look at an example. In this example, we are considering the analysis of a tensile test specimen. The specimen is shown here. All dimensions are given in millimeters. The thickness of the specimen is one millimeter, and the Mooney-Rivlin constants are given here. We need to model only one quarter of the specimen because of symmetry conditions. And in the analysis, we used 14 eight-node elements to represent this quarter of the specimen.
Notice that on the right hand side, we constrained all material particles to move together and only move to the right, as shown by delta here. This represents 1/2 the gauge length of the specimen, and one is particularly interested in seeing the deformations in that gauge length.
This view graph now shows the calculated response. The applied node is given here. The extension is given here. This is the gauge response. And that gives us-- this curve gives us the total response. Of course, quite nonlinear, because the material is a nonlinear elastic material.
The deformed mesh is shown on the next view graph here. At the force 4 Newton, we have this deformed mesh here shown in red. And notice here, we have the undeformed mesh in dashed black line. Notice that the displacements are large, and the strains certainly are very large in this problem.
Finally I would like to share with you some analysis results that we obtained in the analysis of a rubber sheet with a hole. And these analysis results are shown on the slides, so let me walk over here so that we can share that information.
Here, we show the rubber sheet. It's a square sheet, 20 inch by 20 inch, that contains a hole of diameter 6-inch. Notice the sheet is subjected to a force on the left hand side, the distributed pressure, uniformly distributed pressure on the left hand side and on the right hand side. The mass density is given over here also, and the thickness off the sheet is given here as well. We will use the mass density later on when we perform also or look at the results of a dynamic analysis.
The next slide now shows the finite element mesh that we used to model one quarter of the sheet. We only needed to model one quarter because of symmetry conditions. Notice we used four-node elements in this mesh, and that the distributive pressure on this right hand side is modeled by concentrated loads at these nodes.
The next slide shows now the solution results that we obtained. Notice down here, we have a little key for the points that we want to look at. Here is point B. Here is point A. And here is point C. The graphs here shows the displacements at point B, at point A, and point C. In each case, the W displacement, which at point B is into this direction, at point A is into this direction, but at point C is downwards.
Notice the very large displacements, W here, going up to 10, 11 inches for the load applied. Remember the sheet is a square sheet of 20 inches, so we have very large deformations here. We are computing our solution results obtained in the total Lagrangian formulation, with the hyperelastic Mooney-Rivlin material law, with the results obtained by Iding.
In fact, there is almost no difference between our results, and the results of Iding, so that we could not plot the difference. And these curves here refer to our results and the Iding results.
The next slide now shows the deformations that the sheet has undergone at maximum load. Here is the original configuration of the sheet, the initial configuration. And here, we show the deformed configuration of the sheet. Notice that the sheet has undergone very large displacements and very large strains, particularly of course here in the region of the hole.
The final slide then shows the dynamic response that the sheet undergoes when we subject the sheet to a step loading, where the static loading was applied before. Notice the step loading is applied over 10 delta-t where delta-a is 0.0015 seconds. It goes up linearly to 75 pounds per square inch, and then it stays constant.
Once again, we're looking at the displacements at point B, point A, W displacements here, and at point C downward displacement, again denoted as W. Here you see the solution response. Point B, point A, and point C. In this particular analysis, we used first the Newmark method with data equal to 1/2 and alpha equal to 1/4. In fact, that is nothing else than the trapezoidal rule, which we discussed in an earlier lecture.
We also performed the analysis once with the Wilson theta method, theta being 1.4. This method we did not discuss in the earlier lecture, but it is described in the textbook. Notice that the results obtained using the Newmark method and the Wilson theta method are really very close to each other.
It's also interesting to see the static response that the sheet undergoes for this load level at the points A, B, and C, and that is given by these horizontal lines. Once again, you have here a large displacement and very large strain problem.
This brings me to the end of what I wanted to discuss with you in this lecture. Thank you very much for your attention.