Large Rotations with NL solvers - observations / questions

Sound_Spinning

code_aster V16.7.9 - all run from the terminal with run_aster *.export.

I've been trying to get my head around NL solvers in CA for large rotations on bodies.
I started 1st with a 1 POU_D_T_GD beam, and that worked up to 2PI rotations; if rot >= 2PI it breaks, but below that it worked.
I then modelled a solid tet part and tried to drive it either with an angle disp or with a moment history. It works fine up to exactly 1 radian; if rot >= 1rad it breaks again. I used the Energy Balance table from static / dyn NL solvers tests; and the values do show when it goes pear shape >= 1rad, values explode. No matter what I tried couldn't get over the 1rad issue...

I then converted my solid tet part little test into a simple pendulum driven by gravity only, so no enforcing of angle or moment histories at pivot. I have at pivot point a soft coupling (LIAISON_RBE3), the ref (pivot) node with 6DoF so I can query / export the rotations.
To my surprise it all works in terms of x,y,z (translation) disps; however, when I plot the angle of rotation at the pivot node I get the graph you see in animation below, all between 1 <--> -1 units angle.
So at -90deg angle I get -1 unit, at -180deg I get 0 units and at -270deg I get 1 unit angle again, doesn't make sense; but the model seems to do the right thing. This is with the dynamics NL solver with numerical damping:

I also read extensively in help docs about NL solvers, with many caveats and (potential) limitations. All examples in docs / forum I found claiming to do large NL rots, appear to be imposed via x,y,z disps on nodes on surfaces, which I consider a cheat as a general method, as it makes that surface fully rigid, which is no good. Plus you need to know where nodes will be, which defeats itself the main purpose of large rots solve with NL solvers, with flexible bodies.

Before I give up on this tho, has anyone dealt with these cases with large rots imposed on bodies; and have you got any workarounds? Since the pendulum works, I tend to think that it is (maybe) the updating of the angle on those nodes, rather than the actual solve? What's your opinion?

Thanks,
Jesus

PS: I also tried rigid couplings (LIAISON_SOLIDE) with and w/out follower loads, as the help docs imply; it is even worse vs soft coupling.

mf

Sound_Spinning

Hello,
I am not sure about what's going on here, but the graph looks more like a coordinate (something like y/r) rather than an angle. Why should the angle be limited to +/- 1 rad (what is the unit, anyway?): it clearly shouldn't. The y-coordinate becomes ambiguous when the pendulum surpasses the horizontal plane on reversing, the two peaks go away with a lower amplitude. Also, when the pendulum is at the bottom, the angle can't be = 0 like it is here, but y is. Just try a large initial angle like 135° here: https://burchill.nz/javascript-simulations/pendulum-angle/ It looks similar,
Mario.

Sound_Spinning

mf
Thanks for sparing your thoughts.
The angle is wrong for sure. I also always output a table of disps/rots at both points you see: pivot and middle one. The latter is a good one for a check, the X, Y disps are perfect in the output table, while the angle RZ is equally wrong vs pivot one in graph.
When I started with the 1 beam model the angle in rads in tables is correct, till it breaks at 2PI. So, it seems to be something to do with updating the angle value correctly on nodes on couplings or similar when solving them as NL. I think that may be the reason why I couldn't drive my 1st tests with an angle or moment history.
It's so close tho... the disps are all correct. It'd be such a good thing to be able to do MBD like assemblies with the NL solver, as it is possible in other codes.

I found this old thread on a similar subject about large rots:
https://forum.code-aster.org/public/d/19093-dynamics-for-rotating-body-dyna-non-line/10

vishwas

Hi,
Maybe you can use paraview functions.
Try plotting arctan(y/x) for some point.

Anirudh

Sound_Spinning

vishwas
Have done that, no need for PV; I use the tables output from code_aster in text results file, it's all there.

My question remains: has anyone achieved large rots with solid (tets, etc) bodies with NL solvers in CA?

jacob

do you mean large rotations with some kind of pin connector?

Sound_Spinning

jacob
Yes, it is either a pin or a hinge, doesn't really matter as long as we only have body rots around the DoF we want. Like in my pendulum example, but that gives the wrong angle at pivot.
It's basically a usual MBD (multi body dynamics) setup, nothing special for a NL solver.

jacob

how did you even managed that?
If I use LIAISON_SOLIDE
_F(CHARGE=RBE2,
TYPE_CHARGE='SUIV')

I get

<EXCEPTION> <MECANONLINE5_40> ║
║ ║
║ La dynamique avec des chargements de Dirichlet de type suiveur. ║
╚═══════════════════════════════════════════════════════════════════

Sound_Spinning

jacob
That's correct, rigid couplings in CA with rots DoF at ref node, do not play ball with NL solvers, have tested it.
As opposed to what the help docs claim...
see table:

jacob

if I use RBE3 without SUIV, I am getting also -1, 1

