Yade

Use PotentialBlock to create polyhedrons of different sizes and let them fall freely

Asked by weijie on 2020-09-16

Dear all,

I use the PotentialBlocks module to create polyhedrons of different sizes. The sizes of the three polyhedrons are 0.0035, 0.0075, and 0.02. Let these polyhedrons accumulate under the action of gravity, but I found that the polyhedrons will have problems when falling under gravity. Specifically, they do not conform to free fall. Some polyhedrons rotate during the fall. The kinetic energy I generated in the script The figure also illustrates this problem. I think it may be related to the time step, but this problem will arise no matter how the time step is adjusted.

Could you please tell me how to adjust so that many polyhedrons of different sizes can fall freely correctly?

Thanks in advance.
Jie

Here's my script
##################
from yade import polyhedra_utils,pack,plot,utils,export,qt
import numpy as np
import math
import random
import os
import errno
try:
 os.mkdir('./vtk/')
except OSError as exc:
 if exc.errno != errno.EEXIST:
  raise
  pass

# Engines
Kn=1e8
Ks=Kn/10

O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.001, avoidSelfInteractionMask=2),
 InteractionLoop(
  [Ig2_PB_PB_ScGeom(calContactArea=True)],
  [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=Kn, ks_i=Ks, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
  [Law2_SCG_KnKsPBPhys_KnKsPBLaw(allowViscousAttraction=True, traceEnergy=False)]
 ),
 NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
]

#-------------------------------------------
#Material
n = PolyhedraMat(young=7.2e6,poisson=.2,density=2.5e3)
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(30.0),density=2.5e3,label='frictionless'))
#-------------------------------------------
#Dimensions

meanSize = 0.05
wallThickness = 0.5*meanSize
distanceToCentre = 0.05
lengthOfBase =0.25
heightOfBase = 0.6

#-------------------------------------------
#Make Cloud
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),5*heightOfBase,0.5*(lengthOfBase-wallThickness))#原来3.87

sp.makeCloud(mn,mx,psdSizes=[0.02,0.03,0.04],psdCumm=(0,0.6,1),num=-1)

for center,radius in sp:
    if radius<0.015 :
        aa= [-0.14777172323583135, -0.25201095498198656, -0.959931610624574, -0.027735242321152678, -0.6738824603591866, 0.5159189498469127, 0.6363558175107363, 0.9788716835354719, -0.4948792144101887, 0.288732126188569, -0.16073867070584955, 0.4733856703887221]
        bb= [0.962716291198134, -0.5546504574135571, 0.27026769822865254, 0.5121779180611, -0.7190455783175763, 0.4156828115587175, -0.18560473779779554, 0.17691649133242643, 0.21248364580864817, -0.9473529667988755, -0.5712974115474287, -0.32002788489055134]
        cc= [-0.22658521680291194, 0.7930027418994205, -0.07407208798121794, 0.8584314396525701, -0.16987020316168303, -0.7490229885414054, 0.7487337008758547, 0.10252210623594532, -0.8425824965002652, -0.13840561984253785, -0.8048492699250739, -0.8206632439454689]
        dd= [ 0.01046852, 0.00727971, 0.003921, 0.00489801, 0.00656607, 0.00419807, 0.00729218, 0.00571139, 0.00905, 0.00527244, 0.00897574, 0.00696368]
        r=min(dd)/2 #Suggested value
        b=Body()
        b.mask=1
        b.aspherical=True
        color=Vector3(random.random(),random.random(),random.random())
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=[0,0,1], wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
        O.bodies.append(b)
    if radius>0.015:
        aa= [0.3011899152145198, -0.4738405289361289, -0.029664438986950585, -0.32315972490601735, 0.11149576241679122, -0.5941065406372341, -0.06506475333790455, -0.09841007069609463, 0.16099368191543195, 0.6157363079855748, 0.08296443872076406, -0.8425530695803065, 0.6266406273784567, -0.6766736747863946]
        bb= [-0.9121619654730617, -0.5644017335157929, 0.9950903462905665, -0.9274141931828142, -0.7277361711583503, 0.31308496793216756, -0.3589323337903661, 0.13885002723804651, 0.3406605045802146, 0.38826368415512436, -0.47439526687484435, 0.07723626359180513, 0.629874387220669, -0.7136524211976543]
        cc= [-0.27793017777382756, 0.6759628956842952, -0.09442046271285841, -0.18833668384501445, 0.6767338916818655, -0.7409556135336125, 0.9310929908622869, -0.985411654041895, 0.9262998731525696, 0.6856530540984873, -0.8763937657665611, 0.5330467939376249, -0.45888972579708165, -0.18114347785609775]
        dd= [ 0.01006424, 0.00638965, 0.00577747, 0.00372272, 0.00676639 , 0.0088564, 0.00358511, 0.00330544, 0.00560061, 0.00452511, 0.00683327, 0.01007156, 0.00563036 , 0.00749891]

        r=min(dd)/2 #Suggested value
        b=Body()
        b.mask=1
        b.aspherical=True
        color=Vector3(random.random(),random.random(),random.random())
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
        O.bodies.append(b)

