Bending beams mass computation

Here the complete code for copy pasta:

# encoding: utf-8
"An example showing various bending beams."

from builtins import range
from yade.gridpfacet import *

#### Parameter ####
L=10. # length of the beam
n=12 # number of nodes used to generate the beam
r=L/50. # radius of the beam element

#### Engines ####
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_GridConnection_Aabb()]),
 InteractionLoop(
  [Ig2_GridNode_GridNode_GridNodeGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 NewtonIntegrator(gravity=(0,0,-10),damping=0.5,label='newton')
]

#### Create materials and set different properties ####
O.materials.append(CohFrictMat(young=1e6,poisson=0.3,density=1e1,frictionAngle=radians(10),normalCohesion=1e7,shearCohesion=1e7,momentRotationLaw=False,label='mat1'))
O.materials.append(CohFrictMat(young=1e6,poisson=0.3,density=1e1,frictionAngle=radians(10),normalCohesion=1e7,shearCohesion=1e7,momentRotationLaw=True,label='mat2'))
O.materials.append(CohFrictMat(young=3e6,poisson=0.3,density=1e1,frictionAngle=radians(10),normalCohesion=1e7,shearCohesion=1e7,momentRotationLaw=True,label='mat3'))
O.materials.append(CohFrictMat(young=1e7,poisson=0.3,density=1e1,frictionAngle=radians(10),normalCohesion=1e7,shearCohesion=1e7,momentRotationLaw=True,label='mat4'))

#### Vertices ####
vertices1=[]
vertices2=[]
vertices3=[]
vertices4=[]
for i in range(0,n):
  vertices1.append( [i*L/n,0,0] )
  vertices2.append( [i*L/n,1,0] )
  vertices3.append( [i*L/n,2,0] )
  vertices4.append( [i*L/n,3,0] )

#### Create cylinder connections ####
nodesIds=[]
cylIds=[]
cylinderConnection(vertices1,r,nodesIds,cylIds,color=[1,0,0],highlight=False,intMaterial='mat1')
cylinderConnection(vertices2,r,nodesIds,cylIds,color=[0,1,0],highlight=False,intMaterial='mat2')
cylinderConnection(vertices3,r,nodesIds,cylIds,color=[0,0,1],highlight=False,intMaterial='mat3')
cylinderConnection(vertices4,r,nodesIds,cylIds,color=[1,1,1],highlight=False,intMaterial='mat4')

#### Set boundary conditions ####
for i in range(0,4):
   O.bodies[nodesIds[i*n]].dynamic=False
# O.bodies[nodesIds[i*n]].state.blockedDOFs='xyzXYZ'
# O.bodies[nodesIds[i*n]].state.blockedDOFs='xyz'

#### For viewing ####
from yade import qt
qt.View()

#### Set a time step ####
O.dt=1e-05

#### Allows to reload the simulation ####
O.saveTmp()

#### Added mass computation ####
print(f"\nBody zero is a {O.bodies[0].shape} with a mass of {round(O.bodies[0].state.mass,3)}.")
ana_node_mass = (4/3)* r**3 * math.pi * O.bodies[0].material.density
print(f"The mass of a Node/sphere is {round(ana_node_mass,3)} if analytically computed.")

sim_total_mass = sum( [b.state.mass for b in O.bodies] )
print(f"\nSimulation total mass: {round(sim_total_mass,3)} kg.")
ana_beam_mass = r**2 * math.pi * L * O.bodies[0].material.density + ana_node_mass # because the cylinders/beam include a half-sphere at the top and bottom
ana_total_mass = 4*ana_beam_mass # because there are 4 beams
print(f"Analytic total mass: {round(ana_total_mass,3)} kg.\n")

Hello,

>Can someone explain to me how the mass is computed in the simulation and why it differs quite a lot to my analytical computation?

I am not an expert in GridConnections. But I'll add some code to help you continue thinking about it. Yade holds the Gridnodes as bodies as well as the GridConnections as bodies. Hence len(O.bodies) > len(cylIds). I would say these volumes overlap based only on the viewport, so would be better to sum with a GridConnection criteria maybe:

sim_total_mass=0
for b in O.bodies:
  if isinstance(b.shape,GridConnection): sim_total_mass+=b.state.mass

But when I do this I get 23. And if I switch it from GridConnection to just summing mass on GridNode, I get 46. It seems the volume is allowed to overlap and the GridConnection mass estimate is incorrect, or it isn't accounting for the overlaps with GridNodes.

Hope it helps.

Robert

Hello,
thanks for you thoughts and sorry for the late reply.
I tried several things to see how the mass is computed. It turns out there might be even more confusion.

If only nodes are added to the simulation, none of them have any mass at all (0.0kg).
As soon as CylinderConnections are added the nodes and the cylinders receive a mass, that is "incorrectly" computed as shown above.

The point that they overlap is a good one. This might be the reason for a higher mass in the simulation, but as you already have shown this doesn't lead to a valid result either.

I found a paper that deals with those connections in yade [1].
In there it is told that "The mass of the cylinder element is lumped into its two nodes."
Which would explain the higher masses of the nodes, but on the other hand would mean that the cylinders don't have any mass at all.
Did I intrepret that correctly?
Is the paper acutally refering to this exact yade-functions or did they use something else?
Since they used nodes and cylinder connections I thought this might be the same approach.

Would you mind sharing your thoughts on that?

Thanks for your help and best regards,
Paul

[1] https://www.researchgate.net/publication/282996677_A_general_method_for_modelling_deformable_structures_in_DEM

Here is a fix for your test. :)

