Angular acceleration problem

Dear Readers,

I have a case including 3D meshes and point masses connected to the meshes with MPC:s. Load is 6DOF acceleration, three linear components and three angular acceleration components around remote origin axles. The angular rotations cause additional linear components related to the distances from axles. Rotational speed is so low, that centrifugal loads are neglectable. Case is totally linear at this stage.

This should be a simple case, but I can’t figure out how to feed the calculated local linear accelerations back to CalculiX from user dload.f since there should either be three components, one for each rotation, or their combination vector with orientation. And can one only feed the additional components from user subroutine and use normal *DLOAD, GRAV for the linear components.

I assume that either loadtype=2 (gravity with vector) might be a solution, or loadtype=3 (generalized gravity). However, dload.f seems not to have either orientation or component option with acceleration, if the other parameters don’t have hidden options for such use. So how could this be done and if not at all, what approach could be used?

Best regards,


Couldn’t you just connect the mass to the reference point of a rigid body or coupling constraint (attached to the rest of the model) and apply a displacement BC with a proper amplitude to this reference point in a *DYNAMIC step ?

Hi, thanks for your comment. Unfortunately that doesn’t represent what’s needed. The meshes (structures) and mass points attached to them are on top of a huge moving rigid body. Meshed parts are objects of interest. There are multiple load cases, having different motion origins and their own accelerations. Fortunately only linear responses are needed. There doesn’t exist enough data from the cases for a proper dynamic analysis, only “snapshots” of demanding acceleration values. That’s why I am trying to get it done this way, something which is possible with an *ANGULAR ACCELERATION BC in Abaqus. Since that load option does not exist in CalculiX, I have to calculate in a user load subroutine element integration points’ and masses’ resulting linear accelerations.

Do you have some scheme of the problem?. Something simple just to understan it.

A simple 2-D scheme could be as such (Z coordinates omitted):
A rigid body has top line from (0,0) (X,Y) to (10,0).
Mesh is a rectangle (0,0) to (1,0) to (1,1) to (0,1), filled with material.
There is a point mass in (2,0.5) attached to the mesh.
Gravity affects the model downwards Y.
At the time moment of of interest the rigid body has no linear horizonal or vertical acceleration, but it has a centre of rotation at (9,0.5), around which it has known angular acceleration aa.
Resulting vertical acceleration is gravity+8.5aa for the mesh and g+7aa for the point mass. They can be smaller or bigger than 1g, depending on the aa sign.

So this is a very simple schema showing how the angular accelerations affects local linear accelerations. In reality structures modelled with large meshes have different accelerations for each elements’ integration points due to different distances from each three rotation axis. Also there are 3 linear+3 angular accelerations present, angular accelerations having separate independent (no common origin) perpendicular rotation axels. That’s the case type I am studying with numerous acceleration datas.

I’m not sure that user dload is the correct approach, user cload could also be usefull. However I haven’t succeeded to find a way to get either integration point masses or volumes and material dencities into user subroutine for calculating the forces. Both should be in some arrays, but how to access them from user subroutine?

Another question is, whether there is a mechanism to return local linear accelerations as a flux from user dload subroutine. The only parameter in user dload that has three components is displacement. Is that an overloaded parameter, which could carry the acceleration components in return?

Sorry for the sudden change of font, I can’t remove it!

Also I have tried to identify the fortran module where accelerations are set to the masses and integration points, that could be another route to get the wanted effect. Could somebody plese give a hint which module that is!

50% of solving a problem is understanding it and it’s hard to get an idea of it.

Sounds scary

What does top line means?.Is it some kind of pendulum.

We can forget about this by the moment. Just to simplify.

You are talking about time and accelerations. That is not a static problem. I guess you are looking for the accelerations to know the instantaneous forces acting on the body.

Don’t you mean it has no horizontal or vertical “velocity”?. You comment there is Gravity downwards.

The rotation of a Rigid body is arround the REF node position. Is your ref node at (9,0.5) or at (10,0).?