countPBs=0
for b in O.bodies:
 if isinstance(b.shape,PotentialBlock):
  countPBs=countPBs+1
print("number of PotentialBlocks = ", countPBs)

color=[0,0.5,1]
r=0.15*wallThickness

bbb=Body()
bbb.mask=3
wire=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], isBoundary=False, color=color, wire=wire, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [2*lengthOfBase,0,0]
O.bodies.append(bbb)

bA=Body()
bA.mask=3
bA.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bA, bA.shape.volume, bA.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bA.state.pos = [0.5*lengthOfBase,0.5*heightOfBase,0]
lid=O.bodies.append(bA)

bB=Body()
bB.mask=3
bB.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[0.5*wallThickness-r,0.5*wallThickness-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*lengthOfBase-r,0.5*lengthOfBase-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bB, bB.shape.volume, bB.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bB.state.pos = [-0.5*lengthOfBase,0.5*heightOfBase,0]
O.bodies.append(bB)

bC=Body()
bC.mask=3
bC.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], isBoundary=False, color=color, wire=True, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bC, bC.shape.volume, bC.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bC.state.pos = [2*lengthOfBase,0.5*heightOfBase,0.5*lengthOfBase]
O.bodies.append(bC)

bD=Body()
bD.mask=3
bD.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*heightOfBase-r,0.5*heightOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bD, bD.shape.volume, bD.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bD.state.pos = [2*lengthOfBase,0.5*heightOfBase,-0.5*lengthOfBase]
O.bodies.append(bD)
O.dt =2e-5
#-------------------------------------------
# Plot results
def myAddPlotData():
    global wallThickness
    global meanSize
    uf=utils.unbalancedForce()
    if isnan(uf):
        uf = 1.0
    KE = utils.kineticEnergy()
    plot.addData(timeStep1=O.iter,timeStep2=O.iter,timeStep3=O.iter,timeStep4=O.iter,time=O.time,unbalancedForce=uf,kineticEn=KE)

from yade import plot
plot.plots={'timeStep1':('unbalancedForce'),'timeStep2':('kineticEn'),}
plot.plot() #Uncomment to view plots
O.engines=O.engines+[PyRunner(iterPeriod=200,command='myAddPlotData()')]

#-------------------------------------------

# exportPotentialBlocks
from yade import export
vtkExporter_PotentialBlocks = export.VTKExporter('vtk/pb')

def vtkExport():
 vtkExporter_PotentialBlocks.exportPotentialBlocks(what=dict(n='b.id'))

O.engines=O.engines+[PyRunner(iterPeriod=5000,command='vtkExport()')]

from yade import qt
v=qt.View()
v.sceneRadius=10.0
v.ortho=True # I activate orthotropic projection, to make visual comparisons easier
O.saveTmp()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-09-21
Last query:
2020-09-21
Last reply:
2020-09-21
Jan Stránský (honzik) said : #1

Hello,
could you please try to create a MWE [1], a minimal example, i.e. with just a few polyhedrons, but still reproducing your problem?
thanks
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

weijie (amandajoe) said : #2

Hi Jan,

I simplified my script and deleted many unnecessary parts, leaving only a few polyhedra. But I found that the more polyhedrons, the more serious the problem.

Best regards
Jie

