# Angular momentum of pfacet clump not conserved after collision

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

# 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
# --------------------------------------------------------------------------------------------------------------------------------------------------------
import sys, os
sys.path.append(".")

# ---------------------------------------------------------------------------- input parameter
# ----------------------------------------------------- target
target_young = 50e9
target_density = 1000
target_poisson = 0.3

# ---------------------------------------------------------------------------- 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',
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
# 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]
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)

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
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:
Expired
For: