DEM Application
WARNING: This page is still under construction.
The DEM Kratos Team
Theory
This application solve the the equations.... Mathematical approach to the problems.
Nothing numerical
Integration Schemes
Forward Euler Scheme
Contact Laws
Concept of indentation HMD, LSD
Normal Force Laws
Linear Repulsive Force
The most simple representation of a repulsive contact force between a sphere and a wall is given by a linear law, where the force acting on the sphere when contacting a plane is a linear function of the indentation, which in turn would bring a quadratic dependence with the contact radius. The next figure shows this simple law:
Hertzian Repulsive Force
On the other hand, Hertz solved in 1882 the non-cohesive normal contact between a sphere and a plane. In 1971 Johnson, Kendall and Roberts presented the solution (JKR-Theory) for the same problem in this case adding cohesive behaviour. Not much later, Derjaguin, Müller and Toporov published similar results (DMT-Theory).
Both theories are very close and correct, the main difference being that the JKR theory is oriented to the study of flexible, large spheres, while the DMT theory is specially suited to represent the behaviour of rigid, small ones.
The previous figure shows the standard representation of a Linear or Hertzian contact between a sphere and a wall. The distribution of contact pressures between both bodies follow a parabolic law.
JKR Cohesive Force
The preceding capture shows the a representation of a JKR contact between a sphere and a wall. In this case, the distribution of pressures between both bodies is more complex due to the formation of a neck at the boundaries of the contact. A later figure shows a detailed view of the pressures involved.
In the previous graphic, it is very interesting to note the existence of two singular values of contact radius: one for which the total forces acting on the contacting sphere is zero, and another for which the maximum value of adhesion is achieved.
In the previous figure, the blue area represents the distribution of pressures acting on the sphere when contacting a wall if a Hertzian Law is followed. On the other hand, the sum of both green and blue areas represents the JKR distribution. Note the larger values and the existence of adhesive behaviour at both sides of the pressures distribution.
An example of granular simulation without cohesive forces:
The same simulation as before, this time with cohesive forces in both sphere-sphere and sphere-plane contacts.
References:
V. L. Popov. Contact Mechanics and Friction (2010).
Tangential Force Laws
Damping Force Laws
(restit. coef)
Numerical approach (implementation)
Structure of the code (Strategy, Scheme, Element, Node, Utilities, functions frequently used like FastGet,...)
DEM strategies
Non-cohesive materials Strategy
Continuum Strategy
Search strategies
The contact search is a very important part of the method in terms of computational cost (range 60%-90% of simulation time in simulations with large number of particles) and it is possibly the most difficult part to treat when dealing with particles that have no spherical/circular shape.
The contact detection basically consists in determining, for every target object, what other objects are in contact with, and then, judge for the correspondent interaction. It is usually not needed to perform a search every time step, which is generally limited for the stability of the explicit integration of the equations of motion. Due to the fact that the search is an expensive step a lower search frequency can be selected without much loss of accuracy.
The most naïve method of search that can be set is the brute search. For every element, the method does a loop for all other elements checking for the contact. The order of the number of operations needed is quadratic: n2 . Normally, the strategy to avoid such an expensive scheme is to divide the contact search in two basic stages, a global search and a subsequent local resolution of the contact. In this case the computation time of the contact search is proportional to n log n: Global Contact Search It consists on locating the list of potential contact objects for each given target body. There are two different basic schemes: the Grid/Cell based algorithms or the Tree based ones. There are numerous different methods and variations for each type.
Figure 3.3: Tree-based algorithm structure
Figure 3.2: Grid/Cell-based algorithm structure
In the Grid based algorithms a general rectangular grid is defined dividing in cells the entire do-
main; a unified bounding box (rectangular or spheric) is adopted to represent the discrete elements
that can be of any shape; the potential contacts are determined by selecting the surrounding cells
where each target body is centred on.
6
In the Tree based algorithms each element is represented by a point. Starting from a centred one,
it splits the domain in two sub domains obtaining points that have larger X coordinate in one sub
domain and points with smaller values of the X in the other sub domain. The method proceeds for
next points alternating every time the splitting dimension and obtaining a tree structure like the
one shown in Figure 3.3 tree structure. Once the tree is constructed, for every particle, the nearest
neighbours have to be determined following the tree in upwind direction.
DEM schemes
DEM elements
Spheric Particle
Spheric Continuum Particle
Spheric Swimming Particle
Benchmarks
The DEM Benchmarks consist of a set of 9 simple tests which are run every night and whose object is to make sure both that the application performs correctly and that the code did not break after the daily changes. They are the following:
Test1: Elastic normal impact of two identical spheres
Check the evolution of the elastic normal contact force between two spheres with time.
If the coefficient of restitution is 1, the module of the initial and final velocities should be the same. Also, by symmetry, velocities should be equal for both spheres.
_
Test2: Elastic normal impact of a sphere with a rigid plane
Check the evolution of the elastic normal contact force between a sphere and a plane.
If the coefficient of restitution is equal to 1, the module of the initial and final velocity should remain unchanged.
_
Test3: Normal contact with different restitution coefficients
Check the effect of different restitution coefficients on the damping ratio.
If total energy is conserved, the restitution coefficient and the damping ratio values should be identical.
_
Test4: Oblique impact of a sphere with a rigid plane with a constant resultant velocity but at different incident angles
Check the tangential restitution coefficient, final angular velocity and rebound angle of the sphere.
_
Test5: Oblique impact of a sphere with a rigid plane with a constant normal velocity but at different tangential velocities
Check the final linear and angular velocities of the sphere.
_
Test6: Impact of a sphere with a rigid plane with a constant normal velocity but at different angular velocities
Check the final linear and angular velocities of the sphere.
_
Test7: Impact of two identical spheres with a constant normal velocity and varying angular velocities
Check the final linear and angular velocities of both spheres.
By symmetry, the tangential final velocity of both spheres should be zero. Additionally, for a coefficient of restitution of 1, there should be no changes in the modules of both linear and angular velocities and their values should conserve symmetry.
_
Test8: Impact of two differently sized spheres with a constant normal velocity and varying angular velocities
Check the final linear and angular velocities of both spheres.
In this case, it is interesting to note that, the bigger and/or denser sphere 2 is, the more this test resembles the sphere versus plane simulation.
_
Test9: Impact of two identical spheres with a constant normal velocity and different restitution coefficients
Check the effect of different restitution coefficients on the damping ratio.
If total energy is conserved, the restitution coefficient and the damping ratio values should be identical.
References:
Y.C.Chung, J.Y.Ooi. Benchmark tests for verifying discrete element modelling codes at particle impact level (2011).
How to analyse using the current application
Pre-Process
GUI's & GiD
D-DEMPack
D-DEMPack is the package that allows a user to create, run and analyze results of a DEM simulation for discontinuum / granular / little-cohesive materials. It is written for GiD. So in order to use this package, you should install GiD first.
You can read the D-DEMPack manual or follow the D-DEMPack Tutorials for fast learning on how to use the GUI.
C-DEMPack
Continuum / Cohesive
F-DEMPack
Fluid coupling
Post-Process
Application Dependencies
The Swimming DEM Application depends on the DEM application
Other Kratos Applications used in current Application
FEM-DEM
Programming Documentation
The source code is accessible through this site.
Problems!
What to do if the Discrete Elements behave strangely
In the case you notice that some discrete elements cross walls, penetrate in them or simply fly away of the domain at high velocity, check the following points:
In the case of excessive penetration:
- Check that the Young Modulus is big enough. A small Young Modulus makes the Elements and the walls behave in a very smooth way. Sometimes they are so soft that total penetration and trespass is possible.
- Check the Density of the material. An excessive density means a big weight and inertia that cannot be stopped by the walls.
- Check the Time Step. If the time step is too big, the Elements can go from one side of the wall to the other with no appearence of a reaction.
- Check the frequency of neighbour search. If the search is not done frequently enough, the new contacts with the walls may not be detected soon enough.
In the case of excessive bounce:
- Check that the Young Modulus is not extremely big. An exaggerated Young Modulus yields extremely large reactions that can make the Elements bounce too fast in just one time step. Also take into account that the stability of explicit methods depends on the Young Modulus (the higher the modulus, the closer to instability).
- Check the Density of the material. A very low density means a very small weight and inertia, so any force exerted by other elements or the walls can induce big accelerations on the element.
- Check the Time Step. If the time step is too big, the method gains more energy, and gets closer to instability.
- Check the restitution coefficient of the material. Explicit integration schemes gain energy noticeably, unless you use a really small time step. In case the time step is chosen to be big (but still stable), use the restitution coefficient to compensate the gain of energy and get more realistic results.
Contact
Contact us for any question regarding this application:
-Miguel Angel Celigueta: maceli@cimne.upc.edu
-Guillermo Casas: gcasas@cimne.upc.edu
-Salva Latorre: latorre@cimne.upc.edu
-Miquel Santasusana: msantasusana@cimne.upc.edu
-Ferran Arrufat: farrufat@cimne.upc.edu
