spheres rotations go unstable

Hi,

I am simulating a test in which material is sandwiched between two plates and the upper plate is rotated with a sinusoidal function.
I ran the model with CohFrictMat and things go fine.
However, when I try CohBurgersMat (a new class of material derived from CohFrictMat) spheres start to rotate even though the interactions are cohesive.
What's missing is these interactions?

This is the model script for CohFrictMat case:

O.reset()
from yade import utils, plot
from yade import pack, qt

O.engines=[
ForceResetter(),
GlobalStiffnessTimeStepper(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True, setCohesionOnNewContacts=True)],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(damping=0.7,gravity=[0,0,-9.81])
]

id_Mat1=O.materials.append(CohFrictMat(normalCohesion= 1e15, shearCohesion= 1e13, isCohesive= True, young=1e5,
density=2600, poisson=0.3, frictionAngle= 0.4))
CohMat=O.materials[id_Mat1]

lower_plate = O.bodies.append(box((0,0,-0.0003),(.006,.006,.0002),fixed=True,material=CohMat))

upper_plate = O.bodies.append(box((0,0,0.0046),(.006,.006,.0002),fixed=False,material=CohMat))

fc=yade.geom.facetCylinder((0,0,0.001), 0.0051, 0.005, orientation=Quaternion((0, 0, 1), 0), segmentsNumber=100, wallMask=6, angleRange=None, closeGap=False, radiusTopInner=-1, radiusBottomInner=-1,material=CohMat)
O.bodies.append(fc)

pred=pack.inCylinder((0,0,0),(0,0,2.5e-3),5e-3)

O.bodies.append(pack.regularOrtho(pred,radius=0.0005,gap=0.0, material=CohMat))

print(len(O.bodies))

for x in range(len(O.bodies)):
 if (O.bodies[x]):
                if isinstance(O.bodies[x].shape,Sphere):
   O.bodies[x].shape.color=(0,0.3,0.7)
  else:
   O.bodies[x].shape.color=(0.3,0.3,0.3)

O.forces.addF(1,(0,0,-0.02),permanent=True)

O.run(40000,True)

O.bodies[0].state.blockedDOFs='xyzXYZ'
O.bodies[upper_plate].state.blockedDOFs='xyzXYZ'

def calcul():
 w=math.sin(O.time*25.12)*1.0
 O.bodies[upper_plate].state.angVel=(0,0,w)

O.engines=O.engines+[PyRunner(iterPeriod=1,command='calcul()',label='Rate')]
O.engines=O.engines+[PyRunner(iterPeriod=1,command='addPlotData()',label='plotDataCollector')]

plot.plots={'t':('eps',),'t ':('pos'),' t':('n')}

def addPlotData():
       yade.plot.addData(t=O.time,pos=O.bodies[upper_plate].state.pos[2],eps=O.bodies[upper_plate].state.angVel[2],
   n=len([i for i in O.interactions]))
plot.plot()

#====================

But if we run the script with CohBurgersMat, it goes crazy:

#===============================================

O.reset()
from yade import utils, plot
from yade import pack, qt

O.engines=[
ForceResetter(),
GlobalStiffnessTimeStepper(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_CohBurgersMat_CohBurgersPhys(),
Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys(setCohesionNow=True,setCohesionOnNewContacts=True)],
[Law2_ScGeom_CohBurgersPhys_CohesiveBurgers(),Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(damping=0.7,gravity=[0,0,-9.81])
]

wallMat=O.materials.append(FrictMat(young=1e4,poisson=0.3,density=2500,frictionAngle=2.5))

bMat=O.materials.append(CohBurgersMat(kmn=1.05e3,kkn=5.0e3,cmn=1.3e3,ckn=1.6e3,
kms=1.05e3,kks=5.0e3,cms=1.3e3,cks=1.6e3,normalCohesion= 1e12, shearCohesion= 1e12, isCohesive= True,
density=1500, poisson=0.3, frictionAngle= 0.4))

bMat1=O.materials.append(CohBurgersMat(kmn=2.05e3,kkn=5.0e3,cmn=1.3e3,ckn=1.6e3,
kms=2.05e3,kks=5.0e3,cms=1.3e3,cks=1.6e3,normalCohesion= 1e12, shearCohesion= 1e12, isCohesive= True,
density=1500, poisson=0.3, frictionAngle= 0.4))

bMat2=O.materials.append(CohBurgersMat(kmn=2.05e2,kkn=5.0e2,cmn=1.3e2,ckn=1.6e2,
kms=2.05e2,kks=5.0e2,cms=1.3e2,cks=1.6e2,normalCohesion= 1e12, shearCohesion= 1e12, isCohesive= True,
density=1500, poisson=0.3, frictionAngle= 0.4))

lower_plate = O.bodies.append(box((0,0,-0.0003),(.006,.006,.0002),fixed=True,material=wallMat))

upper_plate = O.bodies.append(box((0,0,0.0046),(.006,.006,.0002),fixed=False,material=wallMat))

fc=yade.geom.facetCylinder((0,0,0.001), 0.0051, 0.005, orientation=Quaternion((0, 0, 1), 0), segmentsNumber=100, wallMask=6, angleRange=None, closeGap=False, radiusTopInner=-1, radiusBottomInner=-1,material=wallMat)
O.bodies.append(fc)

pred=pack.inCylinder((0,0,0),(0,0,2.5e-3),5e-3)

O.bodies.append(pack.regularOrtho(pred,radius=0.0005,gap=0.0, material=bMat))

print(len(O.bodies))

for x in range(len(O.bodies)):
 if (O.bodies[x]):
                if isinstance(O.bodies[x].shape,Sphere):
   O.bodies[x].shape.color=(0,0.3,0.7)
  else:
   O.bodies[x].shape.color=(0.3,0.3,0.3)

O.forces.addF(1,(0,0,-0.02),permanent=True)

O.run(40000,True)

O.bodies[0].state.blockedDOFs='xyzXYZ'
O.bodies[upper_plate].state.blockedDOFs='xyzXYZ'

def calcul():
 w=math.sin(O.time*25.12)*1.0
 O.bodies[upper_plate].state.angVel=(0,0,w)

O.engines=O.engines+[PyRunner(iterPeriod=1,command='calcul()',label='Rate')]
O.engines=O.engines+[PyRunner(iterPeriod=1,command='addPlotData()',label='plotDataCollector')]

plot.plots={'t':('eps',),'t ':('pos'),' t':('n')}

def addPlotData():
       yade.plot.addData(t=O.time,pos=O.bodies[upper_plate].state.pos[2],eps=O.bodies[upper_plate].state.angVel[2],
   n=len([i for i in O.interactions]))
plot.plot()

#========================================================

You can have access to CohBurgersMat codes here:

https://www.dropbox.com/s/npodv33qultd2xo/CohBurgersMat.zip?dl=0