simulation of uniaxial compression test with precracked

Hi:

I want to simulate the closed pre-crack example in paper[1] . In paper[1], the author use two methods to creat the model. And I have some questions about creating pre-crack model. Here is my question:

1. About method 1 "removing the cohesive feature of all interactions located along the flaw surface". Can i do this using 'i.phys.isCohesive=0' to creat precrack in model. If i can, how can I identify the particles on either side of the crack.

2.About method 2 "using smooth joint method". I tried using SJM in a uniaxial compression simulation. But there was a problem with the simulation, and particles disappeared when the simulation started. Here is my script:

###########
from yade import ymport, utils, plot,pack

readParamsFromTable(noTableOk=True, # unknownOk=True,
 intRadius=1.25,
 dtSafety=.8,
 damping=0.4,
 strainRateTension=.05,
 strainRateCompression=0.5,
 setSpeeds=True,
 # 1=tension, 2=compression (ANDed; 3=both)
 doModes=3,

 specimenLength=.1,
 specimenRadius=0.025,
 sphereRadius=0.0016,

 # isotropic confinement (should be negative)
 isoPrestress=0,
)

#material
sphereMat=O.materials.append(JCFpmMat(type=1,density=2000,young=3.18e9,poisson=0.3,tensileStrength=6e6,cohesion=6e6,frictionAngle=radians(10)))

#body
pred = pack.inCylinder((0,0,0),(0,0,.1),0.025)
sp=pack.randomDensePack(pred,spheresInCell=2000,radius=0.0016,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=sphereMat)

#uniaxial function
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=dtSafety*PWaveTimeStep()

import gts
# joint
v1 = gts.Vertex(0.01,0,0.05)
v2 = gts.Vertex(0,0.01,0.05)
v3 = gts.Vertex(-0.01,0,0.05)
v4 = gts.Vertex(0,-0.01,0.05)

e1 = gts.Edge(v1,v2)
e2 = gts.Edge(v2,v4)
e3 = gts.Edge(v4,v1)
f1 = gts.Face(e1,e2,e3)

e4 = gts.Edge(v4,v3)
e5 = gts.Edge(v3,v2)
f2 = gts.Face(e2,e4,e5)

s1 = gts.Surface()
s1.add(f1)
s1.add(f2)

facet = gtsSurface2Facets(s1,wire = False,material=sphereMat)
O.bodies.append(facet)

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

O.saveTmp()
# identification of spheres onJoint, and so on:
execfile('identifBis.py')

#O.engines
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True)],
 ),
 NewtonIntegrator(damping=damping,label='damper'),
 UniaxialStrainer(dead=1,strainRate=-strainRateCompression,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
 PyRunner(dead=0,virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotdata',initRun=True),

]

#plot stress-strain curve
plot.plots={'eps':('sigma',)}
plot.plot(subPlots=False)
def addPlotData():
 yade.plot.addData(t=O.time,i=O.iter,eps=-strainer.strain,sigma=-strainer.avgStress/1e6)
 plot.saveDataTxt('/home/liuchangdong/MyYade/install/bin/precrack.txt',vars=('eps','sigma'))

O.timingEnabled=False
#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
is2aabb.aabbEnlargeFactor=-1.
#SIMULATION REALLY STARTS HERE
strainer.dead=0
O.run(20000,1)
###########

Thanks for your help

changdong Liu

[1]https://www.sciencedirect.com/science/article/pii/S0013794415007225
doi:10.1016/j.engfracmech.2015.12.034