*CYCLIC SYMMETRY MODEL yields unphysical result for rotating disk segment

[OliverKunc]

I am trying to simulate a rotating disk consisting of C3D20(R) elements. I have tried both with and without R, the results were not affected. The disk is rotated around the x-axis which is outside of the solid elements.

In a nutshell, here are the control cards of the input deck:

*TIE, NAME = TIE_CYCLIC, CYCLIC SYMMETRY
S_CYC_SLAVE,S_CYC_MASTER
*CYCLIC SYMMETRY MODEL , N = 50, TIE = TIE_CYCLIC
0.,0.,0.,1.,0.,0.
*TRANSFORM,TYPE=C,NSET=N_SUPPRESS
0.,0.,0.,1.,0.,0.
*BOUNDARY
N_SUPPRESS,2,3
*STEP, NLGEOM = NO
*STATIC
*DLOAD
Disk,CENTRIF,1000000.,0.,0.,0.,1.,0.,0.

You can download my complete model on filetransfer.io

I have cyclic symmetry conditions on the cut faces of the disk so it radially carries loads itself. Axial and circumferential rigid body modes are prevented by prescribing Dirichlet boundary conditions to the respective DOF of the nodal set N_SUPPRESS. I have tried varying N_SUPPRESS, e.g. a single node or multiple adjacent nodes, but the result didn’t change.

The result is unphysical because the smallest displacement is not on the lower end and also the displacements vary at the top of the holes.

I have falsified this agains commercial FE solvers. However, I am not able to fix this input deck. Does anyone know how to fix this?

I am using calculix 2.19 and have tried both spooles and pardiso solver on multiple and on a single core, which didn’t seem to affect anything (except runtime).

Thank you very much!

[Calc_em]

@OliverKunc

I have falsified this agains commercial FE solvers.

Can you say more about that ? Which software and with what conditions (also using cylic symmetry or not) ? I submitted this input file in Abaqus and the results are pretty much the same so it’s not a fault of CalculiX but rather the assumptions and settings of this analysis.

Try fully fixing the entire bottom face instead of the current BC - see how this affects the results.

[OliverKunc]

I used Permas with MPC Join conditions w.r.t. a cylindrical coordinate system. As far as I understand, these should be equivalent to the CYCLIC SYMMETRY setup in CalculiX. Here is the result:

This is a very reasonable result. The disk is very stiff, only the upper part gets any significant displacement at all. And displacement magnitude never decreases with increasing radius.

I will try the CalculiX model with fixed lower side later, but still the current definition of N_SUPPRESS should lead to expected results. Especially considering how rotors in compressors and turbines are fixed to the shaft.

[Calc_em]

I would try running the simulation with a simplified disk (no cutouts at the outer edge) and comparing the results obtained in Permas and CalculiX with analytical solution for a rotating disk.

[OliverKunc]

Thanks for the suggestion, will try it out.

But I suspect that I got something elementary wrong with CalculiX because I am new to it. Would much appreciate if @dhondt could take a look at this.

[Disla]

@OliverKunc

The result is unphysical because the smallest displacement is not on the lower end and also the displacements vary at the top of the holes.

Hi Oliver,

Keep in mind the effect of Poisson’s ratio. See difference in pictures. Set it to 0 and check displacement. Should be closer to what you intuitively expect.
If you have set up a high Poisson ratio the volume must preserve and the body around the holes is probably deforming (closing). Deformed or not the important point is that your result must be symmetric which now doesn’t look like.
I would review the mesh. Maybe it is slightly offset .

EDIT: ¿Is your N_SUPPRESS a Symmetric collection respect theta and Z coordinate.?
I would constrain the midplane in X. If not the holes deform unequally.

[Disla]

Hi have applied the cyclic symmetry condition and no other BC is required to avoid Rigid body motion. By other hand ¿Why do you apply a TRANSFORM Card ?

*MATERIAL,NAME=Material
*ELASTIC,TYPE=ISOTROPIC
210000000000,0.3
*DENSITY
7850
*SOLID SECTION,ELSET=Eall,MATERIAL=Material
*TIE,NAME=Cyclic_1,POSITION TOLERANCE=0,CYCLIC SYMMETRY
Slave,Master
*CYCLIC SYMMETRY MODEL,N=50,NGRAPH=50,TIE=Cyclic_1
0,0,0,1,0,0
*STEP,NLGEOM=YES,INC=100
*STATIC
1,1,0,0
*DLOAD
Eall,CENTRIF,1000000,0,0,0,1,0,0
*NODE FILE,GLOBAL=YES
U
*EL FILE,GLOBAL=YES
S
*END STEP

[OliverKunc]

Hi Disla,

thank you very much for your responses and your effort.

@Disla

Keep in mind the effect of Poisson’s ratio.

I used Poisson’s ratio of 0.3 which is the value for steel. This is a standard value, anything much different is unrealistic for applications in this field. For completeness sake, I use mm, t, s as units.

@Disla

Deformed or not the important point is that your result must be symmetric which now doesn’t look like.
I would review the mesh. Maybe it is slightly offset .

The results appear to be symmetric by eyesight and by manual checking of some values. Look at the top corner nodes of the segment interfaces:

Here are the corresponding lines of the .frd file’s displacement block:

-1 3602 7.49161E-03-4.97898E-03 7.91584E-02
-1 3740 7.49161E-03 4.98108E-03 7.91583E-02

-1 4795 6.65408E-04-4.98029E-03 7.91762E-02
-1 4781 6.65408E-04 4.98201E-03 7.91761E-02

This is what I expect, although the accuracy is lower than expected with the y-displacement.

Also, CalculiX accepted the mesh as periodic. It would complain if the periodicity was violated too much.

@Disla

EDIT: ¿Is your N_SUPPRESS a Symmetric collection respect theta and Z coordinate.?

At the moment it is a single point, but i have also tried a symmetric collection of points. Both works fine with PERMAS. Constraining the circumferential and the axial DOF of a single point must suffice because all other possible rigid body modes are constrained by the cyclic symmetry condition.

@Disla

I would constrain the midplane in X. If not the holes deform unequally.

This is not possible for realistic models.

@Disla

have applied the cyclic symmetry condition and no other BC is required to avoid Rigid body motion.

When I do that, displacements are on the order of 10 cm! They must be below 1 mm. This is a 20 mm thick solid steel disk.

@Disla

¿Why do you apply a TRANSFORM Card ?

Because one needs to constrain axial and rotational DOF which are not constrained by the cyclic symmetry model. Assuming CalculiX has no advanced automatic suppression of rigid body modes.

[Disla]

Hi Oliver,

I have been able to load you full model.

I’m receiving a warning.
N_SUPPRESS node 4983 belongs to the Slave Surface.

*ERROR in cascade: zero coefficient on the dependent side of an equation dependent node: 4983 direction: 2

I have change your N_SUPRESS to Node 3568 .( I have also adjusted the nodes on the inner faces including middle ones to a perfect cylinder radius 50mm)

BC in Node 3568 is Y=0 to constrain in plane Rigid body rotation and X=0 to constrain Rigid body displacement in the X direction. NO TRANSFORM is required.
The picture shows X displacement scaled 2000. Y and Z displacements representation are suppressed. Result is symmetric as it should.

*SURFACE,NAME=S_CYC_MASTER,TYPE=NODE
S_CYC_MASTER
*SURFACE,NAME=S_CYC_SLAVE,TYPE=NODE
S_CYC_SLAVE
*MATERIAL,NAME=Material
*MATERIAL,NAME=STAHL
*ELASTIC,TYPE=ISOTROPIC
210000000000,0.3
*DENSITY
8000
*SOLID SECTION,ELSET=DISK,MATERIAL=STAHL
*BOUNDARY
3568,1,0
3568,2,0
*TIE,NAME=Cyclic_1,POSITION TOLERANCE=0,CYCLIC SYMMETRY
S_CYC_SLAVE,S_CYC_MASTER
*CYCLIC SYMMETRY MODEL,N=50,NGRAPH=1,TIE=Cyclic_1
0,0,0,1,0,0
*STEP
*STATIC
*DLOAD
DISK,CENTRIF,1096622,0,0,0,1,0,0
*NODE FILE,GLOBAL=YES
U
*EL FILE
S,E
*END STEP

The disc is stretching in the radial direction due to Centrifugal force and in tangential direction due to cyclic symmetry.
As Poisson Ratio is 0.3, the elemental volumes should be preserved up to a certain degree. The only way is shrinking in the X axis . It is more noticeable close the origin as strain is larger there. Displacement magnitude is computed with all three components.

@OliverKunc

The result is unphysical because the smallest displacement is not on the lower end and also the displacements vary at the top of the holes.

