tdem.html - 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.html (17744B)
---
1
2 <!DOCTYPE html>
3
4 <html xmlns="http://www.w3.org/1999/xhtml">
5 <head>
6 <meta charset="utf-8" />
7 <title>Discrete element method — sphere 2.15-beta documentation</title>
8 <link rel="stylesheet" href="_static/classic.css" type="text/css" />
9 <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
10
11 <script type="text/javascript" id="documentation_options" data-url_root="./" src="_static/documentation_options.js"></script>
12 <script type="text/javascript" src="_static/jquery.js"></script>
13 <script type="text/javascript" src="_static/underscore.js"></script>
14 <script type="text/javascript" src="_static/doctools.js"></script>
15 <script type="text/javascript" src="_static/language_data.js"></script>
16
17 <link rel="index" title="Index" href="genindex.html" />
18 <link rel="search" title="Search" href="search.html" />
19 <link rel="next" title="Fluid simulation and particle-fluid interaction" href="cfd.html" />
20 <link rel="prev" title="Introduction and Installation" href="introduction.html" />
21 </head><body>
22 <div class="related" role="navigation" aria-label="related navigation">
23 <h3>Navigation</h3>
24 <ul>
25 <li class="right" style="margin-right: 10px">
26 <a href="genindex.html" title="General Index"
27 accesskey="I">index</a></li>
28 <li class="right" >
29 <a href="py-modindex.html" title="Python Module Index"
30 >modules</a> |</li>
31 <li class="right" >
32 <a href="cfd.html" title="Fluid simulation and particle-fluid interaction"
33 accesskey="N">next</a> |</li>
34 <li class="right" >
35 <a href="introduction.html" title="Introduction and Installation"
36 accesskey="P">previous</a> |</li>
37 <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
38 </ul>
39 </div>
40
41 <div class="document">
42 <div class="documentwrapper">
43 <div class="bodywrapper">
44 <div class="body" role="main">
45
46 <div class="section" id="discrete-element-method">
47 <h1>Discrete element method<a class="headerlink" href="#discrete-element-method" title="Permalink to this headline">¶</a></h1>
48 <p>Granular material is a very common form of matter, both in nature and industry.
49 It can be defined as material consisting of interacting, discrete particles.
50 Common granular materials include gravels, sands and soils, ice bergs,
51 asteroids, powders, seeds, and other foods. Over 75% of the raw materials that
52 pass through industry are granular. This wide occurrence has driven the desire
53 to understand the fundamental mechanics of the material.</p>
54 <p>Contrary to other common materials such as gases, liquids and solids, a general
55 mathematical formulation of it’s behavior hasn’t yet been found. Granular
56 material can, however, display states that somewhat resemble gases, fluids and
57 solids.</p>
58 <p>The <a class="reference external" href="https://en.wikipedia.org/wiki/Discrete_element_method">Discrete Element Method</a> (DEM) is a numerical
59 method that can be used to
60 simulate the interaction of particles. Originally derived from
61 <a class="reference external" href="https://en.wikipedia.org/wiki/Molecular_dynamics">Molecular Dynamics</a>,
62 it simulates particles as separate entities, and calculates their positions,
63 velocities, and accelerations through time. See Cundall and Strack (1979) and
64 <a class="reference external" href="http://anders-dc.github.io/2013/10/16/the-discrete-element-method/">this blog post</a> for
65 general introduction to the DEM. The following sections will highlight the
66 DEM implementation in <code class="docutils literal notranslate"><span class="pre">sphere</span></code>. Some of the details are also described in
67 Damsgaard et al. 2013. In the used notation, a bold symbol denotes a
68 three-dimensional vector, and a dot denotes that the entity is a temporal
69 derivative.</p>
70 <div class="section" id="contact-search">
71 <h2>Contact search<a class="headerlink" href="#contact-search" title="Permalink to this headline">¶</a></h2>
72 <p>Homogeneous cubic grid.</p>
73 <div class="math">
74 <p><img src="_images/math/5c79b34fe936898cd55bfd2dfd9b85d58bd66fad.png" alt="\delta_n^{ij} = ||\boldsymbol{x}^i - \boldsymbol{x}^j|| - (r^i + r^j)"/></p>
75 </div><p>where <img class="math" src="_images/math/79a3d439d28652c547386f39b555d90d3aaf102d.png" alt="r"/> is the particle radius, and <img class="math" src="_images/math/fa44d7cdd2d3f90f81820cfe85818a142149d124.png" alt="\boldsymbol{x}"/> denotes the
76 positional vector of a particle, and <img class="math" src="_images/math/5aa339d4daf45a810dda332e3c80a0698e526e04.png" alt="i"/> and <img class="math" src="_images/math/e3fc28292267f066fee7718c64f4bbfece521f24.png" alt="j"/> denote the indexes
77 of two particles. Negative values of <img class="math" src="_images/math/a74a3c3e364aec16c0c7ebe002cad3214b489460.png" alt="\delta_n"/> denote that the particles
78 are overlapping.</p>
79 </div>
80 <div class="section" id="contact-interaction">
81 <h2>Contact interaction<a class="headerlink" href="#contact-interaction" title="Permalink to this headline">¶</a></h2>
82 <p>Now that the inter-particle contacts have been identified and characterized by
83 their overlap, the resulting forces from the interaction can be resolved. The
84 interaction is decomposed into normal and tangential components, relative to the
85 contact interface orientation. The normal vector to the contact interface is
86 found by:</p>
87 <div class="math">
88 <p><img src="_images/math/1ecc9f3eb9f630270af6c89a67d8e95b5c55014a.png" alt="\boldsymbol{n}^{ij} =
89 \frac{\boldsymbol{x}^i - \boldsymbol{x}^j}
90 {||\boldsymbol{x}^i - \boldsymbol{x}^j||}"/></p>
91 </div><p>The contact velocity <img class="math" src="_images/math/bdf3ad94aa38c155154e3286b79db56b1731afba.png" alt="\dot{\boldsymbol{\delta}}"/> is found by:</p>
92 <div class="math">
93 <p><img src="_images/math/bc86cfac20bf315a60bfd00614d9e5dd192c6dd5.png" alt="\dot{\boldsymbol{\delta}}^{ij} =
94 (\boldsymbol{x}^i - \boldsymbol{x}^j)
95 + (r^i + \frac{\delta_n^{ij}}{2})
96 (\boldsymbol{n}^{ij} \times \boldsymbol{\omega}^{i})
97 + (r^j + \frac{\delta_n^{ij}}{2})
98 (\boldsymbol{n}^{ij} \times \boldsymbol{\omega}^{j})"/></p>
99 </div><p>The contact velocity is decomposed into normal and tangential components,
100 relative to the contact interface. The normal component is:</p>
101 <div class="math">
102 <p><img src="_images/math/3c1524f2029c88a3186b303042ffef4ec5916196.png" alt="\dot{\delta}^{ij}_n =
103 -(\dot{\boldsymbol{\delta}}^{ij} \cdot \boldsymbol{n}^{ij})"/></p>
104 </div><p>and the tangential velocity component is found as:</p>
105 <div class="math">
106 <p><img src="_images/math/df18cfb7be27872e85cfca746a5b65aeec27187a.png" alt="\dot{\boldsymbol{\delta}}^{ij}_t =
107 \dot{\boldsymbol{\delta}}^{ij}
108 - \boldsymbol{n}^{ij}
109 (\boldsymbol{n}^{ij} \cdot \dot{\boldsymbol{\delta}}^{ij})"/></p>
110 </div><p>where <img class="math" src="_images/math/6087536f1e28c80ec705866fefe5d1760e121703.png" alt="\boldsymbol{\omega}"/> is the rotational velocity vector of a
111 particle. The total tangential displacement on the contact plane is found
112 incrementally:</p>
113 <div class="math">
114 <p><img src="_images/math/53f1505c53f0face6b8b68cb4b5919ec39f84592.png" alt="\boldsymbol{\delta}_{t,\text{uncorrected}}^{ij} =
115 \int_0^{t_c}
116 \dot{\boldsymbol{\delta}}^{ij}_t \Delta t"/></p>
117 </div><p>where <img class="math" src="_images/math/e497b0cd9de2666c43ccf2e42b7636cde6238dba.png" alt="t_c"/> is the duration of the contact and <img class="math" src="_images/math/b4ed9c2e208e08edeca8b1550ec0840acd090276.png" alt="\Delta t"/> is the
118 computational time step length. The tangential contact interface displacement is
119 set to zero when a contact pair no longer overlaps. At each time step, the value
120 of <img class="math" src="_images/math/c46baa2a60c45bb8746615b3d8b545cf4a6523f6.png" alt="\boldsymbol{\delta}_t"/> is corrected for rotation of the contact
121 interface:</p>
122 <div class="math">
123 <p><img src="_images/math/a5d04fe388b8b273356584389cb4a1358bb90ad1.png" alt="\boldsymbol{\delta}_t^{ij} = \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij}
124 - (\boldsymbol{n}
125 (\boldsymbol{n} \cdot \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij})"/></p>
126 </div><p>With all the geometrical and kinetic components determined, the resulting forces
127 of the particle interaction can be determined using a contact model. <code class="docutils literal notranslate"><span class="pre">sphere</span></code>
128 features only one contact model in the normal direction to the contact; the
129 linear-elastic-viscous (<em>Hookean</em> with viscous damping, or <em>Kelvin-Voigt</em>)
130 contact model. The resulting force in the normal direction of the contact
131 interface on particle <img class="math" src="_images/math/5aa339d4daf45a810dda332e3c80a0698e526e04.png" alt="i"/> is:</p>
132 <div class="math">
133 <p><img src="_images/math/c831f45c0dbc852540dc9d70bed0d2c034929af8.png" alt="\boldsymbol{f}_n^{ij} = \left(
134 -k_n \delta_n^{ij} -\gamma_n \dot{\delta_n}^{ij}
135 \right) \boldsymbol{n}^{ij}"/></p>
136 </div><p>The parameter <img class="math" src="_images/math/c713414d12f194f3fab98645df441d23d54164ec.png" alt="k_n"/> is the defined <a class="reference external" href="https://en.wikipedia.org/wiki/Hooke's_law">spring coefficient</a> in the normal direction of the
137 contact interface, and <img class="math" src="_images/math/9cc677554b62a84edc8937d747d04d117f715112.png" alt="\gamma_n"/> is the defined contact interface
138 viscosity, also in the normal direction. The loss of energy in this interaction
139 due to the viscous component is for particle <img class="math" src="_images/math/5aa339d4daf45a810dda332e3c80a0698e526e04.png" alt="i"/> calculated as:</p>
140 <div class="math">
141 <p><img src="_images/math/a928c6f408661dd8d6860f7c924635e2ae9bce96.png" alt="\dot{e}^i_v = \gamma_n (\dot{\delta}^{ij}_n)^2"/></p>
142 </div><p>The tangential force is determined by either a viscous-frictional contact model,
143 or a elastic-viscous-frictional contact model. The former contact model is very
144 computationally efficient, but somewhat inaccurate relative to the mechanics of
145 real materials. The latter contact model is therefore the default, even though
146 it results in longer computational times. The tangential force in the
147 visco-frictional contact model:</p>
148 <div class="math">
149 <p><img src="_images/math/ad8a19ed2553a83794270f0a35f29f66b15c0bb7.png" alt="\boldsymbol{f}_t^{ij} = -\gamma_t \dot{\boldsymbol{\delta}_t}^{ij}"/></p>
150 </div><p><img class="math" src="_images/math/9cc677554b62a84edc8937d747d04d117f715112.png" alt="\gamma_n"/> is the defined contact interface viscosity in the tangential
151 direction. The tangential displacement along the contact interface
152 (<img class="math" src="_images/math/c46baa2a60c45bb8746615b3d8b545cf4a6523f6.png" alt="\boldsymbol{\delta}_t"/>) is not calculated and stored for this contact
153 model. The tangential force in the more realistic elastic-viscous-frictional
154 contact model:</p>
155 <div class="math">
156 <p><img src="_images/math/a3471ecd78ccfe91e0db3d5a607624380afb04d1.png" alt="\boldsymbol{f}_t^{ij} =
157 -k_t \boldsymbol{\delta}_t^{ij} -\gamma_t \dot{\boldsymbol{\delta}_t}^{ij}"/></p>
158 </div><p>The parameter <img class="math" src="_images/math/c713414d12f194f3fab98645df441d23d54164ec.png" alt="k_n"/> is the defined spring coefficient in the tangential
159 direction of the contact interface. Note that the tangential force is only
160 found if the tangential displacement (<img class="math" src="_images/math/0362550f715837c65115e59846809a42dfba49d1.png" alt="\delta_t"/>) or the tangential
161 velocity (<img class="math" src="_images/math/c17d332a01b1635dd5d3b5c297e9d44c5b742c04.png" alt="\dot{\delta}_t"/>) is non-zero, in order to avoid division by
162 zero. Otherwise it is defined as being <img class="math" src="_images/math/dc895ec2bd6380cb07f0f5d723025eee92f7c09c.png" alt="[0,0,0]"/>.</p>
163 <p>For both types of contact model, the tangential force is limited by the Coulomb
164 criterion of static and dynamic friction:</p>
165 <div class="math">
166 <p><img src="_images/math/4397f45ed8f18be640b692425a8419f5cb72ada3.png" alt="||\boldsymbol{f}^{ij}_t|| \leq
167 \begin{cases}
168 \mu_s ||\boldsymbol{f}^{ij}_n|| &
169 \text{if} \quad ||\boldsymbol{f}_t^{ij}|| = 0 \\
170 \mu_d ||\boldsymbol{f}^{ij}_n|| &
171 \text{if} \quad ||\boldsymbol{f}_t^{ij}|| > 0
172 \end{cases}"/></p>
173 </div><p>If the elastic-viscous-frictional contact model is used and the Coulomb limit is
174 reached, the tangential displacement along the contact interface is limited to
175 this value:</p>
176 <div class="math">
177 <p><img src="_images/math/c19adb8434d3b9c150be8c8ab3939958c8edaba9.png" alt="\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)"/></p>
182 </div><p>If the tangential force reaches the Coulomb limit, the energy lost due to
183 frictional dissipation is calculated as:</p>
184 <div class="math">
185 <p><img src="_images/math/2fc641eba6a1d6e4cd948e7e48341561e69fbf10.png" alt="\dot{e}^i_s = \frac{||\boldsymbol{f}^{ij}_t
186 \dot{\boldsymbol{\delta}}_t^{ij} \Delta t||}{\Delta t}"/></p>
187 </div><p>The loss of energy by viscous dissipation in the tangential direction is not
188 found.</p>
189 </div>
190 <div class="section" id="temporal-integration">
191 <h2>Temporal integration<a class="headerlink" href="#temporal-integration" title="Permalink to this headline">¶</a></h2>
192 <p>In the DEM, the time is discretized into small steps (<img class="math" src="_images/math/b4ed9c2e208e08edeca8b1550ec0840acd090276.png" alt="\Delta t"/>). For each time
193 step, the entire network of contacts is resolved, and the resulting forces and
194 torques for each particle are found. With these values at hand, the new
195 linear and rotational accelerations can be found using
196 <a class="reference external" href="https://en.wikipedia.org/wiki/Newton%27s_laws_of_motion">Newton’s second law</a>
197 of the motion of solid bodies. If a particle with mass <img class="math" src="_images/math/e9bc7da808d33a16a8347f27a519bd067186aa66.png" alt="m"/> at a point in time
198 experiences a sum of forces denoted <img class="math" src="_images/math/fb7d5a70e84450d796cf4499cf923fa360b0d35b.png" alt="\boldsymbol{F}"/>, the resultant acceleration
199 (<img class="math" src="_images/math/be9a4fc94d921be9480d891b5c38b77e187b630a.png" alt="\boldsymbol{a}"/>) can be found by rearranging Newton’s second law:</p>
200 <div class="math">
201 <p><img src="_images/math/fb9c8157d9c1c87673619af0faf8dc55f256353a.png" alt="\boldsymbol{F} = m \boldsymbol{a} \Rightarrow \boldsymbol{a} = \frac{\boldsymbol{F}}{m}"/></p>
202 </div><p>The new velocity and position is found by integrating the above equation
203 with regards to time. The simplest integration scheme in this regard is the
204 <a class="reference external" href="https://en.wikipedia.org/wiki/Euler_method">Euler method</a>:</p>
205 <div class="math">
206 <p><img src="_images/math/3daf2a856e4e741d27444b532b2fff5a9d58b1db.png" alt="\boldsymbol{v} = \boldsymbol{v}_{old} + \boldsymbol{a} \Delta t"/></p>
207 </div><div class="math">
208 <p><img src="_images/math/9d7283875d4be9fca20b3d8c2289a89763a3b38e.png" alt="\boldsymbol{p} = \boldsymbol{p}_{old} + \boldsymbol{v} \Delta t"/></p>
209 </div></div>
210 </div>
211
212
213 </div>
214 </div>
215 </div>
216 <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
217 <div class="sphinxsidebarwrapper">
218 <h3><a href="index.html">Table of Contents</a></h3>
219 <ul>
220 <li><a class="reference internal" href="#">Discrete element method</a><ul>
221 <li><a class="reference internal" href="#contact-search">Contact search</a></li>
222 <li><a class="reference internal" href="#contact-interaction">Contact interaction</a></li>
223 <li><a class="reference internal" href="#temporal-integration">Temporal integration</a></li>
224 </ul>
225 </li>
226 </ul>
227
228 <h4>Previous topic</h4>
229 <p class="topless"><a href="introduction.html"
230 title="previous chapter">Introduction and Installation</a></p>
231 <h4>Next topic</h4>
232 <p class="topless"><a href="cfd.html"
233 title="next chapter">Fluid simulation and particle-fluid interaction</a></p>
234 <div role="note" aria-label="source link">
235 <h3>This Page</h3>
236 <ul class="this-page-menu">
237 <li><a href="_sources/dem.rst.txt"
238 rel="nofollow">Show Source</a></li>
239 </ul>
240 </div>
241 <div id="searchbox" style="display: none" role="search">
242 <h3 id="searchlabel">Quick search</h3>
243 <div class="searchformwrapper">
244 <form class="search" action="search.html" method="get">
245 <input type="text" name="q" aria-labelledby="searchlabel" />
246 <input type="submit" value="Go" />
247 </form>
248 </div>
249 </div>
250 <script type="text/javascript">$('#searchbox').show(0);</script>
251 </div>
252 </div>
253 <div class="clearer"></div>
254 </div>
255 <div class="related" role="navigation" aria-label="related navigation">
256 <h3>Navigation</h3>
257 <ul>
258 <li class="right" style="margin-right: 10px">
259 <a href="genindex.html" title="General Index"
260 >index</a></li>
261 <li class="right" >
262 <a href="py-modindex.html" title="Python Module Index"
263 >modules</a> |</li>
264 <li class="right" >
265 <a href="cfd.html" title="Fluid simulation and particle-fluid interaction"
266 >next</a> |</li>
267 <li class="right" >
268 <a href="introduction.html" title="Introduction and Installation"
269 >previous</a> |</li>
270 <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
271 </ul>
272 </div>
273 <div class="footer" role="contentinfo">
274 © Copyright 2014, Anders Damsgaard.
275 Created using <a href="http://sphinx-doc.org/">Sphinx</a> 2.2.0.
276 </div>
277 </body>
278