Contact models for discrete element method simulations of non-spherical particles

Abstract

Computational modelling of granular systems has increased significantly in recent years, particularly using the Discrete-Element Method (DEM), a simulation tool proposed by Cundall & Strack (1979). In the most common implementation of DEM, particles are modelled as rigid bodies, for which the deformations at the contact points are captured by permitting overlaps between the interacting particles. These correlations between force and overlap are referred to as 'contact models'.

Particles in 3D DEM simulations are generally modelled as spheres because of their computational simplicity. Particle shape though is an important factor, causing a major increase in scientific interest in DEM modelling of non-spherical particle systems. One of the greatest limitations of using non-spherical particles in DEM modelling is the lack of appropriate mathematical relationships between the contact force and interparticle overlap. For this reason, contact models suitable for spheres are usually adopted which is not physically justifiable.

One of the most common methods of emulating non-spherical granules is the Multi-Spherical (MS) approach, where a number of spheres are 'glued' together in a cluster in order to construct a more complex shaped granule. Such an approach, however, leads to two sources of errors: (i) the contact model error arising from adopting a contact model suitable for individual spheres and (ii) the shape approximation error arising from modelling a specific desired particle shape using a finite number of spheres. Taking spherocylinders as an exemplar shape, a Finite Element Analysis (FEA) study was conducted to quantify both sources of errors. The study focused on the force responses obtained from the uniaxial compression of spherocylinders and the collisions between two identical spherocylinders. It should be noted that only the normal contact forces are considered in this thesis and any tangential interactions were disregarded.

For the contact model error, the rods used in the tests consisted of three spheres and had a variable aspect ratio by adjusting the sphere overlaps between the constituent spheres. The force response comparison was performed between the FEA results of MS rods and the results that would have been obtained in a DEM setting when different contact models were applied. A number of contour plots were generated showcasing the errors of 3 different contact models: Hertz, Hooke, Tatara. For the shape approximation error all granules tested had a fixed aspect ratio of 3:1. The force response comparison was performed entirely using FEA where the force responses of MS rods with a variable number of spheres was compared against that of a perfect spherocylinder. A number of contour plots were generated showcasing the errors as the number of constituent spheres increases. The results indicate that the MS approach with spherical contact models should generally be avoided for contact overlaps exceeding 1.5% of the constituent spheres' diameters. Additionally, while greater sphere overlaps mean more accurate shape representation, sphere overlaps of 60% should not be exceeded if contact models designed for spheres are to be used.

Going beyond distance-overlap based contact models, two volume-overlap based contact models were evaluated by performing indentation tests on a spherical granule with different shaped indenters using FEA. The energy-conserving model proposed by Feng (2021c) [Feng, Y. T. (2021c). An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles: Contact volume based model and computational issues. Computer Methods in Applied Mechanics and Engineering 373, 113493, doi: 10.1016/j.cma.2020.113493] cannot be recommended for practicality reasons. A linear volume-based model though seems like a better alternative since the stiffness index constant k is independent of the contacting bodies' shapes according to the FEA conducted. However, both models rely on calculating an accurate contact volume overlap which is difficult to obtain for arbitrary shapes and not in good agreement with the deformed volume obtainable from FEM.

Finally, the influence of the Poisson effect on the stress dependence of the elastic moduli of soil was examined. Hertzian spheres neglect the Poisson effect, so it was examined whether this omission contributes to the inability of smooth-sphere DEM simulations to correctly capture the stress dependence of the elastic moduli of soil. This was done by isotropically compressing a spherical granule of silica sand using FEM. At low-to-moderate confining stresses the Poisson effect had a measurable but very limited influence. It was found that the Poisson effect became significant only at confining stresses on the order of 100 MPa. Thus, at lower stresses, rough-surface contact models remain the most justifiable way to match the stress–stiffness response measured in laboratory testing using DEM simulations.