Axisymmetric analysis - Toroidal shell under internal pressure

Hello everyone,

I’m Luca and I’m very new to CalculiX. I have to create a FE model to simulate a spiral case (hydraulic turbine component) using the latter solver. As a first step, I would like to solve a torus under internal pressure via a 2D axisymmetric analysis and then compare the results with analytical ones. The toroidal shell made of steel is subjected to an internal pressure of 10 MPa and it has the following dimensions:

I have created the following .fbd file that generates the geometry and the mesh. The type of elements is CAX8. The displacements along y axis are blocked for the nodes a,b,c and d (cf figure above).

! Geometric tolerance for CGX
GTOL 0.2 

! Center of the torus section
PNT PC 500.0 0.0 0.0

! Points of the torus section aligned with X axis
PNT P01 400.0 0.0 0.0
PNT P02 500.0 100.0 0.0
PNT P03 600.0 0.0 0.0
PNT P04 500.0 -100.0 0.0

PNT P05 405.0 0.0 0.0
PNT P06 500.0 95.0 0.0
PNT P07 595.0 0.0 0.0
PNT P08 500.0 -95.0 0.0

! Lines of the torus section
LINE L01 P01 P02 PC 40.0
LINE L02 P02 P06 2.0
LINE L03 P06 P05 PC 40.0
LINE L04 P05 P01 2.0

LINE L05 P02 P03 PC 40.0
LINE L06 P03 P07 2.0
LINE L07 P07 P06 PC 40.0
LINE L08 P06 P02 2.0

LINE L09 P03 P04 PC 40.0
LINE L10 P04 P08 2.0
LINE L11 P08 P07 PC 40.0
LINE L12 P07 P03 2.0

LINE L13 P04 P01 PC 40.0
LINE L14 P01 P05 2.0
LINE L15 P05 P08 PC 40.0
LINE L16 P08 P04 2.0

! Surfaces of the torus section
GSUR A01 + BLEND + L01 + L02 + L03 + L04
GSUR A02 + BLEND + L05 + L06 + L07 + L08
GSUR A03 + BLEND + L09 + L10 + L11 + L12
GSUR A04 + BLEND + L13 + L14 + L15 + L16

SETA SECTION_SURF s A01 A02 A03 A04
COMP SECTION_SURF up

! Merg points
MERG p SECTION_SURF
MERG l SECTION_SURF
MERG s SECTION_SURF
REP SECTION_SURF

! Points for boundary conditions
SETA PTS_BC p P01 P03 P05 P07

! Line for pressure
SETA L_PRESSURE l L03 L07 L11 L15
COMP L_PRESSURE do


! Element definition
SETA TORUS_SECTION all
ELTY TORUS_SECTION qu8c

! Mesh the model
mesh TORUS_SECTION

MERG n SECTION_SURF

send TORUS_SECTION abq


! Apply boundary conditions
send PTS_BC abq spc 2
send PTS_BC abq nam

! Apply pressure
send L_PRESSURE abq pres 10

Everything seems OK to me when I run the command cgx -b … .fbd.

I have also created the following input file for ccx (.inp) where I define the material properties, BCs and load.

** Nodes and elements
*INCLUDE,INPUT=TORUS_SECTION.msh

** Material parameters
*MATERIAL,NAME=EL 
*ELASTIC 
 200000.0, 0.3
*DENSITY 
 0.000000007850
 
*SOLID SECTION,ELSET=ETORUS_SECTION,MATERIAL=EL
 
** Boundary conditions
*INCLUDE,INPUT=PTS_BC_2.bou

** Analysis
*STEP
*STATIC 
*DLOAD
** Apply pressure
*INCLUDE,INPUT=L_PRESSURE.dlo 

** Results
*NODE FILE,OUTPUT=2D
U
*EL FILE 
S 
*END STEP

when I run the command ccx … .inp , it fails giving me the following error message (for all elements):

*ERROR in e_c3d: nonpositive jacobian determinant in element 1

Could anyone tell me why I’m getting this error and how to solve it to get a solution please? As a beginner I’d appreciate any comments/suggestions/advices regarding my model.

Thank you,

Luca

The modeling approach itself is correct (I recreated this model in Abaqus and it works fine). The problem lies in the definition of elements. They might be somehow inverted - Abaqus gives “zero, small or negative volume” error for this mesh.

you have to switch your surface

add these line before meshing:

