yade-users team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

Hi,

One of the reviews of my next DEM paper was quite critical, and the
reviewer asked for all equations used in calculation to be included
in the paper. That makes a hefty lot of equations, I wonder if he
will change his mind when he will see 42 equations on 6 pages.

Anyway I had to write down everything that is in 

CohesiveFrictionalRelationships.cpp 
CohesiveFrictionalContactLaw.cpp

and other classes. In the process I saw that nearly 70% of code can
be plain deleted and is unused in real calculation. Perhaps it's a
leftover from some older version, or some sub-results used for
debugging are never actually used. I should be removing that too.

anyway I'm sending you all those equations. Maybe someone will find
time to put them on wiki. WARNING: But which one wiki? Now we have
two clones that are independently edited, and are getting out of sync!

Also, if you care to read it carefully, maybe you will find a
mistake! I could forget a \vec{} here or a "+" sign there... If you
find a mistake, please let me know.

best regards
-- 
Janek Kozicki                                                         |

Attachment: DEM_equations.pdf
Description: Adobe PDF document

\documentclass[11pt,a4paper]{article}

\usepackage{graphicx}
\usepackage{a4wide}
\def\quat#1{{\mathring #1}}%
\newcommand{\D}{\displaystyle}%
\def\unit#1{\textrm{\hskip 1pt#1}}

\title{Discrete Element Method}
%\author{The authors name should be here}
\begin{document}
%\maketitle
\section*{Discrete Element Method}

Symbols used in this chapter:

\begin{description}
\item[$\vec{X_A}$, $\vec{X_B}$] --- current position of sphere A and sphere B (a three dimensional vector referring to global coordinate system),
\item[$\vec{V_A}$, $\vec{V_B}$] --- current velocity of sphere A and sphere B,
\item[$\vec{\omega_A}$, $\vec{\omega_B}$] --- current angular velocity of sphere A and sphere B,
\item[$\quat{A}$, $\quat{B}$] --- current orientation of sphere A and sphere B (a unit quaternion referring to global coordinate system),
\item[$\quat{A'}$, $\quat{B'}$] --- orientation of sphere A and sphere B when the contact was established (an initial orientation of contacting spheres, a unit quaternion),
\item[$R_A$, $R_B$] -- radius of sphere A and sphere B,
\item[$E_A$, $E_B$] -- Young's modulus for sphere A and B,
\item[$\nu_A$, $\nu_B$] --- Poisson's ratio for sphere A and B,
\item[$\mu_A$, $\mu_B$] --- friction angle for sphere A and B,
\item[$K_n$] --- normal stiffness     of the contact between two spheres,
\item[$K_s$] --- tangential stiffness of the contact between two spheres,
\item[$K_r$] --- rotational stiffness of the contact between two spheres,
\item[$\beta$] --- dimensionless coefficient for rotational stiffness,
\item[$\mu$] --- friction angle between two spheres,
\item[$\eta$] --- elastic limit of rolling
\item[$\Delta\vec{X}$] --- distance between centers of two spheres,
\item[$\Delta\vec{V}$] --- relative velocity between two spheres,
\item[$\Delta\vec{V_s}$] --- tangential velocity between two spheres,
\item[$\Delta\vec{X_s}$] --- tangential displacement between two spheres,
\item[$\Delta\quat{Q}$] --- orientation displacement between two spheres (a unit quaternion),
\item[$\Delta\vec{\omega}$] --- angular rotation between two spheres (a three dimensional vector),
\item[$U$] --- penetration depth between two spheres,
\item[$\vec{N}$] --- normal of the contact (a vector of unit length),
\item[$\vec{C}$] --- position of the contact point between two spheres (a three dimensional vector referring to global coordinate system),
\item[$\vec{F_n}$] --- normal contact force,
\item[$\vec{F_s}$] --- tangential contact force,
\item[$\vec{M}$] --- contact moment,
\item[$\vec{F_A}$, $\vec{F_B}$] --- force acting on spheres A and B,
\item[$\vec{M_A}$, $\vec{M_B}$] --- moment acting on spheres A and B,
\item[$\alpha$] --- numerical damping coefficient,
\item[$\Delta t$] --- time increment between two iterations.
\end{description}

