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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
