Impact of two differently sized spheres with a constant normal velocity and varying angular velocities

Hi,

I'm currently following along the paper "Benchmark tests for verifying discrete element modelling codes at particle impact level" by Y. C. Chung · J. Y. Ooi to verify yade as asked by my supervisor. In the benchmark test 8, it seeks to simulate the impact of two differently sized spheres with a constant normal velocity and varying angular velocities. However, from the code i had, the recoil vs. incident value seems to be quite off. I'm quite new to both yade and python, any help is appreciated!

the following is my script:

#benchmarks 8
from yade import plot
import math
import sys
import numpy as np

al_alloy=FrictMat(
 young=7.00e10,
 poisson=0.33,
 frictionAngle=np.arctan(0.40),
 density=2700,
 label='al_Mat')
O.materials.append(al_alloy)

dense_al=FrictMat(
 young=7.00e10,
 poisson=0.33,
 frictionAngle=np.arctan(0.40),
 density=2700*1000,
 label='al_dense')
O.materials.append(dense_al)

#nylon=FrictMat(
 #young=2.5e9,
 #poisson=0.40,
 #frictionAngle=np.arctan(0.4),
 #density=1000,
 #label'nylon_Mat')
#O.materials.append(nylon)

r=0.1
ang=4
vst=r*ang
incident=vst/0.2
O.bodies.append([
 sphere(center=(0,0,0),radius=r,material='al_Mat'),
 sphere(center=(0,0.65,0),radius=r*5,material='al_dense')
 ])

O.bodies[0].state.vel=(0,0.2,0)
O.bodies[0].state.angVel[0]=ang

#vy=O.bodies[0].state.vel[1]
#vz=O.bodies[0].state.vel[2]
#recoil=vz/vy

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom()],
 [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.5)],
 [Law2_ScGeom_MindlinPhys_Mindlin()]
 ),
NewtonIntegrator(damping=0.0,gravity=(0,0,0)),
PyRunner(command='recoilangle()',iterPeriod=1)
]

O.dt=0.1*utils.PWaveTimeStep()

def recoilangle():
 vy=O.bodies[0].state.vel[1]
 vz=O.bodies[0].state.vel[2]
 angular=O.bodies[0].state.angVel[0]
 vst=vz-0.2*angular
 recoil=vst/vy
 print 'Vcn:',O.bodies[0].state.vel[1],'Vct:',O.bodies[0].state.vel[2],'Angular:',O.bodies[0].state.angVel[0],'incident:',incident,'recoil angle:',recoil