Inconsistent result with CalculiX on a shell model of a plate in bending

[tcoff27]

Hello! I am using calculiX with prepomax on a calculation of a thin plate (thickness 0.7411mm) in bending, subjected to gravity, supported on 2 edges and with a plane of symmetry on the other two edges. At 2 nodes at one extremity, I block UX movements to avoid rigid body movements.

Note: I directly imported my mesh from external software via the .inp format. These are linear interpolation elements.

In my 1st case, the case of a flat plate, I am OK with the result (I compared it analytically). See attached “flat_plate.pmx” I get a U2 max displacement of 396.8mm. I noticed that the default solver (pastix) gives completely wrong results (the calculation diverges it seems to me). But, with other solvers, it’s OK, for example “Pardiso”. flat_plate.zip

In my 2nd case, all conditions are identical, except the shape of the plate. This plate has a wave shape. The problem is that I have a result that seems completely wrong to me whatever the solver. The maximum deflection U2 obtained is 0.06647mm (with pardiso). I use another software at my disposal, I obtain a result close to 386 mm. Which seems much more realistic to me. See the attached “wave_plate.pmx” file. On prepomax forum, someone has exported the input file from wave_plate.pmx and submitted it in Abaqus. The result is 384.4 mm. So the setup seems correct and it seems to be an issue with CalculiX. https://we.tl/t-ZWAWhirzYv

For information, my ultimate aim is to activate nlgeom to consider large transformations. Thank you in advance.

PS: here are some screenshots of my 1st case :

and of my second case :

[victorkemp]

@tcoff27 I wonder if you could reduce it to a simpler example and an .inp file? There are some known bugs with rotation constraints on shells but it's not clear if they're what's happening here.

Just because a .inp file works in Abaqus doesn't mean CCX is wrong. They may have different requirements for mesh refinement and interpret some cards differently (eg. *ORIENTATION).

[victorkemp]

I can't quite see the BCs from the screenshots but do you have any rotation BCs that don't meet this requirement (from CCX manual):

"... no drilling rotation should be prescribed, unless all rotational degrees of freedom are set to zero in the node. If the shell surface is not aligned along the global coordinate directions, prescribing a moment or rotation aboun an axis perpendicular to the drilling direction may require the definition of a local coordinate system"

[tcoff27]

red_wave_plate.zip Hi @victorkemp ! Thanks for your message and your help! (and sorry for my late answer). Following your suggestion, I reduced my example. You can find it attached to the .inp and .pmx formats. I'm not sure to well understand the requirement about rotations constraints. What is a drilling rotation for a shell element?

I explain deeper my calculation:

• The waved plate has an angle to the Z axis and the modelled plate is a cyclic pattern among Z. So the constraints are UZ=RX=RY=0 on edges 1 and 3.
• The waved plate is simply supported on edges 2 and 4. So the constraints are UY=0.
• To remove rigid body motion, constraints on nodes 5 and 6 are UX=UY=UZ=RX=RY=0.

About the results, I think there is an error in CalculiX. Whatever the solver used, I have a max displacement of -0.001036 mm at the center and not enough displacement at the edges 1 and 3. (Does the issue come from the RX and RY constraints?)

In another software, here is the result: (max displacement of -0.0335mm from edge 1 to edge 3)

[victorkemp]

Thanks for the clarification.

The drilling DOF is the rotation DOF normal to the shell's mid-surface. It looks like all your rotation constraints have some component in this direction so they are definitely at risk of making it go wrong. You should put a *TRANSFORM on each node along the edge and only have a single rotation constraint on each node about the edge tangent, as well as the existing displacement constraints.

I'm not sure if there's another problem too but the displacement of -0.001036 mm is more than I get with fully constrained symmetry edges (-0.00055 mm), so it's consistent with this being the only problem.

[xyont]

hi,

is it possible to provide CAD files, e.g Step or Iges of reduced examples? Interesting to know and study further about the problem.

from my first look in the model, it seems there's over-constrained at some edge. Try to modify as below, Displacement_Rotation-1 : U3, UR1 (remove UR2) Displacement_Rotation-2 : U3, UR1 (remove UR2)

Also, it's required to improved by quadratic element S8R with layered and eliminate using linear element S4, when bending is dominant in the case.

result by removing rotational dof at corner nodes.

[tcoff27]

DONN10.zip Thank you very much for your messages (sorry again for my response time)! @xyont I attached the geometry file (DONN10.igs).

By removing the UR2 at edges 1 and 3 and at nodes 5 and 6, I indeed obtained good results compared to -0.0335mm with my professional software.

Boundary conditions updated:

Results with lineic mesh (S4 elements):

Results with quadratic mesh (S8 elements):

Results with quadratic mesh (S8R elements):

Generally speaking, S8R elements are recommended by default rather than S8?

I'm just coming back to the subject of drilling DOF @victorkemp
I am disturbed by this point. I often work with shell elements. I am used to model an X-Y symmetry plane by constraints UZ=RX=RY=0.

If we consider the case of a flat plate in bending, RY is the drilling rotation, we would therefore have to remove this constraint RY if I understand correctly?

