## 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 Newtons 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:

• The soft matter model which corresponds to a time-driven simulation approach. Time is devided into timesteps $\Delta t$ of equal length. For each time step Newtons equations are solved. If a particle does not overlap with any other particle nor with any boundary it freely propagates. If there is some overlap, a short-range potential is added which then yields a repulsive force. Also friction forces can be added during the overlap. Multi-particle collisions are possible.
• The hard matter model which allows a simplification of the above MD. In this approach it is assumed that collisions between particles occur almost immediately. In particular the duration of a collision is much smaller than the time between two collisions. So one can restrict the particle interactions to interactions of pairs of particles which allows noticeable simplifications. This model is suitable for the simulation of dilute gases with few collisions.
• Osther advanced approaches like MD Monte-Carlo, rigid body models, etc. (see T. Pöschel, T. Schwager, 2005) which are not considered here.

#### 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.

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:

• A Wall is a rectangular plane with finite sidelengths that is orthogonal to one of the coordinate axes. By setting a zero position on the plane, its position in space is then completely determined by the value on the orthogonal axis. It will be allowed to move in a sinusoidal way, therefore we assign a frequency and a amplitude.
typedef struct Wall_ {
// 0,1,2 related to the axis the normal vector is parallel to.
int orientation;
// dimensions in each space direction (from, to)
// pos[orientation] = {value on orientation axis, ...}
double pos[3][2];
// frequency
double omega;
// amplitude
double A;
// thickness
double d;
// coefficient of restitution
double Res;
} Wall;

• A Particle is completely determined by position, velocity, acceleration and its mass and radius.
typedef struct Particle_ {
// position of the particle
double r[3];
// velocity
double v[3];
// acceleration
double a[3];
// mass
double mass;
} Particle;


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:

• propagateWalls: For all walls with $\texttt{omega} > 0$ the current value on the span orientation-axis is set to $\text{starting pos.} + \texttt{A} \cdot \sin(\texttt{omega} \cdot t)$.
• calcAcceleration: For each particle the acceleration $\vec a_i$ is calculated as the sum of gravitation and the above LSDM potentials in the case of collision.
Therefore for each particle the distance \texttt{d} of the middelpoint to the walls in the system and to any other particle is computed (except for the distances that have already been processed because of symmetry). Now the acceleration follows directly from the LSD model. The spring stiffness $k$ is a constant with different values for particle - particle collisions and for particle - wall collisions. (s. a.) The coefficient of restitution that determines the dissipative force is experimental input and also a constant for particle-particle and particle-wall collisions respectively.
• propagateParticles: Here the position-Verlet algorithm is implemented for the accelerations calculated in the last paragraph.
• Output: The current particle positions, velocities and accelerations are written to a file.
• updateObservables: Some physically interesting numbers are derived from the current state of the system, especially the asymmetry parameter $\mathfrak{a}$ and the total energy of all particles. \newline These values are also written to a file.
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. \newline 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. \newline During the development process we faced some problems which we will discuss here shortly:
• Starting positions: With a naive random starting positions approach some particle spawn in each other, which leads to very high starting accelerations and starting velocities at the beginning. To prevent that, we iterate the choosing of random starting position until one particle starts without contact to any other.
• Energy too high: For certain configurations (for instance with high restitution coefficients) some particles developed velocities and thereby energies that were so high that these particles were able to pass the LSDM potentials of the walls. Firstly we introduced a wall-thickness to increase the maximal energy of these potentials. Indeed later it turned out that typical partical energies of our problem were much below this critical value.
• Problems with real-world parameters: Although the setting of the simulation is one to one designed after the real experiment, we had problems to observe the asymmetric effect with the experimental set of parameters like $\epsilon_\text{pp} \approx 0.95$ for particle - particle collisions. This may be due to the absence of some pertubation effects (like the margin of the bottom plate etc.), but until now we do not have a substantive explanation for that. We now use configuration parameters for the simulation that produce the wanted qualitative behaviour.
• Programming errors: Of course, Segfaults, logical errors, ...