Question #693898 on Yade changed:
https://answers.launchpad.net/yade/+question/693898
Status: Open => Answered
Luc Scholtès proposed the following answer:
Hi David,
"What constitutive law would you recommend to achieve these cohesive
bonds?"
You can use the JCFPM. All interactions (onJoint or notOnJoint) can be
defined with or without cohesion/tensile strength.
The tensile strength defines a maximum admissible force in tension (normal
direction to the contact).
The cohesion defines a maximum admissible force in shear (tangential direction
to the contact).
Please find below an example for simulating a uniaxial compression test
on an "intact" sample (no joints) with JCFPM.
Luc
---
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from yade import pack, plot
################# SIMULATIONS DEFINED HERE
#### packing (previously constructed)
OUT='compressionTest_JCFPM'
#### Simulation Control
rate=-0.01 #deformation rate
iterMax=10000 # maximum number of iterations
saveVTK=2000 # saving output files for paraview
#### Material microproperties
intR=1.1 # allows near neighbour interaction (can be adjusted for every packing)
DENS=2500 # could be adapted to match material density:
dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> packing
porosity as to be computed?
YOUNG=20e9
FRICT=7
ALPHA=0.1
TENS=1e6
COH=1e6
#### material definition
def sphereMat(): return
JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
#### create the specimen
pred=pack.inCylinder((0,0,0),(0,1,0),0.25)
O.bodies.append(pack.regularHexa(pred,radius=0.025,gap=0.,material=sphereMat))
# up to you to define another sample here, e.g., with randomDensePack or
anything else.
R=0
Rmax=0
nbSpheres=0.
for o in O.bodies:
if isinstance(o.shape,Sphere):
nbSpheres+=1
R+=o.shape.radius
if o.shape.radius>Rmax:
Rmax=o.shape.radius
Rmean=R/nbSpheres
print 'nbSpheres=',nbSpheres,' | Rmean=',Rmean
#### boundary condition (see 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.
#### coordination number verification
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
if isinstance(O.bodies[i.id1].shape,Sphere) and
isinstance(O.bodies[i.id2].shape,Sphere):
numSSlinks+=1
if i.phys.isCohesive :
numCohesivelinks+=1
print "nblinks=", numSSlinks, " | nbCohesivelinks=", numCohesivelinks, " ||
Kcohesive=", 2.0*numCohesivelinks/nbSpheres
################# SIMULATION REALLY STARTS HERE
strainer.dead=0
O.run(iterMax)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