jacob

well, the docs is talking only about STAT_NON_LINE, so maybe it works for static.

However, even though the results looks ok, I don't think the reactions etc. will be ok for RBE3.
Maybe the way how to do it, is to create some spokes made of POU_D_T_GD from the center node. Maybe that's the way how commercial solvers do that.
Other way would be to update the constrains each time step, but it would be time-step dependent.

Sound_Spinning

jacob
I tried, it doesn't. I did a lot of testing (both statics & dyna NL) before posting the pendulum. I don't think CA can handle real large rots, coz they are historically tricky. Large angles in FEA are not trivial; e.g. a 3D compound rotation must follow an order. Once you've done a 1st turn say around the global X axis, the other 2 axes have now moved... etc.
Also, we have the maths challenge with e.g. SIN(PI) == SIN(-PI) kind of thing. I just don't think it is trivial. It took well known commercial solvers a while to get it right; e.g. Abaqus.
But, I'll give it one last shot in CA, just in case.

PS: POU_D_T_GD is limited to < 2PI, as I mentioned in my 1st post. After that, it breaks.

jacob

what do you mean it breaks? Non convergence may be general issue. When I tried the pendulum with one POU_D_T_GD I needed to add

                               _F(EVENEMENT='DELTA_GRANDEUR',

                                  NOM_CHAM='DEPL',

                                  NOM_CMP='DRZ',

                                  VALE_REF=0.1)),

otherwise, the simulation blowed up even for smaller angle than 90°.

I am trying to do it with beams, like I mentioned, but I am not able to converge after 50°.

Sound_Spinning

jacob
breaks = it doesn't solve when it gets to 2PI angle, below it no prob. That is with 1 beam only.
When I tried to rotate solids (via angle disp or moment ramp) the limit seems to be 1rad, just below 60deg.
Can you pls expand where those commands you posted go? Never seen that before, and CA syntax is not trivial.
Also, I was gonna read about 'Piloting', does it have anything to do with NL loads?

vishwas

Hello,
Long back I was able to simulate this.
http s://youtu.be/1vuVxv7OGKI?si=_rkPf78ioOOML_qD
My youtube channel.

Anirudh

Sound_Spinning

vishwas
Nice. That is with an external load, like the pendulum. But what I'm after is (large) rotating a body via an angle or moment history on single ref node on a soft coupling (*_RBE3). Can your model do that for angles > 60deg with that solid you show?

jacob

yes, with rbe3 it looks somehow correct, but somehow the rotation in the master node is not.

jacob

But on the gif above you have solids with more than 180 deg.

To defi list instant

Sound_Spinning

jacob
That's correct, coz the solid is driven by an external load, gravity. Then the x,y,z disps appear correct but the angle reported at pivot is wrong (the graph). But when you try to drive the solid with an angle disp or a moment history is when it 'breaks' @ >= 1rad. I think both are related, the wrong rot angle reported in the pendulum and the inability to let the pivot being driven by a large rot/moment BC/load history.

I found in the solver code, in file calc_bt_ops.py this function, which seems to shift angle values, but I'm yet to digest what it really does and how/where is used.

    # ==============================================================================
    def red_ang_disp(ang):

        """ """
        ang = ang[np.isnan(ang) == False]
        ad_1 = np.max(ang) - np.min(ang)

        ang_2 = np.zeros((len(ang), 1))
        ang_2 = ang
        ang_2[ang > pi] = ang_2[ang > pi] - pi
        ad_2 = np.max(ang_2) - np.min(ang_2)

        ang_3 = np.zeros((len(ang), 1))
        ang_3 = ang
        ang_3[ang > 3 * pi / 2] = ang_3[ang > 3 * pi / 2] - pi
        ang_3[ang < pi / 2] = ang_3[ang < pi / 2] + pi
        ad_3 = np.max(ang_3) - np.min(ang_3)

        if ad_1 == min(ad_1, ad_2, ad_3):
            return ang
        elif ad_2 == min(ad_1, ad_2, ad_3):
            return ang_2
        elif ad_3 == min(ad_1, ad_2, ad_3):
            return ang_3

mf

Sound_Spinning

Wow, a textbook example of how NOT to write code in Python, this is making me angry. No docstring, meaningless variable names, no type hints, mixing of libraries (np is in there but pi is called from somewhere else, maybe math?), no info about the units...

It takes a list or an np.array (which one should it be??) of angles and applies some math on these values. From -4pi to 4pi the output looks like this:

What's its purpose? We don't know, because nobody filled the docstring! Anybody learning Python right now seeing this? Don't do it that way, undocumented code is the source of all problems,

Mario.

Sound_Spinning

mf
Great work! Thanks for looking into it.
The question I have: where / how is this angle function used? I am not familiar with the code_aster code, yet. Can you help to find where this is used, and see if we can uncover this mystery with large rots? I'd like to dig in more into the subject.