Integration of angular velocity of non-spherical particles is wrong
As explained in "Particle-scale computational approaches to model dry and saturated granular flows of non-Brownian, non-cohesive, and non-spherical rigid bodies" by Wachs (2019), the equation of rotational rigid particle motion, formulated in the body frame, generally reads (Eq. 27):
\frac{d (J \omega)}{d t} = J \frac{d (\omega)}{d t} + \omega \times J \omega = T
The second term featuring the cross product is non-zero if the particle is non-spherical. However, this term is missing in all integrator kernels of PE and mesa-pd (e.g. used in combination with DEM-like collision models).
On the other hand, the HCSITS uses own integration procedures and should also be checked to ensure correctness. There, a difference in the pe-HCSITS (pe/cr/HCSITS.impl.h, line 308):
bodyCache.w_[j] = body->getAngularVel() + dt * ( body->getInvInertia() * ( ( body->getInertia() * body->getAngularVel() ) % body->getAngularVel() ) );
and the mesapd-HCSITS (mesa_pd/kernel/InitParticlesForHCSITS.h)
ac.getAngularVelocityRef(j) = ac.getAngularVelocity(j) + dt * ( ac.getInvInertiaBF(j) * ( ( ac.getInertiaBF(j) * ac.getAngularVelocity(j) ) % ac.getAngularVelocity(j) ) );
has been noted. While the pe uses the worldframe (WF) consistently, mesa_pd mixes WF and body frame (BF) which seems wrong.
Possible solution (assuming that pe HCSITS variant is correct)
-
add getInertia
andgetInvInertia
utility functions to mesa_pd (adapt generation files). See pe/rigidbody/RigidBody.h for implementation. -
Use in XWithShape.h integrators (adapt generation files) according to pe formulation. For VelocityVerlet, it could be that more advanced algorithm needs to be use, see e.g. https://yade-dem.org/doc/formulation.html#orientation-aspherical -
Fix InitParticlesForHCSITS.h accordingly.