about cohFrictMat model

Dear all,
I've been trying to model the cohesive frictional model under uniaxial force. The problems are:
1. when the normalCohesion is low (eg.1e7), the stress-strain relationship behaves wired, going up and down rapidly. At first, I thought it was due to the repulsive force at the beginning, so I make i.phys.unp = i.geom.penetrationDepth as #266828[1] and set dt=0 when create the interaction as # 630910[2]. but still cannot resolve this issue.
2. when I set both normalCohesion and shearCohesion to 1e9 and start counting the crack number, I find that the number of crack fluctuate a lot. this made me confused because I've set the setCohesionOnNewContacts=False, how would the number of cracks decrease if there was no new bond forming? By the way, the pattern of the crack is also a bit strange, what I am expecting is that near the peak strength, there supposed to be a surge in the crack number and remain steady after the peak.
3. also the model seems have no difference whether I set fragile as True or False, any suggestions is appreciated!
Thanks,
Libby

[1].https://answers.launchpad.net/yade/+question/266828
[2].https://answers.launchpad.net/yade/+question/630910

#######code####
from yade import pack,qt,plot
isoPrestress=0
sphereRadius=.015
normalCohesion=1e9
shearCohesion=1e9
RollingStiffness=1

idConcrete=O.materials.append(CohFrictMat(young=23.1e9,poisson=0.4,frictionAngle=atan(0.5),density=2036,normalCohesion=normalCohesion,shearCohesion=shearCohesion,fragile=False))
spp=randomDensePack(inAlignedBox((0,0,0),(.5,.5,2)),radius=sphereRadius,rRelFuzz=1/3,returnSpherePack=True,spheresInCell=2000,material=idConcrete) # memoizeDb='/tmp/bending2DPackCache.sqlite'
spp.toSimulation()
for b in O.bodies:
  b.shape.color=(0,1,0)
### define constant variables
factor=1.5
damping=.4
strainRateTension=50
## render the scale of the displacement
renderer=qt.Renderer()
renderer.dispScale=(10,10,10) ## set the render scale
bb=uniaxialTestFeatures() ## extract the information for the uniaxial test
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
#O.dt=0.5*PWaveTimeStep() ## set critical time step
O.engines=[
  ForceResetter(),
  InsertionSortCollider([
    Bo1_Sphere_Aabb(aabbEnlargeFactor=factor,label='bo1s'),],verletDist=.05*sphereRadius), ## expand the collision detection to creat initial interactions
      ##
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=factor,label='ig2ss'), ## create collision geometry
    ],
    [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)], ## creat collision phys
    [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False)],
   ),
  #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
  NewtonIntegrator(damping=damping,label='damper'), #UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer'),## apply strain on the both end of the body along the axis and the strain speed is set at the beginning directly
  PyRunner(iterPeriod=1,command='addPlotData()'),
  PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'), ## set the stop criteria for the simulation
  #qt.SnapshotEngine(fileBase='3d',iterPeriod=1000,label='snapshot'),
  #PyRunner(command='finish()',iterPeriod=20000
]
#sigma=strainer.avgStress+isoPrestress

## set equlibrum for overlapping spheres
#for b in O.bodies:
  #b.dynamic=False
O.dt=0
O.step() ## go one step to creat interactions
bo1s.aabbEnlargeFactor=1 ## reset the interaction radius
ig2ss.interactionDetectionFactor=1
O.dt=0.3*PWaveTimeStep()
for i in O.interactions:
  i.phys.cohesionDisablesFriction=True
  i.phys.unp = i.geom.penetrationDepth
O.engines=O.engines[:4]+[UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer')]+O.engines[4:]
#for b in O.bodies:
 #b.dynamic=True
### reinforcing the boundary
dim=utils.aabbExtrema()
if strainRateTension>0:
  layerSize=.05
  height=dim[1][axis]-dim[0][axis]
  for b in O.bodies:
    if isinstance(b.shape,Sphere):
      if (b.state.pos[axis]<(dim[0][axis]+layerSize*height)) or (b.state.pos[axis]> (dim[1][axis]-layerSize*height)):
        b.shape.color=(1,1,1)
  for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      if O.bodies[i.id1].shape.color==(1,1,1) or O.bodies[i.id2].shape.color==(1,1,1):
        i.phys.normalAdhesion*=1000
        i.phys.shearAdhesion*=1000

def stopIfDamaged():
  if O.iter<2 or 'sigma' not in plot.data : return
  sigma=plot.data['sigma']
  extremum=max(sigma)
  if abs(strainer.strain)>5e2 : ## stop the test if the sigma is too small or strain is too large ##abs(sigma[-1]/extremum)<0.5 or
   print("damaged")
   O.pause()
#####count crack number
def addPlotData():
  crackCount=0
  for i in O.interactions:
    if i.phys.cohesionBroken==True:
      crackCount+=1
  print('crack number=',crackCount)
  plot.addData(i=O.iter,sigma=strainer.avgStress,crackNumber=crackCount,eps=strainer.strain)
plot.plots={'eps':('sigma'),'eps ':('crackNumber')}
plot.plot()
O.saveTmp()