Discrete Element Method (DEM)

In this guide, we summarize the theory behind DEM. For further details, we refer the reader to the article by Blais et al. [1] and the one by Golshan et al. [2]

\[\begin{split}m_i\frac{d\mathbf{v_i}}{dt} &= \sum_{j\in \mathcal C_i} (\mathbf{F}_{ij}^\mathrm{n} + \mathbf{F}_{ij}^\mathrm{t}) + m_i\mathbf{g} + \mathbf{F}_i^\mathrm{ext} \\ I_i\frac{d\mathbf{\omega_i}}{dt} &= \sum_{j\in \mathcal C_i} (\mathbf{M}_{ij}^\mathrm{t} + \mathbf{M}_{ij}^\mathrm{r}) + \mathbf{M}_i^\mathrm{ext}\end{split}\]

Where:

• \(m_i\) mass of the particule i;

• \(\mathbf{v_i}\) velocity of the particule i;

• \(\mathcal C_i\) particles in contact list;

• \(\mathbf{F}_{ij}^\mathrm{n}\) normal contact force due to the contact between particles i and j;

• \(\mathbf{F}_{ij}^\mathrm{t}\) tangential contact force due to the contact between particles i and j;

• \(\mathbf{g}\) gravitationnal acceleration;

• \(\mathbf{F}_i^\mathrm{ext}\) external forces;

• \(I_i\) moment of inertia of particle i;

• \(\mathbf{\omega_i}\) angular velocity of particle i;

• \(\mathbf{M}_{ij}^\mathrm{t}\) tangential friction torque due to the contact between particles i and j;

• \(\mathbf{M}_{ij}^\mathrm{r}\) rolling friction torque due to the contact between particles i and j;

• \(\mathbf{M}_i^\mathrm{ext}\) external torques;

Contact Force and Torque Models

The normal and tangential contact forces use linear or nonlinear viscoelastic models and are calculated as followed:

\[\begin{split}\mathbf{F}_{ij}^\mathrm{n} &= -(k_\mathrm{n}\delta_{\mathrm{n}})\mathbf{n}_{ij}-(\eta_\mathrm{n}\mathbf{v}_{rn}) \\ \mathbf{F}_{ij}^\mathrm{t} &= -(k_\mathrm{t}\mathbf{\delta}_\mathrm{t})-(\eta_\mathrm{t}\mathbf{v}_{rt})\end{split}\]

Where:

• \(\delta_{\mathrm{n}}\) normal overlap;

• \(\mathbf{\delta_\mathrm{t}}\) tangential overlap vector;

• \(\mathbf{n}_{ij}\) contact normal vector;

• \(k_\mathrm{n}, k_\mathrm{t}\) spring constants;

• \(\eta_\mathrm{n}, \eta_\mathrm{t}\) damping model constants;

• \(\mathbf{v}_{rn}\) relative velocity in the normal direction;

• \(\mathbf{v}_{tn}\) relative velocity in the tangential direction;

Representation of a typical particle-particle contact. [2]

The contact normal vector \(\mathbf{n}_{ij}\) is computed as:

\[\mathbf{n}_{ij}=\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{\left|\mathbf{x}_{j}-\mathbf{x}_{i}\right|}\]

The normal overlap (\(\delta_{\mathrm{n}}\)) is the contact distance between the particles i and j. In the case of a collision between a particle and a wall, the wall is considered as j. The tangential overlap (\(\delta_\mathrm{t}\)) depends on the contact history and is updated during a contact. The normal and tangential overlaps are calculated as follow:

\[\begin{split}\delta_{\mathrm{n}} =& \:R_i + R_j - |\mathbf{x}_{j} - \mathbf{x}_{i}| \\ \mathbf{\delta}_{ij}^{\mathrm{t,new}} &= \mathbf{\delta}_{ij}^{\mathrm{t,old}}+\mathbf{v}_{ij,\mathrm{t}}dt\end{split}\]

Relative Velocities

The relative velocities are calculated to update the tangential overlap and for their contribution in the force models:

