Egor Demidov

libgran: a C++ template library for DEM simulations

libgran is a C++ template library for Discrete Element Method (DEM) simulations. DEM is a technique for simulation of granular media consisting of rigid spherical particles.

Drawing of DEM simulations


DEM is a technique for simulation of granular media consisting of rigid spherical particles. The resultant force and torque acting on each particle are computed and used with Newton’s second law to compute the motion of particles:

mx¨=fIω¨=τ m\ddot{\bf{x}}=\bf{f}\quad\quad\quad I\ddot{\boldsymbol{\omega}}=\boldsymbol{\tau}

The forces that particles experience arise from friction at inter-particle contacts, bonding between particles, inter-particle attraction, field forces, etc. libgran contains a bonded and a non-bonded contact model, a Van der Waals attraction model and is designed to be easily extensible with custom models.

Included force models

Contact force

The contact model is based on constraining the four degrees of freedom of motion of two particles relative to each other: normal translation, tangential translation, torsion, and rolling. Let us consider two particles i and j. We can begin by defining a unit normal vector:


The relative velocity at the point of contact is:

vij=vjvi+ωj×an+ωi×an\mathbf{v}_{ij}=\mathbf{v}_j-\mathbf{v}_i+\boldsymbol{\omega}_{j}\times a\mathbf{n}+\boldsymbol{\omega}_{i}\times a\mathbf{n}

where aa is particle radius corrected for inter-particle overlap/separation, v\mathbf{v} is particle translational velocity, and ω\boldsymbol{\omega} is particle angular velocity. Relative velocity can be decomposed into normal and residual (tangential) components:

vij,n=(vijn)n\mathbf{v}_{ij,\rm n}=\left(\mathbf{v}_{ij}\cdot \mathbf{n}\right)\mathbf{n} vij,t=vijvij,n\mathbf{v}_{ij,\rm t}=\mathbf{v}_{ij}-\mathbf{v}_{ij,\rm n}

Similarly, relative angular velocity:


can be decomposed into normal (torsional) and residual (rolling) components:

ωij,o=(ωijn)n\boldsymbol{\omega}_{ij,\rm o}=\left(\boldsymbol{\omega}_{ij}\cdot\mathbf{n}\right)\mathbf{n} ωij,r=ωijωij,o\boldsymbol{\omega}_{ij,\rm r}=\boldsymbol{\omega}_{ij}-\boldsymbol{\omega}_{ij,\rm o}

To constrain the four degrees of freedom, we insert four springs, as illustrated in the featured image. The length of the normal spring, δ\delta, can be computed directly at any point in the simulation from positions of the particles, x\mathbf{x}, and their radius, rr:


The remaining three springs have zero length at the time the contact is formed and have their lengths accumulated throughout the duration of the contact. Let a spring vector be ξ\boldsymbol\xi. Then the rate of stretching / contraction of a spring, ξ˙\dot{\boldsymbol\xi}, is given by:

ξ˙t=vij,t\dot{\boldsymbol\xi}_{\rm t}=\mathbf{v}_{ij,\rm t} ξ˙o=rωij,o\dot{\boldsymbol\xi}_{\rm o}=r\boldsymbol{\omega}_{ij,\rm o} ξ˙r=ωij,r×an\dot{\boldsymbol\xi}_{\rm r}=\boldsymbol{\omega}_{ij,\rm r}\times a\mathbf{n}

Then, as long as the contact lasts, spring ξ\boldsymbol{\xi} is incremented at each time step to obtain a new spring, ξ\boldsymbol{\xi}', to be used at the next time step:

ξ=ξ+ξ˙Δt\boldsymbol{\xi}'=\boldsymbol{\xi}+\dot{\boldsymbol{\xi}}\Delta t

Bonded contact force