If you rotate the model 90º and look at it in vertical position it seems like a building with discrette masses under seismic loading. Small w—> Small Normal component and horizontal proportional to the height. Alpha *r in your case.

Regarding the multiple rotations all acting at the same time,…If you can only consider a linear aproach some sort of principle of superposition would be required.

Ok, let’s take it again. I simplified a 3D-case down to 2D by taking a slice perpendicular to one rotation axis through the model. Other axels, being parallel to the slicing surface, disappear in 2D case. Since the model is attached on top of a rigid body, the only meaningful item of it in this case is it’s top surface, which reduces to top line in the 2D slice. Rotation axis reduces to a rotation point (you called it REF node), which I located to (9,0.5) so that there are no cos and sin terms present in rotation lever arms, also didnt take any other linear accelerations than gravity in account. As said, in real world rigid bodies have up to six active DOF:s, as is the situation in my case, the motions and accelerations are spectral, from which we have only “snapshots” from some instantaneous time moments, e.g. totally static conditions. Time is not present, and superposition is valid. Please don’t get confused with possible rotations, they only rotate the gravity vector and the world doesn’t turn upside down in the analysis. Constant velocity does not generate mass forces, so velocities are omitted.

Loads are forces caused by the accelerations affecting on any masses: element integration points, point masses etc. Since there is no angular acceleration load option present in CalculiX, one has either to generate equivalent forces to the masses somehow or manipulate the acceleration routine.

For that purpose I’m studying different options how that could be done. Both user cload and user dload have some necessary parameters and returning options missing from their calls. The more I study them, the more acceleration routine manipulation sounds tempting. What’s your opinion?

With user cload for each node
-How to get the necessary mass information into user subroutine.
-How to get list of all elements connected to a node.
-Calculating fractions of joining element masses. Sounds to be too big job.
-Calculating spatial accelerations and mass forces for the node.
-User cload seems to have idof in call, is it called from CalculiX automatically for all DOF:s?

With user dload for each integration point and point mass
-How to get integration point masses into the routine? Are they in an array in Common area? Haven’t found them. For efficiency, I’d like to use the CalculiX internal mass information instead of recalculating it. By default it must have been done before the steps.
-Calculating spatial acceleration and mass force to the integration point.
-Defining loadtype to acceleration or mass force (if they are available).
-Returning acceleration or mass force.
-How to return acceleration direction/componets? This is not explained in the manual.

Manipulating acceleration routine
One could write own acceleration interpretation routine to feed spatial accelerations into mass force calculation. So far I have consumed about ten pots of coffee browing thrue the fortran codes. Any hint of modules involved would be of great help!

I think now I see the problem.

You have already decomposed your rotational acceleration into Normal and Tangential.
As w is small you are disregarding the Normal (centripetal acceleration) but not tangential that depends on the distance to the rotating point that can be large. At(R)=alpha*R.

So, the forces on your body have two sources. One is proportional to the mass (g) and the other is proportional to the distance to the rotation point at the moment of interest.
But we only have available g in Calculix and it applies to the whole system without the possibility to make it spatial dependent.


What about using the concept of apparent weight?

If you are only interested in Fy that formula could also be read as:

Fi=mi*(g+Aplhari)= mi(g+Aplharig/g)= mi*(1+Aplha*ri /g)*g=mi’*g

Being mi’= mi*(1+Aplha*ri/g).

That is a linearly increasing density with the distance to the rotation point.
Extending that to 3D and diferent axis would need about ten pots of coffee more.

Thanks for your reply. Yes, you are right, in my case the centripetal acceleration can be neglected, in general case not. I have been browsing the code a bit, and it seems possible to manipulate gravity forces bodyf(1…3) when they are given values. So far I have found that to happen in subroutine mafillsm. Is that the correct place for all mechanical cases? If I leave the summation intact, centrifugal forces should also be present, right?

I have three angular accelerations with their own rotation axels. At this stage I try to get it working to some level with hard coded data in the acceleration routine. For this I would like to get some confirmations for my variable findings:
-bodyf(3) is gravity vector for one item at a time in mafillsm.f. Same routine in rhs.f.
-bodyfx(3) is also acceleration vector into which bodyf is copied in e_c3d.f
-bodyf in replaced with its value multiplied with density in e_c3d.f
-ff is bodyf*xmass in e_c3d, likely to be mass force caused by acceleration, but where is it’s value stored? If every element has it’s own value in every integration point, that’s exactly where I would like to get mass forces stored.

So far I have no succeeded to follow the entire chain of acceleration inputs and storage, looks like there is only a fixed value for it. Is that correct finding?
So to get done, what I would like, the subroutine call chain should be able to provide individual accelerations or mass forces. The latter does probably not work with ff, because of the array size ff(60). But that’s probably not the general store for the values, I still need to find that.

So any hints are most wellcome, and if there are any flow charts and/or variable meaning descriptions, they are more than wellcome! And if you see any reasons why the case couldn’t be solved at all this way, please let me know so I don’t loose our time and mine. There are other options available.

Best regards,

1 Like

I am very interested in this problem for multibody dynamics. At any instant, any body is subjected to joint forces and torques, applied forces and torques, inertial forces from linear acceleration, angular acceleration, centripetal acceleration and Coriolis acceleration, and gravitational force. I would like to use these forces and torques on the body to calculate stresses and deformations in the body. I would like this to be automated as much as possible. Is there any work being done for this?
Aik-Siong Koh

Are you the dev of the new assembly solver for FreeCAD ? Nice to see you here. Have you seen the very recent work of Jose Egas done with MBDyn and CalculiX ?

Yes, I developed OndselSolver for FreeCAD integrated assembly.
Yes, Jose and I have been collaborating on MbD and CalculiX.
Of course, the video is incomplete because of the lack of angular acceleration and Coriolis acceleration in CalculiX.
Hence, I am here to ask and help to make that capability available in CalculiX.
I am glad to see discussions on angular accelerations already.
What is the status? How can we finish it?
Aik-Siong Koh

You would have to talk with the (only) CalculiX dev - Guido Dhondt. His e-mail is here: About the authors

CalculiX has centrifugal and gravity loads meant mostly for static analyses but no acceleration boundary conditions for dynamic analyses (apart from the acceleration base motion for linear dynamics).

Thanks. I emailed Guido.

No response from Guido. Can someone prompt him to look in his spam folder? ( Thanks.

Can CalculiX allow user to specify distributed body loads that are a function of x,y,z and element mass? Thanks.

Some ideas:

  • Apply a different gravity load to each element.
  • Write a user subroutine dload.f
  • Wild nutty idea - use the NEWTON label in *DLOAD and cleverly position fixed mass elements around the place to define the nonuniform gravity field.

Just my personal opinion here - Guido surely has a lot on his plate and I hope he spends time on more fundamental work like bug fixes than convenience features like this that can be done with preprocessing.

1 Like

What I want and CalculiX needs to be complete is to handle angular acceleration.
We already have Centrif to handle centripetal acceleration (R x Omega^2) in the radial direction.
So, it would be very simple to handle acceleration due to angular acceleration which is (R x Alpha) in the perpendicular direction.

I am looking at “Apply a different gravity load to each element.” as a workaround. I would appreciate an example file since I am a total newbie using FreeCAD FEM workbench.


Sure it would be nice but I think it’s not too hard to write this using dload.f or a script to apply a *DLOAD on each element. Yea it’s some extra work but it’s more achievable for a 3rd party user than, say making shells, beams, prescribed rotations, etc. work correctly.

Prepomax makes traction loads using *CLOAD, for instance, and I’ve done similar things like that without having to recompile CCX. Maybe ask the Prepomax author?

PrePoMax doesn’t have a solution for that. Since Aik-Siong Koh wants to combine FreeCAD MBD solution with CalculiX, he should find a way to handle this with Python scripting so that the ccx input deck is modified properly. FreeCAD FEM has a keyword editor and full scripting interface. As shown in the video I linked above, Jose Egas already managed to handle a similar problem, just without modifying the .inp file directly.