movment along z axis whiles blockedDOFs= zXY

Asked by Mahdeyeh on 2021-03-03

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,

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Karol Brzezinski
Solved:
2021-03-04
Last query:
2021-03-04
Last reply:
2021-03-03
Best Karol Brzezinski (kbrzezinski) said : #1

Hi,

I guess that when bodies are clumped you have to constrain the clump, not the individual bodies. Your condition required material to be Dolomite, however, only individual bodies has this property.

Try to change your loop starting in line 80 to:

for b in O.bodies:
    if b.isClump:
        b.state.blockedDOFs='zXY'
        b.state.angMom[2]=0
        b.state.angVel[2]=0
        b.state.vel[2]=0

Cheers,
Karol

Mahdeyeh (mahdiye.sky) said : #2

Thanks Karol Brzezinski, that solved my question.