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.
_______________________________________________
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