First the distance between spheres A and B is calculated:
%Example creating equations, see equation~\ref{eq_sum}:
\begin{equation}
\Delta\vec{X} = \vec{X_A} - \vec{X_B}
%\label{eq_sum}
\end{equation}

The penetration depth of the contact is calculated:

\begin{equation}
U = R_A + R_B - ||\Delta\vec{X}||
\end{equation}

This paper does not treat about cohesive granular material therefore negative penetration depth means that there is no contact. If the penetration depth is positive ($U > 0$) then the contact between two spheres exists and the calculation can proceed further. The normal of the contact is calculated:
\begin{equation}
\vec{N} = \frac{\Delta\vec{X}}{||\Delta\vec{X}||}
\end{equation}

The contact point $\vec{C}$ is calculated the center of the overlapping volume of two spheres:
\begin{equation}
\vec{C} = \vec{X_A} + (R_A-\frac{U}{2})\vec{N}
\end{equation}

If the contact was not present in previous calculation steps, the initial orientation of both contacting spheres is stored in memory:

\begin{equation}
\quat{A'} = \quat{A}
\end{equation}
\begin{equation}
\quat{B'} = \quat{B}
\end{equation}

The normal stiffness $K_n$ of the contact is calculated as a harmonic average between two Young's modulus, using sphere's radius:

\begin{equation}
K_n = 2 \frac{E_A R_A E_B R_B}{E_A R_A + E_B R_B}
\end{equation}

The tangential stiffness $K_s$ of the contact is calculated as a harmonic average using Poisson's ratio for each sphere:

\begin{equation}
K_s = 2 \frac{E_A R_A \nu_A E_B R_B \nu_B}{E_A R_A \nu_A + E_B R_B \nu_B}
\end{equation}

The rotational stiffness $K_r$ is calculated using both sphere's radius and a dimensionless coefficient $\beta$:

\begin{equation}
K_r = \beta R_A R_B K_s
\end{equation}

The friction angle for the contact between two spheres is taken as the smaller one from both spheres:

\begin{equation}
\mu = \min(\mu_A,\mu_B)
\end{equation}

The normal contact force $\vec{F_n}$ is calculated as:

\begin{equation}
\vec{F_n} = K_n U \vec{N}
\end{equation}

If the contact was not present in previous calculation steps, the initial value of tangential force is set to zero:

\begin{equation}
\vec{F_s} = 0
\end{equation}

The increment of tangential force is calculated by accumulating the path of moving contact point between two spheres. First the relative velocity is calculated. This approach was proposed in (F. Alonso-Marroquin, R. Garcia-Rojo, H.J. Herrmann "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena" ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10):

\begin{equation}
\Delta\vec{V} = (\vec{V_B} + \vec{\omega_B} \times (- R_B \vec{N})) - (\vec{V_A} + \vec{\omega_A} \times (R_A \vec{N}))
\end{equation}

Next the tangential velocity between two spheres is calculated:

\begin{equation}
\Delta\vec{V_s} = \Delta\vec{V} - (\vec{N}\cdot\Delta\vec{V})\vec{N}
\end{equation}

The tangential displacement between two spheres is calculated:

\begin{equation}
\Delta\vec{X_s} = \Delta t \Delta\vec{V_s}
\end{equation}

Then the shear force $\vec{F_s}$ is updated using current tangential displacement and tangential stiffness:

\begin{equation}
\vec{F_s} = \vec{F_s} + K_s \Delta\vec{X_s}
\end{equation}

%A simple Coulomb criterion $||F_s|| - ||F_n||\tan{\mu}\le0$ is enforced with following formula:
%\begin{equation}
%\texttt{if not:~} \left(||F_s|| - ||F_n||\tan{\mu}\le0\right) \texttt{~~then:~} \left(\vec{F_s} = \vec{F_s} \frac{||F_n||\tan{\mu}}{||F_s||}\right)
%\end{equation}

If Coulomb criterion:
\begin{equation}
||\vec{F_s}|| - ||\vec{F_n}||\tan{\mu}\le0,
\end{equation}
is not satisfied, then the shearing force $\vec{F_s}$ needs to be corrected down, to satisfy it:
\begin{equation}
\vec{F_s} = \vec{F_s} \frac{||\vec{F_n}||\tan{\mu}}{||\vec{F_s}||}
\end{equation}

The forces acting on spheres A and B are accumulated:

\begin{equation}
\vec{F_A} = \vec{F_A} - (\vec{F_n} + \vec{F_s})
\end{equation}
\begin{equation}
\vec{F_B} = \vec{F_B} + (\vec{F_n} + \vec{F_s})
\end{equation}

The moments acting on spheres A and B, due to $F_s$ and $F_n$ are accumulated:

\begin{equation}
\vec{M_A} = \vec{M_A} - (\vec{C} - \vec{X_A})\times(\vec{F_n} + \vec{F_s})
\end{equation}
\begin{equation}
\vec{M_B} = \vec{M_B} + (\vec{C} - \vec{X_B})\times(\vec{F_n} + \vec{F_s})
\end{equation}

The angular orientation displacement $\Delta\quat{Q}$ between two spheres is calculated from original relative orientation of two spheres and current relative orientation:

\begin{equation}
\Delta\quat{Q} = \quat{A}~\left(\quat{A'}\right)^{-1}~\quat{B'}~\left(\quat{B}\right)^{-1}
\end{equation}

Next the quaternion $\Delta\quat{Q} = a + b \texttt{i} + c \texttt{j} + d \texttt{k}$ is converted into a three dimensional rotation vector $\Delta\vec{\omega}=\{x,y,z\}$ using the following formula:
\begin{equation}
\Delta\vec{\omega}=\left\{
\begin{array}{l}
x = \phi \D\frac{b}{\sin\left(\phi /2\right)}\\[5mm]
y = \phi \D\frac{c}{\sin\left(\phi /2\right)}\\[5mm]
z = \phi \D\frac{d}{\sin\left(\phi /2\right)}
\end{array}
\right.,\textsl{~where:~}\phi = 2\arccos\left(a\right).
\end{equation}
If the rotation angle $\phi=0$, then axis of the rotation can be anything since no rotation occurs. In this case to avoid the division by zero, the zeros are assigned to $x$, $y$ and $z$. Then the contact moment $\vec{M}$ between two spheres is calculated as:

\begin{equation}
\vec{M} = K_r \Delta\vec{\omega}
\end{equation}

In this paper the rotational stiffness for twist and bending is the same value $K_r$, using two different values can be obtained by decomposing $\Delta\vec{\omega}$ into two vectors, one lying on the plane of contact second being in the direction of the normal of the contact.

If following criterion:
\begin{equation}
||\vec{M}|| - \eta \frac{R_A+R_B}{2} ||\vec{F_n}||\le0,
\end{equation}
is not satisfied, then the moment $\vec{M}$ needs to be corrected down, to satisfy it:
\begin{equation}
\vec{M} = \vec{M} \frac{\eta \frac{R_A+R_B}{2} ||\vec{F_n}||}{||\vec{M}||}
\end{equation}


The moments acting on spheres A and B, due to contact moment are accumulated:

\begin{equation}
\vec{M_A} = \vec{M_A} - \vec{M}
\end{equation}
\begin{equation}
\vec{M_B} = \vec{M_B} + \vec{M}
\end{equation}

After all contacts are processed using above equations, the gravity effect is applied on all spheres:
\begin{equation}
\vec{F_i} = \vec{F_i} + m_i g,
\end{equation}
where $g=9.81\unit{m/s}^2$, $m$ mass of $i$-th sphere and $F_i$ is residual force vector for sphere $i$-th. 

In the model, a local non--viscous damping scheme is adopted to dissipate kinetic energy (Cundall and Hart 1992):
\begin{equation}
\vec{F^k_i} = \vec{F^k_i} - \alpha~\textsl{sgn}(\vec{V}^k_i)|F^k_i|,
\end{equation}
where $F^k_i$ is the $k$-th component of the residual force vector for sphere $i$-th, the $\vec{V}^k_i$ is the $k$-th component of velocity $\vec{V}_i$ of $i$-th sphere, $\alpha$ is a positive numerical damping coefficient less than 1 and sgn$(\bullet)$ returns the sign of argument $\bullet$. The same is applied for moment:
\begin{equation}
\vec{M^k_i} = \vec{M^k_i} - \alpha~\textsl{sgn}(\vec{\omega}^k_i)|M^k_i|,
\end{equation}
where $\vec{M^k_i}$ is the $k$-th component of the residual moment vector for sphere $i$-th and $\vec{\omega}^k_i$ is the $k$-th component of angular velocity of $i$-th sphere.

Then a Newton's force law is applied to all spheres to calculate their acceleration:

\begin{equation}
\vec{a_i} = \vec{F_i}/m_i,
\end{equation}

where $\vec{a_i}$ is acceleration of $i$-th sphere, $\vec{F_i}$ is force acting on it, and $m_i$ is its mass. The same is applied for moment to calculate angular acceleration:

\begin{equation}
\vec{\dot{\omega}^k_i} = \vec{M^k_i}/\vec{I^k_i},
\end{equation}

where $\vec{\dot{\omega}^k_i}$ is $k$-th component of angular acceleration of $i$-th sphere, $\vec{M^k_i}$ is $k$-th component of moment acting on $i$-th sphere and $\vec{I^k_i}$ is $k$-th component of inertia of an $i$-th sphere. Inertial matrix here is expressed as a vector containing only diagonal components of inertial matrix, because the remaining cells in the matrix are zero, for a sphere.

Then the position of all spheres is integrated using a leap from integration scheme:

\begin{equation}
\vec{V_i} = \vec{V_i} + \Delta t~\vec{a_i},
\end{equation}
\begin{equation}
\vec{X_i} = \vec{X_i} + \Delta t~\vec{V_i}.
\end{equation}

And the orientation of all spheres is integrated in similar way, using quaternions. First an angular velocity is calculated:

\begin{equation}
\vec{\omega_i} = \vec{\omega_i} + \Delta t~\vec{\dot{\omega}_i}.
\end{equation}

Then a quaternion $\Delta\quat{Q_i}= a + b \texttt{i} + c \texttt{j} + d \texttt{k}$ representing rotation at time step $\Delta t$ is constructed, assuming that $\vec{\omega_i} = \{x,y,z\}$:

\begin{equation}
\Delta\quat{Q_i}=\left\{
\begin{array}{l}
a = \cos\left(\phi /2\right)    \\[2mm]
b = \sin\left(\phi /2\right) \D\frac{x}{||\vec{\omega_i}||}  \\[4mm]
c = \sin\left(\phi /2\right) \D\frac{y}{||\vec{\omega_i}||}  \\[4mm]
d = \sin\left(\phi /2\right) \D\frac{z}{||\vec{\omega_i}||}
\end{array}
\right.,\textsl{~where:~}\phi = ||\vec{\omega_i}||.
\end{equation}

Then the orientation of $i$-th sphere is updated:

\begin{equation}
\quat{Q_i} = \quat{Q_i}~\Delta\quat{Q_i}.
\end{equation}

Before the calculation of an iteration step starts the residual forces and moments acting on spheres are set to zero:

\begin{equation}
\vec{F_i} = 0.
\end{equation}
\begin{equation}
\vec{M_i} = 0.
\end{equation}

Then the time is incremented:

\begin{equation}
t = t+\Delta t,
\end{equation}

and next calculation step starts.

\end{document}