Engineering Mechanics - Questions and Answers

Some Question and Answers on FE Theory and Modeling [Updated August 2000]

Question:

WVP - Our research group doesn't have much experience with the fatigue behaviour of composite materials. Our aim is to acquire some knowledge about this subject within our department. Till now we have developed a rather simple experimental setup to follow the fatigue deterioration of plain woven E-glass/epoxy composites. Since we hope to describe the fatigue behaviour of the material with some kind of law later on, it is important to accurately determine the stresses in the composite material.

Now there are some problems : - an accurate elastic theory to describe the behaviour of plain woven composites is not developed yet. Most of the theories are micromechanical of nature and are only useful for very simple loading configurations. - when talking about deterioration, phenomenons like delaminations and interlaminar stresses are very important in determining the deteriorated state of the material. The classical plate theory does not take into account these effects. - with our experimental setup composite specimens are loaded as a cantilever, so the more common boundary conditions for plates on 4 simply supported edges don't apply to this case. As a first approach we used the following approximations :

- the plain woven material can be modelled as two cross-ply layers (0/90).
- to acquire the interlaminar stresses, I did some reading in the libraries and found about the "Generalized Laminate Plate Theory" in some articles published in "Computers & Structures" and "The International Journal for Numerical Methods in Engineering". I think this theory would be very appropriate for determining the interlaminar stresses in the plain woven composite. This is important because a damage model should need these stresses to find out whether or not deterioration will continue.

Therefore I would like to ask you :

- are there any commercially available implementations of this theory that run on a plain PC computer with Windows 95/98/NT or do any guidelines exist how you can implement it yourself in C or Fortran ? Besides the composite specimens are simple rectangular plates, so there is no need for a preprocessor or a modeller that can handle various shapes of the specimen. I just need a FEM code that implements the theory and solves the equations. And more important : is it possible to retrieve the source as well, because the degradation law for the fatigue modelling should be implemented later on in this FEM code ?

WVP

Answer

I had written a FORTRAN program for an extension of the GLPT to composite shells. My doctoral thesis was: "Layer-wise Shell Theory for Discretely Stiffend Cylindrical Shells". I studied under Dr. Reddy who fist formulated the theory.

If you need it, I am willing to give you a version of my program (source listing). You can try to understand it (with user manual) and use it for your purpose. Let me know if you are interested.

Question:
1) which version of FORTRAN did you use and in which environment (Unix, GNU,Fortran compiler, ... )?
2) is there any user manual or a few guidelines relating to the format of the input file ( if you have one, I would be pleased to receive one of them ) ?
3) did you publish some articles that are treating the plate shell theory you are implementing in this code ? Then I could search them in the library overhere.

Answer
* I used FORTRAN 90 and the program is compilable with FORTRAN PowerStation. I had used dynamic array dimensioning which is supported by this compiler. It is easy to remove the DAD and compile it by any compiler.
* I am glad that you received the source code. I hope it will benefit you.
* I had a paper published at Composite structures: Vol 41, 1998, pp 13-26

* I don't have an on-line version of the user's manual. i will mail it to you later. Regards and good luck...

Question:
I have been working for a while with your program "COSHELL". However I have difficulties with the finite element modelling of two problems. I have read your publication "Local behavior of discretely stiffened composite plates and cylindrical shells" in Composite Structures, but I cannot detect the error in my input files. I hope you could give some advice. Don't be discouraged by the number of text lines. I have written
out all my assumptions to clear out the problem for you. And to further simplify the case for you, I have reduced the problems to a very simple finite element mesh with four four-noded elements. It is a square plate with 200 mm side and 3 mm total thickness (two layers of 1.5 mm and zero angle).

The first problem is the axial tension of an unstiffened orthotropic plate. The attached figure "figures.jpg" shows a schematic drawing of the problem. A uniform tensile stress is applied to one side of the
orthotropic plate. The first input file "axial.txt" that I created, uses the following settings:

