tupdate documentation - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
DIR commit 81d7b44047babc963125026004031b612e6de13d
DIR parent 0d5bba5c16a5ad5a76f908a2f082722e3d869c34
HTML Author: Anders Damsgaard <andersd@riseup.net>
Date: Thu, 7 Sep 2017 17:25:14 -0700
update documentation
Diffstat:
M doc/html/.buildinfo | 2 +-
M doc/html/_sources/introduction.txt | 13 +++++--------
M doc/html/cfd.html | 196 ++++++++++++++++----------------
M doc/html/dem.html | 86 +++++++++++++++---------------
M doc/html/genindex.html | 10 +++++-----
M doc/html/index.html | 10 +++++-----
M doc/html/introduction.html | 28 +++++++++++++---------------
M doc/html/objects.inv | 0
M doc/html/py-modindex.html | 10 +++++-----
M doc/html/python_api.html | 22 ++++++++++++++++------
M doc/html/search.html | 10 +++++-----
M doc/html/searchindex.js | 4 ++--
M doc/html/sphere_internals.html | 14 +++++++-------
M doc/pdf/sphere.pdf | 0
M doc/sphinx/conf.py | 9 ++++++---
M doc/sphinx/introduction.rst | 13 +++++--------
M src/main.cpp | 4 ++--
17 files changed, 218 insertions(+), 213 deletions(-)
---
DIR diff --git a/doc/html/.buildinfo b/doc/html/.buildinfo
t@@ -1,4 +1,4 @@
# Sphinx build info version 1
# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: fe70a33065653dc9f13fbefebe7ac2f6
+config: f6a3bc901e2eb5ac86d94b7ac7732d4c
tags: 645f666f9bcd5a90fca523b33c5a78b7
DIR diff --git a/doc/html/_sources/introduction.txt b/doc/html/_sources/introduction.txt
t@@ -132,9 +132,11 @@ from the root directory::
$ cmake . && make
-NOTE: If your system does not have a GCC compiler (e.g. GCC-5 for CUDA 8)
-compatible with the installed CUDA version, try using ``clang-3.8`` instead::
+NOTE: If your system does not have a GCC compiler compatible with the installed
+CUDA version (e.g. GCC-5 for CUDA 8), you will see errors at the linker stage.
+In that case, try using ``clang-3.8`` as the C and C++ compiler instead::
+ $ rm -rf CMakeCache.txt CMakeFiles/
$ export CC=$(which clang-3.8) && export CXX=$(which clang++-3.8) && cmake . && make
After a successfull installation, the ``sphere`` executable will be located
t@@ -157,10 +159,6 @@ The output should look similar to this:
.. program-output:: ../../sphere --version
-The build can be verified by running a number of automated tests::
-
- $ make test
-
The documentation can be read in the `reStructuredText
<http://docutils.sourceforge.net/docs/ref/rst/restructuredtext.html>`_-format in
the ``doc/sphinx/`` folder, or in the HTML or PDF formats in the folders
t@@ -193,10 +191,9 @@ After compiling the ``sphere`` binary, the procedure of a creating and handling
a simulation is typically arranged in the following order:
* Setup of particle assemblage, physical properties and conditions using the
- Python API.
+ Python API (``python/sphere.py``).
* Execution of ``sphere`` software, which simulates the particle behavior as a
function of time, as a result of the conditions initially specified in the
input file.
* Inspection, analysis, interpretation and visualization of ``sphere`` output
in Python, and/or scene rendering using the built-in ray tracer.
-
DIR diff --git a/doc/html/cfd.html b/doc/html/cfd.html
t@@ -6,7 +6,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>Fluid simulation and particle-fluid interaction — sphere 1.00-alpha documentation</title>
+ <title>Fluid simulation and particle-fluid interaction — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -14,7 +14,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -25,7 +25,7 @@
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="index.html" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="index.html" />
<link rel="next" title="Python API" href="python_api.html" />
<link rel="prev" title="Discrete element method" href="dem.html" />
</head>
t@@ -45,7 +45,7 @@
<li class="right" >
<a href="dem.html" title="Discrete element method"
accesskey="P">previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -68,31 +68,31 @@ solution procedure and the numerical implementation.</p>
<p>Following the outline presented by <a class="reference external" href="http://www.cimec.org.ar/ojs/index.php/mc/article/view/486/464">Limache and Idelsohn (2006)</a>, the
continuity equation for an incompressible fluid material is given by:</p>
<div class="math">
-<p><img src="_images/math/341e857d6da62b73593a73dea371b22a955fd8e9.png" alt="\nabla \cdot \boldsymbol{v} = 0"/></p>
+<p><img src="_images/math/12f7f7910a11a5fc522425b33139287524c6e643.png" alt="\nabla \cdot \boldsymbol{v} = 0"/></p>
</div><p>and the momentum equation:</p>
<div class="math">
-<p><img src="_images/math/f4c62b9c99051d7d2ca2b222af4b3394d02e1b1d.png" alt="\rho \frac{\partial \boldsymbol{v}}{\partial t}
+<p><img src="_images/math/8b081782b6d73a048e968995b8605086e2c9ff1f.png" alt="\rho \frac{\partial \boldsymbol{v}}{\partial t}
+ \rho (\boldsymbol{v} \cdot \nabla \boldsymbol{v})
= \nabla \cdot \boldsymbol{\sigma}
- \boldsymbol{f}^i
+ \rho \boldsymbol{g}"/></p>
-</div><p>Here, <img class="math" src="_images/math/c8e734bf3818506e254ebfa35047a7957546c2f8.png" alt="\boldsymbol{v}"/> is the fluid velocity, <img class="math" src="_images/math/f574498915fa9e02eeb5141c24835d077eba3e75.png" alt="\rho"/> is the
-fluid density, <img class="math" src="_images/math/b4104067dcc504da8eada0697129bcae9f7cf9b6.png" alt="\boldsymbol{\sigma}"/> is the <a class="reference external" href="https://en.wikipedia.org/wiki/Cauchy_stress_tensor">Cauchy stress tensor</a>,
-<img class="math" src="_images/math/3ac346ec0a0efeed60827b2d73c5776cf4605902.png" alt="\boldsymbol{f}^i"/> is the particle-fluid interaction vector and
-<img class="math" src="_images/math/af2f4149a102bb8a3d74e9afbfd0d21b4c2f7084.png" alt="\boldsymbol{g}"/> is the gravitational acceleration. For incompressible
+</div><p>Here, <img class="math" src="_images/math/e7b48cd405b36440e40f970237c832a836ef198c.png" alt="\boldsymbol{v}"/> is the fluid velocity, <img class="math" src="_images/math/9a51ab9a0b521705e1e8762fac6bdd6f11771758.png" alt="\rho"/> is the
+fluid density, <img class="math" src="_images/math/2c6c8b6b1866df63a964f3c153a8db9eb4aec6d7.png" alt="\boldsymbol{\sigma}"/> is the <a class="reference external" href="https://en.wikipedia.org/wiki/Cauchy_stress_tensor">Cauchy stress tensor</a>,
+<img class="math" src="_images/math/af5ed4f2edd188dadef68b3626d9fe7328f44848.png" alt="\boldsymbol{f}^i"/> is the particle-fluid interaction vector and
+<img class="math" src="_images/math/2be557ceeacba6f01feedb2ed907af87442104a1.png" alt="\boldsymbol{g}"/> is the gravitational acceleration. For incompressible
Newtonian fluids, the Cauchy stress is given by:</p>
<div class="math">
-<p><img src="_images/math/9cb510f5a23b8a2cf167578e3197615a3cab7510.png" alt="\boldsymbol{\sigma} = -p \boldsymbol{I} + \boldsymbol{\tau}"/></p>
-</div><p><img class="math" src="_images/math/3eca8557203e86160952e1c0f735f7417f3285b1.png" alt="p"/> is the fluid pressure, <img class="math" src="_images/math/cfe0e96f004950d55bc7f58233eedae9f8b6abc1.png" alt="\boldsymbol{I}"/> is the identity
-tensor, and <img class="math" src="_images/math/e29be0a3b86507644c82abbc0adfcb22f0edb5e3.png" alt="\boldsymbol{\tau}"/> is the deviatoric stress tensor, given
+<p><img src="_images/math/778cdde6b1b06e6e1e0c059db34e1821a0744b5b.png" alt="\boldsymbol{\sigma} = -p \boldsymbol{I} + \boldsymbol{\tau}"/></p>
+</div><p><img class="math" src="_images/math/27d463da4622be5b3ef1d4176ced7d7a323c6425.png" alt="p"/> is the fluid pressure, <img class="math" src="_images/math/a46643c766fffafb430d4298b0f4a91d7ea54ac0.png" alt="\boldsymbol{I}"/> is the identity
+tensor, and <img class="math" src="_images/math/ecff587acb289da44784794a8cbe92f3ef24951e.png" alt="\boldsymbol{\tau}"/> is the deviatoric stress tensor, given
by:</p>
<div class="math">
-<p><img src="_images/math/2ac5fed54f92fbc94d2fefb2aa5f17a1e326282e.png" alt="\boldsymbol{\tau} =
+<p><img src="_images/math/048aa5b140cc86095a635258ba222f4ce562c46c.png" alt="\boldsymbol{\tau} =
\mu_f \nabla \boldsymbol{v}
+ \mu_f (\nabla \boldsymbol{v})^T"/></p>
</div><p>By using the following vector identities:</p>
<div class="math">
-<p><img src="_images/math/ffe6c17031ff1e2333c5bdc57773967dda2bebcc.png" alt="\nabla \cdot (p \boldsymbol{I}) = \nabla p
+<p><img src="_images/math/54c6ba63db4b565c0d31fd9e64caae2e0990706b.png" alt="\nabla \cdot (p \boldsymbol{I}) = \nabla p
\nabla \cdot (\nabla \boldsymbol{v}) = \nabla^2 \boldsymbol{v}
t@@ -101,29 +101,29 @@ by:</p>
</div><p>the deviatoric component of the Cauchy stress tensor simplifies to the
following, assuming that spatial variations in the viscosity can be neglected:</p>
<div class="math">
-<p><img src="_images/math/f220cc08054c66b15757b7073edbbfe33c2736e6.png" alt="= -\nabla p
+<p><img src="_images/math/3b48976b9058ec32db3d413cfccaede0cd1bf31b.png" alt="= -\nabla p
+ \mu_f \nabla^2 \boldsymbol{v}"/></p>
</div><p>Since we are dealing with fluid flow in a porous medium, additional terms are
introduced to the equations for conservation of mass and momentum. In the
following, the equations are derived for the first spatial component. The
solution for the other components is trivial.</p>
<p>The porosity value (in the saturated porous medium the volumetric fraction of
-the fluid phase) denoted <img class="math" src="_images/math/10e009bdb83f96c5f47c58b34d5d4b12ef268d5b.png" alt="\phi"/> is incorporated in the continuity and
+the fluid phase) denoted <img class="math" src="_images/math/c2f31c22645274c375eff7920cfdfdc18d60341f.png" alt="\phi"/> is incorporated in the continuity and
momentum equations. The continuity equation becomes:</p>
<div class="math">
-<p><img src="_images/math/6875433152b8a1a1624c131527286ababd8a8462.png" alt="\frac{\partial \phi}{\partial t}
+<p><img src="_images/math/ccd454c154f50b005e2c673771bf6f7fc20dcdcd.png" alt="\frac{\partial \phi}{\partial t}
+ \nabla \cdot (\phi \boldsymbol{v}) = 0"/></p>
-</div><p>For the <img class="math" src="_images/math/188c175aac0a8a9c22499336711b5d7256407254.png" alt="x"/> component, the Lagrangian formulation of the momentum equation
-with a body force <img class="math" src="_images/math/ce4f41d90f3fea6e5313ad3c0e8aac4f969ef16b.png" alt="\boldsymbol{f}"/> becomes:</p>
+</div><p>For the <img class="math" src="_images/math/a59f68a4202623bb859a7093f0316bf466e6f75d.png" alt="x"/> component, the Lagrangian formulation of the momentum equation
+with a body force <img class="math" src="_images/math/0d2439f09c0757d8e087ec316120959ccb1243b4.png" alt="\boldsymbol{f}"/> becomes:</p>
<div class="math">
-<p><img src="_images/math/60290a944b15183b06237a1483946560e640e1bf.png" alt="\frac{D (\phi v_x)}{D t}
+<p><img src="_images/math/94cb90a2611b5dcea79f236933ff4f759fddbe71.png" alt="\frac{D (\phi v_x)}{D t}
= \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\sigma}) \right]_x
- \frac{1}{\rho} f^i_x
+ \phi g"/></p>
</div><p>In the Eulerian formulation, an advection term is added, and the Cauchy stress
tensor is represented as isotropic and deviatoric components individually:</p>
<div class="math">
-<p><img src="_images/math/b991196699db1cde269ad006f427d8168fcfd679.png" alt="\frac{\partial (\phi v_x)}{\partial t}
+<p><img src="_images/math/30db4313208df84c2ad59bb33fb8d7cda1ea18ac.png" alt="\frac{\partial (\phi v_x)}{\partial t}
+ \boldsymbol{v} \cdot \nabla (\phi v_x)
= \frac{1}{\rho} \left[ \nabla \cdot (-\phi p \boldsymbol{I})
+ \phi \boldsymbol{\tau}) \right]_x
t@@ -132,7 +132,7 @@ tensor is represented as isotropic and deviatoric components individually:</p>
</div><p>Using vector identities to rewrite the advection term, and expanding the fluid
stress tensor term:</p>
<div class="math">
-<p><img src="_images/math/a85779399f71d0d2d7fd373c5c5f2b572cf3bfd7.png" alt="\frac{\partial (\phi v_x)}{\partial t}
+<p><img src="_images/math/dfc197cc1f3569837c89d6fd4eda57f135c1941a.png" alt="\frac{\partial (\phi v_x)}{\partial t}
+ \nabla \cdot (\phi v_x \boldsymbol{v})
- \phi v_x (\nabla \cdot \boldsymbol{v})
= \frac{1}{\rho} \left[ -\nabla \phi p \right]_x
t@@ -141,24 +141,24 @@ stress tensor term:</p>
+ \phi g_x"/></p>
</div><p>Spatial variations in the porosity are neglected,</p>
<div class="math">
-<p><img src="_images/math/249fc7f1e5fa5aa23b21d1fdf32ed10f0fcef400.png" alt="\nabla \phi := 0"/></p>
+<p><img src="_images/math/f986f0b44e081477adec302df83930cdd31e5a42.png" alt="\nabla \phi := 0"/></p>
</div><p>and the pressure is attributed to the fluid phase alone (model B in Zhu et al.
2007 and Zhou et al. 2010). The divergence of fluid velocities is defined to be
zero:</p>
<div class="math">
-<p><img src="_images/math/96df96ab11be34dbe17cb0837626b7cc2772013e.png" alt="\nabla \cdot \boldsymbol{v} := 0"/></p>
+<p><img src="_images/math/bb89d2c54178c6063bbddbc2ef5cc575a6ad20d6.png" alt="\nabla \cdot \boldsymbol{v} := 0"/></p>
</div><p>With these assumptions, the momentum equation simplifies to:</p>
<div class="math">
-<p><img src="_images/math/e58aec0a9dc37c78b74f5c34874ff852d8839f43.png" alt="\frac{\partial (\phi v_x)}{\partial t}
+<p><img src="_images/math/0c167e8743c59960431a47f16ee7244f2a13028f.png" alt="\frac{\partial (\phi v_x)}{\partial t}
+ \nabla \cdot (\phi v_x \boldsymbol{v})
= -\frac{1}{\rho} \frac{\partial p}{\partial x}
+ \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
- \frac{1}{\rho} f^i_x
+ \phi g_x"/></p>
-</div><p>The remaining part of the advection term is for the <img class="math" src="_images/math/188c175aac0a8a9c22499336711b5d7256407254.png" alt="x"/> component
+</div><p>The remaining part of the advection term is for the <img class="math" src="_images/math/a59f68a4202623bb859a7093f0316bf466e6f75d.png" alt="x"/> component
found as:</p>
<div class="math">
-<p><img src="_images/math/8bd85c906906e83dd3471973969959e831d271d7.png" alt="\nabla \cdot (\phi v_x \boldsymbol{v}) =
+<p><img src="_images/math/145ab6190ce7c0c86da9eda1e2fac70f23b3a5fd.png" alt="\nabla \cdot (\phi v_x \boldsymbol{v}) =
\left[
\frac{\partial}{\partial x},
\frac{\partial}{\partial y},
t@@ -175,10 +175,10 @@ found as:</p>
\frac{\partial (\phi v_x v_x)}{\partial x} +
\frac{\partial (\phi v_x v_y)}{\partial y} +
\frac{\partial (\phi v_x v_z)}{\partial z}"/></p>
-</div><p>The deviatoric stress tensor is in this case symmetrical, i.e. <img class="math" src="_images/math/3a94e291e91082dc48a18009bdd3fd70ebedec76.png" alt="\tau_{ij}
+</div><p>The deviatoric stress tensor is in this case symmetrical, i.e. <img class="math" src="_images/math/3cf4b8afb4d4ab4fd67c0576c1a06e642a0a9b7f.png" alt="\tau_{ij}
= \tau_{ji}"/>, and is found by:</p>
<div class="math">
-<p><img src="_images/math/66dc70aa131124c6f87d7f2cbb96350d8b0397b9.png" alt="\frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
+<p><img src="_images/math/0e293ff1406db52958abdce04458871a8a4c966f.png" alt="\frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
= \frac{1}{\rho}
\left[
\left[
t@@ -217,17 +217,17 @@ found as:</p>
+ \frac{\partial (\phi \tau_{xz})}{\partial z}
\right)"/></p>
</div><p>In a linear viscous fluid, the stress and strain rate
-(<img class="math" src="_images/math/1ee5b5b2fe535711338cab9fde51ba86968186e2.png" alt="\dot{\boldsymbol{\epsilon}}"/>) is linearly dependent, scaled by the
-viscosity parameter <img class="math" src="_images/math/a1a25b9ae9dc5911e30e15d5ca2fc5ee8234f1f1.png" alt="\mu_f"/>:</p>
+(<img class="math" src="_images/math/bdf1f5aad7dd65124011111a75c2b8b6c3252f0f.png" alt="\dot{\boldsymbol{\epsilon}}"/>) is linearly dependent, scaled by the
+viscosity parameter <img class="math" src="_images/math/d24e1c9504f0a7305bbf80f4d758858cd67a566a.png" alt="\mu_f"/>:</p>
<div class="math">
-<p><img src="_images/math/2133ec642b2fa2a0bdc0612dc9080df18f64da53.png" alt="\tau_{ij} = 2 \mu_f \dot{\epsilon}_{ij}
+<p><img src="_images/math/f08a1f9bfdf41e9a3f9852b0aeae4a16500bb993.png" alt="\tau_{ij} = 2 \mu_f \dot{\epsilon}_{ij}
= \mu_f \left(
\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}
\right)"/></p>
</div><p>With this relationship, the deviatoric stress tensor components can be
calculated as:</p>
<div class="math">
-<p><img src="_images/math/4d94b19a340ed5d6cfc3285e3c0678be914c3834.png" alt="\tau_{xx} = 2 \mu_f \frac{\partial v_x}{\partial x} \qquad
+<p><img src="_images/math/5fe663ae852d5c4d779d8d7827ca30a28852936a.png" alt="\tau_{xx} = 2 \mu_f \frac{\partial v_x}{\partial x} \qquad
\tau_{yy} = 2 \mu_f \frac{\partial v_y}{\partial y} \qquad
\tau_{zz} = 2 \mu_f \frac{\partial v_z}{\partial z}
t@@ -239,26 +239,26 @@ calculated as:</p>
\tau_{yz} = \mu_f \left(
\frac{\partial v_y}{\partial z} + \frac{\partial v_z}{\partial y} \right)"/></p>
-</div><p>where <img class="math" src="_images/math/a1a25b9ae9dc5911e30e15d5ca2fc5ee8234f1f1.png" alt="\mu_f"/> is the dynamic viscosity. The above formulation of the
+</div><p>where <img class="math" src="_images/math/d24e1c9504f0a7305bbf80f4d758858cd67a566a.png" alt="\mu_f"/> is the dynamic viscosity. The above formulation of the
fluid rheology assumes identical bulk and shear viscosities. The derivation of
the equations for the other spatial components is trivial.</p>
</div>
<div class="section" id="porosity-estimation">
<h2>Porosity estimation<a class="headerlink" href="#porosity-estimation" title="Permalink to this headline">¶</a></h2>
<p>The solid volume in each fluid cell is determined by the ratio of the
-a cell-centered spherical cell volume (<img class="math" src="_images/math/db94af9f89e3d01956eca6e0dc1be147b0bd4246.png" alt="V_c"/>) and the sum of intersecting
-particle volumes (<img class="math" src="_images/math/754decd655584a824e0979fc2df27ae24379558e.png" alt="V_s"/>). The spherical cell volume has a center at
-<img class="math" src="_images/math/13e4be57cf47de5302fc2e1746d965a6ee0f332c.png" alt="\boldsymbol{x}_i"/>, and a radius of <img class="math" src="_images/math/e8522a3d29a2f0f886fe156c30351ba0489992b9.png" alt="R_i"/>, which is equal to half
+a cell-centered spherical cell volume (<img class="math" src="_images/math/f998254de8fa53a31e3799a264d5fbe702f0f4d0.png" alt="V_c"/>) and the sum of intersecting
+particle volumes (<img class="math" src="_images/math/5bd22601fda700e6f36ba7e93252b6dd66efbef3.png" alt="V_s"/>). The spherical cell volume has a center at
+<img class="math" src="_images/math/5e8dd49b345fe138a2dacb54663b30b69335b6e6.png" alt="\boldsymbol{x}_i"/>, and a radius of <img class="math" src="_images/math/2b2b8f7c354490409ac10cf84c3fff2611531449.png" alt="R_i"/>, which is equal to half
the fluid cell width. The nearby particles are characterized by position
-<img class="math" src="_images/math/ec1ce56da06483e62a9456229761ff166b345ed2.png" alt="\boldsymbol{x}_j"/> and radius <img class="math" src="_images/math/b51a1f74e439e8b3b84c8fe25870d53a571391ef.png" alt="r_j"/>. The center distance is defined
+<img class="math" src="_images/math/1165d3ed99998c0f1ab1f9f25b42757182923a8b.png" alt="\boldsymbol{x}_j"/> and radius <img class="math" src="_images/math/e0d457958f773008171b4fb55f90c1c7a1a64073.png" alt="r_j"/>. The center distance is defined
as:</p>
<div class="math">
-<p><img src="_images/math/68f7c3ca449a28047e11af58f55b159b3156685c.png" alt="d_{ij} = ||\boldsymbol{x}_i - \boldsymbol{x}_j||"/></p>
+<p><img src="_images/math/3c39d664a1a005333c2e2e6088465ca5aabbc861.png" alt="d_{ij} = ||\boldsymbol{x}_i - \boldsymbol{x}_j||"/></p>
</div><p>The common volume of the two intersecting spheres is zero if the volumes aren’t
intersecting, lens shaped if they are intersecting, and spherical if the
particle is fully contained by the spherical cell volume:</p>
<div class="math">
-<p><img src="_images/math/2b8b39498e0fa8ab6324f088a2330d47ef18cbd5.png" alt="V^s_{i} = \sum_j
+<p><img src="_images/math/53e472051fcf13e3695e33c0f473e8d2c64f6ca7.png" alt="V^s_{i} = \sum_j
\begin{cases}
0 & \textit{if } R_i + r_j \leq d_{ij} \\
\frac{1}{12d_{ij}} \left[ \pi (R_i + r_j - d_{ij})^2
t@@ -268,7 +268,7 @@ particle is fully contained by the spherical cell volume:</p>
\end{cases}"/></p>
</div><p>Using this method, the cell porosity values are continuous through time as
particles enter and exit the cell volume. The rate of porosity change
-(<img class="math" src="_images/math/de27eb4f973d27763dc4d176bf2453607ecaae42.png" alt="d\phi/dt"/>) is estimated by the backwards Euler method
+(<img class="math" src="_images/math/778cd9f6a1fa43976de0f124e625d935cf4c5f43.png" alt="d\phi/dt"/>) is estimated by the backwards Euler method
by considering the previous and current porosity.</p>
</div>
<div class="section" id="particle-fluid-interaction">
t@@ -280,53 +280,53 @@ semi-empirical relationships. The drag force scales linearly with the relative
difference in velocity between the fluid and particle phase. On the base of
Newton’s third law, the resulting drag force is applied with opposite signs to
the particle and fluid.</p>
-<p>For fluid cells with porosities (<img class="math" src="_images/math/10e009bdb83f96c5f47c58b34d5d4b12ef268d5b.png" alt="\phi"/>) less or equal to 0.8, the drag
+<p>For fluid cells with porosities (<img class="math" src="_images/math/c2f31c22645274c375eff7920cfdfdc18d60341f.png" alt="\phi"/>) less or equal to 0.8, the drag
force is based on the Ergun (1952) equation:</p>
<div class="math">
-<p><img src="_images/math/c65c93ecc36954688972785ecbae82a6d7a14a3a.png" alt="\bar{\boldsymbol{f}}_d = \left(
+<p><img src="_images/math/86705dc922462e20080973641e4bcf0246476b94.png" alt="\bar{\boldsymbol{f}}_d = \left(
150 \frac{\mu_f (1-\phi)^2}{\phi\bar{d}^2}
+ 1.75 \frac{(1-\phi)\rho_f
||\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p||}{\bar{d}}
\right)
(\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p)"/></p>
-</div><p>here, <img class="math" src="_images/math/082bced08237667a772227c2e3eb8359ff9d4af8.png" alt="\bar{d}"/> denotes the average particle diameter in the cell,
-<img class="math" src="_images/math/22b6d9dab340bfe13993da20589e4720449fe391.png" alt="\boldsymbol{v}_f"/> is the fluid flow velocity, and
-<img class="math" src="_images/math/5b7072008fcf4dedaeed4d537fe5c808adedb826.png" alt="\bar{\boldsymbol{v}}_p"/> is the average particle velocity in the cell. All
+</div><p>here, <img class="math" src="_images/math/22474828ed1ac7525699de27914be48129437c04.png" alt="\bar{d}"/> denotes the average particle diameter in the cell,
+<img class="math" src="_images/math/35a936dcbe9d44868273ba50e463cc18afa0f108.png" alt="\boldsymbol{v}_f"/> is the fluid flow velocity, and
+<img class="math" src="_images/math/ee125b2b5874be1ce504888dd7d7fa42b65e9261.png" alt="\bar{\boldsymbol{v}}_p"/> is the average particle velocity in the cell. All
particles in contact with the previously mentioned cell-centered sphere for
porosity estimation contribute to the average particle velocity and diameter in
the fluid cell.</p>
<p>If the porosity is greater than 0.8, the cell-averaged drag force
-(<img class="math" src="_images/math/fa8a7b2fb896adf1c9f17077664691e511aeb8e5.png" alt="\bar{\boldsymbol{f}}_d"/> is found from the Wen and Yu (1966) equation,
+(<img class="math" src="_images/math/7dcc7d011290597e0833cc160b66996c5d9d4595.png" alt="\bar{\boldsymbol{f}}_d"/> is found from the Wen and Yu (1966) equation,
which considers the fluid flow situation:</p>
<div class="math">
-<p><img src="_images/math/0949063a26c6869766d571d13e61cf78f3a38c06.png" alt="\bar{\boldsymbol{f}}_d = \left(
+<p><img src="_images/math/75a974365e0604c0e8675dc2543166a1216a324e.png" alt="\bar{\boldsymbol{f}}_d = \left(
\frac{3}{4}
\frac{C_d (1-\phi) \phi^{-2.65} \mu_f \rho_f
||\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p||}{\bar{d}}
\right)
(\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p)"/></p>
-</div><p>The drag coefficient <img class="math" src="_images/math/1beccd9a659ce57a6a8083a589fe91ad33fa3b06.png" alt="C_d"/> is evaluated depending on the magnitude of the
-Reynolds number <img class="math" src="_images/math/0e0f0c7e69e0f315aed6d1762a66bf27dd9a155a.png" alt="Re"/>:</p>
+</div><p>The drag coefficient <img class="math" src="_images/math/15746e47c01eb1c4d15e0409dfd2ae36896638ab.png" alt="C_d"/> is evaluated depending on the magnitude of the
+Reynolds number <img class="math" src="_images/math/e7d45090e3dbe0db821c4e5fdfb275a0e0fdf533.png" alt="Re"/>:</p>
<div class="math">
-<p><img src="_images/math/e2872fb63d995fa1112c468dcb10669c69795aef.png" alt="C_d =
+<p><img src="_images/math/5962c6c4e4e193ec5e6b17320744252e30a86315.png" alt="C_d =
\begin{cases}
\frac{24}{Re} (1+0.15 (Re)^{0.687} & \textit{if } Re < 1,000 \\
0.44 & \textit{if } Re \geq 1,000
\end{cases}"/></p>
</div><p>where the Reynold’s number is found by:</p>
<div class="math">
-<p><img src="_images/math/f4bf201ecead6e3486688493ab1b1e82c39d93df.png" alt="Re = \frac{\phi\rho_f\bar{d}}{\mu_f}
+<p><img src="_images/math/8be75291f64c7dd21b34613d8a5346f3b9aff877.png" alt="Re = \frac{\phi\rho_f\bar{d}}{\mu_f}
||\boldsymbol{v}_f - \bar{\boldsymbol{v}}_p||"/></p>
</div><p>The interaction force is applied to the fluid with negative sign as a
-contribution to the body force <img class="math" src="_images/math/ce4f41d90f3fea6e5313ad3c0e8aac4f969ef16b.png" alt="\boldsymbol{f}"/>. The fluid interaction
+contribution to the body force <img class="math" src="_images/math/0d2439f09c0757d8e087ec316120959ccb1243b4.png" alt="\boldsymbol{f}"/>. The fluid interaction
force applied particles in the fluid cell is:</p>
<div class="math">
-<p><img src="_images/math/c24e3d875b09e3973624f97802672222ca57c3be.png" alt="\boldsymbol{f}_i = \frac{\bar{\boldsymbol{f}}_d V_p}{1-\phi}"/></p>
-</div><p>where <img class="math" src="_images/math/4d3f01b8ca8bdab3543eda695ff7e565bcb66381.png" alt="V_p"/> denotes the particle volume. Optionally, the above
+<p><img src="_images/math/4eefba7b01497c1c37254ba0938852a69ca402da.png" alt="\boldsymbol{f}_i = \frac{\bar{\boldsymbol{f}}_d V_p}{1-\phi}"/></p>
+</div><p>where <img class="math" src="_images/math/0e157574dc52c16b7aad0f260830b27b523706f4.png" alt="V_p"/> denotes the particle volume. Optionally, the above
interaction force could be expanded to include the force induced by the fluid
pressure gradient:</p>
<div class="math">
-<p><img src="_images/math/f8484aaf38ecdf73f7489fe6b43dec92ae505398.png" alt="\boldsymbol{f}_i = \left(
+<p><img src="_images/math/5fc35741bb341c4a7bac2b261074cee511b162b9.png" alt="\boldsymbol{f}_i = \left(
-\nabla p +
\frac{\bar{\boldsymbol{f}}_d}{1-\phi}
\right) V_p"/></p>
t@@ -336,10 +336,10 @@ pressure gradient:</p>
<p>The partial differential terms in the previously described equations are found
using finite central differences. Modifying the operator splitting methodology
presented by Langtangen et al. (2002), the predicted velocity
-<img class="math" src="_images/math/f5d499efbd23360f6e9229bdccf14e0951954730.png" alt="\boldsymbol{v}^*"/> after a finite time step
-<img class="math" src="_images/math/5b4271afe7fc7c0ee84172e0ad19b82caf450c00.png" alt="\Delta t"/> is found by explicit integration of the momentum equation.</p>
+<img class="math" src="_images/math/8615f3b3cbbe1dc55cb6b255dad751d458a4ed3c.png" alt="\boldsymbol{v}^*"/> after a finite time step
+<img class="math" src="_images/math/ec002955bdf95ee9869878fbad4f80fc98539359.png" alt="\Delta t"/> is found by explicit integration of the momentum equation.</p>
<div class="math">
-<p><img src="_images/math/4414da5fcbe9d66f08400308688cad041ca9dd9c.png" alt="\frac{\Delta (\phi v_x)}{\Delta t}
+<p><img src="_images/math/3f8cfba42aa84819f05cca10369f3913e4833b5a.png" alt="\frac{\Delta (\phi v_x)}{\Delta t}
+ \nabla \cdot (\phi v_x \boldsymbol{v})
= - \frac{1}{\rho} \frac{\Delta p}{\Delta x}
+ \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
t@@ -355,10 +355,10 @@ presented by Langtangen et al. (2002), the predicted velocity
+ \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
- \frac{1}{\rho} f^i_x
+ \phi g_x"/></p>
-</div><p>We want to isolate <img class="math" src="_images/math/eee8ed92791a9c140bb435f75178642c8bc57146.png" alt="\Delta v_x"/> in the above equation in order to project
+</div><p>We want to isolate <img class="math" src="_images/math/154f421588a5329f4f7216120a639b998a94cade.png" alt="\Delta v_x"/> in the above equation in order to project
the new velocity.</p>
<div class="math">
-<p><img src="_images/math/609c6937a183637c7356d1ee9e2800015d73a5f4.png" alt="\phi \frac{\Delta v_x}{\Delta t}
+<p><img src="_images/math/f8e4fccaabfa73731975fa62731bb79fed02f428.png" alt="\phi \frac{\Delta v_x}{\Delta t}
= - \frac{1}{\rho} \frac{\Delta p}{\Delta x}
+ \frac{1}{\rho} \left[ \nabla \cdot (\phi \boldsymbol{\tau}) \right]_x
- \frac{1}{\rho} f^i_x
t@@ -374,13 +374,13 @@ the new velocity.</p>
+ \Delta t g_x
- v_x \frac{\Delta \phi}{\phi}
- \nabla \cdot (\phi v_x \boldsymbol{v}) \frac{\Delta t}{\phi}"/></p>
-</div><p>The term <img class="math" src="_images/math/8ce03f78ed945f2ef3dac87c8799b55b393527e7.png" alt="\beta"/> is introduced as an adjustable, dimensionless parameter
-in the range <img class="math" src="_images/math/30d3be82eb9bf820400a3dc3e46e84fe58e9f2c9.png" alt="[0;1]"/>, and determines the importance of the old pressure
+</div><p>The term <img class="math" src="_images/math/410a9d0df9c135dd73b269cba7ef04dcfd932b1f.png" alt="\beta"/> is introduced as an adjustable, dimensionless parameter
+in the range <img class="math" src="_images/math/b00c7b17d75d91818f2f9060562323eb06789226.png" alt="[0;1]"/>, and determines the importance of the old pressure
values in the solution procedure (Langtangen et al. 2002). A value of 0
corresponds to <a class="reference external" href="https://en.wikipedia.org/wiki/Projection_method_(fluid_dynamics)#Chorin.27s_projection_method">Chorin’s projection method</a> originally described
in <a class="reference external" href="http://www.ams.org/journals/mcom/1968-22-104/S0025-5718-1968-0242392-2/S0025-5718-1968-0242392-2.pdf">Chorin (1968)</a>.</p>
<div class="math">
-<p><img src="_images/math/7ac8e5caca3e02a9a9b6c9576ed9bd5effddb2fb.png" alt="v_x^* = v_x^t + \Delta v_x
+<p><img src="_images/math/30bcc3939e4e8f039a0879219ae7d4488c349831.png" alt="v_x^* = v_x^t + \Delta v_x
v_x^* = v_x^t
- \frac{\beta}{\rho} \frac{\Delta p^t}{\Delta x} \frac{\Delta t}{\phi^t}
t@@ -390,22 +390,22 @@ v_x^* = v_x^t
+ \Delta t g_x
- v^t_x \frac{\Delta \phi}{\phi^t}
- \nabla \cdot (\phi^t v_x^t \boldsymbol{v}^t) \frac{\Delta t}{\phi^t}"/></p>
-</div><p>Here, <img class="math" src="_images/math/98a432864561550105d4de621bc0a9b61068c043.png" alt="\Delta x"/> denotes the cell spacing. The velocity found
-(<img class="math" src="_images/math/e6eadc1093d6e4f5a41f86b8301c0d7559d2c609.png" alt="v_x^*"/>) is only a prediction of the fluid velocity at time
-<img class="math" src="_images/math/e35ec9f12692dad5164f15f68fc16f969099e241.png" alt="t+\Delta t"/>, since the estimate isn’t constrained by the continuity
+</div><p>Here, <img class="math" src="_images/math/05899b70d35eb61a85aaaba3affb9c7363f0d633.png" alt="\Delta x"/> denotes the cell spacing. The velocity found
+(<img class="math" src="_images/math/69274c46c55ea9c2723428b533616806381a2b5f.png" alt="v_x^*"/>) is only a prediction of the fluid velocity at time
+<img class="math" src="_images/math/acbab363f98c011e1a927c8ff8b4f565421f5def.png" alt="t+\Delta t"/>, since the estimate isn’t constrained by the continuity
equation:</p>
<div class="math">
-<p><img src="_images/math/5a4c725fb6cd73a97b137a6ddc218bd457e95ab1.png" alt="\frac{\Delta \phi^t}{\Delta t} + \nabla \cdot (\phi^t
+<p><img src="_images/math/4fc3cb57fe6e8e2710c92acae1dae02dfb464b2c.png" alt="\frac{\Delta \phi^t}{\Delta t} + \nabla \cdot (\phi^t
\boldsymbol{v}^{t+\Delta t}) = 0"/></p>
</div><p>The divergence of a scalar and vector can be <a class="reference external" href="http://www.wolframalpha.com/input/?i=div(p+v)">split</a>:</p>
<div class="math">
-<p><img src="_images/math/6b14a983094b158179fe53e0a59ebe0d53891195.png" alt="\phi^t \nabla \cdot \boldsymbol{v}^{t+\Delta t} +
+<p><img src="_images/math/dd69e8d47250a0448f829f2e1d42222cbda38058.png" alt="\phi^t \nabla \cdot \boldsymbol{v}^{t+\Delta t} +
\boldsymbol{v}^{t+\Delta t} \cdot \nabla \phi^t
+ \frac{\Delta \phi^t}{\Delta t} = 0"/></p>
</div><p>The predicted velocity is corrected using the new pressure (Langtangen et al.
2002):</p>
<div class="math">
-<p><img src="_images/math/2fcd737952b690517150c61a5981d5e810d76268.png" alt="\boldsymbol{v}^{t+\Delta t} = \boldsymbol{v}^*
+<p><img src="_images/math/849751204b923a6bcea83d33614df1ce13444a71.png" alt="\boldsymbol{v}^{t+\Delta t} = \boldsymbol{v}^*
%- \frac{\Delta t}{\rho} \nabla \epsilon
- \frac{\Delta t}{\rho \phi^t} \nabla \epsilon
\quad \text{where} \quad
t@@ -413,21 +413,21 @@ equation:</p>
</div><p>The above formulation of the future velocity is put into the continuity
equation:</p>
<div class="math">
-<p><img src="_images/math/915ee93516c4e9e5e9df4b596291692f2b283c00.png" alt="\Rightarrow
+<p><img src="_images/math/ba57616deddc047e5626df1a0e2e2c4f49c4e95a.png" alt="\Rightarrow
\phi^t \nabla \cdot
\left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \right)
+
\left( \boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \nabla \epsilon \right)
\cdot \nabla \phi^t + \frac{\Delta \phi^t}{\Delta t} = 0"/></p>
</div><div class="math">
-<p><img src="_images/math/a1bc6e0933f817147b43d6662c901ffc0b2b4151.png" alt="\Rightarrow
+<p><img src="_images/math/5695b27fddc4f8e8a8d1194a6e74fd72eacc1a4d.png" alt="\Rightarrow
\phi^t \nabla \cdot
\boldsymbol{v}^* - \frac{\Delta t}{\rho \phi^t} \phi^t \nabla^2 \epsilon
+ \nabla \phi^t \cdot \boldsymbol{v}^*
- \nabla \phi^t \cdot \nabla \epsilon \frac{\Delta t}{\rho \phi^t}
+ \frac{\Delta \phi^t}{\Delta t} = 0"/></p>
</div><div class="math">
-<p><img src="_images/math/03b863159155c1ff0a1478a3a14f9ceecb42221e.png" alt="\Rightarrow
+<p><img src="_images/math/bf12e075e2a1ffed21117ae6c3b5a57228a35e48.png" alt="\Rightarrow
\frac{\Delta t}{\rho} \nabla^2 \epsilon
= \phi^t \nabla \cdot \boldsymbol{v}^*
+ \nabla \phi^t \cdot \boldsymbol{v}^*
t@@ -435,38 +435,38 @@ equation:</p>
+ \frac{\Delta \phi^t}{\Delta t}"/></p>
</div><p>The pressure difference in time becomes a <a class="reference external" href="https://en.wikipedia.org/wiki/Poisson's_equation">Poisson equation</a> with added terms:</p>
<div class="math">
-<p><img src="_images/math/49cc970cd820d4022dd996e66da2e453cf4aeeda.png" alt="\Rightarrow
+<p><img src="_images/math/af973ef88ea161ebaf13a396d18b9fd419b48392.png" alt="\Rightarrow
\nabla^2 \epsilon
= \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t}
+ \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t}
- \frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t}
+ \frac{\Delta \phi^t \rho}{\Delta t^2}"/></p>
</div><p>The right hand side of the above equation is termed the <em>forcing function</em>
-<img class="math" src="_images/math/0001d02b63ede2fe3219e05a7cd09c82ae6298b6.png" alt="f"/>, which is decomposed into two terms, <img class="math" src="_images/math/aad9832dd45598090e4be6aaeeaa038dfc093751.png" alt="f_1"/> and <img class="math" src="_images/math/ff9078c258266635e46b9767b0c7e248c865e4ab.png" alt="f_2"/>:</p>
+<img class="math" src="_images/math/875eb40014526135383caa89fd500ae40a835f56.png" alt="f"/>, which is decomposed into two terms, <img class="math" src="_images/math/09968cb620b9b121ddcc33bb90eab038355375ab.png" alt="f_1"/> and <img class="math" src="_images/math/a23addf1349384bc83f74ce7c550fe98bbdc2b48.png" alt="f_2"/>:</p>
<div class="math">
-<p><img src="_images/math/ba756dc27a3f5d7f8ff7d7801d8c444393c13e6a.png" alt="f_1
+<p><img src="_images/math/6d2f443abfc2651609cd8dd745ff62139f7d51c4.png" alt="f_1
= \frac{\nabla \cdot \boldsymbol{v}^* \phi^t \rho}{\Delta t}
+ \frac{\nabla \phi^t \cdot \boldsymbol{v}^* \rho}{\Delta t}
+ \frac{\Delta \phi^t \rho}{\Delta t^2}
f_2 =
\frac{\nabla \phi^t \cdot \nabla \epsilon}{\phi^t}"/></p>
-</div><p>During the <a class="reference external" href="http://www.rsmas.miami.edu/personal/miskandarani/Courses/MSC321/Projects/prjpoisson.pdf">Jacobi iterative solution procedure</a> <img class="math" src="_images/math/aad9832dd45598090e4be6aaeeaa038dfc093751.png" alt="f_1"/> remains constant,
-while <img class="math" src="_images/math/ff9078c258266635e46b9767b0c7e248c865e4ab.png" alt="f_2"/> changes value. For this reason, <img class="math" src="_images/math/aad9832dd45598090e4be6aaeeaa038dfc093751.png" alt="f_1"/> is found only
-during the first iteration, while <img class="math" src="_images/math/ff9078c258266635e46b9767b0c7e248c865e4ab.png" alt="f_2"/> is updated every time. The value
+</div><p>During the <a class="reference external" href="http://www.rsmas.miami.edu/personal/miskandarani/Courses/MSC321/Projects/prjpoisson.pdf">Jacobi iterative solution procedure</a> <img class="math" src="_images/math/09968cb620b9b121ddcc33bb90eab038355375ab.png" alt="f_1"/> remains constant,
+while <img class="math" src="_images/math/a23addf1349384bc83f74ce7c550fe98bbdc2b48.png" alt="f_2"/> changes value. For this reason, <img class="math" src="_images/math/09968cb620b9b121ddcc33bb90eab038355375ab.png" alt="f_1"/> is found only
+during the first iteration, while <img class="math" src="_images/math/a23addf1349384bc83f74ce7c550fe98bbdc2b48.png" alt="f_2"/> is updated every time. The value
of the forcing function is found as:</p>
<div class="math">
-<p><img src="_images/math/2a6ab5768c27f4033a07a2bd556e89a39cdff6c0.png" alt="f = f_1 - f_2"/></p>
+<p><img src="_images/math/7145387c0de7584800cc45b61520c507b70038cb.png" alt="f = f_1 - f_2"/></p>
</div><p>Using second-order finite difference approximations of the Laplace operator
second-order partial derivatives, the differential equations become a system of
equations that is solved using <a class="reference external" href="https://en.wikipedia.org/wiki/Relaxation_(iterative_method)">iteratively</a> using Jacobi updates. The total
-number of unknowns is <img class="math" src="_images/math/7774004c6a3397305cc7616b705accbbbfb28632.png" alt="(n_x - 1)(n_y - 1)(n_z - 1)"/>.</p>
+number of unknowns is <img class="math" src="_images/math/046284126cf107c26a2999e7949ec80c95f65457.png" alt="(n_x - 1)(n_y - 1)(n_z - 1)"/>.</p>
<p>The discrete Laplacian (approximation of the Laplace operator) can be obtained
by a finite-difference seven-point stencil in a three-dimensional, cubic
-grid with cell spacing <img class="math" src="_images/math/43e42baf6cef7e9620082be300b7622d524b3a65.png" alt="\Delta x, \Delta y, \Delta z"/>, considering the six
+grid with cell spacing <img class="math" src="_images/math/fffb1eace0724fc01f1ea20e385d9f356b44f629.png" alt="\Delta x, \Delta y, \Delta z"/>, considering the six
face neighbors:</p>
<div class="math">
-<p><img src="_images/math/ae70bc18ed864da95e96ab31ccce2f21f06d6979.png" alt="\nabla^2 \epsilon_{i_x,i_y,i_z} \approx
+<p><img src="_images/math/976d1f7bdf9faa8c8c428dbcc25b1540eedaa139.png" alt="\nabla^2 \epsilon_{i_x,i_y,i_z} \approx
\frac{\epsilon_{i_x-1,i_y,i_z} - 2 \epsilon_{i_x,i_y,i_z}
+ \epsilon_{i_x+1,i_y,i_z}}{\Delta x^2}
+ \frac{\epsilon_{i_x,i_y-1,i_z} - 2 \epsilon_{i_x,i_y,i_z}
t@@ -475,11 +475,11 @@ face neighbors:</p>
+ \frac{\epsilon_{i_x,i_y,i_z-1} - 2 \epsilon_{i_x,i_y,i_z}
+ \epsilon_{i_x,i_y,i_z+1}}{\Delta z^2}
\approx f_{i_x,i_y,i_z}"/></p>
-</div><p>Within a Jacobi iteration, the value of the unknowns (<img class="math" src="_images/math/458d3955e23ef7ba378c27f284a2e426f756bd3e.png" alt="\epsilon^n"/>) is
-used to find an updated solution estimate (<img class="math" src="_images/math/29bbcc61ff8f993f8eca2f33c8faf0c7ad80defc.png" alt="\epsilon^{n+1}"/>).
+</div><p>Within a Jacobi iteration, the value of the unknowns (<img class="math" src="_images/math/d0d1e04f7f59d8b2990d0e765f452a9ae94e4096.png" alt="\epsilon^n"/>) is
+used to find an updated solution estimate (<img class="math" src="_images/math/a87cadafacc1e85e0e69da8c11478ec665517312.png" alt="\epsilon^{n+1}"/>).
The solution for the updated value takes the form:</p>
<div class="math">
-<p><img src="_images/math/11c52ddf6d9f9157b67c1334ccb39807a8ee777c.png" alt="\epsilon^{n+1}_{i_x,i_y,i_z}
+<p><img src="_images/math/52f1699b7bdfcd7da418c1f5639cfdfa40ed99a7.png" alt="\epsilon^{n+1}_{i_x,i_y,i_z}
= \frac{-\Delta x^2 \Delta y^2 \Delta z^2 f_{i_x,i_y,i_z}
+ \Delta y^2 \Delta z^2 (\epsilon^n_{i_x-1,i_y,i_z} +
\epsilon^n_{i_x+1,i_y,i_z})
t@@ -493,42 +493,42 @@ The solution for the updated value takes the form:</p>
</div><p>The difference between the current and updated value is termed the <em>normalized
residual</em>:</p>
<div class="math">
-<p><img src="_images/math/a6f6238418bca503a9616e884c92c10d86d34fbc.png" alt="r_{i_x,i_y,i_z} = \frac{(\epsilon^{n+1}_{i_x,i_y,i_z}
+<p><img src="_images/math/b843ab93a9cd31bd2f99132fe09390b3b2e69b96.png" alt="r_{i_x,i_y,i_z} = \frac{(\epsilon^{n+1}_{i_x,i_y,i_z}
- \epsilon^n_{i_x,i_y,i_z})^2}{(\epsilon^{n+1}_{i_x,i_y,i_z})^2}"/></p>
-</div><p>Note that the <img class="math" src="_images/math/19bc0073dde1bcd1a8e6a32b251e80cced668f04.png" alt="\epsilon"/> values cannot be 0 due to the above normalization
+</div><p>Note that the <img class="math" src="_images/math/65d19c66c148d5016c6a89d26486bf6d1966ded1.png" alt="\epsilon"/> values cannot be 0 due to the above normalization
of the residual.</p>
<p>The updated values are at the end of the iteration stored as the current values,
and the maximal value of the normalized residual is found. If this value is
larger than a tolerance criteria, the procedure is repeated. The iterative
procedure is ended if the number of iterations exceeds a defined limit.</p>
-<p>After the values of <img class="math" src="_images/math/19bc0073dde1bcd1a8e6a32b251e80cced668f04.png" alt="\epsilon"/> are found, they are used to find the new
+<p>After the values of <img class="math" src="_images/math/65d19c66c148d5016c6a89d26486bf6d1966ded1.png" alt="\epsilon"/> are found, they are used to find the new
pressures and velocities:</p>
<div class="math">
-<p><img src="_images/math/e80afabbf8aabeec85a867c9c73d1f2df69743fa.png" alt="\bar{p}^{t+\Delta t} = \beta \bar{p}^t + \epsilon"/></p>
+<p><img src="_images/math/34ae3abc95078f1fbc7d24815a091daacc0293d6.png" alt="\bar{p}^{t+\Delta t} = \beta \bar{p}^t + \epsilon"/></p>
</div><div class="math">
-<p><img src="_images/math/26446b608e4ea66bd2ac56325be3b7dad3ec4261.png" alt="\bar{\boldsymbol{v}}^{t+\Delta t} =
+<p><img src="_images/math/d653d13e7a5e33414bfff96b3015fb597c5355f5.png" alt="\bar{\boldsymbol{v}}^{t+\Delta t} =
\bar{\boldsymbol{v}}^* - \frac{\Delta t}{\rho\phi} \nabla \epsilon"/></p>
</div></div>
<div class="section" id="boundary-conditions">
<h2>Boundary conditions<a class="headerlink" href="#boundary-conditions" title="Permalink to this headline">¶</a></h2>
<p>The lateral boundaries are periodic. This cannot be changed in the current
version of <code class="docutils literal"><span class="pre">sphere</span></code>. This means that the fluid properties at the paired,
-parallel lateral (<img class="math" src="_images/math/188c175aac0a8a9c22499336711b5d7256407254.png" alt="x"/> and <img class="math" src="_images/math/b124ff74afb0914bb434e8fb849eb56d734412f8.png" alt="y"/>) boundaries are identical. A flow
+parallel lateral (<img class="math" src="_images/math/a59f68a4202623bb859a7093f0316bf466e6f75d.png" alt="x"/> and <img class="math" src="_images/math/276f7e256cbddeb81eee42e1efc348f3cb4ab5f8.png" alt="y"/>) boundaries are identical. A flow
leaving through one side reappears on the opposite side.</p>
<p>The top and bottom boundary conditions of the fluid grid can be either:
prescribed pressure (Dirichlet), or prescribed velocity (Neumann). The
(horizontal) velocities parallel to the boundaries are free to attain other
values (free slip). The Dirichlet boundary condition is enforced by keeping the
-value of <img class="math" src="_images/math/19bc0073dde1bcd1a8e6a32b251e80cced668f04.png" alt="\epsilon"/> constant at the boundaries, e.g.:</p>
+value of <img class="math" src="_images/math/65d19c66c148d5016c6a89d26486bf6d1966ded1.png" alt="\epsilon"/> constant at the boundaries, e.g.:</p>
<div class="math">
-<p><img src="_images/math/d04d640e7e6ee0fbcd4ba0e9fc6cb43f3e536d12.png" alt="\epsilon^{n+1}_{i_x,i_y,i_z = 1 \vee n_z}
+<p><img src="_images/math/9b35bf83c6ef9ec4ea134990256960d8e1c204d0.png" alt="\epsilon^{n+1}_{i_x,i_y,i_z = 1 \vee n_z}
=
\epsilon^{n}_{i_x,i_y,i_z = 1 \vee n_z}"/></p>
</div><p>The Neumann boundary condition of no flow across the boundary is enforced by
-setting the gradient of <img class="math" src="_images/math/19bc0073dde1bcd1a8e6a32b251e80cced668f04.png" alt="\epsilon"/> perpendicular to the boundary to zero,
+setting the gradient of <img class="math" src="_images/math/65d19c66c148d5016c6a89d26486bf6d1966ded1.png" alt="\epsilon"/> perpendicular to the boundary to zero,
e.g.:</p>
<div class="math">
-<p><img src="_images/math/f82d0eba51cfeff3da50b567ffa5126a8ed5b15f.png" alt="\nabla_z \epsilon^{n+1}_{i_x,i_y,i_z = 1 \vee n_z} = 0"/></p>
+<p><img src="_images/math/e244cc08185d93b146a75a8fa9618b0c83b80db3.png" alt="\nabla_z \epsilon^{n+1}_{i_x,i_y,i_z = 1 \vee n_z} = 0"/></p>
</div></div>
<div class="section" id="numerical-implementation">
<h2>Numerical implementation<a class="headerlink" href="#numerical-implementation" title="Permalink to this headline">¶</a></h2>
t@@ -598,7 +598,7 @@ e.g.:</p>
<li class="right" >
<a href="dem.html" title="Discrete element method"
>previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/dem.html b/doc/html/dem.html
t@@ -6,7 +6,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>Discrete element method — sphere 1.00-alpha documentation</title>
+ <title>Discrete element method — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -14,7 +14,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -25,7 +25,7 @@
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="index.html" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="index.html" />
<link rel="next" title="Fluid simulation and particle-fluid interaction" href="cfd.html" />
<link rel="prev" title="Introduction and Installation" href="introduction.html" />
</head>
t@@ -45,7 +45,7 @@
<li class="right" >
<a href="introduction.html" title="Introduction and Installation"
accesskey="P">previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -82,10 +82,10 @@ derivative.</p>
<h2>Contact search<a class="headerlink" href="#contact-search" title="Permalink to this headline">¶</a></h2>
<p>Homogeneous cubic grid.</p>
<div class="math">
-<p><img src="_images/math/7a1586ce67a2a8d709a374bb524625ca12559fd0.png" alt="\delta_n^{ij} = ||\boldsymbol{x}^i - \boldsymbol{x}^j|| - (r^i + r^j)"/></p>
-</div><p>where <img class="math" src="_images/math/2ede365ad144ab396916ec60458da03860803078.png" alt="r"/> is the particle radius, and <img class="math" src="_images/math/80ba696b882e52c7e4d833039710ff8b4005e777.png" alt="\boldsymbol{x}"/> denotes the
-positional vector of a particle, and <img class="math" src="_images/math/a581f053bbfa5115f42c13094857cdd12a37ec49.png" alt="i"/> and <img class="math" src="_images/math/d32c78b759903e3f4bd4fd2ce0b86358f7500c5d.png" alt="j"/> denote the indexes
-of two particles. Negative values of <img class="math" src="_images/math/2822d1531935f7a072ebaa08533f4c5695731465.png" alt="\delta_n"/> denote that the particles
+<p><img src="_images/math/9a32e316af25b770f61240f0c87520e61a20170e.png" alt="\delta_n^{ij} = ||\boldsymbol{x}^i - \boldsymbol{x}^j|| - (r^i + r^j)"/></p>
+</div><p>where <img class="math" src="_images/math/eaa6ad49a7f78fe5a13b486690163bf2dc7e3e60.png" alt="r"/> is the particle radius, and <img class="math" src="_images/math/1fb53a6e43d7ce3b8f6b9680713299ee42199ead.png" alt="\boldsymbol{x}"/> denotes the
+positional vector of a particle, and <img class="math" src="_images/math/df0deb143e5ac127f00bd248ee8001ecae572adc.png" alt="i"/> and <img class="math" src="_images/math/6b21e0b0899a0d2879d3b8019087fa630bab4ea2.png" alt="j"/> denote the indexes
+of two particles. Negative values of <img class="math" src="_images/math/f20d537e82e4ab3d4da02f8cfd843220049121ec.png" alt="\delta_n"/> denote that the particles
are overlapping.</p>
</div>
<div class="section" id="contact-interaction">
t@@ -96,12 +96,12 @@ interaction is decomposed into normal and tangential components, relative to the
contact interface orientation. The normal vector to the contact interface is
found by:</p>
<div class="math">
-<p><img src="_images/math/b2f647d3f64c336d2581f03558393e1dbb92c919.png" alt="\boldsymbol{n}^{ij} =
+<p><img src="_images/math/edf49df0f6a5c2f2624ad844873e377c86bd3e50.png" alt="\boldsymbol{n}^{ij} =
\frac{\boldsymbol{x}^i - \boldsymbol{x}^j}
{||\boldsymbol{x}^i - \boldsymbol{x}^j||}"/></p>
-</div><p>The contact velocity <img class="math" src="_images/math/522c9c96bc63142fa2b83abe960622dcafd75875.png" alt="\dot{\boldsymbol{\delta}}"/> is found by:</p>
+</div><p>The contact velocity <img class="math" src="_images/math/18c0426322089124abc473b5cbb17f99ee7aed73.png" alt="\dot{\boldsymbol{\delta}}"/> is found by:</p>
<div class="math">
-<p><img src="_images/math/3e0f4950bf4cf691a83dd643eddd240f024b9806.png" alt="\dot{\boldsymbol{\delta}}^{ij} =
+<p><img src="_images/math/dc9ba965d2e46bab75a9783346180b4e90516118.png" alt="\dot{\boldsymbol{\delta}}^{ij} =
(\boldsymbol{x}^i - \boldsymbol{x}^j)
+ (r^i + \frac{\delta_n^{ij}}{2})
(\boldsymbol{n}^{ij} \times \boldsymbol{\omega}^{i})
t@@ -110,28 +110,28 @@ found by:</p>
</div><p>The contact velocity is decomposed into normal and tangential components,
relative to the contact interface. The normal component is:</p>
<div class="math">
-<p><img src="_images/math/d5f7df718c75314abef3cc68396ac19229c2b55d.png" alt="\dot{\delta}^{ij}_n =
+<p><img src="_images/math/7064634d7b87e0dd07b7dcfee5d1ab549b5964fa.png" alt="\dot{\delta}^{ij}_n =
-(\dot{\boldsymbol{\delta}}^{ij} \cdot \boldsymbol{n}^{ij})"/></p>
</div><p>and the tangential velocity component is found as:</p>
<div class="math">
-<p><img src="_images/math/10eaed38fefc0732b59e674c14396ff74ea89c93.png" alt="\dot{\boldsymbol{\delta}}^{ij}_t =
+<p><img src="_images/math/e9905af08d85a280a4c2996689972016f7833fef.png" alt="\dot{\boldsymbol{\delta}}^{ij}_t =
\dot{\boldsymbol{\delta}}^{ij}
- \boldsymbol{n}^{ij}
(\boldsymbol{n}^{ij} \cdot \dot{\boldsymbol{\delta}}^{ij})"/></p>
-</div><p>where <img class="math" src="_images/math/b95363838d60e863c4cb530ba409fa56dd735728.png" alt="\boldsymbol{\omega}"/> is the rotational velocity vector of a
+</div><p>where <img class="math" src="_images/math/cd1d64597c760ff1fc8d90d0d08209541947f9bd.png" alt="\boldsymbol{\omega}"/> is the rotational velocity vector of a
particle. The total tangential displacement on the contact plane is found
incrementally:</p>
<div class="math">
-<p><img src="_images/math/904bfd360afed4f050929e89e2cc4f7c6f673e97.png" alt="\boldsymbol{\delta}_{t,\text{uncorrected}}^{ij} =
+<p><img src="_images/math/d3e56488ee50f821a999470459f170e3b0cac3d3.png" alt="\boldsymbol{\delta}_{t,\text{uncorrected}}^{ij} =
\int_0^{t_c}
\dot{\boldsymbol{\delta}}^{ij}_t \Delta t"/></p>
-</div><p>where <img class="math" src="_images/math/42ea00932c36940fb8965e6e7f4fb2de5d616ee7.png" alt="t_c"/> is the duration of the contact and <img class="math" src="_images/math/5b4271afe7fc7c0ee84172e0ad19b82caf450c00.png" alt="\Delta t"/> is the
+</div><p>where <img class="math" src="_images/math/b8cd14d20ff38e8c9e305878c37a283e655e25f9.png" alt="t_c"/> is the duration of the contact and <img class="math" src="_images/math/ec002955bdf95ee9869878fbad4f80fc98539359.png" alt="\Delta t"/> is the
computational time step length. The tangential contact interface displacement is
set to zero when a contact pair no longer overlaps. At each time step, the value
-of <img class="math" src="_images/math/8cfa62046b3a57f37068f027941acea702d9bd8a.png" alt="\boldsymbol{\delta}_t"/> is corrected for rotation of the contact
+of <img class="math" src="_images/math/359991b87d5889ece1cd8b1c84cf17fbb37b815d.png" alt="\boldsymbol{\delta}_t"/> is corrected for rotation of the contact
interface:</p>
<div class="math">
-<p><img src="_images/math/ef61e2c4dfb9fcb1053c8d2557e00bcf3d871b80.png" alt="\boldsymbol{\delta}_t^{ij} = \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij}
+<p><img src="_images/math/71b3ac73bf6447e6cc7ea151b5664dc723a38ea9.png" alt="\boldsymbol{\delta}_t^{ij} = \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij}
- (\boldsymbol{n}
(\boldsymbol{n} \cdot \boldsymbol{\delta}_{t,\text{uncorrected}}^{ij})"/></p>
</div><p>With all the geometrical and kinetic components determined, the resulting forces
t@@ -139,17 +139,17 @@ of the particle interaction can be determined using a contact model. <code class
features only one contact model in the normal direction to the contact; the
linear-elastic-viscous (<em>Hookean</em> with viscous damping, or <em>Kelvin-Voigt</em>)
contact model. The resulting force in the normal direction of the contact
-interface on particle <img class="math" src="_images/math/a581f053bbfa5115f42c13094857cdd12a37ec49.png" alt="i"/> is:</p>
+interface on particle <img class="math" src="_images/math/df0deb143e5ac127f00bd248ee8001ecae572adc.png" alt="i"/> is:</p>
<div class="math">
-<p><img src="_images/math/7d72e011d7f06b5df7ff7decf324f701eae33c3a.png" alt="\boldsymbol{f}_n^{ij} = \left(
+<p><img src="_images/math/2b331ff6c55b9e20fb12706828da440a9ef5b20c.png" alt="\boldsymbol{f}_n^{ij} = \left(
-k_n \delta_n^{ij} -\gamma_n \dot{\delta_n}^{ij}
\right) \boldsymbol{n}^{ij}"/></p>
-</div><p>The parameter <img class="math" src="_images/math/66a1d8738af515286414e18d76a6e7f86b989afc.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
-contact interface, and <img class="math" src="_images/math/3640b174f15a44d33ad8f0622bc44dce8d8f1692.png" alt="\gamma_n"/> is the defined contact interface
+</div><p>The parameter <img class="math" src="_images/math/e8082d4f9b7dfa49370bbf2cebb770d246c66f7d.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
+contact interface, and <img class="math" src="_images/math/ae2c9cf2d489e787ac1340cc4436f6100a3449a3.png" alt="\gamma_n"/> is the defined contact interface
viscosity, also in the normal direction. The loss of energy in this interaction
-due to the viscous component is for particle <img class="math" src="_images/math/a581f053bbfa5115f42c13094857cdd12a37ec49.png" alt="i"/> calculated as:</p>
+due to the viscous component is for particle <img class="math" src="_images/math/df0deb143e5ac127f00bd248ee8001ecae572adc.png" alt="i"/> calculated as:</p>
<div class="math">
-<p><img src="_images/math/4305b4408a39583be7a127c3039d65ff0cd8c954.png" alt="\dot{e}^i_v = \gamma_n (\dot{\delta}^{ij}_n)^2"/></p>
+<p><img src="_images/math/89eab38386c7248f117087bc9f112b43b2a60010.png" alt="\dot{e}^i_v = \gamma_n (\dot{\delta}^{ij}_n)^2"/></p>
</div><p>The tangential force is determined by either a viscous-frictional contact model,
or a elastic-viscous-frictional contact model. The former contact model is very
computationally efficient, but somewhat inaccurate relative to the mechanics of
t@@ -157,24 +157,24 @@ real materials. The latter contact model is therefore the default, even though
it results in longer computational times. The tangential force in the
visco-frictional contact model:</p>
<div class="math">
-<p><img src="_images/math/f808e880e30ba99dc461a8dfdaaca903e8433bc1.png" alt="\boldsymbol{f}_t^{ij} = -\gamma_t \dot{\boldsymbol{\delta}_t}^{ij}"/></p>
-</div><p><img class="math" src="_images/math/3640b174f15a44d33ad8f0622bc44dce8d8f1692.png" alt="\gamma_n"/> is the defined contact interface viscosity in the tangential
+<p><img src="_images/math/6fa02b1c3bc8a438f685b666895f7de1f427b702.png" alt="\boldsymbol{f}_t^{ij} = -\gamma_t \dot{\boldsymbol{\delta}_t}^{ij}"/></p>
+</div><p><img class="math" src="_images/math/ae2c9cf2d489e787ac1340cc4436f6100a3449a3.png" alt="\gamma_n"/> is the defined contact interface viscosity in the tangential
direction. The tangential displacement along the contact interface
-(<img class="math" src="_images/math/8cfa62046b3a57f37068f027941acea702d9bd8a.png" alt="\boldsymbol{\delta}_t"/>) is not calculated and stored for this contact
+(<img class="math" src="_images/math/359991b87d5889ece1cd8b1c84cf17fbb37b815d.png" alt="\boldsymbol{\delta}_t"/>) is not calculated and stored for this contact
model. The tangential force in the more realistic elastic-viscous-frictional
contact model:</p>
<div class="math">
-<p><img src="_images/math/e3d643e0cb551cafec2614062e50e2369d7949ab.png" alt="\boldsymbol{f}_t^{ij} =
+<p><img src="_images/math/365556da97364e0ebe011deb8a6114cc57208b69.png" alt="\boldsymbol{f}_t^{ij} =
-k_t \boldsymbol{\delta}_t^{ij} -\gamma_t \dot{\boldsymbol{\delta}_t}^{ij}"/></p>
-</div><p>The parameter <img class="math" src="_images/math/66a1d8738af515286414e18d76a6e7f86b989afc.png" alt="k_n"/> is the defined spring coefficient in the tangential
+</div><p>The parameter <img class="math" src="_images/math/e8082d4f9b7dfa49370bbf2cebb770d246c66f7d.png" alt="k_n"/> is the defined spring coefficient in the tangential
direction of the contact interface. Note that the tangential force is only
-found if the tangential displacement (<img class="math" src="_images/math/d18b93994f644beb72edf4196d4efed3caf5b19c.png" alt="\delta_t"/>) or the tangential
-velocity (<img class="math" src="_images/math/6b038d0d06f6fa38bcb85df9e6f8bed4155ded47.png" alt="\dot{\delta}_t"/>) is non-zero, in order to avoid division by
-zero. Otherwise it is defined as being <img class="math" src="_images/math/b3cb421b6b325deff2bde26de2fa150f94c73246.png" alt="[0,0,0]"/>.</p>
+found if the tangential displacement (<img class="math" src="_images/math/182017e1dc60c069b4f10c93d4e58d03b6832ce0.png" alt="\delta_t"/>) or the tangential
+velocity (<img class="math" src="_images/math/3f1d577c0a6a09ed5b13001794297b2bd2ee9a28.png" alt="\dot{\delta}_t"/>) is non-zero, in order to avoid division by
+zero. Otherwise it is defined as being <img class="math" src="_images/math/6a4dbf76140ff4b044aef154739dc5b6c21301a2.png" alt="[0,0,0]"/>.</p>
<p>For both types of contact model, the tangential force is limited by the Coulomb
criterion of static and dynamic friction:</p>
<div class="math">
-<p><img src="_images/math/5cc160009025ab2cd3b77343dcea454fb8c93795.png" alt="||\boldsymbol{f}^{ij}_t|| \leq
+<p><img src="_images/math/edcde9b8cbb2394c6bbd12f07b845063be99a62e.png" alt="||\boldsymbol{f}^{ij}_t|| \leq
\begin{cases}
\mu_s ||\boldsymbol{f}^{ij}_n|| &
\text{if} \quad ||\boldsymbol{f}_t^{ij}|| = 0 \\
t@@ -185,7 +185,7 @@ criterion of static and dynamic friction:</p>
reached, the tangential displacement along the contact interface is limited to
this value:</p>
<div class="math">
-<p><img src="_images/math/13ac3fb2c8f5e95101d362b1cb57b5fcb70e9a71.png" alt="\boldsymbol{\delta}_t^{ij} =
+<p><img src="_images/math/092547c763f745333b005b7094cd03c2cfdd92f5.png" alt="\boldsymbol{\delta}_t^{ij} =
\frac{1}{k_t} \left(
\mu_d ||\boldsymbol{f}_n^{ij}||
\frac{\boldsymbol{f}^{ij}_t}{||\boldsymbol{f}^{ij}_t||}
t@@ -193,30 +193,30 @@ this value:</p>
</div><p>If the tangential force reaches the Coulomb limit, the energy lost due to
frictional dissipation is calculated as:</p>
<div class="math">
-<p><img src="_images/math/85bd2d6544d7812f80cc49d2b7e94be81541301c.png" alt="\dot{e}^i_s = \frac{||\boldsymbol{f}^{ij}_t
+<p><img src="_images/math/a699d5b66f6ff1048229d376ae6340ce0ed5ed6b.png" alt="\dot{e}^i_s = \frac{||\boldsymbol{f}^{ij}_t
\dot{\boldsymbol{\delta}}_t^{ij} \Delta t||}{\Delta t}"/></p>
</div><p>The loss of energy by viscous dissipation in the tangential direction is not
found.</p>
</div>
<div class="section" id="temporal-integration">
<h2>Temporal integration<a class="headerlink" href="#temporal-integration" title="Permalink to this headline">¶</a></h2>
-<p>In the DEM, the time is discretized into small steps (<img class="math" src="_images/math/5b4271afe7fc7c0ee84172e0ad19b82caf450c00.png" alt="\Delta t"/>). For each time
+<p>In the DEM, the time is discretized into small steps (<img class="math" src="_images/math/ec002955bdf95ee9869878fbad4f80fc98539359.png" alt="\Delta t"/>). For each time
step, the entire network of contacts is resolved, and the resulting forces and
torques for each particle are found. With these values at hand, the new
linear and rotational accelerations can be found using
<a class="reference external" href="https://en.wikipedia.org/wiki/Newton%27s_laws_of_motion">Newton’s second law</a>
-of the motion of solid bodies. If a particle with mass <img class="math" src="_images/math/c4bb40dd65eae6c11b325989b14e0b8d35e4e3ef.png" alt="m"/> at a point in time
-experiences a sum of forces denoted <img class="math" src="_images/math/8f28cef89dfa1989be5e263ac228a61a92736da9.png" alt="\boldsymbol{F}"/>, the resultant acceleration
-(<img class="math" src="_images/math/701df29c5482cdf1e577052c537c176f5a12caf5.png" alt="\boldsymbol{a}"/>) can be found by rearranging Newton’s second law:</p>
+of the motion of solid bodies. If a particle with mass <img class="math" src="_images/math/edba97b4c0d864d26e92ea7ea73767fa38eef3f7.png" alt="m"/> at a point in time
+experiences a sum of forces denoted <img class="math" src="_images/math/e5861e05d72d214bf2fa78e906e8a73638ee8837.png" alt="\boldsymbol{F}"/>, the resultant acceleration
+(<img class="math" src="_images/math/1a238040fa47ba4651506e01c721eb94a03f11ba.png" alt="\boldsymbol{a}"/>) can be found by rearranging Newton’s second law:</p>
<div class="math">
-<p><img src="_images/math/5bab09e6e7765565f55e2f7119251f199f472a4a.png" alt="\boldsymbol{F} = m \boldsymbol{a} \Rightarrow \boldsymbol{a} = \frac{\boldsymbol{F}}{m}"/></p>
+<p><img src="_images/math/7a16f81eae500f9872de54d88258152bfb0ccbc6.png" alt="\boldsymbol{F} = m \boldsymbol{a} \Rightarrow \boldsymbol{a} = \frac{\boldsymbol{F}}{m}"/></p>
</div><p>The new velocity and position is found by integrating the above equation
with regards to time. The simplest integration scheme in this regard is the
<a class="reference external" href="https://en.wikipedia.org/wiki/Euler_method">Euler method</a>:</p>
<div class="math">
-<p><img src="_images/math/a038fb9a65b9138665b4a5c9d0bbba6bdf3b6dbc.png" alt="\boldsymbol{v} = \boldsymbol{v}_{old} + \boldsymbol{a} \Delta t"/></p>
+<p><img src="_images/math/3d4d7c98c0c7511ac5480967994253f526e13af6.png" alt="\boldsymbol{v} = \boldsymbol{v}_{old} + \boldsymbol{a} \Delta t"/></p>
</div><div class="math">
-<p><img src="_images/math/ef4b452d768af650e6ddc2fd95aff430284f74f2.png" alt="\boldsymbol{p} = \boldsymbol{p}_{old} + \boldsymbol{v} \Delta t"/></p>
+<p><img src="_images/math/91bab0359c2bff8094d2fed7299dc2ca80a27d74.png" alt="\boldsymbol{p} = \boldsymbol{p}_{old} + \boldsymbol{v} \Delta t"/></p>
</div></div>
</div>
t@@ -278,7 +278,7 @@ with regards to time. The simplest integration scheme in this regard is the
<li class="right" >
<a href="introduction.html" title="Introduction and Installation"
>previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/genindex.html b/doc/html/genindex.html
t@@ -7,7 +7,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>Index — sphere 1.00-alpha documentation</title>
+ <title>Index — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -15,7 +15,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -26,7 +26,7 @@
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="index" title="Index" href="#" />
<link rel="search" title="Search" href="search.html" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="index.html" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="index.html" />
</head>
<body role="document">
<div class="related" role="navigation" aria-label="related navigation">
t@@ -38,7 +38,7 @@
<li class="right" >
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -959,7 +959,7 @@
<li class="right" >
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/index.html b/doc/html/index.html
t@@ -6,7 +6,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>The sphere documentation — sphere 1.00-alpha documentation</title>
+ <title>The sphere documentation — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -14,7 +14,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -25,7 +25,7 @@
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="#" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="#" />
<link rel="next" title="Introduction and Installation" href="introduction.html" />
</head>
<body role="document">
t@@ -41,7 +41,7 @@
<li class="right" >
<a href="introduction.html" title="Introduction and Installation"
accesskey="N">next</a> |</li>
- <li class="nav-item nav-item-0"><a href="#">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="#">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -161,7 +161,7 @@ party developers. This document is a work in progress.</p>
<li class="right" >
<a href="introduction.html" title="Introduction and Installation"
>next</a> |</li>
- <li class="nav-item nav-item-0"><a href="#">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="#">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/introduction.html b/doc/html/introduction.html
t@@ -6,7 +6,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>Introduction and Installation — sphere 1.00-alpha documentation</title>
+ <title>Introduction and Installation — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -14,7 +14,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -25,7 +25,7 @@
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="index.html" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="index.html" />
<link rel="next" title="Discrete element method" href="dem.html" />
<link rel="prev" title="The sphere documentation" href="index.html" />
</head>
t@@ -45,7 +45,7 @@
<li class="right" >
<a href="index.html" title="The sphere documentation"
accesskey="P">previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -185,9 +185,11 @@ from the root directory:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span>$ cmake . && make
</pre></div>
</div>
-<p>NOTE: If your system does not have a GCC compiler (e.g. GCC-5 for CUDA 8)
-compatible with the installed CUDA version, try using <code class="docutils literal"><span class="pre">clang-3.8</span></code> instead:</p>
-<div class="highlight-default"><div class="highlight"><pre><span></span>$ export CC=$(which clang-3.8) && export CXX=$(which clang++-3.8) && cmake . && make
+<p>NOTE: If your system does not have a GCC compiler compatible with the installed
+CUDA version (e.g. GCC-5 for CUDA 8), you will see errors at the linker stage.
+In that case, try using <code class="docutils literal"><span class="pre">clang-3.8</span></code> as the C and C++ compiler instead:</p>
+<div class="highlight-default"><div class="highlight"><pre><span></span>$ rm -rf CMakeCache.txt CMakeFiles/
+$ export CC=$(which clang-3.8) && export CXX=$(which clang++-3.8) && cmake . && make
</pre></div>
</div>
<p>After a successfull installation, the <code class="docutils literal"><span class="pre">sphere</span></code> executable will be located
t@@ -216,13 +218,9 @@ following commands to check the executable:</p>
| | | |
| |_| Version: 2.15 |
`-------------------------------------´
- A discrete element method particle dynamics simulator.
+ A discrete-element method particle dynamics simulator.
Written by Anders Damsgaard, license GPLv3+.
- https://cs.au.dk/~adc/
-</pre></div>
-</div>
-<p>The build can be verified by running a number of automated tests:</p>
-<div class="highlight-default"><div class="highlight"><pre><span></span>$ make test
+ https://adamsgaard.dk
</pre></div>
</div>
<p>The documentation can be read in the <a class="reference external" href="http://docutils.sourceforge.net/docs/ref/rst/restructuredtext.html">reStructuredText</a>-format in
t@@ -254,7 +252,7 @@ a simulation is typically arranged in the following order:</p>
<blockquote>
<div><ul class="simple">
<li>Setup of particle assemblage, physical properties and conditions using the
-Python API.</li>
+Python API (<code class="docutils literal"><span class="pre">python/sphere.py</span></code>).</li>
<li>Execution of <code class="docutils literal"><span class="pre">sphere</span></code> software, which simulates the particle behavior as a
function of time, as a result of the conditions initially specified in the
input file.</li>
t@@ -325,7 +323,7 @@ in Python, and/or scene rendering using the built-in ray tracer.</li>
<li class="right" >
<a href="index.html" title="The sphere documentation"
>previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/objects.inv b/doc/html/objects.inv
Binary files differ.
DIR diff --git a/doc/html/py-modindex.html b/doc/html/py-modindex.html
t@@ -6,7 +6,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>Python Module Index — sphere 1.00-alpha documentation</title>
+ <title>Python Module Index — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -14,7 +14,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -25,7 +25,7 @@
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="index.html" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="index.html" />
<script type="text/javascript">
t@@ -44,7 +44,7 @@
<li class="right" >
<a href="#" title="Python Module Index"
>modules</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -100,7 +100,7 @@
<li class="right" >
<a href="#" title="Python Module Index"
>modules</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/python_api.html b/doc/html/python_api.html
t@@ -6,7 +6,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>Python API — sphere 1.00-alpha documentation</title>
+ <title>Python API — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -14,7 +14,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -25,7 +25,8 @@
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="index.html" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="index.html" />
+ <link rel="next" title="sphere internals" href="sphere_internals.html" />
<link rel="prev" title="Fluid simulation and particle-fluid interaction" href="cfd.html" />
</head>
<body role="document">
t@@ -39,9 +40,12 @@
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
+ <a href="sphere_internals.html" title="sphere internals"
+ accesskey="N">next</a> |</li>
+ <li class="right" >
<a href="cfd.html" title="Fluid simulation and particle-fluid interaction"
accesskey="P">previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -1154,7 +1158,7 @@ parameter <cite>hydrostatic</cite> is set to <cite>True</cite>, this value will
the fluid cells at the top</li>
<li><strong>hydrostatic</strong> (<em>bool</em>) – Initialize the fluid pressures to the hydrostatic
pressure distribution. A pressure gradient with depth is only
-created if a gravitational acceleration along <img class="math" src="_images/math/84d7271dd9e78c1e05d6c3c6ecb60309ef7dfc73.png" alt="z"/> previously
+created if a gravitational acceleration along <img class="math" src="_images/math/683f2dd9129a91d21aaf1c04afa6f78b39d4cb0a.png" alt="z"/> previously
has been specified</li>
<li><strong>cfd_solver</strong> (<em>int</em>) – Solver to use for the computational fluid dynamics.
Accepted values: 0 (Navier Stokes, default) and 1 (Darcy).</li>
t@@ -3376,6 +3380,9 @@ rendered beforehand using <a class="reference internal" href="#sphere.render" ti
<h4>Previous topic</h4>
<p class="topless"><a href="cfd.html"
title="previous chapter">Fluid simulation and particle-fluid interaction</a></p>
+ <h4>Next topic</h4>
+ <p class="topless"><a href="sphere_internals.html"
+ title="next chapter">sphere internals</a></p>
<div role="note" aria-label="source link">
<h3>This Page</h3>
<ul class="this-page-menu">
t@@ -3407,9 +3414,12 @@ rendered beforehand using <a class="reference internal" href="#sphere.render" ti
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
+ <a href="sphere_internals.html" title="sphere internals"
+ >next</a> |</li>
+ <li class="right" >
<a href="cfd.html" title="Fluid simulation and particle-fluid interaction"
>previous</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/search.html b/doc/html/search.html
t@@ -6,7 +6,7 @@
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
- <title>Search — sphere 1.00-alpha documentation</title>
+ <title>Search — sphere 2.15-beta documentation</title>
<link rel="stylesheet" href="_static/classic.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
t@@ -14,7 +14,7 @@
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: './',
- VERSION: '1.00-alpha',
+ VERSION: '2.15-beta',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
t@@ -26,7 +26,7 @@
<script type="text/javascript" src="_static/searchtools.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="#" />
- <link rel="top" title="sphere 1.00-alpha documentation" href="index.html" />
+ <link rel="top" title="sphere 2.15-beta documentation" href="index.html" />
<script type="text/javascript">
jQuery(function() { Search.loadIndex("searchindex.js"); });
</script>
t@@ -45,7 +45,7 @@
<li class="right" >
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
t@@ -96,7 +96,7 @@
<li class="right" >
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
- <li class="nav-item nav-item-0"><a href="index.html">sphere 1.00-alpha documentation</a> »</li>
+ <li class="nav-item nav-item-0"><a href="index.html">sphere 2.15-beta documentation</a> »</li>
</ul>
</div>
<div class="footer" role="contentinfo">
DIR diff --git a/doc/html/searchindex.js b/doc/html/searchindex.js
t@@ -1 +1 @@
mx1.adamsgaard.dk:70 /src/sphere/commit/81d7b44047babc963125026004031b612e6de13d.gph:1177: line too long