yade-users team mailing list archive

[Vasileios Angelidakis <question690546@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Vasileios Angelidakis proposed the following answer:
Hi Jie,

I had a look in your script.

There are several reasons causing instability:
- Your particles pass through the boundaries because the particles forming the boundaries are very thin. A value of wallThickness=0.05 worked well for me.
- frictionAngle must be given in radians, not degrees
- Some particles were erased because you had a statement to erase them inside the myAddPlotData() function

Some further advice: 
- Use larger values for "r". If your dodecahedra have an inradius d=0.107, then I would advice using r=0.053 (around half the min(d) value). Please recheck the size of your dodecahedra, because I used d=array(dd)-r in the shape class and I am not sure if you had subtracted "r" from the dd values already.
- calm() resets the velocities of all bodies; no need to put it in a loop

Last, I replaced for you the PotentialBlockVTKRecorder with
export.exportPotentialBlocks, so you don't have to worry about using
minAabb, maxAabb. Also, the "id" parameter is no longer used in the
PotentialBlock shape class.

You can compare the updated script below with your previous script using
meld [1]

All the best,
Vasileios

[1] https://meldmerge.org/


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.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(30.0),density=3000,label='frictionless'))

#-------------------------------------------
# Dimensions
meanSize = 0.05
wallThickness = 0.5*meanSize
lengthOfBase = 0.250
heightOfBase = 0.600

#-------------------------------------------
# Make cloud
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*(lengthOfBase-wallThickness),0.5*meanSize,-0.5*(lengthOfBase-wallThickness)),Vector3(0.5*(lengthOfBase-wallThickness),0.5*heightOfBase,0.5*(lengthOfBase-wallThickness))#
R=sqrt(3.0)*0.009654 #You can find exact formulae for the circumradius of dodecahedra
sp.makeCloud(mn,mx,R,0,-1,False)

#r=0.01*meanSize
aa=[0.5257311121191336, -5.836787854364519e-17, 0.5257311121191336, -5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 0.85065080835204, -0.5257311121191336, -0.85065080835204, 5.836787854364519e-17, 5.836787854364519e-17, -0.8506508083520399]
bb=[0.85065080835204, -0.5257311121191336, -0.85065080835204, 0.5257311121191336, 5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, -0.85065080835204, 5.836787854364519e-17, 0.5257311121191336, -0.5257311121191336, -0.0]
cc=[-5.836787854364519e-17, 0.85065080835204, 5.836787854364519e-17, 0.85065080835204, -0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, 5.836787854364519e-17, 0.5257311121191336, -0.85065080835204, -0.85065080835204, -0.5257311121191336]
dd=[0.01074978698202965, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.010749786982029651, 0.01074978698202965, 0.01074978698202965, 0.010749786982029651, 0.01074978698202965]

r=min(dd)/2 #Suggested value

for s in sp:
	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=np.array(dd)-r, isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
	utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=s[0], fixed=False)
	b.state.pos = s[0] #s[0] stores center
	b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
	O.bodies.append(b)

#-------------------------------------------
# Bottom plate
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]
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)

#-------------------------------------------
# Plot results
def myAddPlotData():
	global wallThickness
	global meanSize
	uf=utils.unbalancedForce()
	if isnan(uf):
		uf = 1.0
	KE = utils.kineticEnergy()

#	for b in O.bodies:
#		if b.state.pos[1] < -5.0*meanSize and len(b.state.blockedDOFs)==0: #i.e. fixed==False
#			O.bodies.erase(b.id)

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=20,command='myAddPlotData()')]
O.engines=O.engines+[PyRunner(iterPeriod=200,command='still()',firstIterRun=4000)] #enable this engine after 4000 iterations
O.dt = 0.2*sqrt(0.3*O.bodies[0].state.mass/Kn)*100

#-------------------------------------------
# Invoke calm()
def still():
	KE = utils.kineticEnergy()
	if KE<0.01:
		calm()

#-------------------------------------------
# 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=500,command='vtkExport()')]


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.