yade-users team mailing list archive

[Mahdeyeh <question695833@xxxxxxxxxxxxxxxxxxxxx>]

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

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.

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,



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