Question #676639 on Yade changed:
https://answers.launchpad.net/yade/+question/676639
Status: Open => Answered
Luc Scholtès proposed the following answer:
Hi,
For your intention, please find below a script that should do what you
want.
Regarding your point about the coordination number, you have to keep in
mind that it is directly related to the sample. For instance, a loose
sample will present a smaller coordination number than a dense sample if
you define the same interaction range in both cases. In the paper you
mentioned, I generated my samples beforehand by compressing
(hydrostatically) clouds of polydisperse particles. I choose the set of
parameters (contact stiffnesses and confining pressure) so that the
compacity of the sample was "high enough" and the interpenetration
between particles "limited". Up to you to choose a generation procedure
that suits your needs. IMO, for rock like materials, dense samples give
more consistent results. You may need to experience different generation
procedures to make your choice. randomDensePack is an option but I never
found it satisfactory.
Hope it helps
Luc
### JCFPM script
# -*- 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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
