Polyhedra gravity deposit simulation explodes with no damping

Asked by Andrew Jennings on 2021-03-11

I have been trying to create a simple polyhedra gravity deposit simulation without success. Everything looks reasonable until the particles contact the wall. After the particles hit the wall, a few of them rebound at an extremely high speed and exit the simulation. I have tried smaller and smaller time steps and the result is always the same. The only way I've found to eliminate this problem is to introduce a very large amount of damping in the NewtonIntegrator. My understanding is that damping should only be used in steady state simulations, and ultimately I want transient simulations so I think damping must be 0. Here is my code. How can I get the particles to stabilize?

from yade import pack,plot,ymport,export,polyhedra_utils

#define materials
matPoly = PolyhedraMat(
 young=1.0e11,
 frictionAngle=radians(25),
 density=2160,
 poisson=0.3,
)
id_matPoly = O.materials.append(matPoly)

# Add box to hold our particles
O.bodies.append(geom.facetBox((.5,.5,.5),(.6,.6,.6),wallMask=31))

# Create cloud of polys
radMean=.05
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,.1),rMean=radMean,rRelFuzz=.5, seed=2)
# Replace spheres by polys
for pos,radius in sp:
    t = polyhedra_utils.polyhedra(matPoly,(radius,0.4*radius,0.7*radius))
    t.state.pos = pos
    O.bodies.append(t)

#export the initial packing
vtkExporter = export.VTKExporter('Mixed.vtk')
vtkExporter.exportPolyhedra(ids='all', what=dict(speed='b.state.vel.norm()'))
vtkExporter.exportFacets()

O.engines=[
    ForceResetter(),
    InsertionSortCollider([
        Bo1_Polyhedra_Aabb(),
        Bo1_Sphere_Aabb(),
        Bo1_Facet_Aabb()
    ],verletDist=.05*radMean),
    InteractionLoop(
        [
            Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),
            Ig2_Facet_Polyhedra_PolyhedraGeom(),
        ],
        [
            Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(), # collision "physics"
        ],
        [
            Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(), # contact law -- apply forces
        ],
    ),
    NewtonIntegrator(damping=0.0,gravity=[0,0,-9.81]),

    PyRunner(command='ExportResults()', label='exportres', iterPeriod=200000),
]

def ExportResults():
    vtkExporter.exportPolyhedra(ids='all', what=dict(speed='b.state.vel.norm()'))

O.dt=0.0000001

O.stopAtIter=10000000
O.run()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Karol Brzezinski
Solved:
2021-03-13
Last query:
2021-03-13
Last reply:
2021-03-12
Andrew Jennings (andy98225) said : #1

Here's a slight modification on the code I posted earlier. I chose a different particle type. All the particles are exactly the same. The all drop from exactly the same height. What is interesting is when the bounce off the floor, they do not react the same. Some of the particles start rolling at a very high speed. They hit other particles and soon all the particles are rolling at a high speed. This cannot be right.

from yade import pack,plot,ymport,export,polyhedra_utils

#define materials
matPoly = PolyhedraMat(
 young=1.0e09,
 frictionAngle=radians(25),
 density=2160,
 poisson=0.3,
)
id_matPoly = O.materials.append(matPoly)

# Add box to hold our particles
O.bodies.append(geom.facetBox((.5,.5,.5),(.6,.6,.6),wallMask=31))

# Create cloud of polys
radMean=.05
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,.1),rMean=radMean, seed=2)
# Replace spheres by polys
for pos,radius in sp:
    #t = polyhedra_utils.polyhedra(matPoly,(radius,radius,radius))
    #t.state.pos = pos
    t = polyhedra_utils.polyhedralBall(radius, 6, matPoly,pos)
    O.bodies.append(t)

#export the initial packing
vtkExporter = export.VTKExporter('Mixed.vtk')
vtkExporter.exportPolyhedra(ids='all', what=dict(speed='b.state.vel.norm()'))
vtkExporter.exportFacets()

O.engines=[
    ForceResetter(),
    InsertionSortCollider([
        Bo1_Polyhedra_Aabb(),
        Bo1_Facet_Aabb()
    ],verletDist=.05*radMean),
    InteractionLoop(
        [
            Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),
            Ig2_Facet_Polyhedra_PolyhedraGeom(),
        ],
        [
            Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(), # collision "physics"
        ],
        [
            Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(), # contact law -- apply forces
        ],
    ),
    NewtonIntegrator(damping=0.0,gravity=[0,0,-9.81]),

    PyRunner(command='ExportResults()', label='exportres', iterPeriod=20000),
]

def ExportResults():
    vtkExporter.exportPolyhedra(ids='all', what=dict(speed='b.state.vel.norm()'))

O.dt=0.000001

O.stopAtIter=1000000
O.run()

Best Karol Brzezinski (kbrzezinski) said : #2

Hi,

what version of Yade do you use? I had the same problem, that I posted here [1]. Due to Ákos Orosz's experience: 'Ubuntu 16.04 and Yadedaily 2018.02b is the last stable version for Eliáš' model'. I have managed to solve this (or rather work around it), by slight modification in the source code and compiling it by myself. Now it works for me pretty well on Ubuntu 20.04. This change is described in the mentioned topic [1].

Best wishes,
Karol

[1] https://answers.launchpad.net/yade/+question/695154
[2] https://gitlab.com/yade-dev/trunk/-/issues/192

Andrew Jennings (andy98225) said : #3

Thanks Karol Brzezinski, that solved my question.

Andrew Jennings (andy98225) said : #4

Hi Karol Brzezinski,

I have a more complicated simulation that mixed spheres and polyhedra. You fix only seems to work for pure polyhedra simulations. When polyhedra start contacting spheres, the popcorn effect reappears. Any advice on where to look in the code to fix this issue when spheres are present?

Karol Brzezinski (kbrzezinski) said : #5

Hi Andrew,

I wouldn't start by searching for another bug in the code. Maybe it is a matter of proper timestep or interaction functor. There is an example of polyhedra-sphere interaction in the repository. I would start by checking whether this example works on my computer and slowly modify it to get my simulation. If the problem reappears it should become clear what is the reason (if not, you will have at least MWE to post another question).

Cheers,
Karol

[1] https://gitlab.com/yade-dev/trunk/blob/master/examples/polyhedra/sphere-interaction.py

Andrew Jennings (andy98225) said : #6

Never mind, I've got a sphere poly mix working. One thing I noticed, the polys bounce much higher than real rocks. Is there a way to set a coefficient of restitution on them?

Karol Brzezinski (kbrzezinski) said : #7

I believe there is not such a thing for polyhedra (but I may be wrong, I am only a user). I wouldn't be so strict about the Newton damping. If you want to model a pendulum, damping should be definitely zero, but if you are dropping some rocks into a box... I would give it a try.

Newton damping is not physical, but the coefficient of restitution isn't either. The difference is that the first one slows down also the particles that are not in contact. Some viscosity or plasticity would have better physical meaning. A skilled developer could probably implement it into the source code. Otherwise, one can try to simply write a function in Python and put it into PyRunner. It would slow down the simulation, but polyhedra simulations are very slow anyway, so maybe the difference wouldn't be so big.

I am working on a similar simulation right now. I use default damping (0.2), and I will probably stick to this. After deposition, you can change the damping and check how much your final simulation depends on it.

Best wishes,
Karol