Stable integration of rotations for non-spherical particles in the Discrete Element Method

Efficient integration of motion is crucial for simulating a wide range of physical systems. Despite better algorithms for integrating the motion of non-spherical particles being available for over 25 years, the current State-of-the-art discrete element method (DEM) codes still rely on inaccurate rotation integration algorithms. This issue is particularly noticeable with the increasing popularity of simulations featuring non-spherical objects. We aim to address this issue, highlight the advantages and limitations of existing algorithms, and propose solutions. We have developed a new third-order algorithm that does not require quaternion normalization for each timestep and works for leapfrog and non-leapfrog schemes. The algorithm provides significant improvements over existing methods with only minor increases in computational cost. Our work includes the implementation and a comparison with those currently used in various DEM codes. Our results show that this approach improves accuracy, stability, and overall simulation performance over existing methods. We believe this algorithm can potentially become the new standard in the field. Furthermore, the outcomes of this algorithm are not limited to the discrete element method but can also be valuable for other particle-based techniques such as molecular dynamics (MD).