New question #683855 on Yade:
https://answers.launchpad.net/yade/+question/683855
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=1 # 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 (