For a pair of particles that is connected with a rigid bond, we would like to approximate a rigid-body motion. In other words, the common reference frame of particles i and j can rotate and translate, but any translation or rotation of particle i relative to particle j should be restricted. The distance between all points in the pair of particles should be approximately preserved over time. It can be shown that when using the springs defined in the contact force section to restrict the motion of particles i and j, then the union of particles i and j will undergo rigid body motion as the stiffness of inserted springs approaches infinity. In the simulation we need to use a finite stiffness value, but as long as the amplitude of oscillations is much smaller than the length scale of particles in the simulation, the motion will, approximately, be rigid.

To stabilize the system over time and dissipate any vibrational kinetic energy in the bonds, each spring is supplemented by a dashpot element. The force, f\mathbf{f}, exerted on particle i by each spring in the i-j bond is given by:


where kk is stiffness and γ\gamma is the damping coefficient of the respective spring. And forces exerted on particle j are equal in magnitude and opposite in direction. Forces arising from the normal and tangential springs are applied to the particles. The force associated with the tangential spring will also give rise to torques because the tangential force is not collinear with n\mathbf{n}. Forces computed from the torsion and rolling resistance springs are quasi-forces that are not applied to particles i and j, but are only used to compute torques that will be applied to the particles.

Frictional contact force

The model described in Luding 2008 is used to simulate frictional contacts between non-bonded particles. A brief description is provided here and the reader is referred to Luding 2008 for a more detailed description. Luding’s model uses the same four springs described in an earlier section to compute normal and tangential forces, rolling and torsion resistance torques. Instead of directly setting force proportional to spring elongation, a certain degree of slip is allowed between particles in contact. That is done by computing a test force, f0\mathbf{f}_{0}:


and deciding whether static or dynamic friction should be used based on Coulomb’s law of friction:

fC,s=μsfnf_{C,s}=\mu_s f_{\rm n} fC,d=μdfnf_{C,d}=\mu_d f_{\rm n}

where μs\mu_s is the static friction coefficient, μd\mu_d is the dynamic friction coefficient, and fnf_n is the magnitude of the normal force between the two particles. Then a choice is made whether static or dynamic friction should be used based on the following condition:

if f0fC,s use static friction\text{if}\ \lVert\mathbf{f}_0\rVert\leq f_{C,s}\ \text{use static friction} if f0>fC,s use dynamic friction\text{if}\ \lVert\mathbf{f}_0\rVert> f_{C,s}\ \text{use dynamic friction}

In case the contact is determined to be in the state of static friction, the tangential spring ξ\boldsymbol\xi is incremented as described in contact force section and the test force is used to compute torques and, in the case of the tangential force, is also applied to the contacting particles. In case the contact is determined to be in the state of dynamic friction, the spring is not allowed to stretch anymore. Instead, a certain extent of slipping is allowed. The spring to be used at the next iteration, ξ\boldsymbol\xi', is set to

ξ=1k(fC,df0f0+γξ˙)\boldsymbol\xi'=-\frac{1}{k}\left(f_{C,d}\frac{\mathbf{f}_0}{\lVert \mathbf{f}_0\rVert}+\gamma\dot{\boldsymbol\xi}\right)

and the magnitude of the Coulomb’s force, fC,df_{C,d}, is used to compute torques and, if applicable, accelerate the particles in contact.

The model is only enabled when the normal force is repulsive. Once particles are not overlapping, the frictional force is set to zero and accumulated springs are reset.

Van der Waals attraction force

Hamaker 1937 derived the potential energy, UU, due to Van der Waals attraction between two particles of radius rr whose surfaces are separated by distance δ\delta from each other to be:


where AA is the Hamaker constant - a material property. The magnitude of force acting on the particles can be derived by differentiating potential energy, UU, with respect to separation distance, δ\delta, and the direction of the force will coincide with the normal unit vector n\bf n defined earlier:


Since in the limit as δ\delta approaches 00 the magnitude of force (and potential energy) becomes infinite, a saturation distance δ0\delta_0 is introduced. Attractive force does not increase past the saturation distance. The magnitude of δ0\delta_0 varies between 0.4 and 1 nm (Ranade 1987).

Available software

Source code is available on GitHub.


Comprehensive documentation, including tutorials and class reference, is available externally.