← Back to team overview

yade-users team mailing list archive

[Question #657943]: load and displacement

 

New question #657943 on Yade:
https://answers.launchpad.net/yade/+question/657943

HI guys,

I have written up a script to simulate compression test for sand particles. But I faced some questions:

1-  I can not plot he sigma result and when I run my model just I can see a blank graph.

2- I really lost my way, because I need to get result based on Load per displacement  (mm)

3- How can I check my specimens size in the script? diameter and height of the sample?


​I would be most grateful if you could inform me your feedback.

Thank you in advance for all your kind attention.

Yours truly,
HeSam

Here is my script:
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from yade import pack, plot

################# SIMULATIONS DEFINED HERE

#### packing 
OUT='compression of WR sand_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
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))

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.