Below is my script:
########
from yade import polyhedra_utils,pack,plot,utils,export,qt
import numpy as np
import math
import random
#-------------------------------------------
# Engines
Kn=1e8
Ks=Kn/10
O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.005, avoidSelfInteractionMask=2),
 InteractionLoop(
  [Ig2_PB_PB_ScGeom(calContactArea=True)],
  [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=Kn, ks_i=Ks, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
  [Law2_SCG_KnKsPBPhys_KnKsPBLaw(allowViscousAttraction=True, traceEnergy=False)]
 ),
 NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
]
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(40.0),density=2500,label='frictionless'))
#-------------------------------------------
# Dimensions
lengthOfBase = 0.250
heightOfBase = 0.600
#-------------------------------------------
# Make cloud
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*lengthOfBase,0.4*heightOfBase ,-0.5*lengthOfBase),Vector3(0.5*lengthOfBase,0.6*heightOfBase,0.5*lengthOfBase)
sp.makeCloud(mn,mx,psdSizes=[0.0035,0.005,0.01,0.02,0.04],psdCumm=(0,0.16,0.32,0.84,1),num=-1,distributeMass=False)
#Generate polyhedrons of different sizes
for center,radius in sp:
    if radius<0.0025 :
        r=0.01*0.0025
        b=Body()
        b.mask=1
        b.aspherical=True
        color=[0.3,0.66,0.78]
        distanceToCentre=0.00175
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
        O.bodies.append(b)
    if radius>0.0025 and radius<0.005 :
        r=0.01*0.0025
        b=Body()
        b.mask=1
        b.aspherical=True
        color=[0.1,0.2,0.3]
        distanceToCentre=0.00375
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random())
        O.bodies.append(b)
    if radius>0.01 and radius<0.02 :
        aa= [0.08376984881458294, -0.47781815584193316, -0.6272525512350887, 0.619009994702633, 0.03464399632141613, -0.9832172384126057, -0.2783328604495522, 0.764868485410703, -0.07639578575476973, -0.5670674455594352, 0.621538863800698, -0.7670822373612994, 0.766206355713753, -0.579033908027536, 0.14707077495765938]
        bb= [-0.8413821074662718, -0.8774642832791362, -0.34697163535772957, 0.10830905836646608, 0.9770674225401074, -0.13517146281570994, 0.6998398693326321, 0.6393484283187056, -0.7572080749137317, -0.19334270015810612, -0.5555042677443186, 0.6224892308049091, -0.516288957458027, 0.5202427575278886, 0.25146364935432863]
        cc= [0.533908945106932, -0.04178805471957184, -0.6972552769440008, 0.7778790229425239, -0.21009294450245802, 0.12252566151035581, 0.6578411480642227, 0.07880220321753764, -0.6486906930150549, 0.8006579247607992, 0.5523626067206088, 0.15521597423169936, -0.3825879413556488, -0.6277477252799387, -0.956627524278262]
        dd= [0.012044189,0.008907045,0.012620139,0.012629451,0.012161151,0.005857936,0.00904349,0.007525031,0.012808652,0.014389324,0.006000359,0.01150185, 0.009394797,0.0084963,0.008248162]
        r=min(dd)/2 #Suggested value
        b=Body()
        b.mask=1
        b.aspherical=True
        color=[0.5,0.4,0.3]
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random())
        O.bodies.append(b)
#------------------------------------------
# Bottom plate
color=[0,0.5,1]
wallThickness=0.025
r=0.15*wallThickness
bbb=Body()
bbb.mask=3
wire=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], isBoundary=False, color=color, wire=wire, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [1*lengthOfBase,0,0]
O.bodies.append(bbb)
#-------------------------------------------
# Plot results
def myAddPlotData():
    KE = utils.kineticEnergy()
    plot.addData(timeStep1=O.iter,time=O.time,kineticEn=KE)
from yade import plot
plot.plots={'timeStep1':('kineticEn'),}
plot.plot()
O.engines=O.engines+[PyRunner(iterPeriod=20,command='myAddPlotData()')]
#Timestep
O.dt = 1e-5
#-------------------------------------------
from yade import qt
v=qt.View()
v.sceneRadius=3.0
O.saveTmp()

Jan Stránský (honzik) said : #3

> I simplified my script and deleted many unnecessary parts, leaving only a few polyhedra. But I found that the more polyhedrons, the more serious the problem.

Are you able to reproduce the problem with single / really only a few polyhedrons?
If a single particle does not have this problem, I see no reason why many **non-interacting** polyhedrons would have a problem.

> Specifically, they do not conform to free fall.
> Some polyhedrons rotate during the fall

Are you sure your particles has no interaction? this would be the symptoms of initial interactions..

> I think it may be related to the time step

free fall should be independent of time step

cheers
Jan

weijie (amandajoe) said : #4

Hi Jan,

>Are you able to reproduce the problem with single / really only a few polyhedrons?

This is what makes me feel strange. I only used three polyhedrons in the following script, and the size of the polyhedron is the same as #2. I found that when there are only three polyhedrons, it can fall freely normally. If the number of polyhedrons is increased to the number of #2, it will not conform to free fall, which makes me very confused.

>Are you sure your particles has no interaction? this would be the symptoms of initial interactions..

During the initial generation of particles, enough space is left so that the particles do not touch each other.

>free fall should be independent of time step

This is another point that puzzles me. When there are many polyhedrons, when I use a smaller time step, the situation improves a little, for example, the rotation is not so great.

Best regards
Jie

