yade-users team mailing list archive

[Xue <question684944@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello, everyone！
  I want to compact the particles under a plate by applying a vertical force on the plate. For this purpose, I limit the degree of freedom of the plate except for the z-axis motion. But in the calculation, the plate still rotates.
  Here's my code：
###################################################
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
import math
from math import *

gravel = PolyhedraMat()
gravel.density = 2600 #kg/m^3 
gravel.young = 5.5E7 #Pa 1E7
gravel.poisson = 0.21 # 20000/1E7
gravel.frictionAngle = 0.65 #rad
O.materials.append(gravel)

steel = PolyhedraMat()
steel.density = 2600 #kg/m^3 
steel.young = 5.5E7 #inital steel was 
steel.poisson = gravel.poisson
steel.frictionAngle = 0.3 #rad
O.materials.append(steel)

sphereballast = FrictMat()
sphereballast.density = 2600
sphereballast.young = 5.5E7
sphereballast.possion = 0.21
sphereballast.frictionAngle = 0.65 #rad
O.materials.append(sphereballast)



global a
a=1

O.bodies.append(geom.facetCylinder(center=(0.0,0.0,0.3),radius=0.15,height=0.6,orientation=Quaternion((0,0,1),0),segmentsNumber=24,wallMask=4,angleRange=(0.0,2*pi),closeGap=False,radiusTopInner=0.0,radiusBottomInner=0.0,material=steel))
O.bodies.append(geom.facetCylinder(center=(0.0,0.0,1.0),radius=0.150,height=0.8,orientation=Quaternion((0,0,1),0),segmentsNumber=12,wallMask=4,angleRange=(0.0,2*pi),closeGap=False,radiusTopInner=0.0,radiusBottomInner=0.0,material=steel))
#orientation=Quaternion((1,0,0),0)##Law2_ScGeom_FrictPhys_CundallStrack()


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

psdSizes = [0.025,0.0355,0.045,0.056,0.065]
psdCumm = [0.0,0.15,0.35,0.7,0.97,1.0]
mn = Vector3(-0.11,-0.11,0)
mx = Vector3(0.11,0.11,3.7)
sp=pack.SpherePack()
sp.makeCloud(minCorner=(-0.16,-0.16,0),maxCorner=(0.16,0.16,0.8),porosity=0.38,psdSizes=[0.025,0.0355,0.045,0.056,0.065],psdCumm=[0.15,0.35,0.7,0.97,1.0])
pred=pack.inCylinder((0,0,0),(0,0,0.7),radius=0.15)
sp2 = pack.filterSpherePack(pred,sp,returnSpherePack=True)
sp2.toSimulation(material=sphereballast)

for b in O.bodies:
	if isinstance(b.shape,Sphere):
		if b.shape.radius < 0.01:
			O.bodies.erase(b.id)



O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Sphere_Aabb()]),
	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,-9.8)),
	PyRunner(command='cycle()',iterPeriod=300,label='step'),
	PyRunner(command='Compact()',iterPeriod=1,label='compact'),
]

def cycle():
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			if b.state.pos[0] < -0.3 or b.state.pos[0] > 0.3:
				O.bodies.erase(b.id)
			elif b.state.pos[1] < -0.3 or b.state.pos[1] > 0.3:
				O.bodies.erase(b.id)
			elif b.state.pos[2] < -0.1 or b.state.pos[2] > 4:
				O.bodies.erase(b.id)
			else:
				pass
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			b.state.vel=(0,0,0)
			b.state.angVel=(0,0,0)
			b.state.angMom=(0,0,0)

def Compact():
	global a
	if a == 1:
		ldpltheight=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
		print  ("pos=%.5f"%(ldpltheight))
		if ldpltheight > 0.5: return
		ldpltheight=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
		ldplt=polyhedra_utils.polyhedra(steel,v=((-0.148,0,ldpltheight),(-0.073*sqrt(3),-0.073,ldpltheight),(-0.073,-0.073*sqrt(3),ldpltheight),(0,-0.148,ldpltheight),(0.073,-0.073*sqrt(3),ldpltheight),(0.073*sqrt(3),-0.073,ldpltheight),(0.148,0,ldpltheight),(0.073*sqrt(3),0.073,ldpltheight),(0.073,0.073*sqrt(3),ldpltheight),(0,0.148,ldpltheight),(-0.073,0.073*sqrt(3),ldpltheight),(-0.073*sqrt(3),0.073,ldpltheight),(-0.148,0,ldpltheight+0.02),(-0.073*sqrt(3),-0.073,ldpltheight+0.02),(-0.073,-0.073*sqrt(3),ldpltheight+0.02),(0,-0.148,ldpltheight+0.02),(0.073,-0.073*sqrt(3),ldpltheight+0.02),(0.073*sqrt(3),-0.073,ldpltheight+0.02),(0.148,0,ldpltheight+0.02),(0.073*sqrt(3),0.073,ldpltheight+0.02),(0.073,0.073*sqrt(3),ldpltheight+0.02),(0,0.148,ldpltheight+0.02),(-0.073,0.073*sqrt(3),ldpltheight+0.02),(-0.073*sqrt(3),0.073,ldpltheight+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,-500))
		a=a+1
	elif a == 3:
		ldpltheight=O.bodies[-1].state.pos[2]-0.01
		print ("pos=%.5f"%(ldpltheight))

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