Marcus N. Bannerman
School of Engineering, University of Aberdeen
m.campbellbannerman@abdn.ac.uk
I distribute an Event-Driven Particle Dynamics (EDPD) package, called
$\mathrm{EDPD} \supset \mathrm{hard\ spheres}$!
It is a general technique for simulating particles which can be compared to the more common Time-Driven (TD) techniques:
Why consider another DEM approach?
Event-Driven (ED) versus Time-Driven (TD) dynamics
Event-Driven (ED) versus Time-Driven (TS) dynamics
ED is not just for spheres
(Parallel hard cubes have a visible but continuous "freezing" transition which remains diffusive^{1}.)
Event-Driven (ED) versus Time-Driven (TD) dynamics
(The thin hard rod model has no excluded volume (ideal gas) yet it has complex transport coefficients^{1}.)
Event-Driven (ED) versus Time-Stepping (TS) dynamics
ED is not just for repulsive/instantaneous-contact forces.
(Square-well homopolymers are attractive, have multibody contacts, and are used in fundamental bio-polymer models.)
Event-Driven (ED) versus Time-Stepping (TS) dynamics
ED is not just for impulsive forces.
The granular damper container itself has a finite mass and is forced by a spring.
Rather than step through time, we calculate the time of "events" exactly (using an approximate model).
This is the automatic adaptive stepping of ED. It jumps to times where significant force/impulse is exchanged.
Assuming no drag (our approximate model), the exact solution to the projectile motion is a parabola:
When it hits the ground at a time we have:
The time of this "event" is a root of the overlap equation, : (note $f<0$ is an overlap, $f=0$ is contact, and $f>0$ is not overlapping)
In this case, this is a simple quadratic with two roots:
Event-detection is driven by root finding.
Consider two duelling catapults:
We calculate for all possible events: 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 and particle , which has an energy change of .
Conservation of momentum requires:
The impulse, , is a solution to
the conservation of energy balance:
Event execution is driven by root finding for an impulse.
Note, we do not decide the time of the simulation! We decide how many events to run it for and the time is calculated.
Modern EDPD algorithms are quite complex.
Initially, the simulation is run in synchronous mode, then time-warp is enabled (with the same simulation speed). Only the required updates to the particles are shown.
Much of what is presented so far is well established,
what is unique about
DynamO?
Stability and
generality.
Even though EDPD is analytic, computers are not exact. Round-off error occurs at the limits of the machine precision.
Although this appears insignificant, catasrophic cancellation can occur and the simulation can fail completely^{1}.
Consider repeated bounces of the projectile: Using a constant coefficient of restitution we have:
The velocity on impact/event decays to zero, and an infinite number of events happens in a finite time.
This "inelastic collapse" is not an instability, but a phenomena of the inelastic hard-sphere model.
Trying to simulate an "inelastic collapse" highlights a catastrophic numerical issue, the particle falls through the floor.
Although should equal after an event, limitations of the precision of the computer leave it at .
As , numerical precision is lost and the argument of the square-root can turn negative!
This gives no real-roots (or events) and the projectile falls through the ground!
This is not just an edge case, this problem is systemic.
It even plauges elastic/molecular systems for certain
sequences of rare ($P=10^{-10}$) events:
Machine precision is not enough to stop particles entering forbidden/invalid states, which may get rapidly worse.
Previous approaches typically attempt to prevent overlaps from forming using small adjustments, but round-off error is inherent and cannot be easily overcome.
Events, by their nature, move a system approaching an invalid state towards a valid state.
An algorithm that always moves towards a valid state is: Events occur at time, $t$, if particles are in contact, $f(t)=0$, or in an invalid state, $f(t)<0$, and this state is deterioriating $\dot{f}(t)<0$.
This appears to be stable in all cases, and correctly simulates inelastic collapse^{1}.
Applying this stable algorithm requires finding the roots of the overlap function $f(t)$, and its time derivative $\dot{f}(t)$
These functions become increasingly complex. For example, consider spheres in gravity...
In this system the overlap function is a quartic:
In rotating systems, the functions are transcendental, and yet we must reliably determine all roots (to remain stable).
DynamO has a compile-time Computer Algebra System (CAS) library. It allows simple definitions of overlap functions:
Variable<'t'> t;
Vector<> rij, vij, aij;
double diam2_ij;
/*...*/
delta_t = nextroot(pow<2>(rij+t*(vij+aij*t/2))-diam2_ij, t);
The compiler will transform the overlap function to create derivatives, and select different classes of root finding.
The library can find all roots of -order polynomials and certain classes of other functions in a stable way, allowing a wide range of models/events to be detected.
For up to 3rd order polynomials, solution by radicals is used. Solution by radicals is too numerically unstable for 4th order and impossible for 5th order.
For higher order polynomials, roots are bounded using Local Max Quadratic estimates^{1}, then bisected using Sturm's method to avoid repeated transformations of the polynomial.
For transcendental overlap functions, the roots are bounded using a bounding object with a polynomial overlap function. The bounds are then updated using a "worst case" polynomial approximation constructed using Taylor's theorem. This is accelerated using standard numerical root finding techniques.
Can we simulate "real" contact forces using EDPD?
We recently demonstrated elastic forces can be "stepped".
Conversion of a potential first requires the steps to be placed.
Optimal placement was determined to be at fixed changes in the potential energy, $\Delta U$ (analogue of timestep, $\Delta t$!).
Optimal energies are given by equating the second virial contribution, which is a volumetric average at $T\to\infty$:
These stepped approximations can closely reproduce the phase behaviour of the LJ system:
And the dynamics, for example, the shear viscosity:
Molecular forces are too soft and EDPD is inefficient at high density: harder, granular potentials are more efficient.
In the past few months, we have worked on a spring-dashpot force as an example implementation of viscous forces.
Previous work^{1} in this direction focussed on single-event collision dynamics and required storing "history". It also demonstrated the importance of considering individual particle trajectories.
A pair of colliding particles under with an elastic spring is simulated to test the trajectory.
For a single collision, the key parameter is where it lands on the final step (inset). $\tau = KE_{final}/(U_{final+1}-U_{final})$
There is an issue in the last step: for vanishing remaining energy the particles become stuck! ("elastic collapse"?).
By
enforcing a minimum kinetic energy (immediately bouncing
$\tau<\alpha$, where $\alpha=1/9$ in the last step), the
average collision time is correct.
Published dissipative forces in EDPD have so far been limited to inelastic coefficients.
To implement arbitrary dissipative forces cheaply, it should be approximated using an impulse/energy loss at each step:
Where $\Delta r$ is the width of the step the particle is leaving. No step-change = no energy loss (quasi-static elastic limit/no inelastic collapse).
Even for high inelasticiy ($\gamma$ set such that $e_{exact}=0.4$) this simple approach
is accurate with minor time scattering.
Relatively small numbers of steps are acceptable, with
collision time showing the most scatter.
Maximum distance is rigidly enforced.
And the model demonstrates the required constant
coefficient of restitution for spring-dashpot particles.
"Sleeping particles" algorithm allows free simulation of stationary particles. Possible massive speed-ups in heaps, rolling drums, etc.
Thanks go to Liang Xiao, Christopher Thomson, Severin Strobl, Leo Lue, and Thorsten Pöschel.
Thanks to you for your attention.