Not sure about this quote. I would agree for Poisson ratio =0 and no hole but here…mmmh. I have some doubts.
¿Could you increase the scale factor of displacemnts on PERMAS to see if the disc is centered?. Let’s say, 2000 factor.

Your result is not symmetric. Check the waves in the picture.

[OliverKunc]

@Disla

I have been able to load you full model.

Dear Disla, thank you very much for your effort! I really appreciate it.

@Disla

I’m receiving a warning.
N_SUPPRESS node 4983 belongs to the Slave Surface.

*ERROR in cascade: zero coefficient on the dependent side of an equation dependent node: 4983 direction: 2

N_SUPPRESS consists only of node 4983 which also belongs to S_CYC_MASTER and not to S_CYC_SLAVE. I don’t receive any warning. What version are you using?

@Disla

I have change your N_SUPRESS to Node 3568

All rotor disks that I have seen have a flange on one of the sides of the disk, so in the current example on either face x=+/-10mm. But for this simplified example we can assume it to be bolted directly to the shaft.

@Disla

BC in Node 3568 is Y=0 to constrain in plane Rigid body rotation and X=0 to constrain Rigid body displacement in the X direction.

Could you please help me out, I don’t understand CalculiX syntax at this point: why do you write

*BOUNDARY
3568,1,0
3568,2,0

as opposed to

*BOUNDARY
3568,1
3568,2

? According to the manual, the first number in the data lines is the node number, the second number is the first DOF to be constrained and the third number is:

Last degree of freedom constrained. This field may be left blank if only
one degree of freedom is constrained

So why do you set it to 0?

To add to my confusion, it also seems to make a difference whether you set the third number to 0 or don’t define it. When I set it to 0, I get this result:

When I don’t define it, I get

Please note that I have changed the material properties and the rotational speed back to my original values in order to keep the result comparable.

@Disla

¿Could you increase the scale factor of displacemnts on PERMAS to see if the disc is centered?. Let’s say, 2000 factor.

Here is a scale factor 5000, undeformed mesh showed for comparison:

@Disla

Your result is not symmetric. Check the waves in the picture.

I thought you meant the cyclicity. The displacements are periodic at the slave and master faces. Yes, the waves are unexpected.

Thanks again.

[Calc_em]

@OliverKunc

Could you please help me out, I don’t understand CalculiX syntax at this point: why do you write

*BOUNDARY
3568,1,0
3568,2,0

as opposed to

*BOUNDARY
3568,1
3568,2

BCs defined before the *Step card are homogenous conditions (no prescribed displacement values, just fixing selected DOFs) and their data lines should contain only the node number or nset name as well as first and last constrained DOF with the last one being optional. So that zero specified there should be ignored by the solver.

[Disla]

Hi Oliver,

4983 belongs to MASTER you are right. Doesn’t change the error. Seems previous Transform on that node together with belonging to a Cyclic condition could be interfering.

I,m using ccx v2.19 MKL Pardiso with MECWAY and the errors are resumed in a separate window. I think if you run ccx from the command line you should review the analysis details in the cmd window.

As there is no previous Transform Card , constraining 1 and 2 corresponds to displacements in X and Y Global coordinates.

As Calc_em says this is probably the most correct form.

*BOUNDARY

3568,1,2

Mecway writes it explicitly for each DOF adding a second coma to warranty both DOF’s are constrained. It has been eaten in one of the copies paste operations sorry. I’m not English speaker and I translate and review grammar before posting. This is the correct optional way of constraining the Homogeneous condition.

*BOUNDARY

3568,1(Two comas here)0

3568,2(Two comas here)0

I can’t explain Permas opening holes and not Shrinking in the X axis. Those are probably other BC’s.

EDIT: The website doesn’t show the two comas inbetween the 1 and 0. . You have to imagine them.

[Disla]

@OliverKunc

we can assume it to be bolted directly to the shaft

I would say that constraining one point do not reflect this condition.
With only One point you get this kind of flower shape at the disc.

¿Could you post a simple sketch to see how do you fix that disc.?

[OliverKunc]

@Disla

I can’t explain Permas opening holes and not Shrinking in the X axis. Those are probably other BC’s.

That was the problem. I rotated the disk around the z-axis in PERMAS. Changing that to x-axis, the results are very similar. Sometimes there are just too many lines of code to see the obvious… sorry for that.

Thank you very much for your support, Disla and Calc_em!