Department of Mathematics
Matematiska institutionen
Master’s Thesis
Analysis and implementation of the Smooth
Discrete Element Method in AgX
Thomas Pettersson
LiTH-MAT-EX--2015/07--SELinköping 2015
Department of Mathematics Linköpings universitet SE-581 83 Linköping, Sweden
Analysis and implementation of the Smooth Discrete
Element Method in AgX
Master thesis in Computational Physics at
Algoryx Simulations AB
Thomas Pettersson LiTH-MAT-EX--2015/07--SE
Master’s Thesis: 30 hp
Thesis Level: A
Supervisor: Tomas Berglund
Algoryx Simulation AB
Examiner: Fredrik Berntsson
mai, Linköpings universitet
Division, Department Avdelning, Institution Computational Mathematics Department of Mathematics SE-581 83 Linköping Date Datum 2015-06-17 Language Språk Swedish/Svenska English/Engelska Report category Rapporttyp Licentiate Master Bachelor
URL for electronic version
http://liu.diva-portal.org/smash/record.jsf?pid=diva2% 3A856766dswid=-674 ISBN — ISRN LiTH-MAT-EX--2015/07--SE Serietitel och serienummer Title of series, numbering
ISSN —
Title Titel
Analysis and implementation of the Smooth Discrete Element Method in AgX Analy och implementation av den släta Diskreta Element Metoden i AgX
Author Författare
Thomas Pettersson
Abstract Sammanfattning
We encounter granular materials on a daily basis. We walk up a gravel path or we eat our breakfast cereals. When handling granular materials on an industrial scale it is important to do so efficiently, to avoid unnecessary energy losses, wear and tear. To help designing efficient tools for handling these materials engineers uses numerical simulations.
This project investigates the difference between the two main approaches to simulation of granular materials, the Smooth- and Non-smooth Discrete Element Methods by imple-menting the Smooth method into AgX dynamics were the Non-smooth method already is implemented, and then setup and execute a range of experiments to investigate their differ-ences.
The investigation shows both advantages and weaknesses for both methods. The result of simulations with smooth discrete element method are more consistent than with the non-smooth discrete element method with respect to choice of time step and other parameters that can be chosen for the simulation. Smooth discrete element method have problems when it comes to extreme situations.
The relative simulation time for system as large as treated by this project (more than 1000) can not be shown to depend on the size of the system.
Keywords
Abstract
We encounter granular materials on a daily basis. We walk up a gravel path or we eat our breakfast cereals. When handling granular materials on an industrial scale it is important to do so efficiently, to avoid unnecessary energy losses, wear and tear. To help designing efficient tools for handling these materials engineers uses numerical simulations.
This project investigates the difference between the two main approaches to simulation of granular materials, the Smooth- and Non-smooth Discrete Element Methods by implementing the Smooth method into AgX dynamics were the Non-smooth method already is implemented, and then setup and execute a range of experiments to investigate their differences.
The investigation shows both advantages and weaknesses for both methods. The result of simulations with smooth discrete element method are more consis-tent than with the non-smooth discrete element method with respect to choice of time step and other parameters that can be chosen for the simulation. Smooth discrete element method have problems when it comes to extreme situations.
The relative simulation time for system as large as treated by this project (more than 1000) can not be shown to depend on the size of the system.
Acknowledgments
I want to thank all the people at Algoryx for their friendly treatment. A special thanks to Tomas Berglund, my supervisor, for his patience with my stupid ques-tions and to Nils Hjelte, who has been invaluable when my programing skills haven’t been enough.
I also like to thank my sister, Sofia Viklinder, for helping me revise the lan-guage of the report, my opponent, Yi Chieh Chung, for her valuable comments on my report and my fellow master thesis worker, Simon Agvik, for helping me keeping up the spirit in a otherwise abandoned office during the vacation weeks. Finally I want to thank Fredrik Berndtsson, my examiner, for his excellent advice on everything from planing the project to writing mathematical equations during the entire project.
Linköping, September 2015 Thomas Pettersson
Contents
1 Introduction 1
1.1 Background . . . 1
1.2 Objective . . . 2
1.3 The project partner, Algoryx . . . 2
1.4 AgX . . . 2
2 Theory 5 2.1 Smooth DEM and Non-Smooth DEM . . . 5
2.1.1 SDEM . . . 5 2.1.2 NDEM . . . 7 2.2 Contact models . . . 7 2.2.1 Normal forces . . . 8 2.2.2 Tangential forces . . . 9 2.2.3 Rolling resistance . . . 11 3 Implementation 15 3.1 Models . . . 15 3.2 Model details . . . 16 3.3 Implementation . . . 17 3.4 Parallelism . . . 17 4 Experiments 19 4.1 Silo . . . 20 4.1.1 Results . . . 20 4.1.2 Simulation time . . . 24 4.2 Angle of Repose . . . 25 4.2.1 Results . . . 26 4.3 Drum . . . 27 4.3.1 Results . . . 28 4.3.2 Simulation time . . . 30 4.4 Blade . . . 31 4.4.1 Results . . . 32
4.5 High pressure test . . . 33 vii
viii Contents 4.5.1 Results . . . 34 5 Discussion 35 5.1 Conclution . . . 36 6 Future work 37 Bibliography 39
1
Introduction
1.1
Background
Granular materials are materials composed of rigid grains, e.g. cereal, sand, pota-toes, pharmaceutical tablets, iron and mineral ore balls etc. These materials are the second-most manipulated and transported materials in the world, only sur-passed by water[14].
A lot of the research about granular materials utilize computer simulations since it has several advantages compared to physical experiments. Simulations are not only cheaper in most cases, but also enables researchers to study the pro-cesses down to individual grains, something that can be extremely difficult with ordinary experiments on large numbers of grains. Investigation of these materi-als and how they interact with their surroundings is an important aspect when designing new equipment for industrial processes[15, 13]. Simulations gives in-sight in how a design solution will affect the process and provide an opportunity to test a larger number of designs than ordinary prototyping does. With these simulations engineers can improve efficiency and lifetime of the equipment[19].
The Discrete Element Methods (henceforth known as DEM) is a family of nu-merical methods that simulates particles in large numbers. DEM evolved from molecular dynamic simulations, and it’s the most common way to simulate gran-ular materials[15].
The two main approaches of DEM is the smooth- and the nonsmooth-DEM (in this project known as SDEM and NDEM respectively). The NDEM considers the stick-slip frictional transitions and collisions as instantaneous events and will thus have velocities that change discontinuously in time, following a contact law. The SDEM on the other hand consider the particles as elastic and bases the forces between the particles on the distance they have penetrated each other with. This demands a much shorter time step in the integration, but the calculations in each
2 1 Introduction
time step is not as time consuming as in NDEM[15].
1.2
Objective
The objective of this thesis is to contribute to the knowledge of the difference between the two different approaches to granular simulation, SDEM and NDEM. A few studies have previously been performed on the subject [15, 9, 6]. These studies have all compared different simulation software with each other. No case of a software capable of running both SDEM and NDEM have been found in the pre-study. To get a more fair comparison between the two models part of the project have been to implement SDEM in AgX, a physics engine that already had NDEM implemented. The investigation concerns the differences in experiments of scale that is interesting from a industrial point of view.
Previous investigation have shown that NDEM is computationally favorable for dense quasi static systems with few particles and with stiff particles, while SDEM gains advantages for simulations with more particles, more kinetic energy and softer materials[15, 6]. This project will mainly focus on qualitative differ-ences with some computational aspects.
1.3
The project partner, Algoryx
Algoryx Simulation AB is a small company with close connections to Umeå Uni-versity that specialize in physics simulations. Algoryx suggested this project. At their office, in Umeå, Algoryx provided a desk and a computer with software suitable to the project, as well as a lot of time to guide the project.
In their physics engine AgX Dynamics, Algoryx offers simulation of granu-lar materials using NDEM. AgX Dynamics has been used successfully by early adopters like LKAB and Atlas Copco. Algoryx would like an SDEM implementa-tion for several reasons:
• SDEM is the industry standard.
• NDEM can be faster than SDEM in certain circumstances and Algoryx wants a in-house reference to show customers.
• To show that NDEM produce the same results as SDEM.
• In some situations it would be favorable to use SDEM instead of NDEM.
1.4
AgX
The program that the development and tests have been conducted is the already mentioned AgX dynamics. The software AgX dynamics is marketed with the slogan
1.4 AgX 3
"AgX Dynamics is a professional multi-purpose physics engine for simulators, engineering, large scale granular simulators and more."[2] As the slogan suggests AgX is designed with a lot more than granular simula-tions in mind. The software is for example widely used in training simulators, according to the product brochure 7 out of 9 major suppliers use it [2].
The basis for the mechanics simulations in AgX is SPOOK [12], a variational time-stepping scheme for rigid multibody systems. The SPOOK formalism is designed for reducing real continuous mechanics to a discrete and simulation friendly formulation where the calculations are presented in matrix form.
AgX can simulate many types of objects such as wires, hydraulics, rigid bod-ies, constraints, etc. For this project two of the object types are interesting, the
GranularBodyand the RigidBody. Objects of GranularBody type have six
degrees of freedom, both directional and rotational, and can fully interact with each other and RidgidBodies through collisions. These collisions are the ones managed with NDEM or SDEM. They are, as of now, all spherical but can have different radius, bulk properties (like density and Young’s modulus) and sur-face properties (like friction coefficient and rolling resistance coefficient). The GranularBodyobject is commonly called particle in this report but should not be confused with the Particle object in AgX. The Particle object is much like the
GranularBodybut have no rotational degrees of freedom.
The RigidBody object has the same type of surface and bulk properties as the GranularBody object. The RigidBodies interact with each other as well as with the Granular Bodies but also with various kinds of dynamic constrains that are not described here since they are irrelevant to scope of the project. The
RigidBodyconsists of one or several Geometries. The Geometries can be of
several different shapes including spheres, boxes, cylinders and also with trimesh. The collision detection operates differently depending on the types of geometries involved to optimize the process. The RigidBody can be in three states:
• dynamic where it responds to all sorts of forces. • kinematic with predetermined movement. • static that are completely stationary.
2
Theory
This chapter describes the theoretical foundations of the project. First the fun-damental differences between SDEM and NDEM are described and then several different contact models are considered. The models are more closely described in [13, 8].
2.1
Smooth DEM and Non-Smooth DEM
The particles in granular mechanics typically only affect each other with short range interactions. More intuitively it can be said that all interaction require contact between the grains. All DEM simulations assumes that grains suffer no plastic deformation. That is true for both the SDEM and NDEM but there are several other fundamental differences between the two models[15].
2.1.1
SDEM
If nothing else is referenced the data in this section is taken from [13].
The main idea of SDEM is to integrate Newtons equations of motion and to calculate the forces explicitly by only using the current information about the system. Newtons equations of motion are
Fi(ri,j, vi,j, ϕi,j, ωi,j) = miai
Mi(ri,j, vi,j, ϕi,j, ωi,j) = Iiαi, j = (1, ..., N ) (2.1)
where Fi is the force on particle i and Mi is the torque on particle i, ri,j, vi,j,
ϕi,j, ωi,j is the relative position, relative velocity, relative angular orientation
and relative angular velocity of particle pair (i, j) respectively, mi is the mass of
particle i, Ii is the inertia tensor of the particle i, ai is the acceleration of center
6 2 Theory
of mass of the particle, αi is the rotational acceleration and N is the number of
particles in the simulation.
For a system of spherical particles the inertia tensor, Ii, can be simplified to
a scalar, Ii, and since spheres are the same in all directions the relative angular
orientation, ϕj, becomes irrelevant since all orientations are equivalent. The sim-plified version of Newtons equations then becomes
Fi(ri,j, vi,j, ωi,j) = miai,
Mi(ri,j, vi,j, ωi,j) = Iiαi, j = (1, ..., N ). (2.2)
For efficient calculations it is assumed that the force and torque is a sum of local contact forces, i.e the contact forces does only depend on ri,j, vi,j, ϕi,j, ωi,j for
that particle pair (i, j). The forces and torques are computed as Fi =
N
X
j=1,j,i
Fi,j(ri,j, vi,j, ϕi,j, ωi,j) + mig,
Mi = N
X
j=1,j,i
Mi,j(ri,j, vi,j, ϕi,j, ωi,j),
(2.3)
where Fi,j and Mi,j are the force and torque from particle j on particle i and g is
the gravitational acceleration.
As mentioned earlier the force between a particle pair only exist if the parti-cles are in contact. This leads to the conditions
Fi,j(ri,j > Ri + Rj) = 0,
Mi,j(ri,j > Ri+ Rj) = 0, (2.4)
with Ri and Rj symbolizing the radius of particle i and j.
In every interesting system of particles the number of pairs where Fi,j = 0 will
vastly outnumber the pairs of Fi,j , 0. For computational efficiency a collision detection algorithm creates a list of particles that are in contact and calculates
Fi =X k Fk(rk, vk, ϕk, ωk) + mig, Mi =X k Mk(rk, vk, ϕk, ωk), (2.5)
where k loops over the contacts that particle i is involved in.
The time step of an explicit method must be shorter than the the shortest time scale of the potential. The time scale, h, of a spring potential is
h =r m kn
. (2.6)
This is also the case in a non-linear spring that is often used in SDEM [15]. Thanks to the local nature of the method, SDEM can be executed in parallel without losing any information.
2.2 Contact models 7
2.1.2
NDEM
NDEM will only be explained conceptually since the implementation of the method isn’t a part of the project. For a more technical approach to NDEM, please con-sult [15]. To be able to use a longer time step, NDEM disregard the smaller time scale and replace it with constraints. The constraints usually are on the form particles shouldn’t overlap, they shouldn’t slip, shouldn’t roll and after a bounce the particle should have the correct velocity with respect to the restitution. The algorithm estimates how much force is needed to achieve these conditions during the time frame decided by the damping constant. During the force calculations limitations like the coulomb friction law and the hertzian contact force are con-sidered.
The drawback of this method is that everything has to be considered globally. This is done by a matrix operation. To solve the matrix equation several different approaches can be used. In the project a Projected Gauss-Seidel solver have been used. [15]
The global nature of NDEM makes parallel execution of the force calculation impossible. There are parallel implementations of NDEM but they aren’t fully equivalent [18]. To eliminate sources of error the standard implementation has been used in the project.
The time step, h, in NDEM is determined by how much overlap you are will-ing to permit in your simulation and by the typical normal velocity of a collision,
vn[15]. The formula is
h ≤ 2R vn
, (2.7)
where is the proportion of the diameter that can be permitted to overlap and R are the radius of the particle.
2.2
Contact models
The contact forces can be separated into several different parts. The normal force are the component that makes the particles separate when they overlap acting in the direction of the normal to the contact surface. Friction is added to the sim-ulation in the form of a force tangential to the contact surface. Finally a rolling resistance is needed for the model to give a realistic results, this part comes from a variety of physical properties, irregularities in the shape stands for most of it but even a perfectly spherical particle would have a resistance to rolling due to the energy losses that comes from the compression and expansion that occurs when the particle rolls on a surface[1]. There have been several models proposed for all of the different components of the contact forces. This section gives a brief description of the most commonly used models.
8 2 Theory
2.2.1
Normal forces
As long as only spheres are considered the normal force Fn depends only on the
geometrical overlap g
g = Ri+ Rj− |ri−rj|, (2.8)
where Ri and Rjare the radius of the two particles and ri and rjare the positions
of their centers. If g ≤ 0 there is no force between the particles[13].
Linear Spring–Dashpot model
It is quite trivial to make the assumption that all contact forces consist of an elastic part (because the particles bounce) and a dampened part (because the bounces does not conserve energy). The simplest ansatz and one often used in SDEM simulation is to model the conservative part as a linear spring and the dissipative as a dashpot, i.e.
Fn = kg + c ˙g, (2.9)
where k and c are spring and dampening constant respectively and ˙g is the time derivative of g, also known as the impact normal velocity.
The coefficient of restitution, defined as ≡ |˙g0|
|˙g| where ˙g 0
is the impact velocity after the impact normal and ˙g is the impact normal velocity after the impact. Integration of Newton’s equation of motion gives
= exp − πc 2mef f q k mef f − c 2mef f 2 , (2.10) with mef f ≡ mimj
mi+mj called the effective mass of the colliding particles.
As seen in the equation does not depend on the velocity, that is a problem since experiments tells us that usually is the case[13].
Hertzian contact law
This force law was first derived by Heinrich Hertz theoretically when he was considering elastic spheres. It has later been improved to include dampening from viscous contact. The improved version by Kuwabara and Kono [11] reads for spheres of the same materials
Fn =2Y √ Ref f 3(1 − ν2)(g 3 2 + A √ g ˙g), (2.11)
where Y denotes Young’s modulus, ν Poisson ratio, A is a function of material vis-cosity and Ref f ≡ RiRj
Ri+Rj. InComputational Granular Dynamics an expression for
2.2 Contact models 9
for the elastic part and a reasonable estimation for the dampening part) can be found on page 20[13]. The parameter A can be calculated analytically [4]
A = 4(1 − ν
2)(1 − 2ν)η
15Y ν2 , (2.12)
where η is the material viscosity constant.
The model preforms quite well in reproducing the coefficient of restitution and the impact time. A problem is that it under certain conditions produces an attractive force when A√g ˙g is the dominating term just before the end of the
impact[10].
In [10] an extended version of the model is suggested where the square root sign working on g in the damping part is exchanged for an arbitrary exponent that can be tuned to each element. This tuning will also change A.
Hysteretic models
In hysteretic models the normal contact force is modeled by linear springs with different spring constants during different phases off the impact. The faces are usually loading, ˙g > 0, and unloading, ˙g < 0, but can sometimes include oth-ers. There are also nonlinear versions with elastic descriptions as in the hertzian contact law. For details about the models, see [10].
Some of the models perform very well in reproducing the impact time and coefficient of restitution for different impact speeds and different materials[10].
2.2.2
Tangential forces
The tangential forces between particles are much harder to describe. One rea-son is that the surface interaction between materials aren’t as well known as the bulk behavior. However, the main reason is that the particles that are modeled are not the perfect spheres that we use as a model. In the tangential force this surface roughness is the main source, compared to the normal force where it was reasonable to omit it completely[13].
Most models for tangential friction incorporates Coulomb’s law of friction that states
Ft≤µsFn when vt= 0, (2.13)
Ft= µdFn when vt, 0, (2.14)
where vtis the relative velocities of the contact points of the particles and µsand
µdare the static and dynamical coefficient of friction. In every applied case that
have been studied the formulation have been simplified to
Ft≤µFn, (2.15)
with µs = µd= µ both due to that the coefficients rarely differ that much and due
10 2 Theory
Coulomb’s law of dynamic friction
The simplest way to formulate a model for the the tangential contact is to take Coulomb’s law of friction literally and let
Ft= −µFn
vt
vt
. (2.16)
This approach gives a behavior that have clear similarities with experimental data especially for high impact angles [17]. At low impact angles the model fails to reproduce the reversing of vtthat has been shown experimentally.
Furthermore the force has a discontinuity at vt = 0. That is unimportant in
a dynamic collision since the resulting vibrations will be limited in time but in a static case the vibrations will cause a collapse over time of all kind of unsup-ported structure.
Viscous Friction Force
A model that intuitively feels good and have been used is to say that Ftis
propor-tional to vt, i.e.
Ft= −γtvt. (2.17)
Unfortunately the proportionality constant γ lacks physical interpretation and the physical properties do not look much like experimental data either[17]. It can not form piles since vt ≈0 ⇒ Ft ≈0 and thus no forces support the heap
when the particles are stationary. At least there is no discontinuity.
Haff–Werner Model
In this model you combine the previous two, to get rid of the discontinuity from the first one. The force is calculated
Ft= − min(γ vt, µFn)
vt
vt. (2.18)
The resulting formula still can not form a heap for the same reason as the previ-ous model and do not reproduce the reversing of vt. Despite that it "was
success-fully used in many simulations"[13]. It comes closer to the experimental values for higher γ as it more and more looks like (2.16), see [17].
Tangential spring
To overcome the problem of not being able to form heaps, more complicated models are suggested. One of them is the tangential spring first proposed by Cundall and Strack.
2.2 Contact models 11
Upon colliding with each other two particles initiates a string where they first make contact. Ft= − min(ktxt, µFn) xt xt (2.19) xt = t Z t0 vtdt (2.20)
where ktis the spring constant that only can be determined after comparison with
experiments, t0the point in time when the particles first came into contact and t
is the present time. The model also take Coulomb’s friction law into account as seen in (2.19). The model predicts the experimental data quite well, including the reversal of relative speed for small impact angles[17]. The model is also capable of piling grains without support since the force between two stationary particles is not zero. If the particles leave each other the starting point of the spring is erased. It also seems to be a little less predictive together with the Hertzian con-tact model than with the Linear spring-dashpot model. A viscous friction force is often added as dampening. The resulting model is
Ft = − min(ktxt, µFn)
xt
xt
−γtvt. (2.21)
Walton and Braun model
The Walter-Braun model is based on the assumption that the forces will not change too much from one time step to the next. This leads to an iterative model
Ft= Ft 0
+ ∆Ft (2.22)
∆Ft = B∆xt (2.23)
Where0 marks that this is the value from the previous time step and ∆ marks
that this is the change since the previous time step. The B is not a constant but depends on previous forces,
B = B0 µFn−Ft µFn−F∗ γ , if slipping increases Ft, −B0µFn+Ft µFn+F∗ γ , if slipping decreases Ft, (2.24)
with γ a numerical parameter with a typical value of 13 and B0 the initial value
of B. F∗ is the value of Ft the last time the time the sign of vt changed[13]. The
force has good correspondence with experimental values [17]
2.2.3
Rolling resistance
The concept of rolling resistance [1] is introduced to account for several mecha-nism such as viscous hysteresis and energy loss plastic deformation (two different types of energy loss due to deformation), surface adhesion, shape deviation from
12 2 Theory spheres and other effects. In [1] four different classes of methods which will be described here are discussed.
A rolling friction coefficient, µr, is used in analogy with the friction coefficient,
µ. µr can be calculated from a simple experiment, as
µr= tan(θcrit), (2.25)
with θcritas the minimum angle of inclination that causes a particle to roll down
an inclined plane[1].
Directional constant torque models
This model is very analog to Coulomb’s law of dynamic friction. It adds a constant torque in the opposite of the rolling motion. The resulting rolling resistant torque thus become
Mr = −ωrel
ωrel
µrRef fFn, (2.26)
where ωrelis the relative angular velocity of the two particles. Since it is
mathe-matically equivalent to Coulomb’s law of dynamic friction, it has the same prob-lem i.e it has a discontinuity when passing through ωrel = 0.
It reproduced the stopping of a particle on a plane of no inclination save for the vibrations from the discontinuity. It could not reproduce a particle stopping on a inclined plane nor could it form an unsupported structure.
Viscous model
Viscous models have a resisting torque dependent on the angular velocity. There seems to be several different dependencies, the one deemed representative by [1] looks like
Mr = −µrRef fFn(ωiRi−ωjRj) (2.27)
a formulation that is indistinguishable from the viscous model of tangential fric-tion if the centers of the particles moves the same way, thus it share most of it’s properties, see Subsection 2.2.2. It reproduced the stopping of a particle on a plane of no inclination. It couldn’t reproduce a particle stopping on a inclined plain nor could it form an unsupported structure
Elastic-plastic spring-dashpot model
This model is related with previous models for both tangential and normal forces, the approach adapted by [1] is similar to the Walton and Braun model, in the meaning that the torque is calculated iteratively.
In this model the final torque is the sum of an elastic spring torque Msand a
viscus dampening torque Mdi.e.
2.2 Contact models 13
Starting with Ms it is, as suggested earlier calculated by iteration in each time
step,
Ms= M
0
s+ ∆Ms. (2.29)
If the absolute value exceeds Mmax = µrRef fFnit is set to that absolute value, in
analogy to Coulomb’s friction law. ∆Ms depends on the the incremental relative
rotation, ∆θr, so
∆Ms= −kr∆θr. (2.30)
The constant kris known as the rotational stiffness and there exists several
differ-ent proposed formulas for calculating it. The one used by [1] is derived for a 2D case, since the 2D formulation contains a free parameter, it will be assumed to be free here.
The viscous torque is not calculated incrementally instead
Md = −cr˙θ if Ms < Mmax, −f cr˙θ if Ms = Mmax, (2.31)
where cr is the rolling dampening constant and f is a number that decides how
much the viscosity is allowed to break the analog to Coulomb’s friction law. In [1] f = 0 is set to emphasize the difference between the models.
This model is deemed to be the most versatile, since it reproduces heap build-ing and stops rollbuild-ing particles on an inclined plane with less inclination than
θcrit[1]
Contact independent models
Some studies have used a model where the torque depends directly on the angu-lar velocity of the particle. Meaning that it does not depend on the other parti-cles. An extreme case has been used where the rotation of particles is completely locked. These types of models are ignored by [1] for their violation of physical principles.
3
Implementation
During the project, code for SDEM has been designed within the existing AgX software. The design is, as mentioned before, constructed to allow different con-tact models to calculate forces. The design also allows parallel computation. The collision detection is done by the existing collision detector in AgX which pro-vides contact data such as normals and contact overlaps etc. Likewise the integra-tion is done by the existing integrator.
The force calculation is implemented for parallel execution on an arbitrary number of threads.
3.1
Models
For the investigation the Hertzian contact law were chosen for calculation of the normal force since it is used as a penalty function in NDEM and it’s more based in physics[13]. For the tangential force and rolling resistance elastic spring models with damping have been used, despite the fact that they produce artifacts like oscillations for deaccelerating particles on a plane. The reasons to choose them is:
1. They can produce static assemblies of particles.
2. Academic papers find it more versatile than other models[16, 13].
3. It seems to be the standard in other commercial and freeware SDEM simu-lations.1
1For example of freeware program visit http://www.cfdem.com/
liggghtsr-open-source-discrete-element-method-particle-simulation-code
16 3 Implementation To increase efficiency by avoiding cache misses, all force calculations for a contact are executed together in the program. That happens at the expense of some flexibility and makes changing between different contact models difficult. However, that has not been an issue in the project since only one of models is used for each of the interactions.
3.2
Model details
After some testing a modified version of (2.11) was used, i.e
Fn= 2Y (1−ν2) √ Ref f(g 3 2 + Ag1/4˙g), if Fn≥0 0, otherwise (3.1) with A =q β Y (1−ν2) √ Ref fmef f . (3.2)
The benefit of (3.1) is that the restitution becomes constant. Although that might not be the most physically correct behavior [10], it is consistent of that of the NDEM implementation. This approach also works with the experimental values that are used in the simulation[19]. This study only has one value of restitution and does not specify at what speed the impact occurred.
In equation 3.2 β is a numerical parameter fitted with a second degree poly-nomial to experimental data of restitution.
The tangential friction is based of equation 2.21 with
kt = 8Gef f q Ref fg, (3.3) and γt= −β q ktmef f, (3.4)
where Gef f is the effective shear modulus. Under the assumption that the mate-rial is homogeneous and isotropic Gef f can be calculated from the Young’s mod-ulus Y , and the Poisson’s ratio ν, see [3]. This means that
Gef f = Y
4(2 − ν)(1 + ν). (3.5)
When calculating the elastic part of tangential force, Coulumb’s friction law is taken into account by adjusting the length of x,
x= x x µFn kt , if x > µFn kt , (3.6)
after adding the dampening part of the force the final result is again adjusted to Coulomb’s friction law
Ft= Ft
Ft
3.3 Implementation 17
For the rolling resistance equation (2.30) is used with
kr = ktRef f 2, (3.8)
and
cr = 0. (3.9)
3.3
Implementation
In Subsection 1.4 the fact that AgX in general work with matrices and in Subsec-tion 2.1.2 the same is said about NDEM. On the other hand, SDEM is not usually implemented using a matrix formulation and thus the project could not use the preexisting calculation pipeline to calculate the forces. Fortunately the collision detection and the integration of positions and velocities could be reused.
In AgX different Tasks executes Kernels in parallel. To implement SDEM in AgX the Task for NDEM was copied and the Kernels relevant to the force calculations were replaced with the force calculation Kernels and force appli-cation Kernels. The force calculation Kernel have to have different versions for the interactions between tow GranularBodies and the interaction between a GranularBody and a RigidBody. Similarly the force application Kernel have to be in two different versions, one for the GranularBodies and one for the RigidBodies.
During the project the plan to have separate Kernels for normal force, tan-gential friction force and rolling resistance was abandoned for a composite Ker-nelto make the calculations more effective and avoid some cache misses.
3.4
Parallelism
The implementation is intended to work fully in parallel. It does so by calculating the forces in one step where the the loop is over the contacts (light blue and brown blocks in Figure 3.1) and applying them in another step (purple blocks) where the loop circles all the particles and geometries.
As can be seen in Figure 3.1 most threads are active most of the time. The collision detection including hash table construction to save the values needed for Elastic-plastic spring-dashpot model for rolling resistance and Tangential spring model for tangential friction (the blue box) takes about 50% of the time and about 40% of the time is spent on force calculation and applying forces to the particles (red box).
18 3 Implementation
Figure 3.1: Profiling of a typical time step for the Drum experiment with 4 threads represented by the 4 lines, see Section 4.3 for a riguros explanation of the Drum experiment, the experiment simulates 95052 particles
4
Experiments
All the experiments have been conducted with the same particle parameters, see Table 4.1. The are the parameters validated by Wang et al. [19]. The parameters corresponds to iron ore pellets in the making. For simplicity the same parameters are adopted for the rigid bodies and the interactions between the rigid bodies and the particles if not otherwise specified. The same radius is used for all particles for simplicity.
Property Abbreviation Value
Radius r 0.00636665m
Young’s modulus Y 6.1684 · 106P a
Density ρ 3700kg/m3
Tangential friction coefficient µt 0.91
Rolling resistance coefficient µr 0.3211
Coefficient of restitution e 0.18
Table 4.1:Particle parameter values, common for all experiments. As described in Equation 2.6 the time step, h, of SDEM depends on the mass of the particle and the Young’s modulus so the same time step can be used for all the experiments. A time step h = 5 · 10−5s have been used. The time step for
NDEM depends on the expected maximal particle velocity and is thus individual for the different experiments.
The reason for not using the parallel implementation was that it is not 100% equivalent to the serial one and there was not time to investigate both cases. The numerical Gauss-Seidel solver was set to run 200 iterations, a number that was deemed adequate. The number of iterations are fix due to that AgX also is used for real time simulations[2].
Every experimental setup have been tested five times, each test with a differ-19
20 4 Experiments ent scaling parameter a. In the different experiments a takes the values 2, 4, 6, 8, 10. How the setup depends on a is explained in each experiment description. The number in the legends of all figures symbolize the number of a corresponding to that line.
Old experimental setups for AgX such as the ones used in [15] and in [19] could not be used since they relies heavily on spawning particles at random loca-tions.
When particles being spawn at overlapping positions, NDEM calculates a force big enough to separate the particles in the dampening time span (up to the force dictated by the overlap size and the Hertz model) see Subsection 2.1.2.
SDEM on the other hand has no such considerations relying on the short time step to never put it in this kind of situation and uses the full force of the Hertz formulation at once. Resulting in an explosive separation of the particles.
4.1
Silo
In the silo experiments the behaviour of stationary particles is investigated. The experiment setup consists of a hollow trimesh cylinder with dimensions as in Ta-ble 4.2. The cylinder is placed with the height axis vertically and a geometric cylinder of the same radius * 1.075 placed at the bottom of the hollow cylinder. The particles spawn in a cuboid inscribed in the cylinder and it’s vertical exten-sion to 2.5 times the height of the cylinder. In the cuboid the particles are placed in a quadratic pattern, 2.5 times the particle radius apart. The particles are then given an offset of up to ±10% of the spacing.
Property Abbreviation Value
Radius Rsilo 0.2m
Height Hsilo a ∗ 0.225m
Table 4.2:Silo, geometric measures
Since the goal of this experiment is to investigate the behaviour of stationary particles the time step of NDEM could in principle be very large. But trial and er-ror showed that the time step could not be larger than h = 2 · 10−3s if the particle
was not to tunnel through the walls of the silo during the initial formation of the stable configuration.
4.1.1
Results
One of the characteristic features of granular materials is that force chains appear [7]. A Force chain is when a chain of particles have significantly stronger forces between them than tho other sorounding particles. As can be seen in the force network visualization in Figure 4.1 both SDEM and NDEM shows these force chains. However the force chains of SDEM are much more pronounced than those of NDEM. It appears that NDEM has more force chains that do not stand out as much from the non chain contact as the SDEM force chains do.
4.1 Silo 21
Figure 4.1: Force network visualization for SDEM (left) and NDEM (right)
after 3 s of simulation, a=4, contacts marked with lines, yellow for weak forces and darkening until forces of 8N that are black.
Another key characteristic of a granular material is that the pressure does not increase linearly with the depth but instead stop increasing at some point. In Figure 4.2 the magnitude of the forces between the particles and the geometries are plotted against the height. The first thing to notice is that the dispersion in force is much bigger for SDEM. But that is to be expected considering Figure 4.1 that clearly shows that larger forces are more common in SDEM.
From Figure 4.2 it seems obvious that the pressure of NDEM keeps constant below a certain level. SDEM also looks like it keeps the pressure limited.
22 4 Experiments
Figure 4.2: Plotting particle-geometry normal force against height of the contact, a=10, blue for NDEM and red for SDEM
Figure 4.2 shows that the highest particles are higher in the SDEM experiment than in the NDEM experiment. Figure 4.3 shows the height of the highest particle in the different experiments. The most interesting part of the graph might be that the lines for a = 2. It is barely visible but the NDEM line is slightly above the SDEM line for quite a long time, until around 3s.
Figure 4.3:Height of the highest particle in the silo experiment
The NDEM simulations have more contacts active in all stages of all simula-tions, as seen in Figure 4.4, the difference in the tests with a = 2 is most
note-4.1 Silo 23
worthy since the particle level for that experiment was practically the same and still it shows a clear difference in its number of contacts. This is probably due to the stronger force chains does not have as many contacts with each other as the weaker and more evenly distributed NDEM force chains seen in Figure 4.1. In the larger experiments the effects are increased by the fact that the NDEM have a more compressed configuration.
Figure 4.4:Number of contacts during the simulation
In Figure 4.5 it’s pretty clear that the sum of the absolute values of the normal forces in NDEM is less than those of SDEM. Considering Figure 4.4 the difference in force per contact is even bigger.
24 4 Experiments
Figure 4.5:Sum of normal forces between particles
4.1.2
Simulation time
The times to finish the 5 experiments are plotted in Figure 4.6. The time evolution seems to be linear in this case, both for SDEM and NDEM. Since the geometry of the Silo experiment is such that the radius of the silo is kept constant, the total number of contacts should grow approximately linearly with number of particles. Since the main source of time consumption should come from the number of contacts in NDEM and the collision detection, that also seems like it should be related to the number of contacts. Having that in mind, Figure 4.6 suggests that a stationary system has a linear evolution with respect to number of contacts.
4.2 Angle of Repose 25
Figure 4.6:Time used to preform the silo simulation per particle in the sim-ulation
4.2
Angle of Repose
In the Angle of Repose experiments the behaviour of a large group of particles forming a heap on a plane is investigated.
The experiment setup consists of a hollow trimesh cylinder with measure-ments as described in Table 4.3. The cylinder is placed with the height axis
ver-Property Abbreviation Value
Radius Rsilos a ∗ 0.087m
Height Hsilo a ∗ 0.027m
Table 4.3:Silo in the angle of repose experiment, geometric configuration tically in the middle of a plane with sides 16 times the silo radius. The particles spawn in a cuboid inscribed in the cylinder and it’s vertical extension to 3 times the height of the cylinder. In the cuboid the particles are placed in a quadratic pattern, 2.5 times the particle radius apart. The particles are then given an off-set of up to ±10% of the spacing. After an initial stabilizing phase of 2.0 s the cylinder is lifted away from the plane at a high velocity. To avoid that the parti-cles are disturbed by the motion of the cylinder the friction and rolling resistance coefficients between the particles and the cylinder are set to zero.
The goal of the experiment is to measure the shape after the collapse. It is hard to estimate the maximum collision speed in the experiment but a time step, h = 0.001s, demonstrated good results in initial tests so it was used. The simulated time of the experiment is 6.1 s.
26 4 Experiments
4.2.1
Results
Comparing picture (a) and (b) in Figure 4.7 that is at the same scale1, shows that the pile in the NDEM experiment is slightly higher than in the SDEM experiment. The angle of the repose is lower than the 34◦observed by D. Wang et al[19]. The NDEM pile is steeper than the SDEM pile and is thus closer to the expected angle. The same pattern can be seen in all sizes of the simulation.
(a)NDEM (b)SDEM
(c)35◦protractor for comparison[5]
Figure 4.7:Shape of pile a=10 at the end of the simulation
In Figure 4.8 the part of the simulation that the actual experiment created to investigate does not say very much2.The forces are larger in the NDEM case,
something that is not surprising given the shape of the piles. More surprising is the results before the removal of the cylinder. The setup is very similar to the one in the silo experiments, but still the graph looks vastly different. The proportions of the cylinders are different but that doesn’t explain such a huge difference, instead the lack of friction against the walls are the main difference. It seems like SDEM depends a lot on the friction of the walls to stabilize where NDEM can do that on it’s own. The damping of the normal force is probably the main issue here. To base the damping on the relative velocity of the particles does not seem to be enough to dampen the system when many particles stack on top of each other. This behavior has also been observed for 1D towers of particles during implementation tests.
1At least approximately, movement to get the scene in a good position was not automated, but at
least systematic
4.3 Drum 27
Figure 4.8: Sum of normal forces between particles in Angle of Repose ex-periment
4.3
Drum
In the Drum test a more dynamic situation than in the previous experiments is investigated. The experiment setup consists of a cylinder made out of 200 partially overlapping boxes (dimensions 0.2m x 0.1m x 0.1 m). The boxes are placed in a circle in the X-Z plane and rotated so that one of the longer edges points towards the center of the circle. To prevent particles from falling out of the drum a thin disk is added at either side of the drum and the friction and rolling resistance coefficients are set to 0 in the contacts between the disks and the particles. The drum has a fix angular velocity ωdrum= 0.5rad/s close to the
0.53rad/s observed in [19] The drum on the other hand was significantly smaller than that of D.Wangs paper and is not constructed the same way.
Property Abbreviation Value
Radius Rdrum 1.0m
Width l a ∗ 0.027m
Angular velocity ωdrum 0.5rad/s
Sides of particle spawn area L
√ 2 ∗ a/10 Table 4.4:Geometric data for the drum
The particles spawn in a cuboid as wide as the drum and with sides that mea-sure
√
2 ∗ a/10. The cuboid is centered around the center of the drum. In the cuboid the particles are placed in a quadratic pattern, 2.5 times the particle ra-dius apart. The particles are then given an offset of up to ±10% of the spacing.
28 4 Experiments
Based on (2.7) and the assumption that the maximum impacts speed is equal to the maximal speed of the kinematic objects during the parts of the simulation that are of interest. vn = rdrum ∗ ωdrum = 0.5m/s and = 0.05 gives the time
step, h = 0.001s.
4.3.1
Results
as seen in Figure 4.9 both NDEM and SDEM shows the ability to recreate the observed behaviour that the top layer of the particles in a drum are in a flowing state while the bottom layer are in a stable configuration co-moving with the drum[7]. The shapes of the slopes are slightly different SDEM have a straighter slope than NDEM while the NDEM slope follows the shape of the drum slightly more. According to D. Wang et al. [19] the slope should be at an angle of 35◦
.
Both of the simulations display approximately 35◦
below the midpoint of the cylinder. The smaller radius of the simulated drum makes the importance of the observation questionable.
The flowing particles are moving in about the same velocity in both pictures since they have the same colour.
(a)NDEM (b)SDEM
(c)35◦protractor for comparison[5]
Figure 4.9:Particles in the drum, a=8, colored after velocity at the end of the simulation, blue = slow green=faster
The results in Figure 4.10 shows that the lines follow each other quite well throughout the simulation. In the stabilized part (3.5s and forward) that was the main interest of the experiment the two smallest SDEM experiments stands
4.3 Drum 29
out by having a slightly smaller velocity. Under the assumption that the flowing particles have about the same velocity (Figure 4.9) the conclusion is that a smaller proportion of the particles are flowing in the small SDEM systems but that the difference disappears when there are more particles in the drum.
Figure 4.10 also shows that the SDEM simulations start the transition into the stable flowing region slightly before the NDEM counterpart. That indicate that the NDEM particles should follow the edge of the drum slightly higher. The phenomenon that can be seen in Figure 4.9.
Figure 4.10:Average velocity of particles
The sum of the normal forces in the drum, Figure 4.11, follows the pattern established in the Silo, Figure 4.5, but to a lesser extent than in the silo. The differences are bigger before the particles have started to flow. This observation indicates that the forces between the flowing particles are equal between NDEM and SDEM or potentially that the SDEM produces smaller forces in the flowing region than NDEM do.
30 4 Experiments
Figure 4.11:Sum of Normal forces
4.3.2
Simulation time
Just as in the most static experiment, the Silo, the most dynamic experiment, the Drum, has a time dependence that seems close to linear with the number of par-ticles, Figures 4.12 and 4.6. The increase might be a little more than linear but there is not enough data points to say for sure. The contacts in the simulation once again increase approximately linearly with the number of particles so the linear behavior in number of contacts still seems like a valid approximation for this many particles. From the experiments there is no reason to believe that there is a difference in time evolution between the static and the dynamic case. How-ever NDEM has the potential to be faster in a static experiment since the time step can be changes after the velocity in the collisions.
4.4 Blade 31
Figure 4.12:Time consumption for running the simulation of the drum
4.4
Blade
The goal of the blade experiment is to investigate the interaction between mov-ing rigid bodies and the particles. The experiment setup consists of a cuboid container constructed from AgX’s box geometries. The inner measures of the con-tainer can be readed in Table 4.5. The particles spawn in a cuboid in the concon-tainer and the containers vertical extension to 1.5 times the height of the container. In the cuboid the particles are placed in a quadratic pattern, 2.5 times the particle radius apart. The particles are then given an offset of up to ±10% of the spac-ing. A second after the simulation starts a second rigid body, known as the blade, with the measurements from Table 4.5 , starts to move downward into the parti-cles. After 3 seconds the bottom of the blade has reached 15 of the height of the container. The blade then starts moving to the right at a fix velocity of 0.1 m/s.
Property Value Container width 0.3m Container length 0.5m Container height 0.3m Blade velocity 0.1m/s Blade width 0.01m Blade length 0.1m Blade height 1m
Table 4.5:Values for blade experiment
Equation 2.7 and the assumption that the maximum impacts speed is equal to the maximal speed of the kinematic objects during the parts of the simulation
32 4 Experiments
that are of interest give vn = 0.1m/s and = 0.05. The time step then becomes,
h = 0.005s but since time steps that long resulted in particles tunneling out of
the silo in the silo experiment the time step of the blade experiment was adjusted to h = 0.002s. The experiment simulated 5.1s.
4.4.1
Results
The images from the experiments are shown in Figure 4.13. The Figure shows a slightly smaller angle of repose For SDEM. That is consistent with the earlier images, Figure 4.7, and to some degree Figure 4.9.
(a)NDEM (b)SDEM
Figure 4.13:Blade experiment, image of the particle configuration at the end of the experiment, seen form the left side in the above description
Figure 4.14 shows the sum of the normal forces directed from the blade to the particles. Forces in all directions are added up but the forces directed forward should make up the majority of the sum. In the main area of the interest in this experiment, the section of horizontal movement, the 2 methods produce very similar results. In the period of vertical movement there is a clear difference, SDEM has more forces during this phase. There could be, at least, two different explanations for this. The first possible explanation is that the particles simulated by SDEM are more reluctant to move out of the way of the blade when it affect them vertically. That would result in a higher force acting on the bottom of the blade. The other possible explanation is that the generally higher forces in SDEM observed in Figure 4.5 is affecting the sides of the blade. The data saved in the experiment is not on a form that permits to see on what edge the particle touches in an easy way so that mystery will have to wait for another investigation.
4.5 High pressure test 33
Figure 4.14:Sum of normal forces between the blade and the particles.
4.5
High pressure test
In the high pressure test the behaviour of particles in extreme conditions are investigated. The experiment setup consists of a hollow trimesh with inner mea-surements as in Table 4.6 The cylinder is placed with the height axis in the verti-cal plane and a geometric cylinder of the same radius * 1.075 placed at the bottom of the hollow cylinder. The particles spawn in a cuboid inscribed in the cylinder and the cylinders vertical extension to 2.0 times the height of the cylinder. In the cuboid the particles are placed in a quadratic pattern, 2.5 times the particle radius apart. The particles are then given an offset of up to ±10% of the spacing.
Property Abbreviation Value
Radius Rsilos 0.2m
Height Hsilo a2∗0.028m
Table 4.6:Silo in the High pressure test, geometric properties
Over the particles a cylindrical plate of the same radius as the trimesh cylin-der is placed. At first the plate is stationary for one second. After a second the plate moves downwards at a constant velocity for 4 seconds and reach a height of 0.2 times the start height.
Based on (2.7) and the assumption that the maximum impacts speed is equal to the maximal speed of the kinematic objects during the parts of the simulation that are of interest vn = 0.1m/s and = 0.05 gives the time step, h = 0.0005s
34 4 Experiments
4.5.1
Results
As seen in Figure 4.15 something is broken in the SDEM simulation. The main explanation to this behaviour would be that the dampening is not enough to keep the small vibration in the particles from expanding into fast movement when the particles are too compressed. Another possible explanation is that the time step might be too large. I have not been able to find a derivation of Equation 2.6 for a herzian spring. The use of it might come from a linearization that is not valid at too high pressure. However, if that was the case the expectation would be to see particles escape from the cylinder since that is what have been observed for too high time steps during the implementation. The NDEM simulation keeps stable.
5
Discussion
In most cases both NDEM and SDEM produce results that seems plausible given the known behaviour of granular material. The exception is that SDEM produces results that are strange in two experiments, in the build up section of the An-gle of Repose experiment, see Section 4.2, and in the High Pressure experiment, see Section 4.5. A common denominator of both these situations is that they are completely impossible to recreate in a physical experiment. The Angle of Repose experiments have friction free walls of the silo. In the High pressure test the particles are compressed at such a degree that they could not be expected to re-turn to their original form, but would rather be broken or at least permanently deformed. The basic assumption of DEM, that the particles will return to their original shape is thus not realistic.
The quantitative aspect of the simulations varies significantly between the two models. The only experiment where there have been relevant experimen-tal data to compare with is the Angle of Repose experiment, section 4.2. Here both of the models fail to reproduce the expected results. In the article by D. Wang et al[19], on the other hand, the angle of repose is correctly reproduced by NDEM. That lead us to believe that experiment setup disappointingly enough might not be appropriate to achieve the angle of repose. The reproduced angle is lower for SDEM not only in the Angle of Repose experiment but also in the Drum experiment, section 4.3, and the Blade experiment, section 4.4. Since nei-ther of the results can be verified by experiments it is uncertain which results are the most correct. That NDEM previously reproduces correct angles of repose does not seem favorable for SDEM. Another experimental setup could of course change the angle for SDEM to. All other results can only be described by their difference without evaluating which one is most correct until they can be verified by experiments.
When evaluating the time usage of the experiments it is important to notice 35
36 5 Discussion
that the particles used in these experiments are extremely soft. Harder materi-als would keep the simulation time for NDEM essentially the same but the time used for SDEM would increase dramatically. Important is that it seems like both models have a close to linear evolution with the size of the system. NDEM can also be sped up considerably using the parallel implementation[18].
The results are affected by several factors. Two important factors are the choice of time step and number of iterations of the solver. While the choice of time step is motivated by equation 2.7 there is nothing that dictates what num-ber of iterations that should be chosen in any given situation. At Algoryx they have seen a connection between tall systems and a need for higher numbers of it-erations. Thus the silo experiments with large multiplier are the most vulnerable for these kind of errors. That is probably the reason that the silo experiment have some strange results.
The relative simulation time depends on the stiffness of the material and how long time steps the user are willing to risk in the NDEM simulation. The number of particles does not seem to affect the relative simulation time for systems as large as the ones investigated in this project.
5.1
Conclution
From the project was concluded that:
• SDEM gives a more constant result since there are fewer parameters to chose from and the parameters that you can choose are more restricted. • NDEM gives the user more chooses so the user can compromise between
fast simulations and exact ones.
• SDEM isn’t as stable as NDEM and give weird results in extreme conditions. • It is difficult to say which of the models that gives the best results without
6
Future work
To further investigate which of the models that hold most predictive power, a series of experiments have to be conducted to compare with simulated results. Several kinds of particles with different friction coefficients, Young’s modulus etc. should also be tested.
To compare the computational efficiency more justly there are several things that can be done for both of the models. NDEM can be performed with the par-allel Projected Gauss-Seidel solver [18] and the solver can be fed with the results of the previous time step to start of closer to the correct answer in the solver and thus limit the need for iterations.
SDEM’s good properties for parallel computing makes it interesting to create an implementation for GPU computation. Optimizations in the force calculations for SDEM would probably not be extremely efficient since it only takes about 40% of the time in the loop, see Figure 3.1. To get major increase in speed a new collision detection, maybe GPU based, would also have to be included.
On the most basic level of the area there is still work to do regarding the friction model for SDEM. A new model for friction that does not give rise to the artifacts that the current one does would be welcome. But since friction is an admittedly hard problem to simulate some work on how to minimize the artifacts with the existing spring model would be welcome.
Similarly I do not find it obvious that rolling resistance should be based on the relative rotation. This is mainly an issue while the particles are sliding on each others surfaces. A sensible expression for the mechanisms that gives the rolling resistance behaviour for objects that are sliding against each other would be much appreciated. Some experiments should be designed to find out how the rolling resistance works between two rotating bodies. It might not make a big difference in the end result but it would be important for the understanding.
The unknown dependency of the number of iterations of the solver for NDEM 37
38 6 Future work
is an issue that make the results doubtful since it can never be sure that the results will not be better with more iterations. A guide to what number of iterations that is "good enough" depending on the situation would be most helpful.
Bibliography
[1] J. Ai, J.-F. Chen, J. M. Rotting, and J. Y.Ooi. Assessment of rolling resistance models in discrete element simulations. Powder technology, 206:269–282, 2011. Cited on pages 7, 11, 12, and 13.
[2] Algoryx Simulation AB. Agx product brochure.
http://www.algoryx.se/mainpage/wp-content/uploads/2014/09/AgX-product-brochure.pdf, 9 2014. Cited on pages 3 and 19.
[3] D. R. Askeland and P. P. Phulé. The science and engineering of materials. Springer, 2003. Cited on page 16.
[4] N. V. Brilliantov, F. Spahn, J.-M. Hertzsch, and T. Poschel. Model for colli-sions in granular gases. Phys. Rev. E, 53:5382––5392, 1996. Cited on page 9.
[5] F. C. for instructional technology. Protractor, 35 degrees, 2004. Cited on pages 26 and 28.
[6] S. Hedman. Smooth and non-smooth approaches to simulation of granular matter. Master’s thesis, 2011. Cited on page 2.
[7] H. M. Jaeger, S. R. Nagel, and R. P. Beringer. The physics of granular materi-als. Physics today, 49, 1996. Cited on pages 20 and 28.
[8] M. Jean. The non-smooth contact dynamics method. Computer methods in applied mechanics and engineering, 177(3):235–257, 1999. Cited on page 5.
[9] J. Kleinert, M. Obermayr, and M. Balzer. Modeling of large scale granular systems using the discrete element method and the non-smooth contact dy-namics method: a comparison. In Proceedings of the ECCOMAS Multibody Dynamics Conference, Zagreb, 2013. Cited on page 2.
[10] H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, and V. Scherer. Review and extension of normal force models for the discrete element method. Pow-der Technology, 171:157–173, 2007. Cited on pages 9 and 16.
40 Bibliography [11] G. Kuwabara and K. Kono. Restitution coefficient in collision between two spheres. Japanese Journal of Applied Physics, 26:1230–1233, 1987. Cited on page 8.
[12] C. Lacoursiere. Ghosts and machines: regularized variational methods for interactive simulations of multibodies with dry frictional contacts. PhD the-sis, 2007. Cited on page 3.
[13] T. Pöschel and T. Schwager. Computational Granular Dynamics. Springer-Verlag, 2005. Cited on pages 1, 5, 8, 9, 10, 11, and 15.
[14] P. Richard, M. Nicodemi, R. Delannay, and D. Ribiere, P.and Bideau. Slow re-laxation and compaction of granular systems. Nature Materials, 4(2), 2005. Cited on page 1.
[15] M. Severin, D. Wang, C. Lacouresière, and K. Bodin. Examining the smooth and nonsmooth discrete element approaches to granular matter. Interna-tional Journal for Numerical Methods in Engineering, 97(12):878–902, 2010. Cited on pages 1, 2, 5, 6, 7, and 20.
[16] J. Shäfer, S. Dippel, and D. Wolf. Force schemes in simulations of granular materials. Journal de physique I, 6(1):5–20, 1996. Cited on page 15. [17] J. Shäfer, S. Dippel, and D. Wolf. Force schemes in simulations of granular
materials. Journal de Physique, pages 5–20, 1996. Cited on pages 10 and 11. [18] J. Sundberg. Parallel projected gauss-seidelsolver for large-scale granular
matter: Examining the physics of the parallel solver and development of amultigrid solver. Master’s thesis, 2014. Cited on pages 7, 36, and 37. [19] D. Wang, M. Servin, T. Berglund, K.-O. Mickelsson, and S. Rönnbäck.
Parametrization and validation of a nonsmooth discrete element method for simulating flows of iron ore green pellets. Powder technology, 283:475–487, October 2015. Cited on pages 1, 16, 19, 20, 26, 27, 28, and 35.
