Orthotropic Shell/Plate Analysis


I need some help.

I’m working on the several examples described in the book “Mechanics of Laminated Composite Plates and Shells (JN. Reddy)”
Basic input data (E1,E2, E3,G12,G13,G23,v12).

E3 is conditioning the convergence of the nonlinear analysis with error:

“*ERROR: increment size smaller than minimum best solution and residuals are in the frd file”

I found that E3 must be bigger than certain value or the nonlinear problem fails.

1-I suspect the issue might be that CalculiX do not check and warns if the Material Stiffness Matrix is positive definite for the orthotropic input data as well as the constrains between the individual Poisson’s ratios.

2-¿What happens with the associated Poisson values v31 and v32 when we enter E3? . ¿Are they internally computed?.

3-Additionally , everything has run fine up to Example 10.5.4 where I have seen some strange behaviour with reduced integration S8R elements.
¿Are SR8 Elements allowed with orthotropic materials?. I have found they get distorted in the middle nodes. See picture. (S4R works with limited accuracy but well)

4-Additionally , ¿Is the traction load available for orthotropic shell elements in nonlinear analysis?. It is failing systematically where pressure load converges easily.



I have some more information.

If I refine it fails before in time. Then I need to to increase the E3 value again to make it converge.

I have been able to capture this picture. Seems like the upper surface that is in compression ¿is buckling?

¿Is that possible?.

¿Could this be the response?. If Normal stiffness E3 is too low, middle nodes on quadratic elements misalign and buckle.

2-No they aren’t. Because 9 independent parameters are used for a 3-dimensional orthotropic solid. All of them have to be input. In classical lamination theories the laminae are considered to be in plane stress state, and therefore kinematic relations (Kirchoff, Von Karman) are used for transversal deformations; hence the elastic transversal constants aren’t used (like in Reddy’s example 10.5.4). I am not familiar about CCX implementation but I guess is 3D so the 9 constants are needed. Maybe setting v13 and v23 = 0 and E3=10*max(E1,E2) can reproduce the plane stress behaviour but not sure if the formulation will tolerate this input.

3- the manual reads

Therefore, a different method has been coded to treat composites. Right now, it can only be used for 8-node shells with reduced integration (S8R) and 6-node shell elements (S6)

4-no idea, but I think there are no edge loads available for CCX shells, they have to be input as point loads at edge nodes. Search the forum posts.

Thank you very much Juan.

Seems E3=10*max(E1,E2) eliminates that strange middle node deformation.
It has a lot of sense to me. If v13 = v23 = 0, E3 should be very large to respect the condition:

v13/E1=v31/E3=0 and v23/E2=v32/E3=0.

I’m comparing Laminate composites with the equivalent orthotropic element that results from the ABD Matrix.

Note: I think symmetry of the compliance Matrix and positive sum of works done by all the stresses imposes some relationship between the 9 parameters . They are not all independent. I Guess Ccx should warn the user when those conditions aren’t satisfied. It makes the convergence to fail without warning on this sense.

Hello again,

In the same context as in the previous post, according to the book I’m working with, I have seen that there is not only one possible integration scheme for reduced elements.
That’s new to me.

Note from the book:

“R-R , F-R , F-F . (R-Reduced Integration Scheme, F- Full Integration Scheme)

  • The first letter refers to the integration rule used for the nonlinear terms while the second letter refers to the integration rule used for the shear terms.”

I’m looking at the ccx v2.19 manual but I don’t see information regarding this point.
Figure 58: Only shows a number and location of integration points (2x2x2 integration) and a reference to the shape functions.
¿Should I assume that the integration scheme is R-R for reduced elements in Calculix.?

I’m mainly talking about

C3D20RL: expanded 8-node shell element with reduced integration = S8R
C3D20RLC: expanded composite 8-node shell element with reduced integration = S8R

Thanks in advance

1 Like

Looking at the documentation, you are right, it is not selective reduced integration but uniform reduced integration. My experience is that works well in general but if the problem is regular and I don’t fear any type of misbehaviour I would use S8 and not S8R (more flexible element). For more details you should check the theoretical reference: The Finite Element Method for Three-Dimensional Thermomechanical Applications - Guido Dhondt - Google Libros

1 Like

Thanks again Juan. I really appreciate your support.

Hi ,

As I read from the manual v.2.19 pag 106, ccx expands the composite shell element into a unique 3-D Brick element to solve and then, the stresses are computed using each layer property. This is possible because there are more integration points, and they are assigned the material properties appropriate for the layer they belong to.
¿How is the first step of solving the 3-D Brick element done?. ¿Does Calculix transform the Composite Layered into the equivalent Ortothropic element by means of the ABD matrix?.


Not sure since not documented in the manual, but I think is not necessary since in the corresponding equation for the finite element formulation before for computing the integral just need the values of the material tensor in the required locations to represent the stiffness of the layered material. So the stiffness matrix is built directly using the values at the integration points of the anisotropic material tensor.

Thank you Juan for your comment and the referenced book.

Good afternoon,

I have obtained some unexpected results when trying to validate through thickness shear stresses in a Cantilever Plate.

