why I cannot control the consolidation process when I choose a certain PSD?

Hi everyone,

I am a new user for yade. I am trying to use Yade to simulate a consolidation process and cutting the soil with a plate. When i use a particle distribution by a certain psdSize and psdCumm rather using the way by setting rMean and rRelFuzz. The code will skip the consolidation process. How can I fix it and get a reasonable process of consolidation. The following is the code:

idSand=O.materials.append(FrictMat(young=10e6,poisson=.25,frictionAngle=radians(30),density=2650,label='sand'))
readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
from yade.params.table import *
from yade import pack, plot
O.bodies.append(geom.facetBox((2,2,2),(2,2,2),wallMask=31))
sp=pack.SpherePack()
#sp.makeCloud((0,0,0),(4,4,1),rMean=rMean,rRelFuzz=rRelFuzz,num=8000)
sp.makeCloud((0,0,0),(4,4,4),porosity=0.3,psdSizes=[2e-2,4e-2,6e-2,8e-2,10e-2,12e-2,14e-2,16e-2],psdCumm=[0.01,0.23,0.4,0.6,0.75,0.93,0.98,1],num=8000)
sp.toSimulation(material=idSand)
# psdSizes: sieve sizes for determine particle sizes in test.
# psdCumm: cummulative percentage numbers from 0~100%.

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.85),
   #PyRunner(command='servo()',iterPeriod=1),
   #VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1000),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.5*PWaveTimeStep()

def checkUnbalanced():
   if O.iter<5000: return
   if unbalancedForce()>.1: return
   a=facet([(0,0,0.4),(4,0,0.4),(4,0,0.6)])
   b=facet([(0,0,0.4),(4,0,0.6),(0,0,0.6)])
   O.bodies.append(a)
   global plate
   plate=O.bodies[-1]
   plate.state.vel=(0,1,0)
   O.bodies.append(b)
   global plate
   plate=O.bodies[-1]
   plate.state.vel=(0,1,0)
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
   checker.command='unloadPlate()'

def unloadPlate():
   if abs(O.forces.f(plate.id)[1])>maxLoad:
      plate.state.vel*=-1
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(plate.id)[1])<minLoad:
      plot.saveDataTxt(O.tags['d.id']+'.txt')
      O.pause()

def addPlotData():
   if not isinstance(O.bodies[-1].shape,Facet):
      plot.addData(); return
   #Depth=max([b.state.pos[2]])
   F=-O.forces.f(plate.id)[1]
   plot.addData(F=F,displacement=plate.state.pos[1]-plate.state.refPos[1],unbalanced=unbalancedForce(),stepnumber=O.iter)

#O.cell.trsf=Matrix3.Identity
plot.plots={'displacement':('F',)}
#plot.plots={'displacement':('F',),'stepnumber':('unbalanced',)}
plot.plot()

O.run()
waitIfBatch()