Question #699273 on Yade changed:
https://answers.launchpad.net/yade/+question/699273
Description changed to:
Hello all,
I have a simulation consisting of two spheres (S1, S2) and a clump made
of two spheres (C1, C2). The two spheres, S1 and S2, are given the same
initial speed but in the opposite directions along the horizontal. The
clump is aligned along the vertical and the two spheres, S1 and S2, are
made to hit the spheres in the clump at the same time. Could this be a
bug in my code or the NewtonIntegrator?
I expect the angular momentum to be conserved in the interactions, but
it is not. I have a similar problem with clumps of pfacet [1]. It was
reported in a bug report related to pfacet [2]. I am posting again
because this appears to be an error related to clumps in general. I have
put the code below.
I am using Yadedaily.
Kind regards,
Rohit K. John
[1] https://answers.launchpad.net/yade/+question/697465
[2] https://gitlab.com/yade-dev/trunk/-/issues/210
# script
#
from yade.gridpfacet import *
from yade import *
from yade import geom, utils
from yade import plot
import sys, os
#
input parameter
# - clump
sphere_young= 50e9
sphere_poisson = 0.3
sphere_friction = radians(30)
sphere_offset = 0.3
sphere_radius = 5e-2
sphere_init_vel = 1
sphere_mass = 0.5
sphere_density = sphere_mass / (4/3*pi*sphere_radius**3)
#
Materials
sphere_mat = 'sphere_mat'
O.materials.append(
FrictMat(
young = sphere_young,
poisson = sphere_poisson,
density = sphere_density,
label = sphere_mat,
frictionAngle = sphere_friction,
)
)
#
Engines
O.engines = [
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(),
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()],
)
]
#
objects
# - clump
clump_sph_1 = sphere([0, sphere_offset, 0], sphere_radius, material =
sphere_mat)
clump_sph_2 = sphere([0,-sphere_offset, 0], sphere_radius, material =
sphere_mat)
clump_sph_1_ID = O.bodies.append(clump_sph_1)
clump_sph_2_ID = O.bodies.append(clump_sph_2)
clump_clump_ID = O.bodies.clump([clump_sph_1_ID, clump_sph_2_ID])
clump_mmoi = O.bodies[clump_clump_ID].state.inertia
# - spheres
sp1 = sphere([ 4*sphere_radius, sphere_offset, 0], sphere_radius, material =
sphere_mat)
sp2 = sphere([-4*sphere_radius,-sphere_offset, 0], sphere_radius, material =
sphere_mat)
sp1_ID = O.bodies.append(sp1)
sp2_ID = O.bodies.append(sp2)
O.bodies[sp1_ID].state.vel = [-sphere_init_vel,0,0]
O.bodies[sp2_ID].state.vel = [ sphere_init_vel,0,0]
#
Additional engines
O.engines += [
NewtonIntegrator(
gravity = [0,0,0],
damping = 0
),
PyRunner(
command = 'graph()',
iterPeriod = 1000
),
]
# - Calculating state after
collision
m = sphere_mass
v0 = sphere_init_vel
l = sphere_offset
iz = clump_mmoi[2]
calculated_final_sphere_lin_vel = 0.5*(-4.0*iz*l*m*v0/(iz + 2.0*l**2*m) +
2.0*l*m*v0)/(l*m)
calculated_final_clump_ang_vel = 4.0*l*m*v0/(iz + 2.0*l**2*m)
calculated_final_clump_ang_mom = clump_mmoi[2]*calculated_final_clump_ang_vel
calculated_final_sphere_ang_mom =
2*sphere_offset*sphere_mass*calculated_final_sphere_lin_vel
print("clump: ", calculated_final_clump_ang_mom)
print("sphere: ", calculated_final_sphere_ang_mom)
print("total: ", calculated_final_clump_ang_mom +
calculated_final_sphere_ang_mom)
# - plotting
plot.plots = {'t':('clump_Lz', 'calculated_clump_final_Lz'), 't1':('Lz',
'calculated_Lz')}
def graph():
L_sp1 = get_AngMom(sp1_ID)
L_sp2 = get_AngMom(sp2_ID)
angMom_clump = O.bodies[clump_clump_ID].state.angMom
L_tot = angMom_clump + L_sp1 + L_sp2
plot.addData(
t = O.time, t1 = O.time,
Lx = L_tot[0], Ly = L_tot[1], Lz = L_tot[2],
clump_Lz = angMom_clump[2],
calculated= calculated_final_clump_ang_mom +
calculated_final_sphere_ang_mom,
calculated_clump_final_Lz = calculated_final_clump_ang_mom
)
def get_AngMom(id):
pos_sp1 = O.bodies[id].state.pos
vel_sp1 =