#### mass computation ####
realLength= (O.bodies[n-1].state.pos-O.bodies[0].state.pos).norm() # first node to last node of first beam
ana_beam_mass = r**2 * math.pi * realLength * O.bodies[0].material.density
totalActiveMass = sum( [O.bodies[id].state.mass for id in nodesIds] )
print("analytical mass",ana_beam_mass," vs. numerical mass ",totalActiveMass/4)
##########

Output:
analytical mass 11.519173063162576 vs. numerical mass 11.519173063162569

Arguably the function "cylinderConnection()" can be confusing. It produces a node to to node distance which is less than L/N since it includes the spherical caps in L.
OTOH the total mass is defined by neglecting the spherical caps, hence an apparent inconsistency.

Another thing is that the cylinderConnections are never dynamic, they are just folloing the nodes. They turn out to have have a mass because some constructors define one by default, but it is not actually used anywhere. So total mass of the beam is the mass of its nodes.

For very short objects it might not be a good idea to neglect the mass of spherical caps. On the other hand it would be painfull to compare to beam theory if the mass was not uniformly distributed (currently it is, with end nodes having 0.5*mass of center nodes, it wouldn't with spherical caps).

In the end, all this is just pre-processing by some python helper function. You can easily assign mass per-node the way you like it, and you could also bend a modified version of cylinderConnection() to your liking.

HTH
Bruno

Thanks for your detailed answer Bruno!
Nice to hear from one of the authors of the papers that build the foundation on my current research/work.

To make sure I understood it correctly:
The mass of the cylinders can be completely ignored. Only the mass of the nodes is relevant for the simulation and its results (total mass, deflection, deformation, forces etc). Hence, I only need to set the right density for the CohFrictMat for the nodes and I should be good to go in terms of correct body mass (if nodes have a different material than cylinders). To avoid confusion, setting the density of the cylinderConnection material to zero manually won’t do anything. Right?

Or do I need to set the density of both materials (for nodes and connections) to the same value, so the computation does not get messed up and I simply ignore cylinder-mass afterwards? I am reffering here to one of the creator-functions e.g. gmshPFacet().

In other words, what do I need to do in order to set the mass of a body, that is made out of nodes and cylinders (and maybe PFacets), to a certain value?

Seems like I am overthinking this and it might be just a simple comparison between wanted mass and mass of all nodes and setting the density of all influenced bodies to corresponding value, but I just want to be sure and safe about that.

For my application I want to simulate a beam with a rectangular cross-section. So, I will be using many nodes, cylinder/grid connections and PFacets. Hence it will not be a very short object and I should be “safe” as small inaccuracies like the spherical caps won’t bother my results.

Furthermore, PFacets seem to never receive a mass, but also the mass of the nodes does not increase as it did when connections were added.

Thanks again for you answer.

Since I work a lot with deformable structures in yade based on these to papers [1][2] I would have several more questions according that approach.
Can you recommend or provide resources for that topic? Any further papers, examples or theory? That would help me a lot.
As an example, what interests me I will write down some questions. Nevertheless, they have nothing to do with mass computation, so this is off-topic.
- What’s the difference between grid- and cylinder-connections.
- Why is it important to set a cohesive material for the nodes with super high cohesive values?
- How is the bending/mechanical behaviour in a node that connects 2 cylinders computed? (Same thing along a connection that connects 2 PFacets)

Thanks Bruno and I hope I did not flood you with questions.
Your help is and was very much appreciated!

Best regards,
Paul

[1] https://www.researchgate.net/publication/282996677_A_general_method_for_modelling_deformable_structures_in_DEM
[2] https://doi.org/10.1016/j.geotexmem.2015.07.015

> You are welcome, nice to hear that someone is using it.

I can guarantee you that it will not be the last time I post questions on this topic on the yade launchpad. For simple deformable bodies, my team and I will most likely use it a lot in the future as soon as I know exactly how to use it and managed to properly simulate correct beahviour of a simple body. 😊

> Set the mass of its nodes. That's all.

Alright, thanks. Just applying the correct node density to achieve the same total mass that I want it to have. Total node mass = wanted total mass. Got it!

> I would say they are two alias for the same thing. Am I confused? have seen one or the other in different scripts?

As far as I know they are multiple similar approaches in the “More examples” section of the yade documentation. Just to summarize:
- Chained-cylinders
- Cylinders (or cylinder-connections. The example I initially used to open this question thread.)
- Grids
- PFacet (although this uses creators and gridNodes and gridConnections and only adds the PFacet)
 Since I am a newbie, I am not yet able to tell what the exact differences are, but they are more or less used for the same things (I believe).

> if cohesion is small you can get get plastic deformation of the beam

Ah I see. So, it is kind of needed to apply resistance to pulling. As particles only “generate” a force when moving towards, but not when moving away from each other.
Linear elastic behaviour is totally fine for my use for now. Let’s see where the journey goes.

> Bending is defined for interactions, based on the differential rotation of two nodes. It is not defined "per node".

Since I had a look in the paper you linked: Is it correct to imagine the bending behaviour as a rotational spring like in figure 3 of [3]? This what I initially meant when I proposed the question on how the bending behaviour “in a node” is computed. As a node itself can’t bend it is somewhat clear that it needs to happen in the interactions as this links multiple nodes.

> Effeindzourou PhD dissertation

Nevertheless it is a good summary, thanks for the hint!

> There is an old 2D version of the same thing described in [3]

As in the question above this gave some basic insights (probably). So this is still of value. Thanks for the suggestion!

Best regards and thanks again,
Paul

[3] https://www.yade-dem.org/w/images/1/1b/Chareyre%26Villard2005_licensed.pdf

What I forgot:

> Currently the only option is to do that in your script. If you make it in a python function we can include it to source alongside the cylinderConnection().

I am totally happy for now if I am able to set the total mass of the body correctly so I do not need to consider the volume of the PFacets as I already know the total mass that the body should have.
If it becomes important for me I will provide an helper function for that!

Cheers

> Is it correct to imagine the bending behaviour as a rotational spring like in figure 3 of [3]?

Mmmmh... I would say no, it isn't the right picture. The rotational spring should be together with the normal spring, in between two inertial objects (the nodes).

It could be that the curvature was really evaluated per-node, as in fig. 3, at the time. I would have to check in "fish" scripts. But for sure it is not the case now. The rotational and normal/shear springs are together in between two nodes.

Bruno

> The rotational and normal/shear springs are together in between two nodes.

Alright, thanks for your input!

I am starting to get the hang of it.

Cheers,
Paul