← Back to team overview

yade-users team mailing list archive

Re: [Question #705718]: Error in the calculation of polyhedral particles

 

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

wangkaichao gave more information on the question:
The following is the code for this polyhedral particle model:
********************************************************
from yade import pack,qt,plot,utils,polyhedra_utils,ymport,export,pack,timing
from yade import *
import numpy
from pprint import pprint
import random
from random import uniform
from random import randint
from math import *

matP = PolyhedraMat() #the material of wall
matP.density = 7800 #kg/m^3 
matP.young = 5.5E8 #10 time biger
matP.poisson = 0.21
matP.fionAngle = 0.03 #0.65
O.materials.append(matP)

steel = PolyhedraMat() #the material of plate
steel.density = 2000 #kg/m^3  in order to ignore the gravity of plate
steel.young = 1E9    #inital steel was ,try
steel.poisson = 0.21
steel.frictionAngle = 0.3 #rad
O.materials.append(steel)


global a
a=1
#length unit: m; using stander unit

O.bodies.append(utils.wall((0,0,0),0,sense=1,material=steel))
O.bodies.append(utils.wall((0.4,0,0),0,sense=-1,material=steel))
O.bodies.append(utils.wall((0,0,0),1,sense=1,material=steel))
O.bodies.append(utils.wall((0,0.4,0),1,sense=-1,material=steel))
O.bodies.append(utils.wall((0,0,0),2,sense=1,material=steel))

'''sp=pack.SpherePack()
sp.makeCloud(minCorner=(-0.19,-0.24,0),maxCorner=(0.19,0.24,0.385),porosity=0.3,psdSizes=[0.0025,0.003,0.00395,0.00535,0.00805,0.0105],psdCumm=[0,0.0939,0.3012,0.5693,0.8842,1])
sp.toSimulation(material=spheresand)'''

O.bodies.append(ymport.textPolyhedra('/home/hello/wenjian/packtestlayer4-3data.dat',material=matP))

O.timingEnabled=True

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()],verletDist=.000005),
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(),Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),Ig2_Facet_Polyhedra_PolyhedraGeom(),Ig2_Sphere_Sphere_ScGeom(),Ig2_Sphere_Polyhedra_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_FrictMat_PolyhedraMat_FrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),Law2_ScGeom_FrictPhys_CundallStrack()],   # contact law -- apply forces #Bo1_Cylinder_Aabb()
),
NewtonIntegrator(damping=0.3,gravity=(0,0,-39.8)),
PyRunner(command='cycle()',iterPeriod=100,label='step'),
PyRunner(command='Compact()',iterPeriod=1,label='compact'),
PyRunner(command='save()',iterPeriod=1000,label='compact')
]

def cycle():
for b in O.bodies:
if isinstance(b.shape,Polyhedra):
if b.state.pos[0] < 0 or b.state.pos[0] > 0.4:
O.bodies.erase(b.id)
elif b.state.pos[1] <0 or b.state.pos[1] > 0.4:
O.bodies.erase(b.id)
elif b.state.pos[2] < -0.15 or b.state.pos[2] > 0.85:
O.bodies.erase(b.id)
else:
pass


def Compact():
global a
if a == 1:
l=O.bodies[6].state.pos[2]#!
if l>0.8 : return#!
ldpltheight=max([b.state.pos[2] for b in O.bodies if isinstance(b.shape,Polyhedra)])
print ("pos=%.5f"%(ldpltheight))
height=ldpltheight+0.03
ldplt=polyhedra_utils.polyhedra(steel,v=((0.4,0.4,height),(0.4,0,height),(0,0,height),(0,0.4,height),(0.4,0.4,height+0.02),(0.4,0,height+0.02),(0,0,height+0.02),(0,0.4,height+0.02)),fixed=False, color=(0.75,0.65,0.65))
O.bodies.append(ldplt)
a=a+1
elif a == 2:
Lastnum=O.bodies[-1].id
O.bodies[Lastnum].state.blockedDOFs='xyXYZ'
O.forces.setPermF(Lastnum,(0,0,-5000))
a=a+1
elif a == 3:
Lastnum=O.bodies[-1].id
plateF=O.forces.f(Lastnum)[2]
print ("Load= %.5f"%(plateF))
h=O.bodies[-1].state.pos[2]
print ("h=%.8f"%(h))
print ("force")
if h<0:#!
O.pause()
else:
if utils.unbalancedForce()<0.25:
a=a+1
elif a == 4:
Lastnum=O.bodies[-1].id
O.bodies[Lastnum].state.vel=[0,0,-0.25]
plateF=O.forces.f(Lastnum)[2]
print ("Load= %.5f"%(plateF))
h=O.bodies[-1].state.pos[2]
print ("h=%.8f"%(h))
print ("vel") 
if h>0:return#!
O.pause()


def save():
export.textPolyhedra('/home/hello/wenjian/packtestlayer4-4data.dat')
pass

O.dt=5E-6


yade.timing.stats()

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