TinyDEM: Minimal open granular DEM code with sliding, rolling and twisting friction

TL;DR

TinyDEM tackles the need for an accessible, minimal open-source discrete element method (DEM) implementation in 3D that includes full frictional torque exchange and quaternion-based orientations. It presents a compact, ~600-line C++11 solver with OpenMP, implementing the Hertz–Mindlin visco-elastic contact model with Coulomb friction, and a unified particle–mesh collision routine that enables complex rigid geometries. The work contributes all three torque modes (sliding, rolling, twisting), a robust SPIRAL quaternion integration scheme, and a transparent, dependency-free codebase suitable for education and rapid prototyping, backed by performance benchmarks on multi-core CPUs. This lightweight, permissively licensed distribution provides a practical platform for learning, method development, and extending DEM capabilities toward more complex particle dynamics.

Abstract

This article introduces TinyDEM, a lightweight implementation of a full-fledged discrete element method (DEM) solver in 3D. Newton's damped equations of motion are solved explicitly for translations and rotations of a polydisperse ensemble of dry, soft, granular spherical particles, using quaternions to represent their orientation in space without gimbal lock. Particle collisions are modeled as inelastic and frictional, including full exchange of torque. With a general particle-mesh collision routine, complex rigid geometries can be simulated. TinyDEM is designed to be a compact standalone program written in simple C++11, devoid of explicit pointer arithmetics and advanced concepts such as manual memory management or polymorphism. It is parallelized with OpenMP and published freely under the 3-clause BSD license. TinyDEM can serve as an entry point into classical DEM simulations or as a foundation for more complex models of particle dynamics.

Figures (6)

• Figure 1: Volume and level of detail of the contact model of published open-source soft-sphere DEM codes. Code volume is measured in number of lines as detailed in Table \ref{['tab:1']}, and includes features beyond the DEM part, where applicable. Codes represented by a filled dot have all three modes of torque included in their friction model, gray dots lack the twisting friction, and empty dots lack both rolling and twisting friction.
• Figure 2: The four particle-particle contact modes. Black arrows indicate linear motion, gray arrows indicate rotation, red arrows indicate resulting forces, magenta arrows indicate resulting torques. Directions are exemplary and can also be reversed.
• Figure 3: Particle-particle contact. Fundamental particle properties in red, geometric contact parameters (distances) in blue, normal vectors in green.
• Figure 4: Particle-mesh collisions. Collisions between particles (yellow) and four different mesh primitives (white) are implemented: a point (A), line segment (B), triangle (C), and rectangle (D). With an optional mesh "radius" $R_j\geq0$ (red), each mesh element $j$ can be given an individual thickness with rounded-off edges. The contact point is shown in green, the closes point of approach on the mesh (which may also lie on the dashed edges of the mesh primitive) in white. By combining these elements, complex mesh geometries can be formed (E), allowing for concurrent contacts of a particle with multiple mesh elements, each having their own "radius". Mesh vertices are shown as black dots. The triangular and rectangular mesh elements have $R_j=0$ in E, hence the sharp edges.
• Figure 5: Example simulations. Particle numbers $N$ are given in the lower left corners. All six simulations have accompanying supplementary videos as labeled. A: The default simulation as preconfigured in the supplied code, a polydisperse set of particles falling through an hourglass with additional obstacles at the throat. B: Binary mixture of two streams consisting of 10 million grains. C: Marble run to demonstrate particle rotations and a somewhat more complex mesh geometry. D: A quasi-2D pile of strongly frictional grains between two vertical plates. $\theta$ is the angle of repose. E: Jamming in 2D. Frictional particles fall through three layers of funnels with neck sizes of 10, 8 and 6 particle diameters, from top to bottom. F: Particle size segregation through vibration of the base of a cylindrical container. Seven balls (blue, initially placed on the ground) have a ten-fold diameter than the rest of the beads, but equal material properties. Images show the gradual upward convection of the larger balls after 0, 83 and 175 vibration cycles.