Normal stiffness of contact under tension

Hi guys,

I am a new user of Yade and I am trying to simulate a cohesive material consists of sand particles that are connected by some elastic material. I choose the CohfrictMat as material and Ip2_CohFrictMat_CohFrictMat_CohFrictPhys for bonding. The question is when I define the normal stiffness of bonding (or contact), it is calculated by [1] harmonic average like: 1/(E1.D1)+1/(E2.D2). It is reasonable when two particles are compacting to each other. But when the two particles are under tension stress and leave each other, the stiffness of bonding should have nothing to do with the modulus of particle self and it should only depend on the bonding material itself. Also the shear stiffness should be different. As we know[2] the tensile crack is one of the main mechanism of compression. So how to define the stiffness of tensile?

###############################################here is my code########################

from __future__ import division
from __future__ import print_function

from future import standard_library
standard_library.install_aliases()
from yade import plot,pack,timing
import time, sys, os, copy
readParamsFromTable(noTableOk=True, # unknownOk=True,
 young=1e6,
 poisson=.6,

 sigmaT=3.5e6,
 frictionAngle=atan(0.8),
 epsCrackOnset=1e-4,
 relDuctility=30,

 intRadius=1.5,
 dtSafety=.8,
 damping=0.4,
 strainRateTension=-1,
 strainRateCompression=-.5,
 setSpeeds=True,

 specimenLength=.04,
 sphereRadius=0.001,

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

from yade.params.table import *
#alphaKr=30,alphaKtw=30,
O.materials.append(CohFrictMat(young=young,poisson=poisson,density=2600,frictionAngle=radians(30),normalCohesion=1e3,shearCohesion=1e3,momentRotationLaw=True,etaRoll=0.1,label='spheres'))
#concreteId=O.materials.append(FrictMat(young=15e6,poisson=.4,frictionAngle=0,density=0,label='frictionlessWalls'))

sp=pack.randomDensePack(pack.inAlignedBox((0,0,0),(specimenLength,specimenLength,specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation()
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=0.5*PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(label='is2aabb'),],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True)],
  #setCohesionNow = True,setCohesionOnNewContacts = True
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True),
  Law2_ScGeom_FrictPhys_CundallStrack()],
 ),
 NewtonIntegrator(damping=damping,label='damper'),
 #CpmStateUpdater(realPeriod=.5),
UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
PyRunner(virtPeriod=1e-3/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
 #PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
#
]
plot.plots={'eps':('sigma',)}
def addPlotData():
 yade.plot.addData({'t':O.time,'i':O.iter,'eps':-(strainer.strain),'sigma':-(strainer.avgStress+isoPrestress),})

plot.plot(subPlots=False)

O.timingEnabled=False

O.run()
################

Best regards,
Aaron

[1] https://www.yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Ip2_FrictMat_FrictMat_FrictPhys
[2] Potyondy, D. O., & Cundall, P. A. (2004). A bonded-particle model for rock. International journal of rock mechanics and mining sciences, 41(8), 1329-1364.