yade-users team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

Janek Kozicki said:     (by the date of Thu, 8 Oct 2009 21:31:47 +0200)

> (oops forgot attachment)
> 
> I'm sending the latest update, with only few corrections:
> 
> - forgot _i index for mass of i-th sphere
> 
> - added V(t+dt/2), V(t-dt/2) clarifications in the leap frog integration
>   part. Please alert me if I made some stupid mistake, which is of
>   course possible ;) Thanks to Vincent about his remarks about this part :)

- and one more correction: added Q(t+dt) notation to orientation, just below.

sorry for too much posts with this :)
-- 
Janek Kozicki                                                         |

Attachment: DEM_equations_3.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}

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 $\vec{F_s}$ and $\vec{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 $\vec{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)|\vec{F^k_i}|,
\end{equation}
where $\vec{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)|\vec{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}/\left(\vec{I^k_i}~m_i\right),
\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 \textit{vector} of an $i$-th sphere. Inertial matrix here is expressed as a vector containing only diagonal components of inertial matrix (equal to $\frac{2}{5}R_i^2$), because the remaining cells in the inertial matrix are zero, for a sphere.

Then the position of all spheres is integrated using a leap frog integration scheme (where $\vec{V}_{i,t-\frac{\Delta t}{2}} + \Delta t~\vec{a}_{i,t}$ is a $2^{nd}$ order approximation of $\vec{V}_{i,t+\frac{\Delta t}{2}}$):

\begin{equation}
\vec{V}_{i,t+\frac{\Delta t}{2}} = \vec{V}_{i,t-\frac{\Delta t}{2}} + \Delta t~\vec{a}_{i,t},
\end{equation}
\begin{equation}
\vec{X}_{i,t+\Delta t} = \vec{X}_{i,t} + \Delta t~\vec{V}_{i,t+\frac{\Delta t}{2}}.
\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,t+\frac{\Delta t}{2}} = \vec{\omega}_{i,t-\frac{\Delta t}{2}} + \Delta t~\vec{\dot{\omega}}_{i,t}.
\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\}$ (where $\vec{\omega}_{i,t-\frac{\Delta t}{2}} + \Delta t~\vec{\dot{\omega}}_{i,t}$ is a $2^{nd}$ order approximation of $\vec{\omega}_{i,t+\frac{\Delta t}{2}}$):

\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}{\phi}  \\[4mm]
c = \sin\left(\phi /2\right) \D\frac{y}{\phi}  \\[4mm]
d = \sin\left(\phi /2\right) \D\frac{z}{\phi}
\end{array}
\right.,\textsl{~where:~}\phi = ||\vec{\omega}_{i,t+\frac{\Delta t}{2}}||.
\end{equation}

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

\begin{equation}
\quat{Q}_{i,t+\Delta t} = \quat{Q}_{i,t}~\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}