\[\begin{split}\mathbf{v}_{ij} &= \mathbf{v}_i-\mathbf{v}_j+\left(R_i\mathbf{\omega}_i+R_j\mathbf{\omega}_j\right)\times\mathbf{n}_{ij} \\ \mathbf{v}_{ij,\mathrm{n}} &= \left(\mathbf{v}_{ij}.\mathbf{n}_{ij}\right)\mathbf{n}_{ij} \\ \mathbf{v}_{ij,\mathrm{t}} &= \mathbf{v}_{ij}-\mathbf{v}_{ij,\mathrm{n}}\end{split}\]

Spring and damping constants

The spring and damping constants for the linear and nonlinear viscoelastic models are calculated as follows:

Spring and damping Models used in Lethe.

Parameters	Linear model definitions	Nonlinear viscoelastic model definitions [3]
Normal spring constant	\(k_\mathrm{n} = \frac{16}{15}\sqrt{R_{\mathrm{e}}}Y_{e}\left(\frac{15m_{e}V^2}{16\sqrt{R_{\mathrm{e}}}Y_{e}}\right)^{0.2}\)	\(k_\mathrm{n} = \frac{4}{3}Y_{e}\sqrt{R_{\mathrm{e}}\delta_{\mathrm{n}}}\)
Normal damping model constant	\(\eta_\mathrm{n} = -2\beta\sqrt{m_{e} k_\mathrm{n}}\)	\(\eta_\mathrm{n} = -2\sqrt{\frac{5}{6}}\beta\sqrt{S_\mathrm{n}m_{e}}\)
Tangential spring constant	\(k_\mathrm{t} = 0.4 k_\mathrm{n}\)	\(k_\mathrm{t} = 8G_{e}\sqrt{R_{\mathrm{e}}\delta_{\mathrm{n}}}\)
Tangential damping model constant	\(\eta_\mathrm{t} = -2\beta\sqrt{m_{e} k_\mathrm{t}}\)	\(\eta_\mathrm{t} = -2\sqrt{\frac{5}{6}}\beta\sqrt{S_\mathrm{t}m_{e}}\)

Where:

• \(R_{\mathrm{e}}\) effective radius;

• \(Y_\mathrm{e}\) effective Young’s modulus;

• \(m_\mathrm{e}\) effective mass;

• \(V\) characteristic impact velocity, this parameters is set to 1.0;

• \(e\) coefficient of restitution;

• \(G_\mathrm{e}\) effective shear modulus;

These parameters are computed as follows:

\[\begin{split}\frac{1}{m_\mathrm{e}} &= \frac{1}{m_i}+\frac{1}{m_j} \\ \frac{1}{R_{\mathrm{e}}} &= \frac{1}{R_i}+\frac{1}{R_j} \\ \frac{1}{G_\mathrm{e}} &= \frac{2(2-\nu_i)(1+\nu_i)}{Y_i}+\frac{2(2-\nu_j)(1+\nu_j)}{Y_j} \\ \frac{1}{Y_\mathrm{e}} &= \frac{\left(1-\nu_i^2\right)}{Y_i}+\frac{\left(1-\nu_j^2\right)}{Y_j} \\ \beta &= \frac{\ln{e}}{\sqrt{\ln^2{e}+\pi^2}} \\ S_\mathrm{n} &= 2Y_{e}\sqrt{R_{\mathrm{e}}\delta_{\mathrm{n}}} \\ S_\mathrm{t} &= 8G_{e}\sqrt{R_{\mathrm{e}}\delta_{\mathrm{n}}}\end{split}\]

Where:

• \(\nu_i, \nu_j\) poisson coefficient of particle i or j;

Coulomb’s limit

Coulomb’s criterion is breached when the following condition is broken during a collision:

\[|\mathbf{F}_{ij}^{\mathrm{t}}| \geq \mu |\mathbf{F}_{ij}^\mathrm{n}|\]

A breach means the collision is having gross sliding and tangential force needs to be limited to the Coulomb’s limit. To do so, the tangential overlap \(\mathbf{\delta_\mathrm{t}}\) is first limited and then the tangential force is recalculated.

