Polyhedra gravity deposit simulation explodes with no damping

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()

Thanks Karol Brzezinski, that solved my question.

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?

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

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?

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