flip all

and your script runs like pincese charming

! Element definition
SETA TORUS_SECTION all
ELTY TORUS_SECTION qu8c

flip all

! Mesh the model
mesh TORUS_SECTION

sometimes it is helpful to use symmetric boundaries,
like to calc. only a part or a quarter:

if you work with surfaces,
you have to take care how you create a surface.
a surface has two sides.
you will see these after meshing,
change one surface with - BLEND
see also gsur and of course flip command

positive orientation indicated by the ”+” sign after the surface name

maybe you work with a part of your subject?
maybe you can change to cylindrical system with the use of boundary and results
interesting could also be send:

Another useful method are so called ”cyclic symmetry” equations.

if you like i can share one example about a rotation disk and work with cyl system

wbr

Hello,

Thanks a lot for your help and all the information! With your help, I got results matching with analytical solution.

I will keep in mind to be careful with the surfaces creation and take a closer look to the documentation, but it is already clearer with your explanation.

Indeed I could have used symmetry BC with respect to x axis to model only half of the geometry. I’m not so sure about modeling only a quarter as it is a torus.

Regarding the cyclic symmetry BC, how does it apply to a 2D axisymmetric analysis? I would apply them if I model a sector of the torus in 3D, but I don’t see how in 2D and with CAX elements.

Best regards,

Luca

1 Like

That’s right, to use cyclic symmetry you would need a 3D model of a sector of a structure. That’s just another way to model such symmetric objects.

Did you change something else in your script apart from adding the flip all line suggested by dichtstoff ? Your stress results are smooth while I obtained some spurious concentrations.

Hi,

Yes indeed, I forgot to mention it but I corrected line 69 of the .fbd file which was:

MERG n SECTION_SURF

Now it is:

MERG n TORUS_SECTION

Before this change, nodes were duplicated at the junctions between the surfaces I defined to create the geometry, which was leading to some stress concentrations at these locations.

Here is the full code for the final .fbd:

! Geometric tolerance for CGX
GTOL 0.2 

! Center of the torus section
PNT PC 500.0 0.0 0.0

! Points of the torus section aligned with X axis
PNT P01 400.0 0.0 0.0
PNT P02 500.0 100.0 0.0
PNT P03 600.0 0.0 0.0
PNT P04 500.0 -100.0 0.0

PNT P05 405.0 0.0 0.0
PNT P06 500.0 95.0 0.0
PNT P07 595.0 0.0 0.0
PNT P08 500.0 -95.0 0.0

! Lines of the torus section
LINE L01 P01 P02 PC 40.0
LINE L02 P02 P06 2.0
LINE L03 P06 P05 PC 40.0
LINE L04 P05 P01 2.0

LINE L05 P02 P03 PC 40.0
LINE L06 P03 P07 2.0
LINE L07 P07 P06 PC 40.0
LINE L08 P06 P02 2.0

LINE L09 P03 P04 PC 40.0
LINE L10 P04 P08 2.0
LINE L11 P08 P07 PC 40.0
LINE L12 P07 P03 2.0

LINE L13 P04 P01 PC 40.0
LINE L14 P01 P05 2.0
LINE L15 P05 P08 PC 40.0
LINE L16 P08 P04 2.0

! Surfaces of the torus section
GSUR A01 + BLEND + L01 + L02 + L03 + L04
GSUR A02 + BLEND + L05 + L06 + L07 + L08
GSUR A03 + BLEND + L09 + L10 + L11 + L12
GSUR A04 + BLEND + L13 + L14 + L15 + L16

SETA SECTION_SURF s A01 A02 A03 A04
COMP SECTION_SURF up

! Merg points
MERG p SECTION_SURF
MERG l SECTION_SURF
MERG s SECTION_SURF
REP SECTION_SURF

! Points for boundary conditions
SETA PTS_BC p P01 P03 P05 P07

! Line for pressure
SETA L_PRESSURE l L03 L07 L11 L15
COMP L_PRESSURE do

flip all

! Element definition
SETA TORUS_SECTION all
ELTY TORUS_SECTION qu8c

! Mesh the model
mesh TORUS_SECTION

MERG n TORUS_SECTION

send TORUS_SECTION abq


! Apply boundary conditions
send PTS_BC abq spc 2
send PTS_BC abq nam

! Apply pressure
send L_PRESSURE abq pres 10

Regards,

Luca

1 Like