yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #25719
[Question #697465]: angular momentum of pfacet clump not conserved after collision
New question #697465 on Yade:
https://answers.launchpad.net/yade/+question/697465
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()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
