Dear All,
I have a problem in simulating clumps spherical particles' gravity deposition in 2D. I use:
b.state.blockedDOFs='zXY'
b.state.angMom[2]=0
b.state.angVel[2]=0
b.state.vel[2]=0
to fix all the particles in Z plane. But I found that, particle moves along Z axis and when check it in simulation inspection window (in Bodies), I see blockedDOFs is without any value and linear velocity for Z has value.
ubuntu : 18.04
yade 2018.02b
this is minimum of my code:
from yade import export,polyhedra_utils
import os
from yade import plot
import math
from yade import utils
import pylab
import matplotlib; matplotlib.rc('axes',grid=True)
from matplotlib import pyplot
from yade import qt
import numpy as np
from numpy import *
from yade import export as expt
Dolomite = FrictMat()
Dolomite.density = 2870 #kg/m^3
Dolomite.young = 24.36e9 #Pa
Dolomite.poisson = 0.2
Dolomite.frictionAngle = radians(55.12) #rad
Bound = FrictMat()
Bound.density = 2870 #kg/m^3
Bound.young = 60e10 #Pa
Bound.poisson = 0.2
Bound.frictionAngle = radians(55.12) #rad
v1=((-6.4993,-0.01,2),(-2,-0.01,2),(-6.4993,-0.01,-2))
v2=((-2,-0.01,2),(-2,-0.01,-2),(-6.4993,-0.01,-2))
v3=((-1.99,-0.01,2),(-1.99,-0.01,-2),(-1,-0.01,2))
v4=((-1,-0.01,2),(-1,-0.01,-2),(-1.99,-0.01,-2))
v5=((-0.99,-0.01,2),(-0.99,-0.01,-2),(1,-0.01,2))
v6=((1,-0.01,2),(1,-0.01,-2),(-0.99,-0.01,-2))
v7=((1.01,-0.01,2),(1.01,-0.01,-2),(2,-0.01,2))
v8=((2,-0.01,2),(2,-0.01,-2),(1.01,-0.01,-2))
v9=((2.01,-0.01,2),(2.01,-0.01,-2),(3,-0.01,2))
v10=((3,-0.01,2),(3,-0.01,-2),(2.01,-0.01,-2))
v11=((3.01,-0.01,2),(3.01,-0.01,-2),(4.58,-0.01,2))
v12=((4.58,-0.01,2),(4.58,-0.01,-2),(3.01,-0.01,-2))
v13=((4.7,0.0059,2),(5.3,0.93,2),(4.7,0.0059,-2))
v14=((5.3,0.93,2),(5.3,0.93,-2),(4.7,0.0059,-2))
v15=((5.3,0.94,2),(5.3,0.94,-2),(6.47,2.57,2))
v16=((6.47,2.57,2),(6.47,2.57,-2),(5.3,0.94,-2))
v17=((6.47,2.58,2),(6.47,2.58,-2),(6.85,3.01,2))
v18=((6.85,3.01,2),(6.85,3.01,-2),(6.47,2.58,-2))
O.bodies.append(utils.facet(v1,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v2,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v3,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v4,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v5,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v6,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v7,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v8,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v9,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v10,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v11,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v12,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v13,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v14,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v15,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v16,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v17,color=(1,0,1), material=Bound))
O.bodies.append(utils.facet(v18,color=(1,0,1), material=Bound))
rad = 0.2
O.materials.append(FrictMat(young=24.36e9,density=2870))
IDs=O.bodies.append([
utils.sphere(center=(rad,2,3.5*rad),radius=rad, material=Dolomite),
utils.sphere(center=(rad,10.75*rad,3.5*rad),radius=rad, material=Dolomite),
utils.sphere(center=(1.75*rad,2,3.5*rad),radius=rad, material=Dolomite),
utils.sphere(center=(1.75*rad,10.75*rad,3.5*rad),radius=rad, material=Dolomite),
utils.sphere(center=(rad,2,3*rad),radius=rad, material=Dolomite),
utils.sphere(center=(rad,10.75*rad,3*rad),radius=rad, material=Dolomite),
utils.sphere(center=(1.75*rad,2,3*rad),radius=rad, material=Dolomite),
utils.sphere(center=(1.75*rad,10.75*rad,3*rad),radius=rad, material=Dolomite),
])
O.bodies.clump(IDs)
# block z translation and block x, y rotations
for b in O.bodies:
if b.material is Dolomite:
b.state.blockedDOFs='zXY'
b.state.angMom[2]=0
b.state.angVel[2]=0
b.state.vel[2]=0
def checkUnbalanced():
print "iter %d , unbalanced forces %f"%(O.iter, utils.unbalancedForce()) # %[(keyname)][flags][width][.precision]typecode : String Formatting
iter00=O.iter
Unbalanced=open("Unbalanced iter Unbalanced forces.txt","a")
Unbalanced.write(repr(iter00)+' '+repr(utils.unbalancedForce())+' '+"\n")
Unbalanced.close()
def angVel():
for c in O.bodies:
if c.state.angVel!=(0,0,0) or c.state.angMom!=(0,0,0):
s1=c.id
s2=O.iter
s3=O.realtime
s4=c.state.angVel
s5=c.state.angMom
angVel=open("particle iter realtime angVel angMom.txt","a")
angVel.write(repr(s1)+' '+repr(s2)+' '+repr(s3)+' '+repr(s4)+' '+repr(s5)+' '+"\n")
angVel.close()
# Engines:
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],verletDist=0.001),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()], # collision "physics"
[Law2_ScGeom_FrictPhys_CundallStrack()] # contact law -- apply forces
),
NewtonIntegrator(gravity=(0,-9.81,0),damping=0.2),
]
utils.calm()
O.engines=O.engines+[PyRunner(command='checkUnbalanced()',iterPeriod=100000,label="checker")] # call our function defined above 100000 cycles
O.engines=O.engines+[PyRunner(command='angVel()',iterPeriod=100,label="checker2")] # call our function defined above 100 cycles
O.dt=10e-7
O.run(100,True)
O.save('model1.bz2')
O.run()
def ZSpeed():
for b in O.bodies:
if b.material is Dolomite:
b.state.blockedDOFs='zXY'
b.state.angMom[2]=0
b.state.angVel[2]=0
b.state.vel[2]=0
O.engines=O.engines+[PyRunner(virtPeriod=0.02,command='ZSpeed()',label="ZSpeed")]
Why particle along Z axis move? how do I fix this problem?
Many thanks for your help.
Best regards,