When using nonlinear viscoelastic contact model, the tangential overlap is computed from tangential spring force :

\[\begin{split}\mathbf{\delta_\mathrm{t}} &= \frac{\mathbf{\tilde{F}_{ij}}}{-k_\mathrm{t}} \\ \mathbf{\tilde{F}_{ij}} &= \mathbf{\hat{F}_{ij}} + \eta_\mathrm{t}\mathbf{v}_{ij,\mathrm{t}} \\ \mathbf{\hat{F}_{ij}^\mathrm{t}} &= \mu |\mathbf{F}_{ij}^\mathrm{n}| \frac{\mathbf{F}_{ij}^\mathrm{t}}{|\mathbf{F}_{ij}^\mathrm{t}|}\end{split}\]

Regarding the particle-wall contacts, the applied models are the same as for particle-particle contacts.

Note

When using a cohesive force model, Coulomb’s criterion needs to be modified. For further information on cohesive force models, see Cohesive force models .

Tangential torque

Tangential torque is the torque generated by the tangential force. It can be calculated through:

\[\mathbf{M}_{ij}^\mathrm{t} = R_{i}\mathbf{n}_{ij} \times \mathbf{F}_{ij}^\mathrm{t}\]

Note

As of now, the lethe-particles solver only uses spherical particles, thus the normal force does not generate a torque on the particle during a collision.

Rolling friction models

Rolling friction may be computed through a constant torque model or a viscous torque model. It is also possible to ignore the rolling resistance. The corresponding model can be described by the following equations:

Rolling Friction Models used in Lethe.

Constant resistance	Viscous resistance	No resistance
\(\mathbf{M}_{ij}^\mathrm{r} = -\mu_\mathrm{r}R_{\mathrm{e}}|\mathbf{F}_{ij}^\mathrm{n}| \mathbf{\hat{\omega}}_{ij}\)	\(\mathbf{M}_{ij}^\mathrm{r} = -\mu_\mathrm{r}R_{\mathrm{e}}|\mathbf{F}_{ij}^\mathrm{n}||\mathbf{V}_{\omega}| \mathbf{\hat{\omega}}_{ij}\)	\(\mathbf{M}_{ij}^\mathrm{r} = 0\)

Where:

• \(\mu_\mathrm{r}\) rolling friction coefficient;

• \(\hat{\omega}_{ij}\) relative angular velocity;

• \(V_{\omega}\) contact point relative velocity caused by the angular velocities;

The parameters are computed as follows:

\[\begin{split}\mathbf{\hat{\omega}}_{ij} &= \frac{\omega_{i} - \omega_{j}}{|\omega_{i} - \omega_{j}|} \\ \mathbf{V}_{\omega} &= \left( \omega_{i} \times R_{i}\mathbf{n}_{ij}-\omega_{j} \times R_{j}\mathbf{n}_{ji} \right).\end{split}\]

Cohesive force models

Lethe supports two cohesive force models: the Johnson-Kendall-Roberts (JKR) and the Derjaguin-Muller-Toporov (DMT). Both models describe attractive forces due to van der Waals effects. Choosing the right model can be based on the Tabor parameter \(\mathbf{\tau}\) which represents the ratio between the normal elastic deformation caused by adhesion and the distance at which adhesion forces occur. [4]

This parameter can be described as:

\[\mathbf{\tau} = \left( \frac{R_{\mathrm{e}} \gamma_{\mathrm{e}}^2}{Y_\mathrm{e}^2 z_{\mathrm{o}}^3}\right)^{1/3}\]

Where \(\mathbf{z_{\mathrm{o}}}\) is the equilibrium separation of the surfaces and \(\mathbf{\gamma}_\mathrm{e}\) the effective surface energy. The DMT model is applicable for low \(\mathbf{\tau}\) values (\(\mathbf{\tau} < 1\)) while the JKR model is more appropriate for high \(\mathbf{\tau}\) values (\(\mathbf{\tau} > 1\)) . In essence, the DMT model is preferred for small, hard particles (high \(Y\)) and the JKR model for large, soft particles.

