Clump angular momentum not conserved

Asked by Rohit John

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


# 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'

        young = sphere_young,
        poisson = sphere_poisson,
        density = sphere_density,
        label = sphere_mat,
        frictionAngle = sphere_friction,

# ---------------------------------------------------------------------------- Engines
O.engines = [


# ---------------------------------------------------------------------------- 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 += [
        gravity = [0,0,0],
        damping = 0
        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

        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 = O.bodies[id].state.vel
    angMom_sp1 = O.bodies[id].state.angMom
    mass_sp1 = O.bodies[id].state.mass
    return angMom_sp1 + mass_sp1 * pos_sp1.cross(vel_sp1)

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

Question information

English Edit question
Yade Edit question
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Rohit John (rohitkjohn) said :

Hello all,

I derived the equations for getting the state after the collision using the following code

Kind regards,
Rohit K. John

# ------------------------------------------------------------------------ Script
from sympy import *
m, v0, l, iz, w1, v1 = symbols('m, v0, l, iz, w1, v1')

eq1 = Eq(2*m*v1*l + iz*w1 - 2*m*v0*l, 0) # Angular momentum conservation
eq2 = Eq(m*v1**2 + 1/2*iz*w1**2 - m*v0**2, 0) # Energy consevation

sol = solve((eq1, eq2),(v1, w1))


Revision history for this message
Karol Brzezinski (kbrzezinski) said :


what version of Yade and system are you using? I haven't tried your code yet, but it sounds very similar to the problem that I had [1].
Jan Stránský pointed out that there is a bug in the software version that I was using (incorrectly assigned clump inertia). Switching from yade 2020.01a to yadedaily solved the problem.



Revision history for this message
Rohit John (rohitkjohn) said :


I apologise for not mentioning the version of Yade I use. I am using Yadedaily. I will amend the initial description.

Kind regards,
Rohit K. John

Revision history for this message
Karol Brzezinski (kbrzezinski) said :


I noticed that b.state.angMom for spheres is zero even though they have some angular momentum. I think that it doesn't affect the simulation. You can compute momentum for spheres manually. I changed one line in your function get_AngMom():

angMom_sp1 = np.array(O.bodies[id].state.angVel)*np.array(O.bodies[id].state.inertia)

Maybe there is a more elegant way to do this instead of changing Vector3 to NumPy array, but it was just a quick fix.

After these modifications numbers look good to me (lower plot looks crazy but it is just a scale issue).

Best wishes,

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.