= Buckling of a tube under twist (non-linear) calculix vs code-aster =

jlucas · 2022-08-15 18:12:37

Hello,

I'm trying to perform a simulation to identify buckling of a tube under a twist moment. So far I did not manage to achieve this with code-aster, mainly, because I could not get right the formalism to apply large rotations. I present bellow my progress starting with the achievement using calculix and then each of my failed tries with code-aster. Any help/comment is welcome.

* Geometry of the tube:
radius (m): 0.030
thickness (m): 0.002
lenght (m): 0.950
* Material properties (carbon steel):
E (N/m^2): 207000000000.0
nu: 0.290
* Boundary conditions:
Fixed in one side;
Twist moment is applied on the other side (0 - 6. rad).
All other degrees of freedom except twist are supposed fixed on the loaded face.

The mesh is made of solid 1st order tetrahedral, so very accurate results are not to be expected.
The main objective at present stage, is only to check qualitatively that I the simulation is giving some sensible results.

jlucas · 2022-08-15 18:22:31

== CalculiX. ==
Enclosed are the simulation files, load path and an animation of the deformation.

(Note: The input files from CalculiX are compatible with Abaqus, so it should be easy to adapt and run the simulation using this software (in case someone wants to try)).

The loading is applied by connecting all dofs of a reference node with a surface using a *COUPLING/KINEMATIC card (see manual for reference).

****
*COUPLING, REF NODE=7430,SURFACE=sLoad,CONSTRAINT NAME=CN1
*KINEMATIC
1,2,3
****

Output files are converted to vtu in order to be viewed with paraview. Load path shows clearly stability loss associated with buckling and the following animation shows the result of the deformation. Qualitatively as expected.

Attachments:
buckling_tube_calculix.gif, Size: 192.08 KiB, Downloads: 199

jlucas · 2022-08-15 18:23:35

== CalculiX Files ==

Attachments:
calculix.zip, Size: 646.61 KiB, Downloads: 114

jlucas · 2022-08-15 18:25:45

== Code-Aster: Try A ==

The loading is applied via a discrete element connected with LIAISON_SOLID to the nodes on the face where twist is to be applied. Deformation is set as PETIT, so not really expected to give valid results.

Attachments:
buckling_tube_code_aster_tryA.gif, Size: 318.02 KiB, Downloads: 138

jlucas · 2022-08-15 18:27:12

== Code-Aster: Try A. Files ==

Attachments:
tryA.zip, Size: 757.22 KiB, Downloads: 89

jlucas · 2022-08-15 18:28:35

== Code-Aster: Try B ==
Mainly the same as Try A, but setting the deformation as "GROT_GDEP".

Looking at the load path and the deformation obviously something is wrong. Help is welcome and appreciated.

The doc U4.44.01 (AFFE_CHAR_MECA) (section 3.6) specifies that to use LIAISON_SOLIDE with large deformation hypothesis, the loading needs to be specified as TYPE_CHARGE='SUIV'. But by specifying such leads to the following error:

║ Erreur utilisateur :
║ Le chargement contient des relations cinématiques LIAISON_SOLIDE qui sont non-linéaires
║ lorsque l'on utilise EXCIT / TYPE_CHARGE='SUIV'.
║ Mais ce cas n'est pas traité car il y a au moins un noeud qui porte les degrés de liberté
║ DRX, DRY et DRZ.

So TYPE_CHARGE='SUIV' is not compatible with a rotation?

Attachments:
buckling_tube_code_aster_tryB.gif, Size: 124.4 KiB, Downloads: 128

jlucas · 2022-08-15 18:30:06

== Code-Aster: Try B. Files ==

Attachments:
tryB.zip, Size: 567.86 KiB, Downloads: 88

jlucas · 2022-08-15 18:32:27

== Code-Aster: Try C ==