1) PX (uniform axial load) = -10 Newton/millimeter. I assume that PX is applied to the plane x = 0, like in the figure in your publication; that's why PX is a negative value. However I am not pretty sure what the impact is of setting the hypothesis to "uniform axial load". It is not clear to me in which plane the uniform axial load is defined.
2) NBSF = 0 (array for natural (force) boundary conditions). I have set this value to zero because I thought that the value of PX imposed already the force boundary conditions.
3) the essential boundary conditions I used are the boundary conditions that are shown in your publication for the buckling of the orthotropic stiffened plate.
The calculation runs without any errors, but the strains and displacements are far from uniform through the height of the orthotropic plate. Moreover they are extremely large. I have three questions:

1) if I consequently use "Newton" and "millimeters", is PX then the value of the applied axial stress in N/mm (along the midplane of the orthotropic plate) ? This would be the total force of 2 kN (see figure), divided by 200 mm. I ask this because the displacements are extremely large for such a loading case.
The elastic properties should then be expressed in N/mm^2.

2) in your publication you mention that the buckling plate is simply supported and as a consequence you impose w=0 on the boundaries. Don't you constrain the Poisson-effect in the thickness direction then, because w=0 is set for all interface nodes and not only for the midplane?

3) I have set the RADIUS value extremely large for the three interfaces to simulate a flat orthotropic plate, but is there perhaps any parameter to indicate that it is a flat plate ?

Because this model was not good, I tried a second one where I applied the external force myself instead of using PX. I also reduced the force the make the displacements smaller. The settings are:

The settings are:

1) PW (uniform pressure load) = 10.0
2) ILOD = 1 (internal pressure)
3) the boundary conditions are the ones that you used in your publication for a stiffened cylindrical panel under a point load.
4) NBSF = 0, because I assumed that the PW value already imposes the force boundary conditions.
Again I hesitate what the units of PW are ? Is this N/mm^2, when you use the units "Newton" and "millimeter" ? Because the w-displacements of the node-number 5 are 600 millimeters!

I would be very grateful if you could take a look at these problems, because I did a lot of work trying to solve the difficulties, but it did not succeed.

Answer

If you put the integration order parameter, NRED (4th parameter at data line 3) to either

0 - full integration
1 - selective reduced for transverse shear terms
2 - selective reduced for transverse normal terms
3 - selected reduced for transverse shear and normal stresses

then you will get reasonable displacements and stresses. I had analyzed it myself this morning. I will attach the result files in the next e-mail.

i.e.,
change
0.0 10.0E0 1.0E00 1 4 1 0 1 1.0 0
to
0.0 10.0E0 1.0E00 1 0 1 0 1 1.0 0

The complete reduced integration could give erroneous results for cases like this depending on the loading, geometry and lamination scheme.

However, using full integration could also give errors on a very thin plate like yours (t/b = 100) when subsjected to transverse loads.

In general, you need to keep in mind that plate and shell elements are sensitive to t/b ratio, loading, lamination schemes and boundary conditions. "Shear" and "membrane" locking are major problems and you need to keep your eyes open on that. This is common in most plate/shell elements; but the inclusion of the transverse normal terms in layerwise theory introduces one more complicating factor.

Hope this helps..... If u still have questions and the numbers still look funny, let me know. I will take care of it.
PS
some answers to your questions..

1) PX (uniform axial load) = -10 Newton/millimeter. I assume that PX is applied to the plane x = 0, like in the figure in your publication; that's why PX is a negative value. However I am not pretty sure what the impact is of setting the hypothesis to "uniform axial load". It is not clear to me in which plane the uniform axial load is defined.

Px is applied to plane x=0;

2) NBSF = 0 (array for natural (force) boundary conditions). I have set this value to zero because I thought that the value of PX imposed already the force boundary conditions.

good...

3) the essential boundary conditions I used are the boundary conditions that are shown in your publication for the buckling of the orthotropic stiffened plate. The calculation runs without any errors, but the strains and displacements are far from uniform through the height of the orthotropic plate. Moreover they are extremely large.

I have three questions:
1) if I consequently use "Newton" and "millimeters", is PX then the value of the applied axial stress in N/mm (along the midplane of the orthotropic plate) ? This would be the total force of 2 kN (see figure), divided by 200 mm. I ask this because the displacements are extremely large for such a loading case. The elastic properties should then be expressed in N/mm^2.

All interfaces at that node will be constrained. w=0 at node 1 means that w-1, w-2, and w-3 at node 1 are all 0.0. Hence Poisson effect is also eliminated.

