Some particles move out of the boundary

Asked by jamespaul on 2019-07-19

Hi,

I generated some spheres in a cylinder.The following script shows how I generated them.

##################################
#material
sphere1=O.materials.append(FrictMat(young=6e5,poisson=0.3,density=2460,frictionAngle=float(atan(0.3)),label='sphere1'))

#Cylinder
drumradius=0.071
drumheight=0.005

Cylinder=O.bodies.append(geom.facetCylinder(center=(drumheight/2,drumradius*1.5,drumradius*1.5),radius=drumradius,height=drumheight*2,wallMask=4,orientation=Quaternion(Vector3(0,1,0),(pi/2.0)),segmentsNumber=24,color=(1,0,0),material='sphere1'))

#spheres
pred1 = pack.inCylinder((-drumheight/2,drumradius*1.5,drumradius*1.5),drumheight/2,drumradius*1.5,drumradius*1.5),drumradius*0.8)
sphs1 = pack.randomDensePack(pred1,radius=0.0006,material=sphere1,spheresInCell=3000)
O.bodies.append(sphs1)

But after running,some spheres fly out of the boundry.Maybe:
1.young=6e5 is too small
2.My O.dt=1*PWaveTimeStep() is too big

I don't want to increase young and decrease O.dt,because it will increase the simulation time.

How about making the wall of the cylinder a little thicker to keep out the particles inside?

Thanks,

James

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-07-23
Last query:
2019-07-23
Last reply:
2019-07-22
Robert Caulk (rcaulk) said : #1

Hello:

Please refer to [1], if the PWave is used to estimate max time step, it should be factored by 0.3 to ensure stability. Try:

O.dt=0.3*PWaveTimeStep()

>How about making the wall of the cylinder a little thicker to keep out the particles inside?

Does it work?

Cheers,

Robert

[1]https://yade-dem.org/doc/formulation.html#estimation-of-by-wave-propagation-speed

jamespaul (jamespauljames) said : #2

Thanks Robert,

>Does it work?

Is there a function to do that? ^_^

James

Robert Caulk (rcaulk) said : #3

Sorry, I misunderstood.

It is possible that the wall facet stiffness is too low, resulting in spheres passing through it. The stiffness of facets is computed according to [1], where wall radius is simply 2 times the contacting particle radius (IINM) [2]. You can't control the thickness, but you can control the youngs modulus, which has an equal effect for the final stiffness. I noticed you are using the same youngs modulus for the facet and the spheres, try increasing the youngs modulus of the facet by an order of magnitude.

[1]https://yade-dem.org/doc/formulation.html#normal-stiffness
[2]https://answers.launchpad.net/yade/+question/673800

jamespaul (jamespauljames) said : #4

Thanks Robert,

I've changed some parameters.

sphere1=O.materials.append(FrictMat(young=6e6,poisson=0.3,density=2460,frictionAngle=float(atan(0.3)),label='sphere1'))
drum1=O.materials.append(FrictMat(young=6e10,poisson=0.3,density=2460,frictionAngle=float(atan(0.3)),label='drum1'))

O.dt=.1*PWaveTimeStep()

But there are still some particles that can fly out.

So sad...

Jan Stránský (honzik) said : #5

Please provide a MWE [1], otherwise we would just guess..
cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

jamespaul (jamespauljames) said : #6

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from yade import pack,qt
import math

#material
sphere1=O.materials.append(FrictMat(young=6e6,poisson=0.3,density=2460,frictionAngle=float(atan(0.3)),label='sphere1'))
drum1=O.materials.append(FrictMat(young=6e10,poisson=0.3,density=2460,frictionAngle=float(atan(0.3)),label='drum1'))

#drum
drumradius=0.071
drumheight=0.005

Cylinder=O.bodies.append(geom.facetCylinder(center=(drumheight/2,drumradius*1.5,drumradius*1.5),radius=drumradius,height=drumheight*2,wallMask=4,orientation=Quaternion(Vector3(0,1,0),(pi/2.0)),segmentsNumber=24,color=(1,0,0),material='drum1'))

#generate spheres
pred1 = pack.inCylinder((-drumheight/2,drumradius*1.5,drumradius*1.5),(drumheight/2,drumradius*1.5,drumradius*1.5),drumradius*0.8)
sphs1 = pack.randomDensePack(pred1,radius=0.0006,material=sphere1,spheresInCell=3000)
O.bodies.append(sphs1)

#periodic box
O.periodic = True
O.cell.setBox(drumheight,3*drumradius,3*drumradius)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],allowBiggerThanPeriod=True),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
       [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    NewtonIntegrator(gravity=(0,0,-9.81),damping=0.1),
    PyRunner(command='motion()',virtPeriod=.2,nDo=2),
 PyRunner(command='addPlotData()',virtPeriod=.2),
]

O.dt=.1*PWaveTimeStep()

qt.Controller()
O.saveTmp()
O.run()

def motion():
 O.engines += [RotationEngine(ids=Cylinder,rotationAxis=[1,0,0],rotateAroundZero=True,angularVelocity=1.77,zeroPoint=(0,drumradius*1.5,drumradius*1.5))];

#I erase the spheres which are move out of the cylinder and keep counting the total number of particles
def addPlotData():
 num=0
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   num+=1
   if math.pow((b.state.pos[1]-drumradius*1.5),2)+math.pow((b.state.pos[2]-drumradius*1.5),2)>drumradius*drumradius:
    O.bodies.erase(b.id)

 print 'num=',num

Best Jan Stránský (honzik) said : #7

from [2]:

> # (!) Read [1] carefully

have you read [1] carefully? I have little experience, but I have shown you in [2] that the particle "fly out of the boundry" just because it travels more periods and the contact is not detected any more. The same case might be here.

You can try e.g. not to use allowBiggerThanPeriod. In this case, you would need to make the cylinder with the same height as the "axial dimension" of the cell plus divide the cylinder into several segments along axis:
###
from yade import geom
r = 3.0 # cyl radius
h = 3.0 # cyl "height"
n = 4 # cyl segments along axis
#
facets = []
for i in range(n):
 cyl = geom.facetCylinder((r,r,(i+.5)*h/n),r,h/n,segmentsNumber=24,wallMask=4,wire=False)
 facets.extend(cyl)
O.bodies.append(facets)
O.periodic = True
O.cell.setBox(10,10,h)
#
rs = 0.5
s = sphere((1.1*rs,r,1000),rs) # z=1000 -> completely "outside"
s.state.vel = 1e-1*Vector3(1,1,.5)
O.bodies.append(s)
###

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.InsertionSortCollider.allowBiggerThanPeriod
[2] https://answers.launchpad.net/yade/+question/682137

jamespaul (jamespauljames) said : #8

Thanks Jan Stránský, that solved my question.

jamespaul (jamespauljames) said : #9

Jan,it works.You're amazing!

Thank you!

Thanks Robert too!