Johnson-Kendall-Roberts force model

The Johnson-Kendall-Roberts (JKR) model describes attractive forces due to van der Waals effects. [5] This model modifies the Hertz formulation by defining a larger contact path radius (\(\mathbf{a}\)) and by taking into account the effective surface energy (\(\mathbf{\gamma}_{e}\)). The model is defined by:

\[a^{3} = \frac{3 R_{\mathrm{e}}}{4 Y_{\mathrm{e}}} \left[|\mathbf{F_{n}^{JKR}}| + 3\pi\gamma_{\mathrm{e}}R_{\mathrm{e}} + \sqrt{6 |\mathbf{F_{n}^{JKR}}| \pi\gamma_{\mathrm{e}}R_{\mathrm{e}} + (3\pi\gamma_{\mathrm{e}}R_{\mathrm{e}})^2 }\right]\]

Where \(\mathbf{F_{n}^\mathrm{JKR}}\) corresponds to the normal spring force and attractive force combined and \(\mathbf{\gamma_{\mathrm{e}}}\) is the effective surface energy. Note that if the effective surface energy is equal to zero, the JKR model reverts to Hertz model.

The effective surface energy can be computed as:

\[\gamma_{\mathrm{e}} = \gamma_{i} + \gamma_{j} - 2\gamma_{i,j}\]

