[EXTERNAL] Re: Problems with Hertz Mindlin and momentum conservation


I tried the things I mentioned in the previous email, and they did not fix the issue you’re seeing. What does fix it is applying the torque slightly differently to each particle. The current code applies torque on each particle according to r_contact X F_t, where F_t is the tangential force, r_contact is the ‘contact vector’, i.e. the vector to the point of contact. Currently, the point of contact is taken to be at the surface of each (undeformed) particle. Because particles are slightly overlapping, this means the torque is applied in slightly different points to each particle. If I correct this, and instead apply the torque at the center of the overlap region (so same location for each particle), the problem you’ve pointed out disappears. Great catch!

This is an easy fix and may be something we’ll want to correct for all the granular pair styles. I just want to test that it doesn’t do anything else strange, I don’t imagine it would. As a side note, I don’t think this would have any noticeable effect on any sort of bulk behavior (especially once you add tangential damping).

Thanks again!


Thank you Dan for your time and attention!

Driven by your intuition I’ve searched in the source code and with an easy calculation I’ve seen That (in a simple planar impact) the loss of angular momentum in a time step dt is F_tdtcsi where F_t and csi are the instantaneous tangential force and overlap.

If i had well understood the source code It would be sufficient a substitution of both “radi” and “radj” with “r/2” in the following lines to fix the problem:

        torque[i][0] -= radi*tor1;
        torque[i][1] -= radi*tor2;
        torque[i][2] -= radi*tor3;

        if (newton_pair || j < nlocal) {
          f[j][0] -= fx;
          f[j][1] -= fy;
          f[j][2] -= fz;
          torque[j][0] -= radj*tor1;
          torque[j][1] -= radj*tor2;
          torque[j][2] -= radj*tor3;

With “r” defined as the distance between the centers of the particles.
Let me know if you agree with me…

I don’t Know how much time you will need to fix all the granular pair styles so do you think that it would be a good idea for me to try to change the source code on my machine by my own in the meanwhile?

Regarding the possibility of this error to affect some bulk behaviour of a large system of grains it is interesting the case that I’m studying for my PhD where this tiny numerical additional torque could sum up with a physical one (done by vibrating wall) that could be of same order. Thus, With the current interaction is very difficult to understand what of the granular collective behaviour observed is due to the real torque.

Thank you again and tell me when the new pair style is ready because i have other scripts to test it in many different (and more complex) granular simulations.



I think r/2 would work only if radi = radj. Otherwise something like:

torque[i][0] -= (radi-0.5*delta)*tor1;

torque[i][1] -= (radi-0.5*delta)*tor2;

torque[i][2] -= (radi-0.5*delta)*tor3;

torque[j][0] -= (radj-0.5*delta)*tor1;

torque[j][1] -= (radj-0.5*delta)*tor2;

torque[j][2] -= (radj-0.5*delta)*tor3;

where delta = radi + radj - r;

I’ll make this addition to the new pair style, that code should hopefully make it into the next stable LAMMPS release.

Thanks again for finding this.