Another example (taken from a Prepomax tutorial on YouTube):

This model has 3 symmetries planes : edge 1 = symmetry Y-Z, edge 2 = symmetry XY and edge 3 = symmetry X-Z. Given the curvature of the shells on this vessel model, we are obliged to block the drilling rotation at a given moment. For example, at the bottom of edge 1, the drilling rotation is RZ. At the top of edge 1, the drilling rotation is RY. According to the manual, I understood that blocking drilling rotation was prohibited, but then how to deal with this case? Yet the result seems correct with the current DOF.

Maybe I didn't well understand the manual. "... no drilling rotation should be prescribed, unless all rotational degrees of freedom are set to zero in the node."

I am therefore very suspicious of the use of boundary conditions when using shells elements. What do you think ?

Once again thank you very much for your help!

[xyont]

in my understanding, restrained nodes of shell element applied at local axes. In case no specific definition by user orientation, local axis being used is global. However, depending on element by element normal angle transition the definitions can differ at curved edges. Drilling dof's probably not a problem at straight or flat edge, but can affect at curved condition.

Regarding an example from YouTube, it seems to have over-constrained also. Edge-1 is required U1 and UR2 or UR3. Edge-2 is required U3 and UR1 or UR2. Edge-3 is required U2 and UR1 also UR3. Needed me to investigate further, i'm also frequently get unexpected result when user local axes is not explicitly defined. Fortunately, there's classical shell element available and extruded model result is possible to export in true solid element.

*confrimed by solid element model (C3D20R) the model of shell element from YouTube example is over-contained and lead to high stress spot at areas of edge transition from vertical to curved.

[victorkemp]

I am used to model an X-Y symmetry plane by constraints UZ=RX=RY=0.

This is the usual way in FEA but you don't actually need to constraint the drilling DOF because the displacement constraints on the adjacent nodes are what really provides that. In most software, that DOF effectively gets ignored so it's harmless to include it and saves you having to decide the normal direction. In CCX it usually gets ignored but sometimes causes wrong constraints.

That screenshot from the Prepomax tutorial looks like it's at risk of this problem too, and should be done with *TRANSFORM and only tangential UR constraint on each node.

Correct, don't constraint RY on a flat plate in the XZ plane.

[tcoff27]

Thanks a lot for your answers! I think I understand how it works better. This is well noted for the prepomax tutorial.

Consequently, on my waved plate calculation, I also removed the UX on corner nodes 5 et 6. I had the same results. Then I apply the TRANSFORM keyword: TRANSFORM,NSET=Edge-1,TYPE=R 1.,0.,0.,1.,1.,0. TRANSFORM,NSET=Edge-2,TYPE=R 1.,0.,0.,1.,1.,0. TRANSFORM,NSET=Edge-3,TYPE=R 1.,0.,0.,1.,1.,0. *TRANSFORM,NSET=Edge-4,TYPE=R 1.,0.,0.,1.,1.,0.

Because my aim is to have a symmetry plane X-Y (global coordinate) on edges 1 and 3, the defined local coordinate is the same as global, but based on your advice, this allows to be sure to have the right coordinate on the edge. On edges 3 and 4, I just constraint UY=0, I keep the same coordinate. Does it look OK for you?

On my reduce model, result is good compared to the result of my professional software (-0.0335mm) : -0.033 mm.

About the *TRANSFORM command @victorkemp, to get only tangential UR on a curved edge (as the example of the vessel on youtube), it seems to me not easy. Are we obliged to select one new coordinate for each node of the edge? What's more it would depend on the meshing resolution…

To finish, I consider my initial calculation, with the large model (refer to my 1st post). I used the same improvements on constraints and TRANSFORM keyword. In my initial message, I had a displacement of 0.06647mm. Now I have 111.1 mm. Unfortunately, I expect a result close to 386 mm (based on my professional software and the result of a flat plate). Do you have any idea on what is wrong with my calculation setup?

Just in case, I share with you the .igs (surface) and .inp (mesh) files of my model. waved_shells.zip

Expected result:

[xyont]

i'm not in detail, the discrepancy is too large, try to exclude *Transform first. Also, questioning related to large deformation (Nlgeom) is turned it off or not?

[tcoff27]

i'm not in detail, the discrepancy is too large, try to exclude *Transform first. Also, questioning related to large deformation (Nlgeom) is turned it off or not?

I run the calculation after removing *TRANSFORM. Nlgeom is turned off. (I plan to consider large deformation later) Here is the result:

I confirm that boundary conditions are:

The waved plate measures 1400 mm x 150 mm x 0.7411 mm (thickness) Young's modulus = 195300 MPa / Poisson coef = 0.31 Thanks for your help

[victorkemp]

I've applied my previous suggestion of using *TRANSFORM to make the RX constraints aligned with the shells and it's only 7% different from Abaqus (-411 mm vs 384.4 mm). I did this using another preprocessor.

quadratic3.zip

This is with linear elements. Maybe quadratic is better but the midside nodes should ideally follow the curves and they don't in @tcoff27 's mesh.

[tcoff27]

Thanks for this file .inp !

Here is what I did. I import this file in Prepomax, I can clearly see the *TRANSFORM keywords on all the nodes on the edges. I confirm that "Material" and "Section-1" are OK.

The BC are missing after importing the inp file, so I apply the following:

And the result is -151.6 mm, I don't have the -411 mm. Maybe I miss something?

May I ask you what preprocessor did you use to define the *TRANSFORM?

[victorkemp]

I did the BCs differently at the corner nodes (omitting UY and moving UZ in one node) because the transforms affect the direction of UY there but I'd be surprised it makes that much difference. Make sure you get the same -411 solution running the file directly in CCX.

I used Mecway but modified it to do this. In particular, removing the Z components from the *TRANSFORMS and URY constraints.

[tcoff27]

Hi! Thanks for your message. I try to run again my calculation (in prepomax) by removing some constraints at the corner nodes, but I get a very close result at around 151mm...

So as suggested by @victorkemp, I wanted to try to launch your file "quadratic3.inp" directly in CCX but I only had Prepomax (I am a beginner with calculix ^^). So I tried to find a calculix executable for windows.

I first install calculixforwin, then I select the quadratic3.inp file and click on "Run Solver - CCX" but unfortunately I get the message "*ERROR in u_calloc: error allocating memory variable=nodempc, file=ccx_2.13.c, line=261, num=376864419, size=4"

So I tried to install bconverged, but I'm having trouble loading a file or just updating the solver version. With the old solver version 2.10, I tried to run the calculation in command lines with ccx quadratic3

************************************************************

CalculiX Version 2.10, Copyright(C) 1998-2015 Guido Dhondt
CalculiX comes with ABSOLUTELY NO WARRANTY. This is free
software, and you are welcome to redistribute it under
certain conditions, see gpl.htm

************************************************************

You are using an executable made on Mon May 23 13:24:06 2016

  The numbers below are estimated upper bounds

  number of:

   nodes:     10635248
   elements:       163555
   one-dimensional elements:            0
   two-dimensional elements:       163555
   integration points per element:            8
   degrees of freedom per node:            3
   layers per element:            1

   distributed facial loads:            0
   distributed volumetric loads:            3
   concentrated loads:            0
   single point constraints:     11787192
   multiple point constraints:     19629409
   terms in all multiple point constraints:    125621473
   tie constraints:            0
   dependent nodes tied by cyclic constraints:            0
   dependent nodes in pre-tension constraints:            0

   sets:         1251
   terms in all sets:       328387

   materials:            1
   constants per material and temperature:            2
   temperature points per material:            1
   plastic data points per material:            0

   orientations:            0
   amplitudes:            4
   data points in all amplitudes:            4
   print requests:            0
   transformations:         1250
   property cards:            0

 *WARNING reading *NODE FILE or *EL FILE: label not applicable
          or unknown;
 *WARNING reading *NODE FILE/OUTPUT or *EL FILE/OUTPUT or *CONTACT FILE/OUTPUT . Card image:
          S,NOE,E,ENER

 STEP            1

 Static analysis was selected

 Decascading the MPC's

 Determining the structure of the matrix:
 number of equations
 2458005
 number of nonzero lower triangular matrix elements
 67123475

 Using up to 1 cpu(s) for the stress calculation.

 Using up to 1 cpu(s) for the symmetric stiffness/mass contributions.

 Factoring the system of equations using the symmetric spooles solver
 Using up to 1 cpu(s) for spooles.

The calculation seems to have run, but I don't succeed to display the results. I just succed to display the meshing when running cgx quadratic3.inp. Am I on the right track? Really sorry for my beginner questions..

[victorkemp]

That's CCX 2.10 which is very old. It's confusing with all the different information around but it's actually very easy. Download http://dhondt.de/calculix_2.21_4win.zip and run ccx_static.exe directly from the command line. You don't need bconverged, calculixforwin, cgx or anything else just to solve a .inp file. Open the solution's .frd file in Prepomax.

I think there's a problem with my *TRANSFORMS in quadratic3.inp. Maybe they're slightly misaligned. Here's another one that was generated in a cleaner way and the displacement is even closer to the other software (-381 mm) or about 1% error.

quadratic4.zip

If you just need to solve this problem, it's probably best to use solid elements (C3D20). Then all the troubles with constraint directions will disappear.

Since this discussion is moving away from being an "issue" in CCX, I think it might be polite to continue on https://calculix.discourse.group/ . If you post there, I'll see it.

[tcoff27]

Thank you very much for explaining this procedure to me. It’s actually much simpler that way! So I managed to execute the quadratic4.inp calculation directly in ccx. I can well check the result of 381 mm in the frd file on prepomax. Yes, we strayed from the initial topic, sorry. Ok, I will restart my calculation with solid elements, I will post on the forum if needed. Thanks once again for your help!

I summarize the topic: we have to be careful when using curved shells with boundary conditions. In particular, drilling rotation must not be constrained. The inconsistent results I got are due to these overconstraints on the rotations at the edges. This can be remedied by adding *TRANSFORM to the nodes.

The problem is therefore more of the order of calculation setup. So, I don't know if it should be seen as a bug and fixed. I'll let the developers judge.