yade-users team mailing list archive

[Rohit John <question699273@xxxxxxxxxxxxxxxxxxxxx>]

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 = 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()
O.saveTmp()
plot.plot()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.