CohFrictMat fragile mode does not break interaction

Hi everybody !

I am trying to test what exactly do the fragile and plastic mode of CohFrictMat, and for this I have implemented a very simple simulation, appending a fixed sphere, and a second one who is moving the opposite way from initial contact, in order to break the cohesion bond by exceeding the tensile strength limit.

I am plotting Fn versus un displacement, and when I run the simulation using Fragile = False in the simulation I get first elastic increase of fn, and then flat plastic behavior, then the two particles are set apart, and interaction disappears. Which is fine.

However, when I use fragile = True, the same elastic behavior initially occurs, up until the fragile limit, then, at that precise timestep, no interaction is detected by Yade, and then this interaction is reset, with a very high value, and that does not drop to zero until the particles are really apart ( they are still overlapping when the strenght limit is reached, and that is because they are initially overlapping).

So I don't really know what is going on here. Could it be that it is the setting of the equilibrium distance again ? since unp relates to plastic displacement ?

any help would be appreciated !

Here is the script, running already until the strenght limit is reached, pressing run you can see what happens next !

# coding=utf-8
# TEST 8 : Appliance of horizontal disp one ball (rebound)-- > no grav/no coh
# EXPECTED : we want to see breakage of cohesion bond

from __future__ import division
import numpy as np
from yade import plot,pack,timing
from pprint import pprint

# DATA COMPONENTS

nb_iter = 1

# INPUTS
# Geometry
RSph = 0.1
CoeffSph = 0.9
Pos1 = (RSph*CoeffSph,RSph,0.)
Pos2 = (-RSph*CoeffSph,RSph,0.)

# Material
EyCohFrictMat = 1e7
poissonCohFrictMat = 0.25
densityCohFrictMat = 2600.
FrictAngSphere = 30.*pi/180.
NCohesionCohFrictMat = 30000
SCohesionCohFrictMat = 30000

# Calculation
damp=0.5
gz=0.
SphereVel= (-1.5,0.,0.)

#EXPORTATION

#création des différentes listes de forces

lisForce = []
lisForN = []
lisForSh = []
lisTime = []

################################################################################"
# SIMULATION BEGINNING - SIMULATION BEGINNING - SIMULATION BEGINNING

## CREATE MATERIAL - CREATE MATERIAL - CREATE MATERIAL

O.materials.append(
  CohFrictMat(
  young = EyCohFrictMat,
  poisson= poissonCohFrictMat,
  density= densityCohFrictMat ,
  frictionAngle = FrictAngSphere,
  normalCohesion = NCohesionCohFrictMat ,
  shearCohesion = SCohesionCohFrictMat ,
  fragile = True,
  momentRotationLaw=True,
  etaRoll=0.1,
  etaTwist=-1,
  alphaKr=0.,
  alphaKtw=0.,
  isCohesive = True,
  label='MatSpheres'))

## CREATE SPHERES - CREATE SPHERES - CREATE SPHERES -

# they could use the default material (utils.defaultMat)
O.bodies.append([
   # fixed: particle's position in space will not change
   sphere(center=Pos1,radius= RSph,fixed=True, material ='MatSpheres'),
   # this particles is free, subject to dynamics
   sphere(Pos2,RSph, material ='MatSpheres')
])

setBodyVelocity(1,SphereVel, 'xyz')

## FUNCTIONAL COMPONENTS - FUNCTIONAL COMPONENTS - FUNCTIONAL COMPONENTS -

O.engines=[

   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=1.5,label='ss2sc')], # collision geometry
      [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True)], # collision "physics"
      [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm = True)] # contact law -- apply forces
   ),
    #Apply gravity force to particles. damping: numerical dissipation of energy.
   NewtonIntegrator(gravity=(0.,0.,gz)),
   PyRunner(command='addPlotData()',iterPeriod=1, dead= True, label = 'dataPlot')

]

def setEquilibrium():
 O.dt = 0.
 O.step()

 for i in O.interactions :
   i.phys.unp = i.geom.penetrationDepth

 O.step()

 dataPlot.dead = False
 ss2sc.interactionDetectionFactor=1. # now reset the interaction radius and go ahead

 O.dt=5e-4*PWaveTimeStep()

def addPlotData():
 yade.plot.addData({'t':O.time,'i':O.iter,'Fn':O.interactions[0,1].phys.normalForce[0] ,'Fs':O.interactions[0,1].phys.shearForce[0], 'un':SphereVel[0]*O.time
   })
 #pprint (plot.data)

setEquilibrium()

plot.plots={'un':('Fn')}
plot.plot()

O.run(250,True)

yade.qt.Controller(), yade.qt.View()

cheers,
Manon