For the integration of the equations of motion we use a fifth order predictor-corrector MD-scheme, see Ref. [18, 19]. Since we are interested in a static situation in 2D, with almost monodisperse particles, no particle has more than six nearest neighbors. For the simulation we keep the neighbors in memory in order to reduce the computational effort.
There are two forces acting on
particle i when it overlaps with particle j, i.e. when the distance
.
We use an elastic force
with the spring constant k, acting on particle i in normal direction
.
The second force in normal direction is dissipative
accounting for the inelasticity of the contacts.
In Eq. 2 the constant
is a phenomenological dissipation coefficient, and
is the relative velocity of
the particles i and j.
As mentioned above, we neglect tangential forces.
The contact of a particle with a wall or an immobile particle
is mimicked by setting the mass of the immobile contact partner
to infinity. Finally, the influence of gravity
g is readily included into the equations of motion.