Trying a "rotational" boundary condition on all nodes as in example ssnv507 (V6.04.507: Rotation d'une inclusion rigide avec X- FEM)

DTH=6.0;

tournex = FORMULE(VALE='sqrt(X*X+Y*Y)*cos(atan2(Y,X)+DTH*INST)-X',
                  DTH=DTH,
                  NOM_PARA=['X', 'Y', 'INST'],)

tourney = FORMULE(VALE='sqrt(X*X+Y*Y)*sin(atan2(Y,X)+DTH*INST)-Y',
                  DTH=DTH,
                  NOM_PARA=['X', 'Y', 'INST'],)



load=AFFE_CHAR_MECA_F( MODELE=model,
                       DDL_IMPO=(_F(GROUP_MA='sLoad',
                                    DX=tournex,DY=tourney)),
);

Results are still not right. See load path and deformation. Also, I don't know with this formulation how to specify in addtion that all other degrees of freedom except twist are fixed.

Attachments:
buckling_tube_code_aster_tryC.gif, Size: 181.76 KiB, Downloads: 125

jlucas · 2022-08-15 18:33:37

== Code-Aster. Try C. Files. ==

Attachments:
tryC.zip, Size: 548.26 KiB, Downloads: 102

jeanpierreaubry · 2022-08-17 08:35:12

hello

i think you are loosing track with try C

at first lets forget the critical buckling load and concentrate on the reaction forces

if you were using a much coarser mesh, a few tens times less elements !
you could apply nodal forces on the perimeter at the top to mimic the torsion couple
first you could find that the results are similar with a coarser mesh
second you could find (as i did) that the results are the same as with LIAISON_SOLIDE
which then cannot be considered at fault

just simply the calculation of reaction is completely erroneous with GROT_GDEP
and this seems to be the case whatever the element type 1D, 2D, 3D or whatever the load as i found out from other studies

and if the reactions are wrong one can wonder if the stress are not wrong either

i made several post here about that with no reaction from EDF team

are we at fault in our way of putting the problem together?
or is code_aster at fault?
it would be nice to have an answer

jean pierre aubry

jeanpierreaubry · 2022-08-17 13:50:08

attached is a simplified similar problem
same overall geometry
only plate elements with just only 240 elements
torsion couple is imitated with only 4 forces at each quadrant
so do not look too much at result around the load area
the general behavior is very similar to jlucas mesh with 22567 elements
applied torsional couple 20000
results are

with GROT_GDEP
sum of reaction 7000
maximum tangential displacement 0.0088
VMIS near fixed points 1.2e9

with PETIT
sum of reaction 20000
maximum tangential displacement 0.024
VMIS near fixed points 3.43e9

analytical
sum of reaction 20000
maximum tangential displacement 0.0224
VMIS near fixed points 1.89e9

comments
reaction with GROT_GDEP are wrong 0.35 times the applied load
displacement with GROT_GDEP are 0.36 times the one with PETIT
VMIS is not so far out from analytical

what to conclude about that????

Attachments:
plate_model.zip, Size: 7.04 KiB, Downloads: 104

jlucas · 2022-08-19 10:08:55

Thank you Jean Pierre Aubry for looking more deeply into this and for simplifying the problem. It might help to grasp it easily.

I totally agree that the silence on this topic from development team is quite uncomfortable. (I guess most members might be on Holiday?)

I post here some selected links with related topics:

Challenge: Buckling of a tube under twisting moment...
(can not offer any other prize than personal gratification and honor mention in present post!)

Steel tube under twist moment. Large deformations. (code-aster v15.4)

Buckling of a cylinder under a twist moment.

Also might help:
LIAISON_SOLIDE with TYPE_CHARGE='SUIV' and GROT_GDEP

Rotation of a cylinder

DEFORMATION='GROT_GDEP’ and FORC_NODA moment equilibrium

jlucas · 2022-10-19 23:03:08

== Back to the buckling of a steel tube. ==

I must say I am quite disappointed and frustrated for not having any additional reply on this topic. I was hopping to get some justification/hint from someone with the knowledge (including the development team) to why things are not working as expected when using the large deformations formulation, specifically why LIAISON_SOLIDE is not working for rotations or why POST_RELEVE_T operator is not able to compute the reaction moments correctly.

Anyway, by playing with options, I've managed to perform a simulation in code-aster where the such said buckling condition is achieved (see attached Figure). The files for the simulation are also attached (bellow). There is only a small change in the geometry (in relation to tryC) which is meant to simplify the analysis. A new physical group (cLoad) is created to apply the torsion load along the outer line of the driven top surface of the cylinder. (This way the number of loaded nodes decreases substantially). The other more important change (in the command file) was to change the TYPE_CHARGE from 'FIXE_CSTE' to 'DIDI'. i.e:

result=STAT_NON_LINE( ...
                      EXCIT=( _F(CHARGE=fix,),
                              _F(CHARGE=load, TYPE_CHARGE='DIDI'),
                      ),
                      ...
);

I don't really understand what such type of loading is doing. U4.51.03-Opérateur STAT_NON_LINE (pp.15) states:

Si TYPE_CHARGE vaut 'DIDI' alors les conditions de Dirichlet (déplacements imposés,
conditions linéaires) s'appliqueront sur l'incrément de déplacement à partir de l'instant donné
sous ETAT_INIT/NUME_DIDI (par défaut l'instant de reprise du calcul) et non sur le
déplacement total. Par exemple pour un déplacement imposé (mot clé DDL_IMPO de
l'opérateur AFFE_CHAR_MECA) la condition sera de la forme u−u 0 =d où u 0 est le
déplacement défini par NUME_DIDI et non u=d .

The automatic translation into English did not really help!!!.

I would expect that the right option to apply is TYPE_CHARGE='SUIV', as such option (as I understood)
applies the loading on the actualized geometry: But doing so, I get the following error:

╔═══════════════════════════════════════════════════════════════
║ <EXCEPTION> <DVP_1>
║
║ Program error.
║
║ Condition not met:
║ .false.
║ File
║ /home/salome/workspace/SALOME-MECA_V2021_PLATFORM/3/salome/V2021/tools/src/Code_aster_stable-1
║ 540/bibfor/load_util/ischar_iden.F90, line 99
║
║
║ Il y a probablement une erreur dans la programmation.
║ Veuillez contacter votre assistance technique.
╚════════════════════════════════════════════════════════════

Any hint on what DIDI is really doing or why SUIV is not working is very welcome and appreciated.

Attachments:
buckling_tube_caster_tryD.gif, Size: 868.02 KiB, Downloads: 76

jlucas · 2022-10-19 23:05:39

== Files for tryD ==

Attachments:
tryD.zip, Size: 353.16 KiB, Downloads: 67

jeanpierreaubry · 2022-10-20 09:06:13

hello

i am really interested in your problem

but for god's shake consider providing a less refined mesh, many (10 ?) times coarser to speed up the conclusion

for trials we do not care that much about things like the stress refinement near the boundary condition

once we have a conclusion we can revert to a more realistic mesh

and spending well over half a hour on each trial calculation on my modest komputer is for me a waste of time!

two practical questions while i am waiting my first run to finish

why is there no FONC_MULT with CHARGE=load,
why using AFFE_CHAR_CINE

jean pierre aubry

jlucas · 2022-10-20 09:59:58

Hi,

I can not really argue about begin consistent by using such options. But sometimes unknown inconsistencies lead to great discoveries (probably not this time!).

I think the use of these options should not affect much the result.

AFFE_CHAR_CINE instead of AFFE_CHAR_MECA because remember to read somewhere that such type of condition should be computationally more efficient if using homogeneous boundary conditions.

U4.44.03:

Cette commande peut être utilisée avec un modèle mécanique, thermique ou acoustique. Le traitement de ces conditions "cinématiques" se fera sans dualisation et donc sans ajout de degrés de liberté de Lagrange.

The non use of FONC_MULT with CHARGE=load, can not justify other than saying that followed an example (see Rotation.zip) in the following post.

post for Rotation.zip file

The reason for a relative fine mesh is because, I think the buckling condition will not be captured if using a coarser mesh. This is because the simulation will stop with some convergence problem (large distortion of elements) before reaching buckling. Also a large tube is required for buckling to occur which implies a considerable number of elements. I must say, I also tried an adjusted version of the 2d example you provided without success. (But maybe here the problem is the structured mesh?)

However, I understand your point and will try to provide a coarser mesh where we reach the buckling phenomena.

jeanpierreaubry · 2022-10-20 10:07:15

you maybe right with the mesh fineness
i had to interrupt the calculation before the end because of non convergence around INST=7
however the results which i cannot understand seem very strange to me
do you get valid results with these files

jeanpierreaubry · 2022-10-20 12:29:41

are you sure of the FORMULE for tournex and tourney

jlucas · 2022-10-20 13:20:58

Hello,

My calculation does also not reach the end. The last computed time step is about 0.7482. I include a plot for the computed torsional stiffness (load path). In this simulation the instability is occurring for an angle of about 2.86 rad (163 deg) at a twist moment of about 92 kN.m. I guess the strange results you are refering to relate with the disagreement with expected analytical calculation. The curious thing is that the load path computed with CalculiX agrees pretty well (check the load_path_calculix.png also inside the calculix.zip on a previous post of this topic). The other curious thing is, by following your suggestion, and by using coarser mesh (6269 tetrahedral elements) the instability is reached at an even larger twist angle of about 7 rad and moment of 200kN.m with both simulations having a comparable torsional stiffness. That said I think valid results, are much dependent on the mesh and our current mesh is already quite coarse.

PS: The shared files include in addition to the plot for the load path, (i.e torsional_stiffness.png), a plot of the trajectory of a node on the cLoad group (selected randomly), a table with the computed total forces/moments and twist angle (same data as the used in the torsional_stiffness plot) and a python script. The python script might be useful for others to read and merge code_aster tables as such should be done to compute the moments when considering large displacement hypothesis. Note that for this case tables DEPL and REAC_NODE need to be read and merged (i.e. synchronized) on the nodes to give access to compute the actualised position of the node to be used in the computation of the reaction moment on each node (i.e: vec_M(t) = _vec_r(t) x vec_F(t)). The toal reaction moment/forces is obtained by summing up the forces/moments on all nodes of the selected nodeset (in this case the loaded nodes on cLoad)

Attachments:
analysis_tryD_files.zip, Size: 84.69 KiB, Downloads: 72

jlucas · 2022-10-20 13:49:39

No, i am not sure for tournex and tourney, specially when using TYPE_CHARGE='DIDI'.

These formulas were copied from other forum posts and can be found on some official example test (ssnv507) always with TYPE_CHARGE='FIXE_CSTE'.

However, noticed that the formula for large rotations might be other than the assumed:

vec_r = vec_r0 + theta*(vec_n x vec_r0 )

see attached pdf taken from:

Dhondt, G. "The Finite Element Method for Three-dimensional Thermomechanical Applications", John Wiley & Sons Ltd, 2004

Attachments:
Dhondt2004__pp155-159.pdf, Size: 53.51 KiB, Downloads: 81

jeanpierreaubry · 2022-10-20 13:53:43

before dealing with all that i had like to have your opinion about the FORMULE tournex and tourney included in the command file

tournex produces a non null value for FX at X=0

i suppose DTH is the nodal tangential load

do not le the tree hide you the forest

jlucas · 2022-10-20 20:02:29

The formula for tournex/tourney is very similar to the cylindrical to rectangular coordinate transformation.

x = r cos(theta)
y = r sin(theta)

with:

r = sqrt(X^2+Y^2)

the radius of the node on a XY plane and

theta = atan2(Y,X) + DHT*INST

the angular increase on the increment.

That was my interpretation.

Maybe AsterO'dactyle can clarify better as he proposed such formula on the post as a

"rotational" boundary condition

without further explanation.

jeanpierreaubry · 2022-10-21 08:42:16

as far as i can understand
in the attached document disk;comm the formulaes for tournex and tourney seems valid with:

1- DTH= 2*PI radians is the total amount of rotation at the end of the problem (over 2Pi would lead to a nasty behavior)

2- the applied load is an imposed rotation angle, not a force or moment

3- are related to a certain spatial position of the object

try to calculate the result of the formulaes for a given node from INST=0 to INST=1 to find out in your case

jlucas · 2022-12-19 12:08:35

Hello, Sorry this is a test as I was not able to reply on this topic for a while...

jlucas · 2022-12-19 13:00:19

Since I could not really understand how tryD was working (with TYPE_CHARCE='DIDI'). I've tried other things .

I've got the same result by applying twist in a more reasonable way. The procedure is based on test case zzzz364a, where a rotation is applied via rotation matrix (see for example en.wikipedia.org/wiki/Rotation_matrix -- the formula in use is related with the one found under the topic "Rotation matrix from axis and angle" and is the same as the one found in the previously attached file: Dhondt2004__pp155-159.pdf).

#1 2022-08-15 18:12:37

= Buckling of a tube under twist (non-linear) calculix vs code-aster =

#2 2022-08-15 18:22:31

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#3 2022-08-15 18:23:35

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#4 2022-08-15 18:25:45

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#5 2022-08-15 18:27:12

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#6 2022-08-15 18:28:35

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#7 2022-08-15 18:30:06

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#8 2022-08-15 18:32:27

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#9 2022-08-15 18:33:37

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#10 2022-08-17 08:35:12

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#11 2022-08-17 13:50:08

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#12 2022-08-19 10:08:55

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#13 2022-10-19 23:03:08

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#14 2022-10-19 23:05:39

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#15 2022-10-20 09:06:13

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#16 2022-10-20 09:59:58

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#17 2022-10-20 10:07:15

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#18 2022-10-20 12:29:41

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#19 2022-10-20 13:20:58

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#20 2022-10-20 13:49:39

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#21 2022-10-20 13:53:43

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#22 2022-10-20 20:02:29

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#23 2022-10-21 08:42:16

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#24 2022-12-19 12:08:35

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

#25 2022-12-19 13:00:19

Re: = Buckling of a tube under twist (non-linear) calculix vs code-aster =

Board footer