Below is my script:
#############
from yade import polyhedra_utils,pack,plot,utils,export,qt
import numpy as np
import math
import random
#-------------------------------------------
# Engines
Kn=1e8
Ks=Kn/10
O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.005, avoidSelfInteractionMask=2),
 InteractionLoop(
  [Ig2_PB_PB_ScGeom(calContactArea=True)],
  [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=Kn, ks_i=Ks, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
  [Law2_SCG_KnKsPBPhys_KnKsPBLaw(allowViscousAttraction=True, traceEnergy=False)]
 ),
 NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
]
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(40.0),density=2500,label='frictionless'))
#-------------------------------------------
# Dimensions
lengthOfBase = 0.250
heightOfBase = 0.600
#-------------------------------------------
#Generate polyhedrons of different sizes
#polyhedron A
r=0.01*0.0025
b=Body()
b.mask=1
b.aspherical=True
color=[0.3,0.66,0.78]
distanceToCentre=0.00175
b.shape=PotentialBlock(k=0.0, r=r, R=0.0,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=[0,0,0],fixed=False)
b.state.pos =[lengthOfBase,0.5*heightOfBase ,0]
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
O.bodies.append(b)
#polyhedron B
r=0.01*0.0025
b=Body()
b.mask=1
b.aspherical=True
color=[0.1,0.2,0.3]
distanceToCentre=0.00375
b.shape=PotentialBlock(k=0.0, r=r, R=0.0,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=[0,0,0],fixed=False)
b.state.pos =[lengthOfBase+0.02,0.5*heightOfBase ,0]
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random())
O.bodies.append(b)
#polyhedron C
aa= [0.08376984881458294, -0.47781815584193316, -0.6272525512350887, 0.619009994702633, 0.03464399632141613, -0.9832172384126057, -0.2783328604495522, 0.764868485410703, -0.07639578575476973, -0.5670674455594352, 0.621538863800698, -0.7670822373612994, 0.766206355713753, -0.579033908027536, 0.14707077495765938]
bb= [-0.8413821074662718, -0.8774642832791362, -0.34697163535772957, 0.10830905836646608, 0.9770674225401074, -0.13517146281570994, 0.6998398693326321, 0.6393484283187056, -0.7572080749137317, -0.19334270015810612, -0.5555042677443186, 0.6224892308049091, -0.516288957458027, 0.5202427575278886, 0.25146364935432863]
cc= [0.533908945106932, -0.04178805471957184, -0.6972552769440008, 0.7778790229425239, -0.21009294450245802, 0.12252566151035581, 0.6578411480642227, 0.07880220321753764, -0.6486906930150549, 0.8006579247607992, 0.5523626067206088, 0.15521597423169936, -0.3825879413556488, -0.6277477252799387, -0.956627524278262]
dd= [0.012044189,0.008907045,0.012620139,0.012629451,0.012161151,0.005857936,0.00904349,0.007525031,0.012808652,0.014389324,0.006000359,0.01150185, 0.009394797,0.0084963,0.008248162]
r=min(dd)/2 #Suggested value
b=Body()
b.mask=1
b.aspherical=True
color=[0.5,0.4,0.3]
b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=[0,0,0],fixed=False)
b.state.pos =[lengthOfBase+0.05,0.5*heightOfBase ,0]
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random())
O.bodies.append(b)
#------------------------------------------
# Bottom plate
color=[0,0.5,1]
wallThickness=0.025
r=0.15*wallThickness
bbb=Body()
bbb.mask=3
wire=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], isBoundary=False, color=color, wire=wire, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [lengthOfBase,0,0]
O.bodies.append(bbb)
#-------------------------------------------
# Plot results
def myAddPlotData():
    KE = utils.kineticEnergy()
    plot.addData(timeStep1=O.iter,time=O.time,kineticEn=KE)
from yade import plot
plot.plots={'timeStep1':('kineticEn'),}
plot.plot()
O.engines=O.engines+[PyRunner(iterPeriod=20,command='myAddPlotData()')]
#Timestep
O.dt = 1e-5
#-------------------------------------------
from yade import qt
v=qt.View()
v.sceneRadius=3.0
O.saveTmp()

Best Jan Stránský (honzik) said : #5

> I found that when there are only three polyhedrons, it can fall freely normally

then the initial overlaps are most suspicious

> During the initial generation of particles, enough space is left so that the particles do not touch each other.

In your original code, if I add this at the end
###
O.step()
for i in O.interactions:
    pd = i.geom.penetrationDepth
    b1,b2 = [O.bodies[id] for id in (i.id1,i.id2)]
    c1,c2 = [b.shape.__class__.__name__ for b in (b1,b2)]
    print(pd,c1,c2)
###
I got something like (depending on the run and makeCloud mood, sometimes there is no interaction, sometimes there are more interactions)
###
3.365690143044668e-05 PotentialBlock PotentialBlock
0.00021375937087687118 PotentialBlock PotentialBlock
###
meaning that you DO have some non-negligible initial interactions (and therefore the space left is NOT enough), disturbing the perfect free fall.
Under the presence of interactions, time step of course does have some influence.

cheers
Jan

weijie (amandajoe) said : #6

Thanks Jan Stránský, that solved my question.