Header-Logo: genug-davon

Simulation

Motivation

Since the experimental approach to granular matter effects is labour-intensive and in our case restricted to a low number of particles ($\leq 80$) involved, we use computer simulations as a second possibility to discuss the phenomenon described in the introdutory section.

Possible Simulation Models

For a given system of $N$ particles with positions $\vec r_i$ and velocities $\vec v_i$ the fundamental problem is to solve Newton`s equations of motion simultaneously for all particles $i$ (see T. Pöschel and T. Schwager. Computational Granular Dynamics. Springer, 2005.): $$\ddot{ \vec r}_i = \frac{1}{m_i} \vec F_i(\vec r_i, \vec v_i, \dots, \vec r_N, \vec v_N)$$ For a high number of involved particles and complicated interaction forces an analytic solution to this system of ODE is not achievable. Hence numerical methods need to be applied. A computer simulation solving these equations is called molecular dynamics (MD). In the literature many possible realisations of granular matter MD are discussed.

One can distinguish between the following models:

Our Model

In our simulation we use the above described soft matter approach because it is easy to implement and as the number of particles in our system is small, efficient computation is not of highest importance for us. Furthermore we neglect any tangential forces that could occur from particle interactions. These would only lead to torque and rotation of the particles which we assume has no impact on the wanted effect.

The important question now is how to model the forces between particles and between wall and particles. The formulas previously described correspond to the so called Hertz model (ebd.).

Since this model is not necessary to see the asymmetry effect which our project aims at, we use a simplier model following S. Lüding, the linear spring dashpot model (LSDM): We assume linear springs with spring konstant $k$ (harmonic potentials) between two particles that induce a repulsive force $$F_\text{el} = - k \xi$$ where $\xi$ is the overlap of two particles. A similar assumption holds for the dissipative force: $$F_\text{diss} = - \nu \dot \xi$$ with the viscous drag coefficient $\nu$.

One can then derive the equation of motion for the overlap: (Luding, S. 10f) $$m_\text{eff} \ddot \xi = - \nu \dot \xi - k \xi$$ Here $m_\text{eff} = m_1 m_2 / (m_1 + m_2)$ is the reduced mass of two colliding particles with masses $m_1$ and $m_2$. Let us now introduce the variables $\eta = \nu / (2 m_\text{eff})$ and $\omega_0 = \sqrt{k / m_\text{eff}}$.
Then we receive the well-known equation of motion of the damped harmonic oscillator $$\ddot \xi + 2 \eta \dot \xi + \omega_0^2 \xi = 0$$ with the general solution $$\xi(t) = \frac{\dot \xi(0)}{\omega} \exp(- \eta t) \left( -\eta \sin(\omega t) + \omega \cos(\omega t) \right) \\ \text{where} \qquad \omega = \sqrt{\omega_0^2 - \eta^2}$$.

We can then derive the coefficient of restitution in this model (Luding, S. 11), $$\epsilon = \exp \left( \frac{- \pi \eta}{\omega} \right)$$ and from that the drag coefficient $$\nu = \frac{\sqrt{4 \cdot m_\text{eff} \cdot k}}{1 + \left( \frac{\pi}{\log(\epsilon)} \right)^2}$$ So by the measuring of the coefficient of restitution this coefficient is determined. The coefficient $k$ has no influence on the energy dissipation, just on the duration of a collision. Since this time span is small compared to the time between two collisions, we choose any value $k$ that produces reasonable behaviour.

Comments on the Implementation

We implemented the above model in a C program. Therefore we followed the described time-driven approach.

We use mainly two types of objects (C structs), Particle and Wall:

We define the total simulation time $t_\text{max}$ (typically 30min), that correponds to the real time and the time step $\Delta t$ (around $10^{-3}$s). These two values determine the number of iterations $N = t_\text{max}/\Delta t$.

As integration algorithm for the Newton equations we use the position-Verlet algorithm. Let $\vec r_i^n$ denote the position of particle $i$ at the $n$-th time step. $\vec v_i^n$ and $\vec a_i^n$ are the corresponding velocity and acceleration. Position-Verlet then states (see J. Clewett. Emergent Surface Tension in Boiling Granular ./media. Dis- sertation, University of Nottingham, 2013.): $$\vec r_i^{n+1} = 2 \vec r_i^n - \vec r_i^{n-1} + \vec a_i^n \Delta t^2$$ $$\vec v_i^{n+1} = \frac{\vec r_i^{n+1} - \vec r_i^n}{\Delta t}$$

Each iteration step consists of the following operations: 

int i_max = t_max/DT;
for (i = 0; i < i_max; ++i) {
    propagateWalls(i*DT);
    
    // ps referes to the system of all particles
    calcAcceleration(&ps);
    propagateParticles(&ps, &ps_last);
    {
        ... output ...
    }
    // observables is the struct that saves all observables
    updateObservables(&ps, &obs);   
    memcpy(ps_last.p, tmp.p, tmp.n_parts * sizeof(Particle));
}

The Development Process

Before we started the work on the final three dimensional program, we considered some easier test cases to check the correctness of our algorithm, especially the energy disspation from particle collisions. Therefore we started with a one-dimensional test system including a collision of two particles. After this worked out fine, we included all three space dimensions and the walls at the boundaries of the box. During the development process we faced some problems which we will discuss here shortly: