Shear force is lower than MaxShearForce but sliding occurs

Hi,

I am testing on a two-grain case, from which I would like to explore when the particle slides under a given force.

In the MWE below, it demonstrates two grains, the bottom one (s1) is fixed, while a permanent force is applied to the upper one (s2). s2 first moves down and attaches to the bottom one (to build interaction), and then s2 is set to allow free translations. So if the shear force (Fs) is less than the maximum shear force (Fs_max), it should not slide. However, I found that s2 slides when Fs is less than Fs_max.

I must misunderstand something, or there is something wrong with the MWE code. Do you have any idea about this problem?

Thanks,
Leonard

########### MWE
from yade import pack, plot
import math
import numpy as np

sandR=0.01
distance=sandR/1
appliedF=1000

frictionAngle=30

young=1e9

KrKt=0.1
eta=0.1

sand = O.materials.append(CohFrictMat(isCohesive=True,young=young,alphaKr=KrKt,alphaKtw=KrKt,\
                                         poisson=0.3,frictionAngle=radians(frictionAngle),etaRoll=eta,etaTwist=eta,\
                                         density=2650.0,normalCohesion=0, shearCohesion=0,\
                                         momentRotationLaw=True,label='sand'))

s1=utils.sphere((0,0,0),radius=sandR,color=[1,1.,1.],fixed=True,material='sand')
s2=utils.sphere((distance,sandR+sandR,0),radius=sandR,color=[1,0,0],dynamic=True,material='sand')

O.bodies.append(s1)
O.bodies.append(s2)

Gl1_Sphere.quality=3
s1.state.blockedDOFs='xyzXYZ' ## no freedom for s1
s2.state.blockedDOFs='xzXYZ' ## only vertical movement at this state.

O.forces.setPermF(1,(0,-appliedF,0)) ## apply a permanent force on s2

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="Ip2Coh")],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)]
    ),
    NewtonIntegrator(damping=0.6),
]

O.step()

setContactFriction(radians(frictionAngle))

O.run(500000,1) ## run sometime to let s2 move and attach s1

def getForce():
    s2.state.blockedDOFs='XYZ' # allow s2 to slide if any
    result=[]
    for i in O.interactions:
        fn=i.phys.normalForce.norm()
        fs=i.phys.shearForce.norm()

        fsmax=fn*np.tan(np.deg2rad(frictionAngle))

        result.append(fn)
        result.append(fs)

        result.append(fsmax)

    return "Shear force is",result[1],"Max shear force is",result[2],

print "please execute getForce() and then click run, you will see s2 slides"