yade-users team mailing list archive

[weijie <question692957@xxxxxxxxxxxxxxxxxxxxx>]

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

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()

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