yade-users team mailing list archive

[Luc Scholtès <question689153@xxxxxxxxxxxxxxxxxxxxx>]

Question #689153 on Yade changed:
https://answers.launchpad.net/yade/+question/689153

Luc Scholtès proposed the following answer:
Hi Yang,

What do you want to achieve?

It seems that there is no loading applied to your sample. No loading
means that the forces are all null...

Please try the following script and tell me if it works (uniaxial
compression with the JCFPM).

Luc

###

from yade import pack, plot

################# SIMULATIONS DEFINED HERO.bodies.append(pack.randomDensePack(pred,radius=D/30.,rRelFuzz=0.4,spheresInCell=1000,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=False,color=(0.9,0.8,0.6),material=sphereMat)) # RK: radius will determine the nb of particles!
E

#### packing (previously constructed)
OUT='compressionTest'

#### Simulation Control
rate=-0.05 #deformation rate
iterMax=200000 # maximum number of iterations                       
saveVTK=int(iterMax/5.) # saving output files for paraview                   

#### Microproperties (interparticle parameters)
intR=1.4 # allows near neighbour interaction (can be adjusted for every packing / the bigger -> the more brittle / careful when intR is too large -> bonds can be created "over" particles) -> intR can be calibrated to reach a certain coordination number K (see calculation on line 115)
DENS=3000 # this one can be adjusted for different reasons (porosity of packing vs porosity of material / increase time step (no gravity -> no real effect on the result)
YOUNG=10e9 # this one controls the Young's modulus of the material
ALPHA=0.15 # this one controls the material Poisson's ratio of the material
TENS=3e6 # this one controls the tensile strength UTS of the material
COH=30e6 # this one controls the compressive strength UCS of the material, more precisely, the ratio UCS/UTS (from my experience: COH should be >= to TENS, >= 10*TENS for competent materials like concrete)
FRICT=30 # this one controls the slope of the failure envelope (effect mainly visible on triaxial compression tests)

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,tensileStrength=TENS,cohesion=COH,frictionAngle=radians(FRICT))

#### create the specimen
L=0.10
D=0.05
pred=pack.inCylinder((0,0,0),(0,0,L),D/2.)
#O.bodies.append(pack.regularHexa(pred,radius=D/20.,gap=0.,material=sphereMat)) # regular packings should be avoided as failure is controlled by particle arrangement
O.bodies.append(pack.randomDensePack(pred,radius=D/30.,rRelFuzz=0.4,spheresInCell=1000,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=False,color=(0.9,0.8,0.6),material=sphereMat)) # RK: radius will determine the nb of particles!

R=0
Rmax=0
Rmin=1e6
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
   if o.shape.radius<Rmin:
     Rmin=o.shape.radius
Rmean=R/nbSpheres

print('nbSpheres=',nbSpheres,' | Rmean=',Rmean, ' | Rmax/Rmin=',
Rmax/Rmin)

#### help define boundary conditions (see utils.uniaxialTestFeatures)
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# DEM loop + 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.8,defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=1,initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks','bstresses'],Key=OUT,dead=1,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)

################# PREPROCESSING

#### manage interaction detection factor during the first timestep and then set default interaction range
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

#### remove particles with less than Kmin contacts (avoid "floating particles") -> Kmin can be adjusted
Kmin=3
erase=0
for o in O.bodies:
    if (not o) or (o.id in negIds) or (o.id in posIds) : continue 
    cont=0
    for i in O.interactions.withBody(o.id) :
	if not i.isReal : continue
	if i.phys.isCohesive :
	    cont+=1
    if cont<Kmin :
	erase+=1
	O.bodies.erase(o.id)
	#print erase, ' | id=', o.id, ' | cont=', cont
print("we removed ", erase, " particles with less than ", Kmin, " contacts")

#### coordination number calculation
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if not i.isReal : continue
    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 ("Kmean=", 2.0*numCohesivelinks/nbSpheres) 

vtk.dead=0
O.step()

################# SIMULATION REALLY STARTS HERE
strainer.dead=0
vtk.iterPeriod=saveVTK
O.run(iterMax)

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.