TinyDEM: Minimal 3D DEM Solver

• TinyDEM is a minimal yet complete DEM solver that simulates dry, soft granular media using a full suite of contact and friction physics.
• It employs a viscoelastic Hertz–Mindlin contact model with sliding, rolling, and twisting friction, integrated via a semi-implicit Euler method and quaternion dynamics.
• The solver's design features a compact codebase of just two C++11 files with OpenMP parallelization, serving as both an educational tool and a foundation for advanced particle simulations.

Searching arXiv for the specified TinyDEM paper and closely related records. tool call: arxiv_search({"query":"(Vetter, 16 Jul 2025) OR TinyDEM discrete element method sliding rolling twisting friction", "max_results": 5, "sort_by": "relevance"}) TinyDEM is a lightweight implementation of a full-fledged discrete element method (DEM) solver in 3D for a polydisperse ensemble of dry, soft, granular spherical particles. It is presented as a deliberately minimal but complete 3D DEM solver whose main contribution is not a novel physical model, but a compact, readable, standalone implementation of a classical soft-sphere DEM code that still includes inelastic normal contact, sliding friction, rolling friction, twisting friction, rotational dynamics, quaternion-based orientation tracking, and a unified particle–particle / particle–mesh collision pipeline. The code is written in simple C++11, parallelized with OpenMP, published under the 3-clause BSD license, and intended both as an entry point into classical DEM simulations and as a foundation for more advanced particle-dynamics models (Vetter, 16 Jul 2025).

1. Minimal scope and design philosophy

TinyDEM is explicitly designed as a compact standalone program that is easy to inspect and easy to adapt. The implementation consists of only two files, tinydem.hpp and tinydem.cpp, with no external dependencies. The code avoids explicit pointer arithmetic, manual memory management, and polymorphism beyond what the C++ standard library provides. The paper presents this combination of features as part of the main scientific contribution: a small, readable, reusable reference implementation of a full DEM solver (Vetter, 16 Jul 2025).

The solver is deliberately minimal in software structure, but not minimal in contact physics. The paper emphasizes that many small DEM codes omit parts of the contact model, whereas TinyDEM includes the full contact physics needed for dry granular media: inelastic normal contact, tangential friction with contact history, rolling friction, twisting friction, torque exchange, rotational dynamics, and a unified treatment of particle–particle and particle–mesh collisions. This suggests that the term “minimal” refers to implementation compactness rather than to a reduced physical model.

The intended role is twofold. First, TinyDEM serves as an entry point for students and researchers who want a compact code that shows how the essential physics and numerics fit together. Second, it is presented as a base for extensions such as cohesion, more complex particle shapes, moving meshes, or periodic boundaries. A plausible implication is that the code is meant to occupy the niche between pedagogical toy solvers and larger production frameworks.

2. Governing equations and particle kinematics

TinyDEM evolves each particle through damped Newton–Euler equations. For particle $i$, the translational and rotational equations of motion are

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$

$$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$

Here,

$$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$

where $m_i$ and $I_i$ are the mass and moment of inertia of a sphere, ${v}_i$ and ${\omega}_i$ are linear and angular velocity, ${g}$ is an external acceleration such as gravity, and $c_i,\gamma_i$ are optional linear and angular damping rates (Vetter, 16 Jul 2025).

The paper notes that these damping terms can be interpreted as Stokes drag in a viscous fluid, giving

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$0

This interpretation is presented as optional rather than intrinsic to the dry-granular model. In the baseline formulation, the damping terms are simply dissipative contributions in the translational and rotational equations.

The particles are spherical and polydisperse, but the inclusion of full rotational dynamics and explicit orientation tracking prepares the code for future generalization beyond spheres. The paper states this explicitly in connection with quaternion-based orientation updates.

3. Contact law, dissipation, and frictional modes

TinyDEM uses a viscoelastic Hertz–Mindlin contact model with Coulomb friction. For two colliding spheres, the overlap and contact normal are defined by

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$1

The normal and tangential forces are computed as

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$2

The normal displacement is $$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$3, while $$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$4 is a history-dependent tangential spring displacement. For Hertzian contact, the stiffness coefficients are

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$5

with

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$6

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$7

Dissipation is introduced through

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$8

$$m_i\frac{d{v}_i}{dt} = \sum_{j\neq i}\left({f}_\mathrm{n}^{ij}+{f}_\mathrm{s}^{ij}\right) + m_i {g} - c_i m_i {v}_i,$$9

The normal restitution coefficient $$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$0 sets the inelasticity of the collision. The paper emphasizes that the tangential damping is not simply copied from the normal one; TinyDEM uses a heuristic expression chosen to keep the tangential restitution behavior sensible over the whole parameter range (Vetter, 16 Jul 2025).

A notable feature is the explicit inclusion of three frictional modes: sliding, rolling, and twisting. Sliding friction satisfies the Coulomb bound

