There are so many engineering problems where capturing the micromechanics of individual particles plays a crucial role in determining the overall system behavior. Some examples include sand piling, powder metallurgy, mineral processing, particle packing etc. These problems are numerically modelled using Discrete Element Method (DEM).
DEM is a finite difference method used for predicting the contact behavior of enormous number of independently moving discrete particles. It is a versatile tool useful for modeling particle behavior in mining, ceramic, food, pharmaceutical, agricultural, chemical, oil and gas and other industries. This article gives a brief introduction to the formulation theory of DEM, time step requirements and elements used for modeling granular materials.
How Particles Interact In DEM
The real crux of this method comes in defining how particles interact when they come into contact with each other. When two spherical particles come into contact, they can essentially interact in 3 ways. They could be just touching without any deformation, the particles could deform at the contact without any penetration, or they could be rigid with some contact penetration.
There are two simulation methods for DEM based on whether deformation is allowed to occur upon contact between the discrete particles. The first one is hard-sphere approach in which the colliding particles do not deform and only momentum is exchanged between them. The second one is soft-sphere approach in which the colliding particles are assumed to be rigid but are allowed to have small overlaps which represent the deformations. Soft-sphere method is used in the Abaqus DEM methodology.
DEM uses Contact Mechanics and Newton’s laws to determine the position of the particles in the system. The contact detection algorithms are used to check for contacts in the domain and calculate forces based on the amount of overlap between the colliding particles. For example, the contact between spherical particles is detected when distance between the centers is less than or equal to their radii summation. Hence the algorithm checks for the distance between all the system particles. This part of the simulation is computationally expensive step and accounts for about 70 – 80% of the total computational effort. For this purpose, the domain is discretized into small cells generally 3-5 times the radius of smallest particle and active contacts are identified as shown in the figure below.
Once the active contacts are identified, contact force models are used to determine the force magnitudes acting on the particles based on contact stiffness and particle overlap. Contact force models such as Hertz normal contact and Johnson-Kendal-Roberts (JKR) model are available in Abaqus for this purpose.
The Math Behind Particle Interaction
The contact between particles can be represented by a spring – damper system as shown in the figure below. Kn represents the normal contact stiffness which can be either linear or non-linear and Cnrepresents the normal damping. The friction between the particles is represented by tangential spring stiffness Kt, coefficient of friction μ and contact damping Ct.
After calculating the forces, Newton’s second law is used to calculate the translational and rotational accelerations of the particles. The rotational motion is calculated using the equation below,
where I is moment of inertia, w is angular velocity and M is the angular moment acting on the particle.
The translational acceleration is calculated for the equation below,
where m is the mass of the particle, v is the velocity and F is the resultant force acting on the particle.
The velocity and position of the particles are then calculated by numerically integrating the acceleration over the time step. Forward time integration technique is used for this purpose in Abaqus.
where v(t) is the velocity of particle, x(t) is the displacement, a(t) is the acceleration and delta t is the time step. Once the new positions of the particles are calculated, they are repositioned in the grid. Similarly, the rotational position and velocity of the particles are updated, and the process is repeated to calculate the new collisions occurring in the next time step.
Time Step Requirements
Conditionally stable explicit dynamic procedure is used in DEM. The choice of time step plays a crucial role in avoiding excessive overlaps between the colliding particles, numerical instabilities, and preventing the effects of Rayleigh waves in DEM simulations. Rayleigh waves are the disturbance wave propagations coming from contact of granular particles. These waves affect the movement of particles that are far away from the origin point. To prevent the Rayleigh waves from propagating further than its neighboring particles, a small-time step based on wave speed is used. Usually, time increments less than
are suggested, where m is mass of particle and k is the contact stiffness. Very large particle velocities require a smaller time increment than the numerical stability limit.
DEM Elements
Each of the DEM particle is modeled using a rigid element with single node (PD3D). Each element has 6 degrees of freedom and can have translational and rotational motions. These are Lagrangian elements and hence we can apply other features such as connectors and constraints. Although these elements can only have spherical shape, complex granular shapes can be modeled by clustering the spherical particles together as shown in the figure below, which can be really useful when modeling things like pharmaceuticals mixing. The particles below approximate an elliptical shape. To get more accurate shape, we can add particles of different size to this cluster.
The particles in each cluster may overlap and can be either connected rigidly or by using constraint conditions. Unless contact exclusions are specified, contact forces will try to push the overlapping particles apart. For example, the particles in the above cluster have beam-type multi-point constraints between nodes of individual particles creating a rigid cluster. We can also define connector elements between nodes to get a deformable cluster. PD3D elements do not have any element output but all the nodal output variables are available in Abaqus.
Final Thoughts
The Discrete Element Method (DEM) is a numerical tool useful for analyzing discontinuous, granular material behavior. It plays a crucial role in understanding the behavior of individual particles in a system and thus reduces physical testing requirements early in the design process. Hopefully this blog helps you to understand the method, its formulation, and applications.
If you need help with FEA or simulation in general, don’t hesitate to give us a call! We’re passionate about what we do, and we would love to learn more about your unique engineering challenges.