Vector-based model of elastic bonds for DEM simulation of solids

A new model for computer simulation of solids, composed of bonded rigid body particles, is proposed. Vectors rigidly connected with particles are used for description of deformation of a single bond. The expression for potential energy of the bond and corresponding expressions for forces and moments are proposed. Formulas, connecting parameters of the model with longitudinal, shear, bending and torsional stiffnesses of the bond, are derived. It is shown that the model allows to describe any values of the bond stiffnesses exactly. Two different calibration procedures depending on bond length/thickness ratio are proposed. It is shown that parameters of model can be chosen so that under small deformations the bond is equivalent to either Bernoulli-Euler rod or Timoshenko rod or short cylinder connecting particles. Simple expressions, connecting parameters of V-model with geometrical and mechanical characteristics of the bond, are derived. Computer simulation of dynamical buckling of the straight discrete rod and half-spherical shell is carried out.

💡 Research Summary

The paper introduces a novel bond model for discrete element method (DEM) simulations of solid materials, termed the vector‑based model (V‑model). Traditional DEM approaches treat inter‑particle bonds as simple linear springs and dashpots, which cannot simultaneously capture axial extension, shear, bending, and torsional deformations. In the V‑model each rigid particle carries three permanently attached vectors. The relative translation and rotation of these vectors between two neighboring particles fully describe the six‑degree‑of‑freedom deformation of a single bond.

A potential energy function is defined as a sum of four quadratic terms, each associated with a distinct stiffness: axial (k_n), shear (k_s), bending (k_b), and torsional (k_t). The energy expression reads

U = ½ k_n (Δl)² + ½ k_s |Δt|² + ½ k_b θ_b² + ½ k_t θ_t²,

where Δl is the change in bond length, Δt the shear displacement vector, and θ_b, θ_t the bending and torsional rotation angles, respectively. By differentiating U with respect to the relative position and orientation, the authors obtain explicit formulas for the force and moment transmitted through the bond. This derivation guarantees that the force–moment pair is energetically consistent, eliminating non‑physical energy drift that can plague conventional DEM formulations.

A central contribution of the work is the analytical mapping between the four bond stiffnesses and the classical mechanical properties of a cylindrical connector of length L and radius r. Using continuum mechanics, the authors show that

k_n = EA/L, k_s = GA/L, k_b = EI/L, k_t = GJ/L,

with E the Young’s modulus, G the shear modulus, A = πr² the cross‑sectional area, I = πr⁴/4 the second moment of area, and J = 2I the polar moment. Two calibration procedures are proposed depending on the aspect ratio L/r.

1. Long, slender bonds (L/r ≫ 1) – The transverse shear deformation is negligible, and the bond behaves as a Bernoulli‑Euler rod. In this regime only k_n and k_b are needed; k_s can be set to a very large value to suppress shear.

2. Shorter or thick bonds (moderate L/r) – Shear deformation becomes significant, and the Timoshenko beam theory must be employed. Here both k_s and k_b are retained, and the calibration uses the full set of continuum relations.

For bonds whose length is comparable to the particle diameter (short cylinders), the authors derive a generalized stiffness matrix that simultaneously accounts for shear, bending, and torsion, allowing any combination of (k_n, k_s, k_b, k_t) to be reproduced exactly. Consequently the V‑model can represent any desired bond stiffness tensor, making it suitable for heterogeneous or composite materials where axial, shear, bending, and torsional responses differ markedly.

The paper validates the model through two numerical experiments. The first simulates a straight discrete rod subjected to axial compression. By introducing a small geometric imperfection and gradually increasing the load, the rod buckles at a critical force that matches the Euler buckling prediction, confirming that the bending stiffness is correctly captured. The second experiment builds a half‑spherical shell from particles and applies uniform external pressure. The shell exhibits classic membrane‑bending buckling modes; the V‑model reproduces the coupled bending‑shear‑torsion response and conserves total mechanical energy throughout the dynamic simulation.

Overall, the V‑model overcomes the limitations of conventional spring‑dashpot bonds by providing a physically transparent, fully coupled description of bond mechanics. Its parameters are directly linked to material constants, enabling designers to calibrate the model against known rod, beam, or cylinder theories. The ability to reproduce Bernoulli‑Euler, Timoshenko, and short‑cylinder behavior makes the V‑model a versatile tool for multiscale DEM simulations, fracture and damage modeling, and the analysis of complex composite structures. Future work suggested by the authors includes extending the framework to nonlinear material laws, temperature‑dependent behavior, and explicit failure criteria, thereby broadening its applicability to a wider range of engineering problems.