DEM Particle Interactions with Solid and Static Bodies

DEM particles require user-defined models for describing interactions with solid surfaces and other DEM particles. As illustrated below, these models describe how a (repelling) contact force is correlated to the overlap distance between two contacting particles. For stiff/elastic particles with very small overlap distances, these interactions can be modeled using Hertz contact theory. This theory links the contact force between two DEM particles (or between a DEM particle and a flat solid surface) to the physical properties of each material (e.g., elastic modulus, Poisson’s ratio, friction coefficients, etc). For particles with more complex solid properties or interactions, users can provide custom UDF contact force expressions. In all cases, these contact forces are calculated on a pairwise basis.

A Hertzian Contact Model


Within the Hertz contact model, interactions between spherical particles can be modeled as:

\[f_{ij}=f_{n,ij}+f_{t,ij},\]

where \(f_{n,ij}\) and \(f_{t,ij}\) are the normal and tangential forces, defined by:

\[f_{n,ij}=-k_{n,ij}\delta_{n,ij}-\gamma_{n,ij}\dot\delta_{n,ij},\]

and

\[f_{t,ij}=-k_{t,ij}\delta_{t,ij}-\gamma_{t,ij}\dot\delta_{t,ij},\]

where \(k_{n,ij}\) and \(k_{t,ij}\) are the normal and tangential stiffnesses; \(\gamma_{n,ij}\) and \(\gamma_{t,ij}\) are the normal and tangential damping coefficients; \(\delta_{n,ij}\) and \(\delta_{t,ij}\) are the normal and tangential overlaps; and \(\dot\delta_{n,ij}\) and \(\dot\delta_{t,ij}\) are the time derivatives.

The stiffnesses are related to the physical properties of the particles and the overlap distances, such that:

\[k_{n,ij}=\frac{4}{3} Y^*_{ij} \sqrt{R^*_{ij} \delta_{n,ij} },\]

and

\[k_{t,ij}= 8 G^*_{ij} \sqrt{R^*_{ij} \delta_{n,ij} },\]

with the equivalent Young’s modulus, \(Y^*_{ij}\), is defined by:

\[\frac{1}{Y^*_{ij}}= \frac{1-\nu_i^2}{Y_i}+\frac{1-\nu_j^2}{Y_j},\]

and the equivalent shear modulus, \(G^*_{ij}\), is defined by:

\[\frac{1}{G^*_{ij}}= 2\frac{\left(2+\nu_i\right)\left(1-\nu_i\right)}{Y_i}+2\frac{\left(2+\nu_j\right)\left(1-\nu_j\right)}{Y_j},\]

where \(Y\), \(G\), and \(\nu\) are the Young’s modulus, shear modulus, and Poisson’s ratio of particle \(i\), and the equivalent radius, \(R^*_{ij}\), is defined by:

\[\frac{1}{R^*_{ij}}= \frac{1}{R_i}+\frac{1}{R_j}.\]

The damping parameters are likewise related to the physical properties of the particles and the overlap distances:

\[\gamma_{n,ij}=-2\left(\sqrt{\frac{5}{6}}\right)\left(\frac{\ln{e_r}}{\sqrt{\ln^2{e_r}+\pi^2}}\right)\sqrt{\frac{2}{3} k _{n,ij} m^*_{n,ij} },\]

and

\[\gamma_{t,ij}=-2\left(\sqrt{\frac{5}{6}}\right)\left(\frac{\ln{e_r}}{\sqrt{\ln^2{e_r}+\pi^2}}\right)\sqrt{k_{t,ij} m^*_{n,ij} },\]

where \(m^*_{ij}\) (in addition to making for a terrific company name) is the equivalent mass defined by:

\[\frac{1}{m^*_{ij}}= \frac{1}{m_i}+\frac{1}{m_j},\]

where, \(Y\), \(G\), and \(\nu\) are the Young’s modulus, shear modulus, and Poisson’s ratio of particle \(i\).

This particle-particle contact theory can be generalized to handle particle-surface contact mechanics by assuming the second particle has an infinite radius. Custom UDF contact forces may be arbitrarily complex, to include combined attractive/repulsive interparticle interaction and time-dependent particle contact forces. As with the Hertz contact model, the key output from any contact model is the total contact forces between two interacting solids. The sum of all these pairwise interaction forces will then be superimposed on any other applied particle forces (e.g., drag, gravity) when integrating the equations of motion.

For both Hertz and UDF contact forces, the pairwise interactions are subject to an interaction cutoff distance. Pairwise interactions are not considered for particles separated by a distance greater than the cutoff. Appealing to this cutoff value increases the efficiency of the solver by eliminating unnecessary pairwise distance calculations. For the Hertz contact model, the interaction cutoff distance is automatically defined as two times the radius of the largest particle. For UDF contact forces, the cutoff distance is user specified.

Note

Since the cutoff radius is equal to two times the radius of the largest particle, minimizing the range of DEM particle diameters will maximize calculation efficiency. Monodispersed particle sets should be considered to maximize simulation solution speed.

When particle rotation is tracked, the time-evolution of the angular velocity of particle \(i\), \(\omega_i\), follows Newton’s second law:

\[I_i \frac{d \mathbf{\omega}_i}{dt}=\displaystyle \sum_j \left(\mathbf{M}_{t,ij}+\mathbf{M}_{r,ij} \right),\]

where \(j\) are the particles interacting with particle \(i\), and the tangential torque is given by:

\[\mathbf{M}_{t,ij}=\mathbf{r}_{i} \times f_{ct,ij},\]

subject to:

\[\left| f_{ct,ij} \right| \leq \mu_{s,ij} \left| f_{cn,ij} \right|,\]

where \(\mathbf{r}_i\) is the distance to the contact point.

The rolling friction torque is given by:

\[\mathbf{M}_{r,ij}=-\mu_{r,ij} \left| f_{cn,ij} \right| \frac{\mathbf{\omega}_{ij}}{\left|\mathbf{\omega}_{ij}\right|}R^*_{ij},\]

where \(\mu_{s,ij}\) and \(\mu_{r,ij}\) are the sliding and rolling friction coefficients between the particle pair.

Return to Create: Particles | DEM Particles