Unbreakable bonds in Tension Test

Hello everyone,

I had the following short conversation with Luc Scholts about the Uniaxial Tension test he performed in [1], and he thought this will also be of an interest to the forum.

Original message:
I am currently trying to duplicate your tension test from [1]. My only question is about the ubreakable bonds that you apply for particles in the recovery zone. How did you define these bonds? Did you define them as clumped, or did you just enter a bigger number for the tension (and Elastic modulus?) of the particles there?

Thanks in advance,
Yaniv

[1] [1]L. Scholtès and F.-V. Donzé, “A DEM model for soft and hard rocks: Role of grain interlocking on strength,” Journal of the Mechanics and Physics of Solids, vol. 61, no. 2, pp. 352–369, Feb. 2013.

Luc's response:
I just define higher strength on the bonds located close to the boundaries.

Please find attached a script where the procedure is implemented (lines 73 to 85).

Cheers

Luc

His script:

# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from yade import ymport, plot

#### Simulation of a uniaxial test (compression or tension depending on the sign of the loading rate "rate" defined below)
# everything is documented at https://yade-dem.org/doc/index-toctree.html -> use the quick serach bar on the left of the page

################# SIMULATIONS DEFINED HERE

#### Packing (previously built)
PACKING='core_112_04_020.spheres'

#### Simulation Control
rate=0.003 # deformation rate (0.003 for tension, -0.03 for compression)
iterMax=10000 # maximum number of iterations
saveVTK=2000 # saving output files for paraview
OUT='tensionTest_r0.003' # output files

#### Material microproperties for JCFPM (see https://yade-dem.org/doc/yade.wrapper.html?highlight=jcfpm#yade.wrapper.JCFpmMat)
intR=1.2 # allows near neighbour interaction
DENS=2500 # could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM)
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=20e9,poisson=0.1,frictionAngle=radians(7),tensileStrength=1e6,cohesion=1e6)

#### Loading of the packing into the simulation (O: Omega: the simulation)
O.bodies.append(ymport.text(PACKING,material=sphereMat))

#### Set up boundary conditions (see https://yade-dem.org/doc/yade.utils.html?highlight=uniaxialtestfeatures#yade.utils.uniaxialTestFeatures)
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# ENGINES DEFINED HERE

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 UniaxialStrainer(strainRate=rate,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
         'eps':strainer.strain,
         'sigma':strainer.avgStress,
         'tc':interactionLaw.nbTensCracks,
         'sc':interactionLaw.nbShearCracks,
         'te':interactionLaw.totalTensCracksE,
         'se':interactionLaw.totalShearCracksE,
         'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)

# if you want to plot during simulation
plot.plots={'i':('sigma')}
plot.plot()

################# PREPROCESSING

#### 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.
Saabb.aabbEnlargeFactor=-1.

# reinforce bonds at boundaries for tensile tests
dim=utils.aabbExtrema()
if rate>0:
 layerSize=0.15
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   if ( o.state.pos[longerAxis]<(dim[0][longerAxis]+layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) or ( o.state.pos[longerAxis]>(dim[1][longerAxis]-layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) :
    o.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.FnMax*=100
    i.phys.FsMax*=100

################# SIMULATION REALLY STARTS HERE
strainer.dead=0
O.run(iterMax)