Where \(\gamma_{i}\) and \(\gamma_{j}\) are the surface energy of each material (particle or wall) and where \(\gamma_{i,j}\) is the interface energy which is equal to zero when both surfaces are the same material. The interface energy term is approximated using [#israelachvili–289]_:

\[\gamma_{i,j} \approx \left( \sqrt{\gamma_{i}} - \sqrt{\gamma_{j}} \right)^{2}\]

To compute the \(\mathbf{F_{n}^{JKR}}\), the contact patch radius needs to be determined. The contact patch radius can be related to the normal overlap as follows:

\[\delta_{\mathrm{n}} = \frac{ a^{2} }{ R_{\mathrm{e}} } - \sqrt{ \frac{2 \pi \gamma_{\mathrm{e}} a }{ Y_\mathrm{e} }}\]

This equation can be rewritten as a fourth-order polynomial function with two complex and two real roots.

\[0 = a^{4} - 2R_{\mathrm{e}}\delta_{\mathrm{n}}a^{2} - 2\pi\gamma_{\mathrm{e}}R_{\mathrm{e}}^{2}a + R_{\mathrm{e}}^{2}\delta_{\mathrm{n}}^{2}\]

Since we are always solving for the same real root, a straightforward procedure, described by Parteli et al. can be used [7]:

\[\begin{split}c_\mathrm{0} &= R_{\mathrm{e}}^{2}\delta_{\mathrm{n}}^{2} \\ c_\mathrm{1} &= \frac{-2\pi\gamma_{\mathrm{e}}R_{\mathrm{e}}^{2}}{Y_{\mathrm{e}}}\\ c_\mathrm{2} &= -2R_{\mathrm{e}}\delta_{\mathrm{n}}\\ P &= -\frac{c_{\mathrm{2}}^{2}}{12} - c_{\mathrm{0}} \\ Q &= - \frac{c_{\mathrm{2}}^{3}}{108} + \frac{c_{\mathrm{0}}c_{\mathrm{2}}}{3} - \frac{c_{\mathrm{1}}^{2}}{8} \\ U &= \left[ -\frac{ Q }{ 2 } + \sqrt{ \frac{ Q^{2} } {4} + \frac{ P^{3} }{ 27 } } \right]^{ \frac{1}{3} } \\ s &= \begin{cases} -5c_{\mathrm{2}}/6 + U - \frac{P}{3U} &{if}\: P \neq 0 \\ -5c_{\mathrm{2}}/6 + Q^{\frac{1}{3}} &{if}\: P = 0 \end{cases}\\ \omega &= \sqrt{c_{\mathrm{2}} + 2 s} \\ \lambda &= \frac{c_{\mathrm{1}} }{2 \omega}\\ a &= \frac{1}{2}\left(\omega + \sqrt{\omega^{2} - 4(c_{\mathrm{2}} + s + \lambda ) } \right)\end{split}\]

Finally, the \(\mathbf{F_{\mathrm{n}}^{JKR}}\) can be computed as follows:

\[\mathbf{F_{n}^{JKR}} = \left( \frac{4 Y_\mathrm{e} a^{3}}{3 R_{\mathrm{e}}} - \sqrt{8 \pi \gamma_{\mathrm{e}} Y_{\mathrm{e}} a^{3}} \right) \mathbf{n}_{ij}\]

The normal damping, tangential damping and tangential spring constants need to be computed using the same procedure as the nonlinear model.

A simplified version of the JKR model (SJKR-A) is implemented in Lethe. Please refer to C. J. Coetzee and O. C. Scheffler for more information on the different versions of the JKR model and their specific features [5].

A modified Coulomb’s limit, based on the work of C. Thornton [9], is used for the JKR model. Using the usual limit can result in permanent slip since the total normal force can be equal to zero even when there is a substantial overlap between particles.

The modified Coulomb’s criterion is breached when the following condition is broken during a collision:

\[|\mathbf{F}_{ij}^{t}| \geq \mu |\mathbf{F_{n}^{JKR} + 2F_{\mathrm{po}}}|.\]

Where \(\mathbf{F_{\mathrm{po}}}\) is the pull-off force, which can be computed as follows:

\[\mathbf{F_{\mathrm{po}}} = \left(1.5\pi\gamma_{\mathrm{e}}R_{\mathrm{e}}\right) \mathbf{n}_{ij}\]

Derjaguin-Muller-Toporov force model

The Derjaguin-Muller-Toporov (DMT) model describes attractive forces due to van der Waals effects. This model is more suitable for particles with smaller diameter, lower surface energy and higher Young’s modulus. In Lethe, the DMT model is implemented following the work of Meier et al. [10]. This implementation includes non-contact forces between particles. The model is described by the following equations:

\[\begin{split}\mathbf{F_{ad}^{DMT}} = \begin{cases} F_{\mathrm{po}} = -2\pi\gamma_{\mathrm{e}}R_{\mathrm{e}}, & \delta_\mathrm{n} \leq \delta_{\mathrm{o}} \\ \frac{-AR_{\mathrm{e}}}{6 \delta_{n}^2}, & \delta_{\mathrm{o}} < \delta_\mathrm{n} < \delta^* \\ 0, & \delta^* \leq \delta_{\mathrm{n}} \end{cases}\end{split}\]

where \(A\) is the Hamaker constant which is used to quantify the strength of van der Waals forces. \(\delta_{\mathrm{o}}\) represents the distance at which the van der Waals force curve equals the pull-off force \(F_{\mathrm{po}}\) and \(\delta^*\) represents a cut-off radius at which the van der Waals has a relative decline of \(C_{\mathrm{FPO}}\) [10]. They are computed using:

\[\begin{split}\begin{align} \delta_{\mathrm{o}} &= - \sqrt{\frac{ -A R_{\mathrm{e}}}{6 F_{\mathrm{po}}}}\\ \delta^* &= \frac{\delta_{\mathrm{o}}}{ \sqrt{C_{\mathrm{FPO}}}} \end{align}\end{split}\]

where \(C_{\mathrm{FPO}}\) is a user parameter used to determined the cut-off distance at which the non-contact forces are being performed.

The Coulomb’s limit threshold for the DMT model is computed in the same way as for the non-linear viscoelastic model. This means that the adhesion force term is not taken into account when computing the norm of the normal force. For further information, see Coulomb’s limit .

Integration Methods

Two types of integration methods are implemented in Lethe-DEM:

• Explicit Euler method;

• Velocity Verlet method

Explicit Euler method is calculated as:

\[\begin{split}\mathbf{v}_{i}^{n+1} &= \mathbf{v}_{i}^{n} + \mathbf{a}_{i}^{n}dt \\ \mathbf{x}_{i}^{n+1} &= \mathbf{x}_{i}^{n} + \mathbf{v}_{i}^{n}dt\end{split}\]

And velocity Verlet method is calculated with half-step velocity as:

\[\begin{split}\mathbf{v}_{i}^{n+\frac{1}{2}} &= \mathbf{v}_{i}^{n} + \mathbf{a}_{i}^{n}\frac{dt}{2} \\ \mathbf{x}_{i}^{n+1} &= \mathbf{x}_{i}^{n} + \mathbf{v}_{i}^{n+\frac{1}{2}}dt \\ \mathbf{v}_{i}^{n+1} &= \mathbf{v}_{i}^{n+\frac{1}{2}} + \mathbf{a}_{i}^{n+1}\frac{dt}{2}\end{split}\]

References

[1]

B. Blais, D. Vidal, F. Bertrand, G. S. Patience and J. Chaouki, “Experimental Methods in Chemical Engineering: Discrete Element Method—DEM,” Can. J. Chem. Eng., vol. 97, pp. 1964-1973, 2019, doi: 10.1002/cjce.23501.

[2] (1,2)

S. Golshan, P. Munch, R. Gassmöller, M. Kronbichler, and B. Blais, “Lethe-DEM: an open-source parallel discrete element solver with load balancing,” Comput. Part. Mech., vol. 10, no. 1, pp. 77–96, Feb. 2023, doi: 10.1007/s40571-022-00478-6.

[3]

R. Garg, J. Galvin-Carney, T. Li, and S. Pannala, “Documentation of open-source MFIX–DEM software for gas-solids flows,” Tingwen Li Dr., p. 10, Accessed: Sep. 2012, Available: https://mfix.netl.doe.gov/doc/mfix-archive/mfix_current_documentation/dem_doc_2012-1.pdf.

[4]

D. S. Grierson, E. E. Flater, and R. W. Carpick, “Accounting for the JKR–DMT transition in adhesion and friction measurements with atomic force microscopy,” Journal of Adhesion Science and Technology, vol. 19, no. 3–5, pp. 291–311, Jan. 2005, doi: 10.1163/1568561054352685.

[5] (1,2)

C. J. Coetzee and O. C. Scheffler, “Review: The Calibration of DEM Parameters for the Bulk Modelling of Cohesive Materials,” Processes, vol. 11, no. 1, Art. no. 1, Jan. 2023, doi: 10.3390/pr11010005.

[6]

J. N. Israelachvili, “Chapter 13 - Van der Waals Forces between Particles and Surfaces,” in Intermolecular and Surface Forces, 3rd ed., J. N. Israelachvili, Ed., Boston: Academic Press, 2011, pp. 253–289, doi: 10.1016/B978-0-12-391927-4.10013-1.

[7]

E. J. R. Parteli, J. Schmidt, C. Blümel, K.-E. Wirth, W. Peukert, and T. Pöschel, “Attractive particle interaction forces and packing density of fine glass powders,” Sci Rep, vol. 4, no. 1, Art. no. 1, Sep. 2014, doi: 10.1038/srep06227.

[8]

G. Violano, G. P. Demelio, and L. Afferrante, “On the DMT Adhesion Theory: From the First Studies to the Modern Applications in Rough Contacts.” Procedia Structural Integrity, vol. 12, pp. 58–70, Jan. 2018, doi: 0.1016/j.prostr.2018.11.106.

[9]

C. Thornton, “ Interparticle sliding in the presence of adhesion,” Journal of Physics D: Applied Physics, vol. 24, no. 11, pp. 1942–1946, 1991, doi: 10.1088/0022-3727/24/11/007.

[10] (1,2)

C. Meier, R. Weissbach, J. Weinberg, W. A. Wall, and A. John Hart, “Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations,” Powder Technology, vol. 343, pp. 855–866, Feb. 2019, doi: 10.1016/j.powtec.2018.11.072.