Motivation
Since the experimental approach to granular matter effects is labourintensive
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:

The soft matter model which corresponds to a timedriven simulation approach. Time is devided into timesteps $\Delta t$ of equal length.
For each time step Newton`s 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 shortrange potential is added which then yields a repulsive force. Also friction forces can be added
during the overlap. Multiparticle 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 MonteCarlo,
rigid body models, etc. (see T. Pöschel, T. Schwager, 2005) which are not considered here.
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 wellknown 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 timedriven 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;
// radius
double radius;
} 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 positionVerlet 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.
PositionVerlet 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^{n1} + \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 orientationaxis 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 particleparticle and particlewall collisions respectively.
 propagateParticles: Here the positionVerlet 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 onedimensional 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 wallthickness 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 realworld 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, ...