$$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$1

and analogous bounds are used for rolling and twisting:

$$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$2

The exchanged torques are

$$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$3

Rolling and twisting are treated with Kelvin–Voigt forms,

$$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$4

with relative velocities derived from the relative spin:

$$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$5

The paper argues that rolling and twisting friction are often neglected in DEM codes, but matter for realistic granular packings, static equilibria, and phenomena such as heap formation. In TinyDEM, the history variables $$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$6, $$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$7, and $$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$8 are updated every timestep. If the Coulomb limit is not exceeded, the contact remains in the static regime and the spring is advanced by $$I_i\frac{d{\omega}_i}{dt} = \sum_{j\neq i}\left({\tau}_\mathrm{s}^{ij} + {\tau}_\mathrm{r}^{ij} + {\tau}_\mathrm{t}^{ij}\right) - \gamma_i I_i{\omega}_i.$$9. If the limit is exceeded, the force is projected back to the Coulomb threshold and the spring is reset consistently, so the force changes continuously when switching between static and kinetic friction. The paper also notes the simplification that static and dynamic coefficients are taken equal, and for twisting uses $$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$0.

4. Time stepping, semi-implicit Euler integration, and quaternion dynamics

Time integration in TinyDEM is done with the semi-implicit Euler method, chosen for simplicity and because it is symplectic when dissipation is switched off. After accelerations are computed, the updates are

$$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$1

The ordering matters: the position update uses the already-updated velocity. This is one of the explicit numerical details highlighted in the paper (Vetter, 16 Jul 2025).

For orientation, TinyDEM uses unit quaternions to avoid gimbal lock. The paper gives the synchronous SPIRAL update

$$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$2

with

$$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$3

Quaternion multiplication is performed as a Hamilton product,

$$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$4

The paper states that quaternion updates are mainly used for accurate orientation tracking and visualization, while also preparing the code for future extension beyond spheres. In the present solver, this means that rotational state is treated as a first-class quantity even though the particles themselves are spherical.

5. Unified collision handling and rigid mesh geometry

TinyDEM handles particle–particle and particle–mesh contacts through one unified collision routine. For particle–mesh interactions, the mesh can consist of points, line segments, triangles, and rectangles, each with an optional radius $$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$5. This allows edges and vertices to be rounded into spherocylindrical or spherical-like primitives, and permits fairly complex rigid boundaries and obstacles to be built from simple elements (Vetter, 16 Jul 2025).

Collision detection is split into a broad phase and a dense phase. For particles, a linked-cell grid with cell size $$m_i=\frac{4}{3}\pi \rho_i R_i^3, \qquad I_i=\frac{2}{5}m_iR_i^2,$$6 is used, so potential contacts are found only within a local 27-cell neighborhood, giving near-linear scaling in the number of particles. For meshes, each element gets an axis-aligned bounding box, which restricts distance checks to the cells it overlaps. The dense phase then reduces to standard point–point, point–segment, point–triangle, or point–rectangle distance calculations.

The mesh-contact strategy is described as pragmatic rather than perfect. Because contacts are allowed with multiple mesh elements simultaneously, TinyDEM avoids the complexity of continuous collision detection and seam tracking, but it can exhibit “ghost collisions” when a particle moves across shared edges between mesh elements. The paper frames this as an acceptable tradeoff in a minimal code and recommends coarse meshes where possible. Even with this simplification, the particle–mesh routine enables simulations of walls, funnels, hourglasses, marble runs, and other complex rigid environments.

6. Extensibility, licensing, and research role

TinyDEM is released under the permissive 3-clause BSD license and parallelized with OpenMP. The implementation is described as simple C++11 and explicitly devoid of advanced constructs such as manual memory management or polymorphism. This combination is central to the paper’s claim that TinyDEM is both compact enough to understand in one sitting and physically rich enough to serve as a serious starting point for DEM work (Vetter, 16 Jul 2025).

The paper presents the solver as a foundation for more advanced particle-dynamics models rather than as an endpoint. The conclusion explicitly suggests adding cohesion, more complex particle shapes, moving meshes, periodic boundaries, or other extensions. This suggests a modular research role: the existing code already contains the essential mechanics, contact laws, frictional history treatment, rotational updates, and geometry handling needed for a broad class of dry granular simulations, while leaving room for domain-specific elaboration.

A common misconception about minimal codes is that they necessarily sacrifice physical completeness. TinyDEM is presented as a counterexample to that assumption. Its minimality is in source-code size and software complexity; its mechanics remain those of a classical 3D soft-sphere DEM solver with inelastic normal contact, sliding friction, rolling friction, twisting friction, torque exchange, rotational dynamics, quaternion orientation tracking, and a general particle–mesh collision routine. In that sense, TinyDEM occupies a distinctive position as a deliberately small but complete reference implementation of classical granular DEM.