To obtain some detail in the shear stress through thickness, I have created an homogeneous composite (layered) material with two layers with exactly the same properties and half the thickness each one. This forces Ccx to generate additional through thickness integration points. See ccx manual v.2.19 pag 106.
The properties of each layer being E1=E2 =E3 and G12=G23=G31=E/2(1+Nu)

This should also be a comparison test by itself as the composite should theoretically approximate to the isotropic behavior under these conditions. “Pure” Isotropic is also included in the test.

Some preliminary conclusions:

1-CalculiX shell element provides uniform shear stress through thickness compared to other formulations. (Could be non-conservative but that was expected) .

2-The "Isotropic composite” approach behaves as expected.

3-Shear result for the Static Analysis and ISOTROPIC material give wrong results when using full integration S8. (Composite is always S8R and gives good result)

4-Shear under Nonlinear Quasi-Static seems to have some problem as I’m obtaining wrong results for both Isotropic (no matter S8 or S8R) and Composite.

Appreciated any comment if this could be some bug in Through thickness shell shear stress + Nonlinear.

Attached report. I also posted the free run file in the Mecway forum in case it is a Pre/Post processor error.

most shells use Midlin-Reissner formulation (also known as first order shear deformation theory FSDT) so they can only represent constant transverse shear strains. CCX shells are based on hexaedral elements that if parabolic (2nd order) only can represent linear transverse deformations/strains, so in the end are almost equivalent to FSDT shells because the first order approximation of a parabola is a constant value, therefore what you get is completely consistent. If you are interested in transversal shear stresses you have to resort to a 3D model where you have at least 3 elements in the thru the thickness direction, so a parabola can be approximated by straight line segments. This is not changed by using laminate approach since strains are described by a single element in the thickness direction. Non-linearity here only introduces the stress stiffening effect due to large displacements I think making more difficult for the model to accurately represent the actual distribution. I’d say S8 should behave better in this case but not sure.

Hi JuanP74 and thanks for your response.

Through thickness shear stress and Sxx response agrees perfectly with the expected result for Static as shown on the report .

In fact, I would say that they agree too well if that’s possible as all four models shear gives exactly -300Psi all over the body.
That is the same result as for FSDT as shown in Figure 6.1.

What is surprising to me is that going to Nonlinear Quasi-static I get such wrong results in all four models.

It doesn’t even provide a uniform value of shear as it should be for a cantilever plate.


I have added the Solid option (3 elements in the thru) C3D20R and it shows the same issue for Sxy.
Sxx and Displacement looks ok.
Refinement do not change significantly the result.

SXY can’t be positive and we reach a value of 7.794 Psi

linear cantilever has a constant shear force along the whole length, so stress is constant, however if non linear effects (large displacements) are relevant then you don’t get constant shear because of that, since the load is no longer a dead load but changes direction with the beam deformation. I do not know if this is the case, check displacements. If that were the case there is nothing a priori wrong in your results as I said but for the thru the thickness distribution.

Hi Juan,

I have reduced the load a factor of 10 to avoid Nonlinear geometric effects.

Refined the solid element in thickness x 8

Refined the solid length x 4

Static Analysis shows a beautiful “parabolic” profile through the thickness and Shear Stress vanishes at the free surface as it must.

Nonlinear Quasi-static shows a linear increase from -80 Psi to 78 Psi at each opposite surface.

I’m 99% sure Nonlinear Quasi-static result is wrong.

Hello, I think I realized what’s going on…the stresses in nonlinear analysis are given in global axis as well however the beam is rotated precisely because of nonlinear geometrical effects so what you see are the stresses in global axis that are difficult to interpret. In your picture if you look at the beam root you see the parabolic distribution because here rotations are small. Try loading the beam with very small load so deflections and rotations are so small that the angle is almost zero, you will recover the parabolic stress distribution you were looking for. Anyway is a very good thing you come up with this because most of the time I work with invariants and do not know how to get the stress tensor in a different set of axes! Will be very helpful for the forum if you find out.

Hi ,

I have set up the problem in a different way. The element has some initial slope such as the free opposite edge ends up flat.

Both Static and Nonlinear analysis reproduce in a rough way the profile at the oposite extreme. You are right. Shear plot is not following the element curvature.

I understand by your response that there is no option yet for Ccx to plot out of plane Shear Stresses in element local coordinates ¿isn’t it?.

Static solution was so close that has completely confused me and I was expecting the same or better in Nonlinear.

Thank you once again, Juan.

I will post if I found some workaround to represent Shear in the element coordinates. It seems a hard problem.

There is a way to do that. Tomorrow, I will attach a test file I created awhile ago.

These are the commands to get stress in the element axis:


The GLOBAL=NO is the thing you want. Victor, from Mecway, gets all the credit here. He helped me with all my questions. Once I finally figured out how to do composites with Mecway and CCX, it was really easy. But it took me about a month of trying things before I got to that point.

Below is a link to a bunch of my old test files. They are Mecway files that use the cards above. I don’t really remember a whole lot about the various tests. However, there are some things in there that should help you out.