Angular momentum of pfacet clump not conserved after collision

Asked by Rohit John on 2021-06-09

Hello,

In the following simulation, two sphere are fired in opposite direction at a pfacet box clump. Ideally after collision, the box starts rotating and the spheres move in the opposite direction (opposite to their corresponding initial velocity) such that the angular momentum of the whole system is conserved. I want to change the inertia and the mass of the pfacet box clump, so I did that using the lines

O.bodies[target_clump_ID].state.inertia = Vector3([10,10,10])
O.bodies[target_clump_ID].state.mass = 1

When I do this, the angular momentum of the system is not conserved. Am I not seeing something here, have I made a mistake?

Moreover, in the comment #20 of [1], it was found that for pfacet clump to have correct dynamics, we should set the mass of the cylinder connections to 0 (as the mass is supposed to be concentrated at the nodes). But when I do that I am no longer able to access the angular momentum of the clump using O.bodies[target_clump_ID].state.angMom as it returns Vector3(0,0,0). How do I fix this? The only work around I found was to set the mass close to zero (1e-10) instead of 0. But is there a proper fix for this?

Kind regards,
Rohit K. John

https://answers.launchpad.net/yade/+question/695558

# GTS file
# --------------------------------------------------------------------------------------------------------------------------------------------------------
14 36 24 GtsSurface GtsFace GtsEdge GtsVertex
0.5 0.5 0.5
0.5 0.5 -0.5
0.5 -0.5 0.5
0.5 -0.5 -0.5
-0.5 0.5 0.5
-0.5 0.5 -0.5
-0.5 -0.5 0.5
-0.5 -0.5 -0.5
0.5 0.0 0.0
0.0 -0.5 0.0
0.0 0.0 0.5
-0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 -0.5
6 8
2 6
1 2
8 7
3 4
5 6
3 7
1 3
8 4
7 5
5 1
4 2
9 3
2 9
4 9
1 9
10 7
4 10
8 10
3 10
11 3
5 11
7 11
1 11
12 6
7 12
8 12
5 12
13 2
5 13
6 13
1 13
14 8
2 14
6 14
4 14
7 21 23
7 17 20
6 25 28
9 33 36
8 13 16
3 29 32
12 14 15
5 15 13
3 16 14
9 18 19
4 19 17
5 20 18
11 22 24
10 23 22
8 24 21
4 26 27
1 27 25
10 28 26
6 30 31
2 31 29
11 32 30
2 34 35
1 35 33
12 36 34

# script
# --------------------------------------------------------------------------------------------------------------------------------------------------------
from yade.gridpfacet import *
from yade import *
from yade import geom, utils
from yade import plot
import sys, os
sys.path.append(".")

# ---------------------------------------------------------------------------- input parameter
# ----------------------------------------------------- target
target_young = 50e9
target_density = 1000
target_poisson = 0.3
target_friction = radians(30)

p_radius = 5e-2

# ---------------------------------------------------------------------------- Materials
target_int_mat = 'pfacet_int_mat'
target_ext_mat = 'pfacet_ext_mat'

O.materials.append(
    FrictMat(
        young = target_young,
        poisson = target_poisson,
        density = target_density,
        label = target_ext_mat,
        frictionAngle = target_friction,
    )
)

O.materials.append(
    CohFrictMat(
        young = target_young,
        poisson = target_poisson,
        density = target_density,
        label = target_int_mat,

        frictionAngle = target_friction,
        normalCohesion = 3e100,
        shearCohesion = 3e100,
        momentRotationLaw = True,
    ))
# ---------------------------------------------------------------------------- Engines
O.engines = [
                ForceResetter(),

                InsertionSortCollider([
                    Bo1_GridConnection_Aabb(),
                    Bo1_PFacet_Aabb(),
                    Bo1_Sphere_Aabb(),
                ]),

                InteractionLoop(
                    [
                        Ig2_PFacet_PFacet_ScGeom(),
                        Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
                        Ig2_GridNode_GridNode_GridNodeGeom6D(),
                        Ig2_GridConnection_PFacet_ScGeom(),
                        Ig2_Sphere_PFacet_ScGridCoGeom(),
                    ],
                    [
                        Ip2_FrictMat_FrictMat_FrictPhys(),
                        Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(
                            setCohesionNow = True,
                            setCohesionOnNewContacts = False
                            ),
                    ],
                    [
                        Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
                        Law2_ScGeom_FrictPhys_CundallStrack(),
                        Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
                        Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
                    ],
                )
            ]
