Marcus N. Bannerman
[email protected]
Each of the $N$ particles represents some portion of the system's mass.
The motion of each particle is described by Newton's law:
As analytical solutions are not usually available for
$N>2$ systems
(it will be demonstrated that they are
available for some forces)
numerical integration is often used.
In granular (dissipative) systems, standard explicit techniques such as Runge-Kutta (i.e., ode45) and Gear's predictor-corrector (i.e. ode15s) are used.
In molecular (elastic) systems, they are much more picky. Only Symplectic (phase volume conserving) techniques are used.
Consider the simplest particle system, the spring-mass system.
Its exact solution is well known. For $r(t=0)=1$ and $v(t=0)=0$, we have:
\begin{align*} r(t) &= \cos(\omega\,t) & v(t) &= \sin(\omega\,t)\\ \end{align*}
\begin{align*} E(t) = \frac{1}{2}m\,v^2 + \frac{1}{2}k\,r^2 = \frac{1}{2} \end{align*}
Integrating the same system using several numerical schemes of timestep $\Delta t$ gives the following exact results:
\begin{align*} E_\text{Euler}(t) &= \frac{1}{2}e^{\frac{\ln(1+\Delta t^2)}{\Delta t}t}\approx\frac{1}{2}e^{t\,\Delta t} \\ E_\text{RK2}(t) &= \frac{1}{2}e^{\frac{\ln(1+\Delta t^4/4)}{\Delta t}t}\approx\frac{1}{2}e^{t\,\frac{\Delta t^3}{4}} \\ E_\text{RK4}(t) &= \frac{1}{2}e^{\frac{\ln(1-\Delta t^6/72+\Delta t^8/576)}{\Delta t}t}\approx\frac{1}{2}e^{-t\,\frac{\Delta t^5}{72}} \\ E_\text{Vel.Verlet}(t) &= \frac{1}{2}+\frac{\Delta t^2}{4-\Delta t^2}\sin^2\left(\frac{\theta}{\Delta t} t\right) \\ \theta&= \tan^{-1}\left(\Delta t\frac{\sqrt{1-\Delta t^2/4}}{1-\Delta t^2/2}\right) \end{align*}
\begin{align*} E_\text{Vel.Verlet}(t) &= \frac{1}{2}+\frac{\Delta t^2}{4-\Delta t^2}\sin^2\left(\frac{\theta}{\Delta t} t\right) \end{align*}
Interestingly, the last integrator is stable, energy is roughly conserved at long times (as the integrator is symplectic).
It can be shown that the numerical system is in-fact following a new definition of the potential energy: \begin{align*} E_\text{Vel.Verlet}^{(potential)}(r(t)) = \frac{r^2}{2}+\frac{\Delta t^2}{4-\Delta t^2}\left(r^2-1\right) \end{align*}
This is why symplectic integrators are so stable and powerful, they are still conservative, and obey the dynamics of a related Hamiltonian.
\begin{align*} E_\text{Vel.Verlet}^{(potential)}(r(t)) = \frac{r^2}{2}+\frac{\Delta t^2}{4-\Delta t^2}r^2 \end{align*}
As $\Delta t\to2$, the approximate solution diverges from the exact solution, thus it may be stable for $\Delta t<2$ but it is increasingly innaccurate.
Even careful and modern granular studies^{1} use timesteps 80% of the critical limit (for speed), giving 178% more potential energy than the underlying model!
^{1}Hosn, R. A., Sibille, L., Benahmed, N., and Chareyre, B., "Discrete numerical modelling of loose soil with spherical particles and interparticle rolling friction", Gran. Matt. (2017), 19:4,
doi:10.1007/s10035-016-0687-0
This is compounded by the use of artificially soft particles to allow larger timesteps.
The time step is 1/20th of the collision time of an undamped Hertz model with $\nu=0.3$, $R=6.204$ mm, $\rho=7800$ kg m^{-3}, and impact velocity of 1 m s^{-1}.
Modern DEM simulations are largely unrelated to the underlying model parameters, which are indirectly related to the bulk properties.
This is out of necessity to keep simulations stable and fast.
Are there any universally stable approaches?
I.e., Are there any methods which are exact solutions
of their underlying models?
We are used to solving certain dynamics problems via exact solutions of approximate models.
Consider how you would solve this problem:
Is your solution exactly correct? No.
Your model is an approximation but the solution of the model is exact.
Event driven dynamics has a similar approach.
Rather than step through time, we calculate the time
of "events" exactly
(using an approximate model).
"Events" are impulses, which may represent soft conservative or dissipative interactions.
Building an approximate model with no drag, the exact solution to the projectile motion is a parabola: \begin{align*} r_y(t+\Delta t) = \color{aqua}{r_y(t)} + \color{lime}{v_y(t)}\,\Delta t+g_y\,\frac{\Delta t^2}{2} \end{align*}
\begin{align*} r_y(t+\Delta t) = \color{aqua}{r_y(t)} + \color{lime}{v_y(t)}\,\Delta t+g_y\,\frac{\Delta t^2}{2} \end{align*} When it hits the ground at a time $t+\Delta t$ we have: \begin{align*} r_y(t+\Delta t)=r_{ground} \end{align*}
To find the time of this "event", we are looking for the roots to the equation: \begin{align*} r_y - r_{ground} + v_y\,\Delta t+g_y\,\frac{\Delta t^2}{2} = 0 \end{align*}
In this case, this is a simple quadratic with two roots: \begin{align*} \Delta t = \frac{-v_y \pm \sqrt{v_y^2 - 4\,g_y(r_y-r_{ground})}}{2\,g_y} \end{align*}
Event-detection is driven by root finding.
Consider two dueling catapults:
We calculate $\Delta t$ for all possible events: \begin{align*} \vec{\Delta t}=\left\{\Delta t_{1-ground},\,\Delta t_{2-ground},\,\Delta t_{1-2}\right\} \end{align*} The earliest event is the next event!
Event processing is driven by sorting.
What actually happens at an event?
Consider a two-body event between particle $i$ and particle $j$, which has an energy change of $\Delta E$.
Conservation of momentum requires: \begin{align*} m_i\,\Delta \vec{v}_i = - m_i\,\Delta \vec{v}_i = \Delta \vec{P} \end{align*}
Where the impulse, $\Delta \vec{P}$, is a solution to
energy balance:
\begin{align*}\scriptsize
\frac{1}{2}m_i\, {v}_i^2 + \frac{1}{2}m_j\, v_j^2 + \Delta E = \frac{1}{2}m_i\, \left(\vec{v}_i +\frac{\Delta \vec{P}}{m_i}\right)^2 + \frac{1}{2}m_j\, \left(\vec{v}_j-\frac{\Delta \vec{P}}{m_i}\right)^2
\end{align*}
Event execution is driven by root finding for an impulse.
The time of the simulation is calculated from the number of events processed.
Exact solution of the model allows some interesting optimisations.
Initially, the simulation is run in synchronous mode. The time-warp algorithm is then enabled (with the same simulation speed) and only the required updates to the particles are shown. As events only involve small, localized particles, the system can run out of sync. i.e., each particle is at different points in time:
Can we simulate "normal" molecular models using EDMD?
The simulation can reproduce the underlying model up to a specified energy approximation, i.e., $\left|E^{sim}-E^{model}\right|<\Delta E$.
These stepped approximations can closely reproduce the phase behavior of the LJ system:
What about real granular systems?
M. Bannerman, J. E. Kollmer, A. Sack, M. Heckel, P. Müller, T. Pöschel, "Movers and shakers: Granular damping in microgravity," Phys. Rev. E, 84, 011301 (2011)
This led to the development of optimal design equations for dampers (see paper).
\begin{align*} L_{opt}=\pi\, A\, \sqrt{\frac{M}{M+N\,m}} + \sigma_{layer} \end{align*}