3) I have set the RADIUS value extremely large for the three interfaces to simulate a flat orthotropic plate, but is there perhaps any parameter to indicate that it is a flat plate ?

that is perfectly fine...

Because this model was not good, I tried a second one where I applied the external force myself instead of using PX. I also reduced the force the make the displacements smaller. The settings are:

1) NBSF = 9, because the force will be applied in the 3 interface nodes (= 2 layers (N=2) + 1) for each of the three nodes at the tensile side of the plate.
2) the {q_u} force vector is assigned a value for each of the 3 interface nodes. The problem is that I don't know how I have to distribute the force between each of the nine interface nodes. I think this depends on the degree of the interpolation function. However, when the calculation is done, the displacements are not uniform again. Thickness-wise, i use linear interpolation....for 2 layers, teh distribution will be therefore 25%, 50% and 25%..etc, etc...

The second problem is an unstiffened cylindrical panel that is subjected to internal pressure. It is a part of a composite pressure vessel that is shown on the left in the figure. The settings are:
1) PW (uniform pressure load) = 10.0
2) ILOD = 1 (internal pressure)
3) the boundary conditions are the ones that you used in your publication for a stiffened cylindrical panel under a point load.
4) NBSF = 0, because I assumed that the PW value already imposes the force boundary conditions.
Again I hesitate what the units of PW are ? Is this N/mm^2, when you use the units "Newton" and "millimeter" ? Because the w-displacements of the node-number 5 are 600 millimeters !

The integration parameter should take care of that.....

Question:
Thank you very much for your help. I was not aware of the fact that such an integration scheme could make such trouble. I still use version 1.0 of your program COSHELL, and I saw in your output result files, that there is mentioned something about the Qi3, Qi4 and Qi5 terms concerning the type of integration scheme used. What does this "reduced integration" refer to ? To the integration of the number of Gauss points ?

* Yes, "reduced integration" makes huge huge difference in the numerical behavior of models. All plate and shell elements suffer from this numerical behavior, so called-membrane or shear locking. The reason is that interpolation functions used for these elements end up making the shear or membrane components stiff for thin plates and shells. If interested, I can give you references on this. It is an interesting area; but causes a headache for the unsuspecting user.

Question:
I have one more little question for you. I have been reading some publications of Reddy and you and others about transverse shear and normal stresses. It is not always clear to me what the implication is of the assumed displacement fields on the transverse continuity and traction-free boundary conditions.
I have tried to summarize it for myself:

1) if w(x,y,z,t) = w_0(x,y,z,t), then the transverse normal strain is zero. The transverse shear strains are constant in each layer. When the transverse shear stresses are derived from the constitutive equations, the transverse shear stress continuity at layer interfaces and the traction-free boundary conditions are NOT satisfied for the transverse shear stresses. Robbins and Reddy reported that layerwise models that neglect transverse normal strain, satisfy the traction-free boundary conditions for transverse shear stresses at the laminate edge only in the integral sense and not in the local sense, despite the level of refinement through the thickness. A more accurate prediction of the transverse shear stresses can be obtained by using the equilibrium method. However Reddy mentions in a study of the ARALL-laminates "Interlaminar stresses (sigma_xz, sigma_yz, sigma_zz) are easily computed from the equilibrium equations of 3D-elasticity when exact analytical solutions are available. An approximate technique is used here to integrate the equilibrium equations and compute the derivatives sigma_xz,z and sigma_yz,z directly from the finite element computation. The equations (with derivatives) are:

sigma_xz,z = - (sigma_xx,x + sigma_xy,y)
sigma_yz,z = - (sigma_xy,x + sigma_yy,y)
(end of citation)

However I don't understand why the condition of the existence of an exact analytical solution has to be fulfilled. And if you would calculate sigma_zz from the values of sigma_xz and sigma_yz, will the estimate be reasonable ? Finally Chaudhuri and Seide have demonstrated that the traction-free boundary conditions are not satisfied for the equilibrium method in the case of asymmetrically laminated plates.

Answer
======
OK....to start with a good reference is a work done by my friend Don Robbins who was also Dr. Reddy's student. He had published a paper in 92 or 93. You may look for it. If you want, I can make a copy of a report he wrote with Dr. Reddy and mail it to you.