# ---------------------------------------------------------------------------- objects
# ----------------------------------------------------- target
(
pnode,
pcyl,
pfacet
) = gtsPFacet(
    'cube.gts',
    radius = p_radius,
    shift = (0,0,0),
    scale = 1,
    wire = False,
    fixed = False,
    color = [0.1,0.5,0.1],
    materialNodes = 'pfacet_int_mat',
    material = 'pfacet_ext_mat',
)

target_ids = pnode + pcyl + pfacet
# for i in pcyl:
# O.bodies[i].state.mass = 0
# See https://answers.launchpad.net/yade/+question/695558 comment #20 for why I am doing this
# If this is uncommented, then we are not able to access the angular momentum of the clump using
# O.bodies[target_clump_ID].state.angMom

target_clump_ID = O.bodies.clump(target_ids)
O.bodies[target_clump_ID].state.inertia = Vector3([10,10,10])
O.bodies[target_clump_ID].state.mass = 1
print(target_clump_ID)

# ----------------------------------------------------- spheres
sp1 = sphere([ 0.61, 0.3, 0], 5e-2, material = target_ext_mat)
sp2 = sphere([-0.61, -0.3, 0], 5e-2, material = target_ext_mat)

sp1_ID = O.bodies.append(sp1)
sp2_ID = O.bodies.append(sp2)

O.bodies[sp1_ID].state.vel = [-1,0,0]
O.bodies[sp2_ID].state.vel = [ 1,0,0]
# ---------------------------------------------------------------------------- Additional engines
ids = target_clump_ID
O.engines += [
    NewtonIntegrator(gravity = [0,0,0], damping = 0),
    PyRunner(command = 'graph()', iterPeriod = 1000)
]

# ----------------------------------------------------- plotting
plot.plots = {'t':('wx', 'wy', 'wz'), 't1':('Lz')}

def graph():
    w = O.bodies[ids].state.angVel
    pos_sp1 = O.bodies[sp1_ID].state.pos
    vel_sp1 = O.bodies[sp1_ID].state.vel
    angMom_sp1 = O.bodies[sp1_ID].state.angMom
    mass_sp1 = O.bodies[sp1_ID].state.mass
    L_sp1 = angMom_sp1 + mass_sp1 * pos_sp1.cross(vel_sp1)

    pos_sp2 = O.bodies[sp2_ID].state.pos
    vel_sp2 = O.bodies[sp2_ID].state.vel
    angMom_sp2 = O.bodies[sp2_ID].state.angMom
    mass_sp2 = O.bodies[sp2_ID].state.mass
    L_sp2 = angMom_sp2 + mass_sp2 * pos_sp2.cross(vel_sp2)

    angVel_clump = O.bodies[target_clump_ID].state.angVel
    mmoi_clump = O.bodies[target_clump_ID].state.inertia
    mmoi_mat_clump = Matrix3([
        mmoi_clump[0], 0, 0,
        0, mmoi_clump[1], 0,
        0, 0, mmoi_clump[2]
    ])
    L_clump = mmoi_mat_clump * angVel_clump

    L_tot = L_clump + L_sp1 + L_sp2
    angMom_clump = O.bodies[target_clump_ID].state.angMom
    # L_tot = angMom_clump + L_sp1 + L_sp2
    print(angMom_clump)

    plot.addData(
        t = O.time, t1 = O.time,
        wx = w[0], wy = w[1], wz = w[2],
        Lx = L_tot[0], Ly = L_tot[1], Lz = L_tot[2],
        )

plot.plot()
# ---------------------------------------------------------------------------- GUI controller
from yade import qt
qt.Controller()
qtv = qt.View()
qtr = qt.Renderer()
qtr.light2 = True
qtr.lightPos = Vector3(1200,1500,500)
qtr.bgColor = [1,1,1]
qtv.ortho = True # orthographic view in viewport
Gl1_Sphere.stripes = True

O.dt = 2e-7#utils.PWaveTimeStep()
O.saveTmp()

Question information

Language:
English Edit question
Status:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:

Can you help with this problem?

Provide an answer of your own, or ask Rohit John for more information if necessary.

To post a message you must log in.