Understanding the Discrete Element Method
Description
More details
Other editions
Additional editions
Persons
Content
About the Authors xv
Preface xvii
Acknowledgements xix
List of Abbreviations xxi
1 Mechanics 1
1.1 Degrees of freedom 1
1.1.1 Particle mechanics and constraints 1
1.1.2 From point particles to rigid bodies 3
1.1.3 More context and terminology 4
1.2 Dynamics of rectilinear degrees of freedom 5
1.3 Dynamics of angular degrees of freedom 6
1.3.1 Rotation in two dimensions 6
1.3.2 Moment of inertia 7
1.3.3 From two to three dimensions 9
1.3.4 Rotation matrix in three dimensions 12
1.3.5 Three-dimensional moments of inertia 13
1.3.6 Space-fixed and body-fixed coordinate systems and equations of motion 16
1.3.7 Problems with Euler angles 19
1.3.8 Rotations represented using complex numbers 20
1.3.9 Quaternions 21
1.3.10 Derivation of quaternion dynamics 27
1.4 The phase space 29
1.4.1 Qualitative discussion of the time dependence of linear oscillations 31
1.4.2 Resonance 34
1.4.3 The flow in phase space 35
1.5 Nonlinearities 39
1.5.1 Harmonic balance 40
1.5.2 Resonance in nonlinear systems 42
1.5.3 Higher harmonics and frequency mixing 44
1.5.4 The van der Pol oscillator 45
1.6 From higher harmonics to chaos 47
1.6.1 The bifurcation cascade 47
1.6.2 The nonlinear frictional oscillator and Poincaré maps 47
1.6.3 The route to chaos 51
1.6.4 Boundary conditions and many-particle systems 52
1.7 Stability and conservation laws 53
1.7.1 Stability in statics 54
1.7.2 Stability in dynamics 55
1.7.3 Stable axes of rotation around the principal axis 56
1.7.4 Noether's theorem and conservation laws 58
1.8 Further reading 61
Exercises 61
References 63
2 Numerical Integration of Ordinary Differential Equations 65
2.1 Fundamentals of numerical analysis 65
2.1.1 Floating point numbers 65
2.1.2 Big-O notation 67
2.1.3 Relative and absolute error 69
2.1.4 Truncation error 69
2.1.5 Local and global error 71
2.1.6 Stability 74
2.1.7 Stable integrators for unstable problems 74
2.2 Numerical analysis for ordinary differential equations 75
2.2.1 Variable notation and transformation of the order of a differential equation 75
2.2.2 Differences in the simulation of atoms and molecules, as compared to macroscopic particles 76
2.2.3 Truncation error for solutions of ordinary differential equations 76
2.2.4 Fundamental approaches 77
2.2.5 Explicit Euler method 77
2.2.6 Implicit Euler method 78
2.3 Runge-Kutta methods 79
2.3.1 Adaptive step-size control 79
2.3.2 Dense output and event location 81
2.3.3 Partitioned Runge-Kutta methods 82
2.4 Symplectic methods 82
2.4.1 The classical Verlet method 82
2.4.2 Velocity-Verlet methods 83
2.4.3 Higher-order velocity-Verlet methods 85
2.4.4 Pseudo-symplectic methods 88
2.4.5 Order, accuracy and energy conservation 88
2.4.6 Backward error analysis 89
2.4.7 Case study: the harmonic oscillator with and without viscous damping 90
2.5 Stiff problems 92
2.5.1 Evaluating computational costs 93
2.5.2 Stiff solutions and error as noise 94
2.5.3 Order reduction 94
2.6 Backward difference formulae 94
2.6.1 Implicit integrators of the predictor-corrector formulae 94
2.6.2 The corrector step 96
2.6.3 Multiple corrector steps 97
2.6.4 Program flow 98
2.6.5 Variable time-step and variable order 98
2.7 Other methods 98
2.7.1 Why not to use self-written or novel integrators 98
2.7.2 Stochastic differential equations 100
2.7.3 Extrapolation and high-order methods 100
2.7.4 Multi-rate integrators 101
2.7.5 Zero-order algorithms 101
2.8 Differential algebraic equations 103
2.8.1 The pendulum in Cartesian coordinates 103
2.8.2 Initial conditions 106
2.8.3 Drift and stabilization 107
2.9 Selecting an integrator 109
2.9.1 Performance and stability 109
2.9.2 Angular degrees of freedom 109
2.9.3 Force equilibrium 109
2.9.4 Exploring new fields 110
2.9.5 ODE solvers unsuitable for DEM simulations 110
2.10 Further reading 111
Exercises 113
References 125
3 Friction 129
3.1 Sliding Coulomb friction 129
3.1.1 A block on a slope 130
3.1.2 Static and dynamic friction coefficients 132
3.1.3 Apparent and actual contact area 134
3.1.4 Roughness and the friction coefficient 135
3.1.5 Adhesion and chemical bonding 136
3.2 Other contact geometries of Coulomb friction 136
3.2.1 Rolling friction 137
3.2.2 Pivoting friction 138
3.2.3 Sliding and rolling friction: the billiard problem 140
3.2.4 Sliding and rolling friction: cylinder on a slope 143
3.2.5 Pivoting and rolling friction 144
3.3 Exact implementation of friction 144
3.3.1 Establishing the difference between dynamic and static friction 145
3.3.2 Single-particle contact 148
3.3.3 Frictional linear chain 151
3.3.4 Higher dimensions 152
3.4 Modeling and regularizations 153
3.4.1 The Cundall-Strack model 153
3.4.2 Cundall-Strack friction in three dimensions 155
3.5 Unfortunate treatment of Coulomb friction in the literature 155
3.5.1 Insufficient models 156
3.5.2 Misunderstandings concerning surface roughness and friction 158
3.5.3 The Painlevé paradox 158
3.6 Further reading 158
Exercises 159
References 159
4 Phenomenology of Granular Materials 161
4.1 Phenomenology of grains 161
4.1.1 Interaction 161
4.1.2 Friction and dissipation 162
4.1.3 Length and time scales 162
4.1.4 Particle shape, and rolling and sliding 163
4.2 General phenomenology of granular agglomerates 164
4.2.1 Disorder 164
4.2.2 Heap formation 165
4.2.3 Tri-axial compression and shear band formation 166
4.2.4 Arching 168
4.2.5 Clogging 168
4.3 History effects in granular materials 168
4.3.1 Hysteresis 169
4.3.2 Reynolds dilatancy 170
4.3.3 Pressure distribution under heaps 171
4.4 Further reading 173
References 173
5 Condensed Matter and Solid State Physics 175
5.1 Structure and properties of matter 176
5.1.1 Crystal structures in two dimensions 176
5.1.2 Crystal structures in three dimensions 178
5.1.3 From the Wigner-Seitz cell to the Voronoi construction 180
5.1.4 Strength parameters of materials 182
5.1.5 Strength of granular assemblies 185
5.2 From wave numbers to the Fourier transform 186
5.2.1 Wave numbers and the reciprocal lattice 186
5.2.2 The Fourier transform in one dimension 188
5.2.3 Properties of the FFT 189
5.2.4 Other Fourier variables 193
5.2.5 The power spectrum 193
5.3 Waves and dispersion 194
5.3.1 Phase and group velocities 194
5.3.2 Phase and group velocities for particle systems 196
5.3.3 Numerical computation of the dispersion relation 199
5.3.4 Density of states 200
5.3.5 Dispersion relation for disordered systems 202
5.3.6 Solitons 204
5.4 Further reading 206
Exercises 206
References 210
6 Modeling and Simulation 213
6.1 Experiments, theory and simulation 213
6.2 Computability, observables and auxiliary quantities 214
6.3 Experiments, theories and the discrete element method 215
6.4 The discrete element method and other particle simulation methods 217
6.5 Other simulation methods for granular materials 218
6.5.1 Continuum mechanics 218
6.5.2 Lattice models 219
6.5.3 The Monte Carlo method 220
References 221
7 The Discrete Element Method in Two Dimensions 223
7.1 The discrete element method with soft particles 223
7.1.1 The bouncing ball as a prototype for the DEM approach 224
7.1.2 Using two different stiffness constants to model damping 227
7.1.3 Simulation of round DEM particles in one dimension 228
7.1.4 Simulation of round particles in two dimensions 228
7.2 Modeling of polygonal particles 229
7.2.1 Initializing two-dimensional particles 229
7.2.2 Computation of the mass, center of mass and moment of inertia 231
7.2.3 Non-convex polygons 237
7.3 Interaction 237
7.3.1 Shape-dependent elastic force law 238
7.3.2 Computation of the overlap geometry 240
7.3.3 Computation of other dynamic quantities 244
7.3.4 Damping 246
7.3.5 Cohesive forces 248
7.3.6 Penetrating particle overlaps 249
7.4 Initial and boundary conditions 250
7.4.1 Initializing convex polygons 250
7.4.2 General considerations 252
7.4.3 Initial positions 253
7.4.4 Boundary conditions 255
7.5 Neighborhood algorithms 257
7.5.1 Algorithms not recommended for elongated particles 258
7.5.2 'Sort and sweep' 263
7.6 Time integration 271
7.7 Program issues 272
7.7.1 Program restart 272
7.7.2 Program initialization 274
7.7.3 Program flow 274
7.7.4 Proposed stages for the development of programs 276
7.7.5 Modularization 278
7.8 Computing observables 280
7.8.1 Computing averages 280
7.8.2 Homogenization and spatial averages 281
7.8.3 Computing error bars 282
7.8.4 Autocorrelation functions 284
7.9 Further reading 285
Exercises 286
References 286
8 The Discrete Element Method in Three Dimensions 289
8.1 Generalization of the force law to three dimensions 289
8.1.1 The elastic force 290
8.1.2 Contact velocity and related forces 291
8.2 Initialization of particles and their properties 292
8.2.1 Basic concepts and data structures 292
8.2.2 Particle generation and geometry update 294
8.2.3 Decomposition of a polyhedron into tetrahedra 296
8.2.4 Volume, mass and center of mass 299
8.2.5 Moment of inertia 300
8.3 Overlap computation 301
8.3.1 Triangle intersection by using the point-direction form 301
8.3.2 Triangle intersection by using the point-normal form 305
8.3.3 Comparison of the two algorithms 309
8.3.4 Determination of inherited vertices 310
8.3.5 Determination of generated vertices 312
8.3.6 Determination of the faces of the overlap polyhedron 315
8.3.7 Determination of the contact area and normal 320
8.4 Optimization for vertex computation 322
8.4.1 Determination of neighboring features 323
8.4.2 Neighboring features for vertex computation 324
8.5 The neighborhood algorithm for polyhedra 325
8.5.1 'Sort and sweep' in three dimensions 325
8.5.2 Worst-case performance in three dimensions 326
8.5.3 Refinement of the contact list 327
8.6 Programming strategy for the polyhedral simulation 329
8.7 The effect of dimensionality and the choice of boundaries 332
8.7.1 Force networks and dimensionality 332
8.7.2 Quasi-two-dimensional geometries 332
8.7.3 Packings and sound propagation 333
8.8 Further reading 333
References 333
9 Alternative Modeling Approaches 335
9.1 Rigidly connected spheres 335
9.2 Elliptical shapes 336
9.2.1 Elliptical potentials 337
9.2.2 Overlap computation for ellipses 337
9.2.3 Newton-Raphson iteration 339
9.2.4 Ellipse intersection computed with generalized eigenvalues 340
9.2.5 Ellipsoids 344
9.2.6 Superquadrics 344
9.3 Composites of curves 345
9.3.1 Composites of arcs and cylinders 345
9.3.2 Spline curves 345
9.3.3 Level sets 347
9.4 Rigid particles 347
9.4.1 Collision dynamics ('event-driven method') 347
9.4.2 Contact mechanics 348
9.5 Discontinuous deformation analysis 349
9.6 Further reading 349
References 349
10 Running, Debugging and Optimizing Programs 353
10.1 Programming style 353
10.1.1 Literature 354
10.1.2 Choosing a programming language 355
10.1.3 Composite data types, strong typing and object orientation 356
10.1.4 Readability 356
10.1.5 Selecting variable names 357
10.1.6 Comments 359
10.1.7 Particle simulations versus solving ordinary differential equations 361
10.2 Hardware, memory and parallelism 362
10.2.1 Architecture and programming model 362
10.2.2 Memory hierarchy and cache 364
10.2.3 Multiprocessors, multi-core processors and shared memory 365
10.2.4 Peak performance and benchmarks 365
10.2.5 Amdahl's law, speed-up and efficiency 367
10.3 Program writing 369
10.3.1 Editors 370
10.3.2 Compilers 370
10.3.3 Makefiles 371
10.3.4 Writing and testing code 372
10.3.5 Debugging 377
10.4 Measuring load, time and profiles 378
10.4.1 The 'top' command 379
10.4.2 Xload 379
10.4.3 Performance monitor for multi-core processors 380
10.4.4 The 'time' command 380
10.4.5 The Unix profiler 383
10.4.6 Interactive profilers 383
10.5 Speeding up programs 383
10.5.1 Estimating the time consumption of operations 383
10.5.2 Compiler optimization options 384
10.5.3 Optimizations by hand 389
10.5.4 Avoiding unnecessary disk output 390
10.5.5 Look up or compute 390
10.5.6 Shared-memory parallelism and OpenMP 390
10.6 Further reading 391
Exercises 392
References 392
11 Beyond the Scope of This Book 395
11.1 Non-convex particles 395
11.2 Contact dynamics and friction 395
11.3 Impact mechanics 396
11.4 Fragmentation and fracturing 396
11.5 Coupling codes for particles and elastic continua 396
11.6 Coupling of particles and fluid 398
11.6.1 Basic considerations for the fluid simulation 398
11.6.2 Verification of the fluid code 398
11.6.3 Macroscopic simulations 399
11.6.4 Microscopic simulations 399
11.6.5 Particle approach for both particles and fluid 400
11.6.6 Mesh-based modeling approaches 402
11.7 The finite element method for contact problems 402
11.8 Long-range interactions 403
References 403
A MATLAB R¿ as Programming Language 407
A. 1 Getting started with MATLAB R¿ 407
A. 2 Data types and names 408
A. 3 Matrix functions and linear algebra 409
A. 4 Syntax and control structures 413
A. 5 Self-written functions 415
A. 6 Function overwriting and overloading 416
A. 7 Graphics 417
A. 8 Solving ordinary differential equations 418
A. 9 Pitfalls of using MATLAB R¿ 420
A. 10 Profiling and optimization 424
A. 11 Free alternatives to MATLAB R¿ 425
A. 12 Further reading 425
Exercises 426
References 430
B Geometry and Computational Geometry 433
B. 1 Trigonometric functions 433
B. 2 Points, line segments and vectors 435
B. 3 Products of vectors 436
B.3. 1 Inner product (scalar product, dot product) 436
B.3. 2 Orthogonality 437
B.. 3 Outer product 438
B.3. 4 Vector product 438
B.3. 5 Triple product 440
B. 4 Projections and rejections 441
B.4. 1 Projection of a vector onto another vector 441
B.4. 2 Rejection of one vector with respect to another vector 442
B. 5 Lines and planes 442
B.5. 1 Lines and line segments 442
B.5. 2 Planes 444
B. 6 Oriented quantities: distance, area, volume etc. 446
B. 7 Further reading 449
References 449
Index 451
1
Mechanics
We start with an outline of classical mechanics, to provide a framework for the discrete element method (DEM). While most of the material in this chapter can be found scattered in various books on mechanics, no text seems to be available which covers concisely the concepts needed for DEM simulation. This chapter is intended as a crash course in theoretical mechanics, with an emphasis on issues relevant to computer implementation and testing. We give a list of secondary literature that the reader may refer to for further details.
1.1 Degrees of freedom
Before discussing the dynamics of a mechanical system, we need to understand the nature of the variables in the system. There are independent variables on the one hand, usually called ‘degrees of freedom’, and then there are dependent variables which depend on the degrees of freedom, via algebraic relations or derivatives.
1.1.1 Particle mechanics and constraints
The concept of a ‘mass point’ means that we neglect the size of the mass and are interested only in its trajectory. The position of a single mass point moving along the Cartesian x-axis is described by the value of x, which corresponds to a single degree of freedom. A point moving in the xy-plane has two degrees of freedom, r2D = (x, y), and a point moving in three-dimensional real space will have three degrees of freedom, r3D = (x, y, z). Although we can describe the motion of a point in three-dimensional space by four ‘space-time coordinates’ using the tuple (x, y, z, t), in classical mechanics t is not considered a degree of freedom but rather a parameter, i.e. an independent variable which cannot be influenced.
Two mass points moving independently along the x-axis represent two degrees of freedom, r1 and r2 (here and in the following, we assume equal masses). If we ‘glue’ these two particles together at distance d = r1 − r2 as in Figure 1.1, one degree of freedom gets lost, and we are left with only a single degree of freedom; in this case we can use either of r1, r2 or the average (r1+r2)/2 to determine the position of both particles uniquely. This means that one constraint between two position variables eliminates one degree of freedom.
Figure 1.1 In two dimensions, the number of degrees of freedom ndof for 1, 2, 3 or 4 constrained particles with an increasing number of constraints introduced. Newly added constraints are in black; previous constraints are in gray.
In two dimensions, for two point particles at r1 = (x1, y1) and r2 = (x2, y2) we have four degrees of freedom, x1, y1, x2 and y2. If we again fix the distance between the particles at a constant distance d, so that
(1.1)we can choose any three variables from {x1, y1, x2, y2} and the fourth will then be determined from (1.1) by elementary geometry. Alternatively, we can introduce new variables, such as the position of the center of mass, (x, y) = (r1 + r2)/2 for particles of the same mass, the displacement (x, y) = (x2 — x1, y2 — y1) between the particles, and the angle θ that the line segment between the two particles makes with the x-axis. In any case, we end up with three independent variables to describe the positions of the two particles fully. This means that a single constraint (1.1) reduces the number of degrees of freedom, i.e. the number of independent variables in the system, by 1.
In three-dimensional space, for two particles at positions (x1, y1, z1) and (x2, y2, z2) as shown in Figure 1.2, a constraint
(1.2)Figure 1.2 In three dimensions, the number of degrees of freedom ndof for 1, 2, 3, 4 or 5 particles constrained so that the resulting cluster has no internal degrees of freedom. Newly added constraints are in black; previous constraints are in gray.
will again reduce the number of degrees of freedom by 1, so if we want to work with the center of mass
we need two angles, ϕ and θ say, to describe the orientation of the ‘rod’ in space. Rotation around the orientation of the rod is not a degree of freedom, as it does not change the positions of the two points. In principle, it does not matter how one defines the degrees of freedom, whether it is with six variables and one constraint (1.2), with three Cartesian coordinates for the center of mass and two angles, or with three Cartesian coordinates for one endpoint and two angles. In each case the number of degrees of freedom is the same, namely 5.
1.1.2 From point particles to rigid bodies
When we introduce one more point mass at (x3, y3, z3) to our set-up, we have 9 variables in total. If we connect this new point to both ends of our rod with the additional constraints
(1.3) (1.4)we get a triangle, as in the middle diagram of Figure 1.2. Again, we can give an alternative description of its position in space using the center of mass, and use three angles, ϕ, θ and ψ, to describe the orientation. So the formula
again holds. If we connect a fourth particle rigidly to the cluster of three particles so that it does not lie in the plane described by the other three, as shown in the fourth diagram from the left in Figure 1.2, then the three extra constraints exactly compensate for the additional three coordinates (x4, y4, z4) of the new particle. In fact, for four or more spatially connected particles, the total number of degrees of freedom is always 6. Note that the rigid body formed by the connected particles need not be three-dimensional; for example, although a triangle is a two-dimensional shape, if it can rotate in three dimensions, then it also has six degrees of freedom. Through the reasoning above, we have derived that an extended rigid body has six degrees of freedom, irrespective of its size. The angular degrees of freedom ϕ, θ, ψ are obtained from the rectilinear degrees of freedom (x1, y1, z1), (x2, y2, z2), ... of the particles upon introducing constraints of finite length between the particles.
‘Mathematically’ one can define a point particle as an object having ‘zero extension’ and a rigid body as one having ‘zero deformation’. A more pragmatic definition of a point particle is an object whose extent is much smaller than the distances that it covers in the processes under investigation; after all, the Earth is pretty extended, but the point-mass approach to describing its trajectory around the sun works rather well. Likewise, a rigid body is an object for which the deformations are much smaller than the scales that are of interest in the processes being investigated.
1.1.3 More context and terminology
In principle, a ‘continuum’ has infinitely many degrees of freedom; but in order to solve continuum problems with a computer, we have to first discretize the continuum to obtain a finite number of degrees of freedom. We could, for instance, decompose the continuum into representative mass points and model the elasticity by springs between the mass points. The deformation of a spring can be computed from the positions of the bodies, so the springs will not be degrees of freedom, while the coordinates of the mass points will be degrees of freedom. With a finite element discretization, we decompose the elastic continuum into a space-filling partition of elements for which elastic stress relations hold, and the degrees of freedom are the nodes of the elements. Depending on the choice of boundary conditions, there may be as many nodes as there are elements, or more; therefore, from the nodes one can calculate the center of mass of the elements, but not vice versa. Describing the physics via the motion of particles, for example of centers of mass, is called the ‘Lagrangian representation’. This approach is natural for particulate systems, so we will adopt it in this book. Formulating the physics for a reference system in which, e.g., density amplitudes change is called the ‘Eulerian representation’; this representation is preferable for many continuum problems. In a Lagrangian representation, velocities of mechanical bodies are not degrees of freedom: they can be obtained as the time derivatives of the positions on which they depend. On the other hand, when we simulate a fluid volume where velocities are assigned to the nodes of a finite element or finite difference approximation in ‘Eulerian representation’, it is the velocities that are the degrees of freedom.
In the previous two subsections, we introduced constraints as algebraic relations between positions, but we remark here that constraints (whose associated functions are usually denoted by g in formulae) can also be imposed on velocities. For a pendulum of length l swinging around the origin as in Figure 1.3(a), the constraint g(x, y) stating that the bob (whose diameter we will neglect) stays at constant distance from the origin is
(1.5)In § 2.8 we will discuss the numerical solution of a problem where, in addition to constraints on x and y, constraint relations for...
