Large Rotations with NL solvers - observations / questions

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

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.

mf

Sound_Spinning
Hello,
it only seems to be used in the file where you found it: calc_bt_ops.py

The function is defined in line 1736 and used in line 214. The function it is used by is called "angle_mean". At least this is commented, it says:

Description:
Computes the angular mean of a group of angles contained in a cell.
Input:
Mesh_ .- FE nodla list
U .- Cathetus meassured along x axis
V .- Cathetus meassured along y axis
vertices - nodal list of the cells containing the data
regions .- Connectivity matrix of the cells
Output:
mean_CQ .- average slope computed for each cell

I searched the whole codebase of code_aster (src-main) with VSCode, this is the only use of this function, the other hit in the search is its own definition. Beat's me what it does, it is assigned to a var named "reduced_ang".

Mario.

If you want to take a closer look at the values, rename this to *.py and run it in Python....(i don't know why the extensions are still wrong in this forum)

disp-red-ang-disp.bin

1kB

Sound_Spinning

mf
Thank you! I'll take a look.

mf

Sound_Spinning

I forgot to say, you could also modifiy the code to print some useful info in your own test case like:

       if sum(AS[:, i_]) != 0:
            reduced_ang = red_ang_disp(angle_[AS[:, i_]])
            print(f"****** {reduced_ang = }  {angle_[AS[:, i_]] = }  {other_stuff = }")

They are 'only' Python scripts, so you can do what you want to get more info out of this without compiling CA. Don't know if the print really appears in the .mess though... :-(

Mario.

jacob

I don't know what are you trying to achieve.

LIASION_SOLIDE won't work for large rotations, maybe for large displacements in static and rbe3 reports the rotation wrong because of the coupling, otherwise the results seems ok.

The solution with POU_D_T_GD seems to work, but for large angles it diverges, which could be maybe solved by some dynamics solver settings.

Quite simple way how to model a hinge only with solids is to use only solids.

and just constrain the displacement in those nodes. However if you want to report rotation, then I guess it must be done in postprocessing as if you would have a contact there.

Sound_Spinning

jacob
All good, but you are talking about how-to model s'thing, which is not my problem.
With solids only, you have no control on driving / monitoring NL assemblies on rots Dof, which is key for good quality stuff: large angles, moment histories, ang velocities and ang. accelerations all matter in NL mechanisms, and being able to load/query at all relevant points in the model.
I'm not trying to solve one model, but seeing if CA can offer a general method for large disps / rots. Also, I'm not so confident yet that the NL disps (x,y,z) results are bang on, just coz it looks OK by eye. Working on it on a test...

jacob

https://forum.code-aster.org/public/d/23652-weird-behavior-of-pou-elements-with-grot-gdep/2

For my simulation with POU_D_T_GD with R=10, 20 there is no problem, for R = 30 it diverges around 1 rad.

Sound_Spinning

jacob
That sounds like the same issue I found with my 1st test models (not posted in this thread) driven via angle or moment histories.
Also, using beams we don't know what the stiffness change is on cylindrical surfs, it'll over constrain them. Hence the use of soft (RBE3) couplings, they allow curved surfaces to ovalise realistically, w/out adding artificial stiffness.

Thanks for sharing, very useful.

Sound_Spinning

I've now done some more post-processing on the simple pendulum case, see if I could belive the NL results from CA.
I modelled the same pendulum with the ABAQUS Learning Edition, same transient dyn solver: SCHEMA = 'HHT', ALPHA = -0.075. The small alpha value should put only minor numerical damping, which apt for a free pendulum.
I also tweaked the defaults in CA for the auto time stepping setting, so that I could get very similar number of load increments on both solvers, which I did, both close to 350 time incs. This is for a like-to-like comparison, so that the load inc doesn't jump/increase too aggressively.
I also did a hand calc on the angle based on predicted X,Y disps at the mid point shown in my pendulum animation, in the middle hole.

1.- Disps & angles CA vs Abaqus vs angle hand calc:

The 1st thing I noticed is that CA prediction on the Y disp curve shape may not make sense, it seems to follow the (wrong) angle prediction shape, which not sure is correct for circular motion.
Also, the 'drift' observed with CA on mag disp & disp_X curves doesn't feel right. The only reason I can thing is that the numerical damping in CA is applied in a very odd way, maybe ....

2.- Energy Balance:
I turned on in CA the energy calcs: ENERGIE = _F( CALCUL = "OUI", ).
Here is the comparison to Abaqus:

Although I used a small alpha value, CA clearly applies a lot of numerical damping; and we see the 'drifting' upwards again, which I'm not sure what it means in real life on such a simple model; it only spins, no 'travelling in global translations' involved as a rigid body.

I must conclude for now that CA is not up to scratch on large disps / rots as a general method. That doesn't take away that it may work for some simple models, and useful for comparative work as a guide.

Furthermore, SimScale (I believe) have recently dropped code_aster as their NL solver in favour of Marc (Hexagon). There must be a reason for it, since they've been at it for 10+ years I reckon.

I'm all ears if anyone has real examples where the NL solvers in CA have shined, would be good to know.

And thanks for your time on this thread too, appreciated.
Cheers,
Jesus

« Previous Page Next Page »