about cohFrictMat model

Asked by libby on 2020-10-12

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()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2020-10-19
Last query:
2020-10-19
Last reply:
2020-10-16
Best Jérôme Duriez (jduriez) said : #1

Hi,

1. Sound to me like a typical "too high loading rate" issue

2. Your crackCount is the number of currently existing interactions which are now purely frictional after having been initially cohesive-frictional (i.e. after having lost their cohesive nature because of excess forces).

I do not see a reason for this number to be monotonously increasing.

If you want to consider the past (passed history of interactions breakage) and not only the present, you can not rely on "for i in O.interactions:" since this loop only considers present interactions.

In the JCFpm model for instance, there are bodies-related variables for this purpose, like JCFpmState.nbBrokenBonds

libby (libby-eleven) said : #2

Thanks Jerome!
1. That's true. After I decrease the strain rate, the model seems to behave normally
2. thanks for pointing out the mistake. I used 'initial bond - current bond' to count the number of cracks instead and it seems working!
############
for i in O.interactions:
  initialBond+=1
print('number of initial bond=',initialBond)
currentBond=0
  for i in O.interactions:
    if i.phys.cohesionBroken==False:
      currentBond+=1
  print('crack number=',initialBond-currentBond)
#################
is that an appropriate way to track the history of interactions breakage?

3. I prepare a simple script to see the effect of how fragile works, but then error occur:' ValueError: shape mismatch: objects cannot be broadcast to a single shape'. I can't figure out why, can you help me with this ?
###########
from yade import pack,qt,plot
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))
O.bodies.append([sphere((0,0,0),.5,fixed=True,material=idConcrete),sphere((0,0,-1),.5,material=idConcrete)])

### define constant variables
factor=1.5
damping=.4
strainRateTension=5
bb=uniaxialTestFeatures() ## extract the information for the uniaxial test
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=0.00005*PWaveTimeStep() ## set critical time step
O.engines=[
  ForceResetter(),
  InsertionSortCollider([
    Bo1_Sphere_Aabb(aabbEnlargeFactor=factor,label='bo1s'),]), ## 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)],
   ),
  UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer'),
  PyRunner(iterPeriod=10,command='addPlotData()'),
  NewtonIntegrator(gravity=(0.,0.,-1000)),
]
O.step() ## go one step to creat interactions
bo1s.aabbEnlargeFactor=1 ## reset the interaction radius
ig2ss.interactionDetectionFactor=1
O.interactions[0,1].phys.unpMax=-1

def addPlotData():
  plot.addData({'i':O.iter,'d':O.bodies[0].state.pos[2]-O.bodies[1].state.pos[2]-1,'f':O.interactions[0,1].phys.normalForce[2],'e':O.interactions[0,1].phys.unp})
plot.plots={'i':('f'),'i ':('e'),' i':('d')}
plot.plot()
O.saveTmp()
###############################

Jérôme Duriez (jduriez) said : #3

2- Looks like yes, at least under the assumption you're starting from initialBond=0

3- A known issue coming from an external library, see https://gitlab.com/yade-dev/trunk/-/issues/110. You can try modifying your use of plot.plot() / abandon live plotting (?)

libby (libby-eleven) said : #4

thanks, that's really helpful!

libby (libby-eleven) said : #5

Thanks Jérôme Duriez, that solved my question.