Normal contact - gnomeCreative/HYBIRD GitHub Wiki

Normal contact

For the normal collision, two models have been implemented: the damped Hertzian model, and the simpler linear-dashpot model. The model choice depends on the desired level of accuracy, and can be controlled by the user with the command contactModel

The linear dashpot

The linear dashpot is one of the earliest developed models, and idealizes the contact as a parallel connection with a spring of stiffness $k_\textrm{n,L}$ and a damper with viscous damping coefficient $\alpha_\textrm{n,L}$. The normal collision force acts in the direction $n$ normal to the contact surface and is expressed as

$F_\textrm{n,coll}=k_\textrm{n,L} \xi + 2\alpha_\textrm{n,L}\sqrt{k_\textrm{n,L} \tilde{m}} \dot{\xi}$,

with $\tilde{m}$ being the effective mass defined as $\tilde{m}=m_\textrm{p}m_\textrm{q}/(m_\textrm{p}+m_\textrm{q})$.

Representation of the DEM and of the contact law.

The overlap obeys Newton's law in the form of a harmonic oscillator

$\ddot{\xi}+2\gamma_\textrm{L}\dot{\xi}+\omega_\textrm{L}^2 \xi=0$,

where $\gamma_\textrm{L}= \alpha_\textrm{n,L}\sqrt{k_\textrm{n,L} /\tilde{m}}$ is the damping ratio, and $\omega_\textrm{L}=\sqrt{k_\textrm{n,L}/\tilde{m}}$ is the frequency of the undamped oscillator. The restitution coefficient $\zeta$, defined as the ratio between the velocity of the particles before and after contact, is constant in this model and is expressed as

$\zeta=\exp { \left( -\frac{\pi \gamma_\textrm{L}}{\sqrt{\omega_\textrm{L}^2-\gamma_\textrm{L}^2}} \right)}$,

if the system is not overdamped. This behavior is unphysical, since experiments~\cite{Kuwabara1987} showed that the coefficient of restitution changes with the impact velocity, albeit not dramatically. The idealization of a constant coefficient of restitution is however desirable from a practical point of view. The coefficient of restitution can be easily estimated with simple experiments, providing access to the determination of $\gamma_\textrm{n,L}$ and $\alpha_\textrm{n,L}$. The time of contact $t_{\textrm{coll},\textrm{L}}$ is also constant in this model, and is inversely proportional to the oscillator frequency

$t_\textrm{coll,L}=\frac{\pi}{\sqrt{\omega_\textrm{L}^2-\gamma_\textrm{L}^2}}$.

In a DEM simulation the time step is chosen so that contacts are represented with an acceptable degree of accuracy. This is usually satisfied by imposing a time step $\Delta t^\textrm{DEM}\lt \frac{1}{10} t_\textrm{coll,L}$. A constant contact time means that the accuracy of the simulation can be easily controlled, which is a favorable property of this model.

The Hertzian model

The linear-dashpot model presented in the last section does not take correctly into account how the deformation $\xi$ is a function of the shape of the particle. If these are spheres, the contact law should be modified according to the theory of Hertz (1894), which has the advantage of being expressed in terms of the physical properties of the material, namely Young's modulus $E_\textrm{p}$ and Poisson's ratio $\nu_\textrm{p}$. The normal force reads:

$F_\textrm{n,coll}=k_\textrm{n,H}\xi+2 \alpha_\textrm{n,H}\sqrt{k_\textrm{n,H} \tilde{m}} \frac{d\xi}{dt}$,

where the tangent stiffness $k_\textrm{n,H}$ is now a function of the overlap, yielding a non-linear behavior. Its expression is

$k_\textrm{n,H}= \frac{2}{3} \frac{E_\textrm{p}}{1-\nu_\textrm{p}^2} \sqrt{\tilde{r}} \xi^{1/2}$,

where $\tilde{r}=r_\textrm{p1}r_\textrm{p2}/(r_\textrm{p1}+r_\textrm{p2})$ is the effective radius. This equation slightly differs from the canonical Hertz theory. It is a model developed by Tsuji et al. (1992), with the purpose of obtaining a pseudo-Hertzian model characterized by a constant coefficient of restitution. While being not completely sound, this approach is practical since the estimation of the coefficient of restitution $\zeta$ is usually an easy task. The coefficient $\alpha_\textrm{n,H}$ can then be determined from $\zeta$ (Antypov & Elliott, 2011) as

$\alpha_\textrm{n,H}=\frac{-\sqrt{5}\ln{\zeta}}{\sqrt{\ln^2{\zeta}+\pi^2}}$,

In contrast to the linear dashpot, this model however does not provide a simple way to estimate the time step. The collision time, as a matter of fact, is a function of the normal collision velocity, as

$t_\textrm{coll,H} \simeq 1.1 \left(\frac{E_\textrm{p}}{1-\nu_\textrm{p}^2} \rho_\textrm{p}\right)^{2/5} \frac{\tilde{r}}{(u_\textrm{n,coll})^{1/5}}$.

The DEM time step can then be adapted according to the principle $\Delta t^\textrm{DEM}< \frac{1}{10} t_\textrm{coll,H}$ by an evaluation of the maximum particle velocity in the system. Alternatively, the maximum velocity can be predicted in the beginning. This allows the use of a constant time step, a simplification that comes at the expenses of the overall performance of the method.