Now, the main thing you need to know when using the equilibrium method to evaluate transverse stresses in FE solutions is that since they are functions of the in-plane stresses themselves, their accuracy depends very much on how accurately the in-plane stresses were calculated to start with. My experience which is supported by the works of others is that at stress-concentration points like supports, the in-plane stresses are not evaluated very accuractely. This translates to bad results in your transverse shear/normal stresses at these points.

If you, on the other hand, have an exact analytical solution for the in-plane stresses, then your transverse stresses can reasonably be calculated anywhere in your element.

Remember stresses from FE analysis are reasonably accurate at Gauss points only. As you move away from the Gauss points, your stresses are interpolated, thereby introducing approximations.

2) if w(x,y,z,t) = w_0(x,y,z,t) + W(x,y,z,t), like in your model, the transverse normal strain is taken into account as well. When the transverse stresses are derived from the constitutive equations, will transverse stress continuity (both for shear stress and normal stress) and traction free boundary conditions automatically be satisfied ?

Answer

Depends on the refinement through the thickness. Remember in 3D Layerwise theory (which is essentially descritized 3D), your elements are interpolated both surface and depth-wise. Now if you use enough elements through the thickness, then the conditions will be staisfied very well. If the descetization in the thickness direction is inadequate, then the continuity conditions will not be met. This is very much expected.

Actually, that is why Don Robbins is saying below....

According to Reddy: "The transverse shear stresses, as computed from the constistutive relations, represent the average transverse shear stress within a particular thickness subdivision. However continuous transverse shear stress variation through the thickness can be computed via the 3-D equilibrium equations after having obtained the in-plane stresses accurately." According to Robbins and Reddy: "Note that the form of the natural boundary conditions guarantees that the transverse shear stresses will satisfy the traction-free edge conditions on a local basis as the number of subdivisions through the thickness is increased."

So, when deriving the transverse shear stresses from the constitutive equations, the traction-free edge condition is not satisfied, because the shear stresses represent the average transverse shear stress within a particular thickness subdivision, but if you increase the number of elements through the thickness, the average transverse shear stress in the boundary layers will tend to zero. And what when you use the equilibrium method in this case, are then all conditions satisfied ?

Answer

I will qualify your above statement by saying:

"when deriving the transverse shear stresses from the constitutive equations, the traction-free edge condition CAN BE satisfied" instead of saying "condition is not satisfied"

And I have explained the reason above..i.e., more sub-divisions depth-wise satisfies the reqt.

In your program COSHELL the in-plane stresses and the transverse normal stress are written out. Could you use here the constitutive equations or better the equilibrium method to calculate the transverse shear stresses ?

Answer
Yes, indeed...that is one of the reasons I like the layerwise formulations. It is essentially a 3D solution; hence at Gauss points, you could get accurate in-plane and transverse stresses because these stresses are directly calculated from the constitutive equations. Both Robbins and mine humble work support that. I consider this as one of the powers of the layerwise theory....

Question
The reason why I wrote you about transverse shear stress continuity, is that I started thinking of the way that you apply loads. For example consider a cantilever beam, which is clamped at one side and with three free edges. Suppose that the cantilever beam is loaded with a transverse line loading F [N/m] at its free edge, opposed to the clamped edge. I have attached a drawing "loadcases.jpg" to make it clear. When you want to apply this load in your COSHELL program, there are basically two ways:

(1) you define this line loading as a transverse load, associated with the transverse displacement w_0(x,y,t) of the midplane of the plate.
(2) you define the line loading as a transverse shear stress and you distribute the parabolic shear stress distribution over the various integration points. Then transverse shear stress continuity becomes an important matter.

The same applies to a bending moment M [Nm/m] (see also attached drawing):

(1) you can define the bending moment as two opposite forces, applied to the upper and lower surface of the plate.
(2) you can define the bending moment as a triangular normal stress distribution, which is contributed to the various integration points.

In fact, I am not sure which way to choose in both cases. I have the feeling that the virtual work, done by both definitions, will not be the same, since the interpolation functions of the finite elements play their role. Due to the separate surface- and depth-wise integration, the weight factors of the several integration points through the thickness do not seem to be easily defined in case of the parabolic shear stress distribution or the triangular normal stress distribution.

What do you think ?