Unstable Granular Simulation


I am trying to simulate a granular system in Lammps. The particles are supposed to hit the bed from the above, but I am trying to create a bed first.
As I decrease the diameter from 8 mm to 0.8 mm, the system becomes very unstable even while damping the kinetic energy by "Fix Viscous "; it is helpful to some extent.

I’ll attach the input.

Is there any way to increase the damping ability of the floor too much or is it going to damage the impact properties of the bed which I want to simulate after I have created the bed.

I’d appreciate any comment.



in.txt (1.74 KB)


In the future I suggest laying out your problem more clearly. I don’t know what is going on from your description, and I’m probably not going to go through much effort to analyze such cryptic problems. Nor is anyone else.

Please analyze the strength of all your relevant terms. I’ve calculated several, but until you have a handle on that you won’t understand why simulations like this will give you garbage. For example, look at the steady state overlap between your particles and the wall (simple potential energy calc). Do you think that a hard particle (like in granular models) should deform so much under its own weight relative to it’s diameter?(let alone when other particles pile on top) Now compare that to how far a particle travels in a given time step at its terminal velocity. What do you think that tells you about your system dynamics?

Read some literature about choosing parameters.

Hi Eric,

Thanks for the comment. I’ll definitely check the parameters and analyze the situation. Would you please suggest some literature that you found effective regarding.



If you’re asking if you could define different damping constants for
some particles (floor) vs others, that’s a current limitation of
granular models in LAMMPS. It would be possible to define
those params on an atom type basis via the pair coeff command,
like other pair styles do, so the floor particle could be a different type.
But we haven’t done that. You could look at the LIGGGHTS add-on
to LAMMPS, which I think has that generalization.



I think the problem is the time step.
It's only 1E-4, which is pretty small for .8 mm particles.
I agree that they should try using LIGGGHTS instead of LAMMPS, since it's designed specifically for granular stuff.
And they should read up on Hertzian contact time; I can't recommend a book, but there's enough stuff out there on google.

email: [email protected]

Sorry for the delay Mehdi. There are several that you could read, I started with CS Campbell’s work for understanding some of the scaling arguments. But Silbert’s work on the LAMMPS page even gives how to choose the Hookean time-step. Note that this chute flow is from that paper.

As for other advice, changing to LIGGGHTS isn’t going to solve problems of stability when parameters were simply chosen incorrectly. Its a simple problem of stability related to both the time-step (you travel ~10% of the diameter in one time-step at your terminal velocity), and the fight between potential energies governing deformations (again ~10% deformation under its own weight). As for choosing time steps be wary of advice that doesn’t have the whole reasoning - you need a time-scale to compare with - you cannot simply get it from length-scales without also estimating some kinematic/dynamics scale.

LIGGGHTS does have more capabilities than LAMMPS for sure, especially if you’re trying to do device calculations. But this problem is software agnostic.

Sorry for the long e-mail and ignore my previous estimates; missplacing decimals can screw with things (your stiffness is fine for the static case). So I took a closer look at your script. I don’t know if you are already reading in a configuration before you pour. Look at the velocities that you are pouring at: -3.04. In one time-step this particle travels a distance of -3.040.0001=3.04e-4 or ~40% the of its diameter. Given acceleration by the time your particle hits the wall at it is traveling at v=-3.04m/s±sqrt(2(9.8m/s^2)*(0.55m-(-.5m)))=-23.62m/s. In one time step it travels -0.0023m or nearly 3 diameters.

Your particles never experience damping (since these groups are never redefined) before a collision with either the wall or particles if you have some already in the group “bed” or group “aab”. If it crashes into a bed of particles, its already deep into the bed before it feels any forces from the particles, and probably causes a cascade of unrealistically large overlaps and forces.

A possible problem as well, you define the upper wall below your simulation domain. Is this effectively trying to fling all your particles into the cavity between -1.0 and -.5? Simple experimentation could answer that question.