yade-users team mailing list archive

[Luc Scholtès <question676639@xxxxxxxxxxxxxxxxxxxxx>]

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.