yade-users team mailing list archive

[weijie <question692957@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

weijie is still having a problem:
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()

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