tdem.rst - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tdem.rst (9356B)
---
1 Discrete element method
2 =======================
3 Granular material is a very common form of matter, both in nature and industry.
4 It can be defined as material consisting of interacting, discrete particles.
5 Common granular materials include gravels, sands and soils, ice bergs,
6 asteroids, powders, seeds, and other foods. Over 75% of the raw materials that
7 pass through industry are granular. This wide occurrence has driven the desire
8 to understand the fundamental mechanics of the material.
9
10 Contrary to other common materials such as gases, liquids and solids, a general
11 mathematical formulation of it's behavior hasn't yet been found. Granular
12 material can, however, display states that somewhat resemble gases, fluids and
13 solids.
14
15 .. The discrete element method (or distinct element method) was initially
16 formulated by Cundall and Strack (1979). It simulates the physical behavior and
17 interaction of discrete, unbreakable particles, with their own mass and inertia,
18 under the influence of e.g. gravity and boundary conditions such as moving
19 walls. By discretizing time into small time steps, explicit integration of
20 Newton's second law of motion is used to predict the new position and kinematic
21 values for each particle from the previous sums of forces. This Lagrangian
22 approach is ideal for simulating discontinuous materials, such as granular
23 matter.
24 The complexity of the computations is kept low by representing the particles as
25 spheres, which keeps contact-searching algorithms simple.
26
27 The `Discrete Element Method
28 <https://en.wikipedia.org/wiki/Discrete_element_method>`_ (DEM) is a numerical
29 method that can be used to
30 simulate the interaction of particles. Originally derived from
31 `Molecular Dynamics <https://en.wikipedia.org/wiki/Molecular_dynamics>`_,
32 it simulates particles as separate entities, and calculates their positions,
33 velocities, and accelerations through time. See Cundall and Strack (1979) and
34 `this blog post
35 <http://anders-dc.github.io/2013/10/16/the-discrete-element-method/>`_ for
36 general introduction to the DEM. The following sections will highlight the
37 DEM implementation in ``sphere``. Some of the details are also described in
38 Damsgaard et al. 2013. In the used notation, a bold symbol denotes a
39 three-dimensional vector, and a dot denotes that the entity is a temporal
40 derivative.
41
42 Contact search
43 --------------
44 Homogeneous cubic grid.
45
46 .. math::
47 \delta_n^{ij} = ||\boldsymbol{x}^i - \boldsymbol{x}^j|| - (r^i + r^j)
48
49 where :math:`r` is the particle radius, and :math:`\boldsymbol{x}` denotes the
50 positional vector of a particle, and :math:`i` and :math:`j` denote the indexes
51 of two particles. Negative values of :math:`\delta_n` denote that the particles
52 are overlapping.
53
54
55 Contact interaction
56 -------------------
57 Now that the inter-particle contacts have been identified and characterized by
58 their overlap, the resulting forces from the interaction can be resolved. The
59 interaction is decomposed into normal and tangential components, relative to the
60 contact interface orientation. The normal vector to the contact interface is
61 found by:
62
63 .. math::
64 \boldsymbol{n}^{ij} =
65 \frac{\boldsymbol{x}^i - \boldsymbol{x}^j}
66 {||\boldsymbol{x}^i - \boldsymbol{x}^j||}
67
68 The contact velocity :math:`\dot{\boldsymbol{\delta}}` is found by:
69
70 .. math::
71 \dot{\boldsymbol{\delta}}^{ij} =
72 (\boldsymbol{x}^i - \boldsymbol{x}^j)
73 + (r^i + \frac{\delta_n^{ij}}{2})
74 (\boldsymbol{n}^{ij} \times \boldsymbol{\omega}^{i})
75 + (r^j + \frac{\delta_n^{ij}}{2})
76 (\boldsymbol{n}^{ij} \times \boldsymbol{\omega}^{j})
77
78 The contact velocity is decomposed into normal and tangential components,
79 relative to the contact interface. The normal component is:
80
81 .. math::
82 \dot{\delta}^{ij}_n =
83 -(\dot{\boldsymbol{\delta}}^{ij} \cdot \boldsymbol{n}^{ij})
84
85 and the tangential velocity component is found as:
86
87 .. math::
88 \dot{\boldsymbol{\delta}}^{ij}_t =
89 \dot{\boldsymbol{\delta}}^{ij}
90 - \boldsymbol{n}^{ij}
91 (\boldsymbol{n}^{ij} \cdot \dot{\boldsymbol{\delta}}^{ij})
92
93 where :math:`\boldsymbol{\omega}` is the rotational velocity vector of a
94 particle. The total tangential displacement on the contact plane is found
95 incrementally:
96
97 .. math::
98 \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij} =
99 \int_0^{t_c}
100 \dot{\boldsymbol{\delta}}^{ij}_t \Delta t
101
102 where :math:`t_c` is the duration of the contact and :math:`\Delta t` is the
103 computational time step length. The tangential contact interface displacement is
104 set to zero when a contact pair no longer overlaps. At each time step, the value
105 of :math:`\boldsymbol{\delta}_t` is corrected for rotation of the contact
106 interface:
107
108 .. math::
109 \boldsymbol{\delta}_t^{ij} = \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij}
110 - (\boldsymbol{n}
111 (\boldsymbol{n} \cdot \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij})
112
113 With all the geometrical and kinetic components determined, the resulting forces
114 of the particle interaction can be determined using a contact model. ``sphere``
115 features only one contact model in the normal direction to the contact; the
116 linear-elastic-viscous (*Hookean* with viscous damping, or *Kelvin-Voigt*)
117 contact model. The resulting force in the normal direction of the contact
118 interface on particle :math:`i` is:
119
120 .. math::
121 \boldsymbol{f}_n^{ij} = \left(
122 -k_n \delta_n^{ij} -\gamma_n \dot{\delta_n}^{ij}
123 \right) \boldsymbol{n}^{ij}
124
125 The parameter :math:`k_n` is the defined `spring coefficient
126 <https://en.wikipedia.org/wiki/Hooke's_law>`_ in the normal direction of the
127 contact interface, and :math:`\gamma_n` is the defined contact interface
128 viscosity, also in the normal direction. The loss of energy in this interaction
129 due to the viscous component is for particle :math:`i` calculated as:
130
131 .. math::
132 \dot{e}^i_v = \gamma_n (\dot{\delta}^{ij}_n)^2
133
134 The tangential force is determined by either a viscous-frictional contact model,
135 or a elastic-viscous-frictional contact model. The former contact model is very
136 computationally efficient, but somewhat inaccurate relative to the mechanics of
137 real materials. The latter contact model is therefore the default, even though
138 it results in longer computational times. The tangential force in the
139 visco-frictional contact model:
140
141 .. math::
142 \boldsymbol{f}_t^{ij} = -\gamma_t \dot{\boldsymbol{\delta}_t}^{ij}
143
144 :math:`\gamma_n` is the defined contact interface viscosity in the tangential
145 direction. The tangential displacement along the contact interface
146 (:math:`\boldsymbol{\delta}_t`) is not calculated and stored for this contact
147 model. The tangential force in the more realistic elastic-viscous-frictional
148 contact model:
149
150 .. math::
151 \boldsymbol{f}_t^{ij} =
152 -k_t \boldsymbol{\delta}_t^{ij} -\gamma_t \dot{\boldsymbol{\delta}_t}^{ij}
153
154 The parameter :math:`k_n` is the defined spring coefficient in the tangential
155 direction of the contact interface. Note that the tangential force is only
156 found if the tangential displacement (:math:`\delta_t`) or the tangential
157 velocity (:math:`\dot{\delta}_t`) is non-zero, in order to avoid division by
158 zero. Otherwise it is defined as being :math:`[0,0,0]`.
159
160 For both types of contact model, the tangential force is limited by the Coulomb
161 criterion of static and dynamic friction:
162
163 .. math::
164 ||\boldsymbol{f}^{ij}_t|| \leq
165 \begin{cases}
166 \mu_s ||\boldsymbol{f}^{ij}_n|| &
167 \text{if} \quad ||\boldsymbol{f}_t^{ij}|| = 0 \\
168 \mu_d ||\boldsymbol{f}^{ij}_n|| &
169 \text{if} \quad ||\boldsymbol{f}_t^{ij}|| > 0
170 \end{cases}
171
172 If the elastic-viscous-frictional contact model is used and the Coulomb limit is
173 reached, the tangential displacement along the contact interface is limited to
174 this value:
175
176 .. math::
177 \boldsymbol{\delta}_t^{ij} =
178 \frac{1}{k_t} \left(
179 \mu_d ||\boldsymbol{f}_n^{ij}||
180 \frac{\boldsymbol{f}^{ij}_t}{||\boldsymbol{f}^{ij}_t||}
181 + \gamma_t \dot{\boldsymbol{\delta}}_t^{ij} \right)
182
183 If the tangential force reaches the Coulomb limit, the energy lost due to
184 frictional dissipation is calculated as:
185
186 .. math::
187 \dot{e}^i_s = \frac{||\boldsymbol{f}^{ij}_t
188 \dot{\boldsymbol{\delta}}_t^{ij} \Delta t||}{\Delta t}
189
190 The loss of energy by viscous dissipation in the tangential direction is not
191 found.
192
193
194 Temporal integration
195 --------------------
196 In the DEM, the time is discretized into small steps (:math:`\Delta t`). For each time
197 step, the entire network of contacts is resolved, and the resulting forces and
198 torques for each particle are found. With these values at hand, the new
199 linear and rotational accelerations can be found using
200 `Newton's second law <https://en.wikipedia.org/wiki/Newton%27s_laws_of_motion>`_
201 of the motion of solid bodies. If a particle with mass :math:`m` at a point in time
202 experiences a sum of forces denoted :math:`\boldsymbol{F}`, the resultant acceleration
203 (:math:`\boldsymbol{a}`) can be found by rearranging Newton's second law:
204
205 .. math::
206 \boldsymbol{F} = m \boldsymbol{a} \Rightarrow \boldsymbol{a} = \frac{\boldsymbol{F}}{m}
207
208 The new velocity and position is found by integrating the above equation
209 with regards to time. The simplest integration scheme in this regard is the
210 `Euler method <https://en.wikipedia.org/wiki/Euler_method>`_:
211
212 .. math::
213 \boldsymbol{v} = \boldsymbol{v}_{old} + \boldsymbol{a} \Delta t
214
215 .. math::
216 \boldsymbol{p} = \boldsymbol{p}_{old} + \boldsymbol{v} \Delta t
217