Fast, flexible particle simulations An introduction to MercuryDPM

We introduce the open-source package MercuryDPM , which we have been developing over the last few years. MercuryDPM is a code for discrete particle simulations. It simulates the motion of particles by applying forces and torques that stem either from external body forces, (gravity, magnetic fields, etc.) or particle interactions. The code has been developed extensively for granular applications, and in this case these are typically (elastic, plastic, viscous, frictional) contact forces or (adhesive) short-range forces. However, it could be adapted to include long-range (molecular, self-gravity) interactions as well. MercuryDPM is an object-oriented algorithm with an easy-to-use user interface and a flexible core, allowing developers to quickly add new features. It is parallelised using MPI and released under the BSD 3-clause licence. Its open-source developers’ community has developed many features, including moving and curved walls; state-of-the-art granular contact models; specialised classes for common geometries; non-spherical particles; general interfaces; restarting; visualisation; a large self-test suite; extensive documentation; and numerous tutorials and demos. In addition, MercuryDPM has three major components that were originally invented and developed by its team: an advanced contact detection method, which allows for the first time large simulations with wide size distributions; curved (non-triangulated) walls; and multicomponent, spatial and temporal coarse-graining, a novel way to extract continuum fields from discrete particle systems. We illustrate these tools and a selection of other MercuryDPM features via various applications, including size-driven segregation down inclined planes, rotating drums, and dosing silos.


Introduction
Granular materials -conglomerations of discrete, macroscopic particles -are ubiquitous in both industry and nature. They range from natural materials like snow, sand, soil, coffee, rice and coal to man-made agglomerates such as medicinal tablets, catalysts or animal feed. Understanding the behaviour of granular media is of paramount importance to the pharmaceutical, mining, food processing and manufacturing industries, and highly relevant to the prediction and prevention of landslides, earthquakes and other geophysical phenomena.
MercuryDPM [1][2][3][4] is an open-source package for simulating granular materials with the discrete particle method (DPM) [5]. It simulates the motion of N particles in a system constrained by N w walls and body forces. It assumes: (i) Particles are unbreakable; however, breakage can be included by forming clusters of 'primary' particles, see Section 6.3 for details. (ii) Particles are undeformable, such that the particle masses m i and inertia tensor I i are constant in the body-based frame. Note that clusters can be used to model deformable particles, see Section 6.3. (iii) All interactions between the particles are binary, i.e. all internal forces/torques are due to particle pair interactions. (iv) Each particle pair i, j has at most a single contact point c ij at which the interaction forces f ij and torques τ ij act.
(v) All external forces/torques acting on a particle i are either body forces f b i or interaction forces f w ik with a wall k. The same is true for torques.
The force and torque acting on each particle i can then be computed as with the branch vector r ij = c ij − r i connecting the particle position r i with the contact point c ij . For given initial conditions, Newton's second law can then be used to evolve the particles' velocities v i , positions r i , angular velocities ω i and orientations For computational stability, the orientation is stored as a quaternion q i ∈ R 4 , which requires the use of a transformation matrix C(q i ); see [6,7] for details. The above differential equations are solved numerically using the Velocity-Verlet algorithm, which is symplectic (thus, energy is conserved in case of elastic forces) and second-order accurate. Using a higher-order accurate time integration scheme would not increase accuracy, because most DPM contact models are non-differentiable.
MercuryDPM is written mainly in object-oriented C++, using many modern features from C++11. We try to follow the latest developments in the C++ language; however, we also guarantee the release version will work on two-year old compilers. The latest version uses parts of the Fortran LAPACK library [8]. However, the parts we use have been incorporated into the MercuryDPM source code, thereby avoiding an external dependency; only a Fortran compiler is required. The code has an easy-to-use user interface and a flexible core, allowing developers to quickly add new features. It is parallelised using MPI and released under the BSD 3-clause licence. Thus, it can be used as part of closed-source derivatives, as long as the derived software acknowledges the MercuryDPM team.
We have tried (and hopefully succeeded) in making Mercury-DPM easy to learn for new users. Its installation process is simple and the package includes many tutorials as well as example codes that demonstrate the package's features. It is also supplemented by a detailed reference manual (docs.mercurydpm.org). Users otherwise unfamiliar with C++ have found the project's coding style intuitive, allowing them to focus on modelling problems.
The code is being developed by a global network of researchers and in the last few years has received contributions from universities such as Cambridge, Stanford, EPFL, Birmingham, Strathclyde, Sydney, Northwestern, Rome, Delft and Manchester, as well as industry, such as MercuryLab in Enschede and RCPE and ECCPM in Graz. We encourage all MercuryDPM users to merge the features they develop into MercuryDPM, thus becoming Mercury-DPM developers.
As the code is fully open-source, all features we develop can be accessed and reused freely for both non-commercial and commercial use. The open-source philosophy allows the code base to grow quickly, and the open-source development reduces the amount of coding errors, as you get near-imminent feedback from other developers. Its open-source community has developed many features, including moving and curved walls; state-of-theart granular contact models; specialised classes for common geometries; non-spherical particles; general interfaces; restarting; visualisation; a large self-test suite; extensive documentation; and numerous tutorials and demos. In the following, we review some of these features.

Coding philosophy
MercuryDPM is written in an object-oriented programming style, i.e. it uses classes to define objects: spherical particles are objects of type SphericalParticle; planar walls are of type InfiniteWall; and periodic boundaries are of type Period-icBoundary. A myriad of other classes have been implemented and ready for use, many of which are described in following sections The user can also derive their own classes by inheriting from existing ones, adding extra functionality.
The clear and structured nature of MercuryDPM means it is quick and easy to develop new features; however, the level of C++ required is still demanding to some users. Therefore, we are developing a graphical interface, opening it up to a whole new set of users, both academic and industrial.

Major components
MercuryDPM has three major components that were originally developed by its team: (i) a contact detection algorithm [10] that can efficiently simulate highly polydisperse packings, which are common in industry; (ii) curved walls [3], making it possible to model real industrial geometries exactly, without triangulation errors; and (iii) MercuryCG [11][12][13][14], a state-of-the-art analysis tool that extracts local continuum fields, providing accurate analytical/rheological information often not available from experiments or pilot plants.

Contact detection
Contact detection -determining which particle pairs are in contact -is one of the most complex parts of any DPM algorithm and can consume the majority of the computational time, if it is not carefully implemented.
The most basic contact detection simply loops through all particle pairs; this algorithm is of quadratic complexity, O(N 2 ), where N is the number of particles in the simulation. Because the rest of the DPM algorithm is of linear complexity, O(N), such a contact detection would make large simulations prohibitively expensive. A more efficient contact detection algorithm is needed.
Most DPM algorithms use the linked-cell algorithm for contact detection [15], illustrated in Fig. 1 left: Particles are placed into a grid whose cell size is the diameter of the largest particle. Thus, particles can only be in contact with particles in the same or in a neighbouring cell, reducing the number of necessary checks. For (nearly) monodispersed simulations, this algorithm is of linear complexity, O(N), and thus sufficiently efficient. However, the order complexity increases to quadratic, O(N 2 ), for highly polydisperse simulations, because the cell size is based on the largest particle diameter.
MercuryDPM uses the hierarchical grid (hGrid) [9,10,16], an advanced contact detection method that uses several grids for different particle sizes, as shown in Fig. 1 middle. This contact method gives MercuryDPM its name: hGridDPM → HgDPM → MercuryDPM. By carefully selecting the number of levels and cell sizes, linear complexity of the algorithm can be guaranteed even for the most challenging particle size distributions.
The effectiveness of the approach has been proven theoretically and in simulations [9]: for highly polydisperse situations, the hierarchical grid is two orders of magnitude quicker than the linked-cell algorithm, see Fig. 1 right. Furthermore, the CPU time varied only minimally when comparing monodisperse, bidisperse and polydisperse systems with the same number of particles and volume fraction [9]. This feature allows for the first time large simulations with wide size-distributions.
The hierarchical grid algorithm is made up of two phases: mapping and contact detection. In the first phase, all particles are mapped onto a grid level. In the second phase, the potential contact partners for every particle in the system are determined. Both phases are of linear complexity, and allow straightforward parallelisation.
The two-or three-dimensional hierarchical grid is a set of L regular grids with different cell sizes s 1 < s 2 < · · · < s L .
Each particle p is mapped into a specific cell in the lowest level grid big enough to contain the particle. A new mapping is done before every time step, as this is cheaper than tracking changes and updating the map. It must be noted that the cell size s h of each level h can be set independently, in contrast to contact detection methods which use a tree structure for partitioning the domain [17], where the cell sizes are taken as double the size of the previous lower level of hierarchy, hence s h+1 = 2s h .
The flexibility of independently choosing s h allows one to select the optimal cell sizes according to the particle size distribution, further improving the performance of the simulations [9]. The contact detection is split into two steps, and the search is done by looping over all particles p and performing the first and second steps consecutively for each p. The first step is a contact search at the level of insertion, using the classical linkedcell method: the search is done in the cell onto which p is mapped, and in its neighbour (surrounding) cells. Only half of the surrounding cells are searched, to avoid testing the same particle pair twice. The second step is the cross-level search. For a given particle p, one searches for potential contacts only at levels of smaller grid size. This implies that the particle p will be checked only against the smaller ones, thus avoiding double checks for the same pair of particles. See [9] for details of the algorithm.

Application: Segregation in rotating drums
Segregation of grains by size is a scientifically interesting and industrially relevant problem. In industrial situations, sizedistributions often range over orders of magnitude and are highly polydispersed; whereas academic studies often consider bidispersity with a factor of only 2-10 in size. One key reason for this discrepancy is computational cost. However, the hierarchical grid, at the heart MercuryDPM, is over three orders of magnitude faster for bidispersed mixtures with a size ratio of 100, and even faster for truly polydisperse packings; see [16] for details. Fig. 2 shows a simulation where each particle is a unique size and the ratio of small to largest radii is 100. This is visualised using ParaView and was run on a single core on a normal desktop computer in a few hours.

Curved walls
A distinguishing feature of MercuryDPM is its support of curved geometric surfaces, or walls. Many types of walls are implemented and ready for use, such as • Flat walls (implemented in InfiniteWall), • Convex polygons (2D) or polyhedra (3D) (IntersectionOfWalls), • Conical and cylindrical shapes, created by rotating a polygon around an axis (AxisymmetricIntersectionOfWalls), Source: Data from [9].

Fig. 2.
Simulations visualised in ParaView of size-based segregation in drum. Colour indicates particle size, from red (small) to green (medium) to blue (large).
Each wall k has a position p k and orientation q k that can be either a fixed value or a prescribed function, to simulate moving and rotating walls.
For contact detection, each wall type has a getDistanceAnd-Normal(particle, distance, normal) function that computes the contact normal n ij and distance d from the wall for any given particle. These values are necessary to compute the contact force.
For many walls, the contact normal and distance is computed analytically (and thus efficiently): For example, for a flat wall k and a spherical particle i, the normal direction n ik can be computed from the wall orientation, and the distance to the particle is given by d = (p i −p k )·n ik . Analytic solutions are also used for In-tersectionOfWalls, AxisymmetricIntersectionOfWalls, TriangleWall, NURBSWall, and LevelSetWall. Note that In-tersectionOfWalls are not simple flat surfaces, but have defined face, edge and vertex contacts; Fig. 4 shows the normal and distance computed for the case of a triangular wall.
More complex wall shapes like the Coil or Screw require an iterative scheme. In both cases, Newton iteration is used to minimise the distance to the wall. Fig. 5 shows the normal and distance for a Screw.
In particular NURBS surfaces are very general and can be used to simulate many types of walls. However, if a particular surface type is not yet implemented, the user can define it by creating a new wall type and writing an appropriate getDistanceAnd-Normal function.
Most other codes approximate curved surfaces via triangulated walls. This can be done in MercuryDPM as well: triangulated walls can be stored as STL or VTK files, and read-in using the readTriangleWalls function. However, it is generally not recommended to use triangulated walls for the following reasons: Firstly, the discretisation error of triangulating surfaces can be significant, especially for surfaces with high local curvature, such as coils or helicoidal shapes, or for moving surfaces that are only separated by a narrow gap. Secondly, as you refine your triangulation, you very quickly get to a large number of triangles, which slows down the contact detection; whereas, with the MercuryDPM curved wall support you have just one wall.

Application: Industrial mixers
All of the above-mentioned triangulation problems occur when studying industrial mixers. One such example is the Nautastyle mixer shown in Fig. 6 right. In MercuryDPM, this mixer is composed of only four curved surfaces: two conical walls for the casing and the base, and a helical screw with a cylindrical shaft, which rotate around their axis as well as along the casing. The high curvature of the helical screw and the narrow gap between the screw and the outer casing are hard to resolve using triangulated surfaces. Thus, triangulated geometries need to be highly refined, with thousands of triangles representing a single surface, which is less efficient and less accurate than using curved walls.

Application: Tunnel boring machine
Thanks to MercuryDPM's support of curved walls, it was possible to simulate a Tunnel Boring Machine (TBM). TBMs are used to excavate tunnels with a circular cross section; they have a rotating wheel, called cutter-head, used to excavate the soil. When the ground is soft, Earth Pressure Balance Machines (EPB) are used. They get this name because they use the excavated material to balance the pressure at the tunnel face and this is obtained using    a screw. The screw allows the maintenance of the prescribed pressure inside the excavation chamber. EPBs have a complex geometry, but with MercuryDPM it is possible to obtain a simplified version using a novel hybrid of complex and triangulated walls. A simplified EPB was obtained using already implemented shapes, like AxisymmetricIntersectionOfWalls, Screw and Tri-angleWall. The EPB's body (Fig. 7a) was created using curved shapes, while the cutter-head, which has a more complex shape, was read in as a triangulated wall from an STL file, using read-TriangleWall. The cutter head (Fig. 7b) has been designed from a physical EPB model used in the Laboratory of Civil Engineering and Building Sciences of ENTPE in Lyon (France). With this model, it was possible to simulate the excavation phase and analyse the behaviour of the tunnelling ground on-site, varying parameters such as the EPB's velocity, the cutterhead's angular velocity and the screw's angular velocity [21].

Coarse graining
Coarse graining (CG) is a micro-macro transition method: it extracts continuum fields (density, momentum, stress, etc.) from discrete particle simulations, allowing the validation and calibration of macroscopic models [11][12][13][14]. Unlike binning methods, which extract average values in small volumes, coarse graining evaluates continuum fields as a function of time and space. Thus, unlike binning, the resulting fields are continuous, satisfy local mass and momentum conservation exactly, and the spatial and temporal averaging scales (w and w t ) are well-defined [12,13].
The approach is flexible and the latest version can model both bulk and mixtures [13,22], boundaries and interfaces [11], non spherical particles [23], time-dependent [13], steady and static situations. It is available in MercuryDPM either as a postprocessing tool or it can be run in real-time, i.e., concurrent with the simulation.
The aim of coarse-graining is to define continuum fields that automatically satisfy the continuum model. Most continuum models are based on mass and momentum conservation; in that case, we need to define density ρ, velocity u and stress σ such that these fields satisfy Spatial coarse-graining defines density by applying a smoothing kernel φ(r) to the statistical-mechanics definition of density, The integral over density should equal mass, density should be non-negative, and the density at a point r should only depend on the particles within the neighbourhood of r. Therefore, we require that the kernel function φ: • has compact support: ∃c ∈ R: φ(r) = 0 for all |r| > c.
A typical coarse-graining function is the cut-off Gaussian, with appropriate prefactor C , or cut-off polynomial functions such as the Lucy kernel, which is smoother (2nd-order differentiable) and more efficient to evaluate than φG; see [12] for details.
To satisfy (1), velocity is defined as Similarly stress has to be defined to satisfy (2).
Note, this integral can usually be computed exactly, without requiring numerical quadrature. For e.g. a Gaussian kernel we obtain other coarse-graining functions (cutoff Gaussian, polynomials), the definitions of ψ ij are more complex, but also explicit.
Output data can be coarse-grained in MercuryDPM using the MercuryCG tool. 1 For example, the following command will apply CG to the output of the FiveParticles application: The result is shown in Fig. 8 (centre). More examples can be found in [24].
MercuryDPM is the only code where coarse-graining can be applied during a simulation. This allows for efficient computation of high-resolution continuum fields without the need for large output files. It further allows for a two-way coupling between the continuum fields and the particle simulation, e.g. for solid-particle coupling (force-controlled walls) and particle-fluid coupling (suspensions); see Section 9 for details.
MercuryCG can also be applied to data from other DPM simulation softwares, and even experiments: it has e.g. been used to analyse experimental data from particle position tracking [25].
There is a lot more to say on the details of coarse-graining, and we will soon publish more details about the algorithm.
Note, the MercuryCG tool is currently being updated, and will soon get a new interface, and more capabilities, including: temporal smoothing kernels; coarse-grained liquid distribution and coarse-grained particle size distribution; Lagrangian-style coarsegraining where the evaluation points move with the flow [26]; and the ability to define your own coarse-grained fields. 1 The tool was called fstatistics in previous publications.

Boundary conditions
Basic particle simulations start with an assembly of particles and walls, and simulate the particle (and wall) motion over time. However, many process simulations require more complex setups, including the insertion or deletion of particles during the simulations; periodic boundary conditions; or stress-, displacement-or temperature-control. Such extensions to the basic DPM simulation setup have the suffix Boundary in MercuryDPM, and are stored in the BoundaryHandler. We now review the most commonly used boundaries.

Insertion boundaries
Insertion boundaries are used to insert new particles during the simulation. Most commonly used is the CubeInsertion-Boundary, which inserts particles in a rectangular region of the domain. The user can define the insertion region, a (variable or constant) insertion rate; a particle size distribution; and a velocity distribution. Other insertion boundaries have been implemented for different geometries of the insertion region (ChuteInsertionBoundary, HopperInsertionBoundary).

Deletion boundaries
Deletion boundaries are used to remove particles during the simulation. The CubeDeletionBoundary removes all particles that enter a rectangular section of the domain.

Periodic boundaries
Periodic boundary conditions are used to simulate small representative volume elements, allowing the study of certain flow or deformation conditions without the influence of walls. Uni-, bi-or triaxial compression, simple shear flow, chute flow, Hele-Shaw or Couette geometries are just some of the many examples. MercuryDPM has three kinds of periodic boundary conditions: • PeriodicBoundary simulates a periodic region between two parallel planes.
• AngledPeriodicBoundary simulates a wedge-shaped periodic region between two non-parallel planes The implementation uses ghost particles, i.e. copies of real particles that are close to the periodic boundary, in order to transfer forces across the periodic boundary. A sketch of the three boundary conditions is given in Fig. 9.

Stress-and strain-controlled periodic boundaries
The StressStrainControlBoundary simulates a wide range of stress-and strain-controlled shear flow and compression tests [27]. It combines a normal periodic boundary in z-direction with a Lees-Edwards boundary in x and y [28][29][30], as shown in Fig. 10. The user can specify a combination of targets for the stress and strain rate tensor. The advantage of this boundary is the freedom of choosing the control parameters freely, allowing the user to specify both a target stress tensor, σ, and a strain-rate tensor,ε, as input parameters. For example, • for constant-stress uniaxial compression, specify σ = (σ xx 0 0, 0 0 0, 0 0 0) and setε to zero; • for constant-rate uniaxial compression, set σ to zero and specifyε = (ε xx 0 0, 0 0 0, 0 0 0); • for triaxial compression, which is mostly used in sample preparation to achieve a homogeneous initial packing, set σ = (σ xx 0 0, 0 σ yy 0, 0 0 σ zz ) andε = 0 for stresscontrol, orε = (ε xx 0 0, 0ε yy 0, 0 0ε zz ) andσ = 0 for volume-control.
Note that the same element in the target stress and strain-rate tensors cannot be set simultaneously, e.g. the user could not set σ = (σ xx 0 0, 0 0 0, 0 0 0) withε = (ε xx 0 0, 0 0 0, 0 0 0) at the same time, or the two control targets will conflict with each other, resulting in an invalid deformation mode. Further, because the Lees-Edwards boundary conditions are applied in the xy-plane, shear can only be applied in the xy-direction (not in xz and yz). However, one can achieve all possible physical constant-rate deformation modes with 3 diagonal elements and one off-diagonal element.

Other boundary conditions
Many other interesting boundary conditions have been implemented in MercuryDPM:  • SubcriticalMaserBoundary and ConstantMassFlow-MaserBoundary insert particles from a periodic system into a larger setup, such as chute flows [31], allowing the simulation of steady-state inflow conditions.
• HeaterBoundary acts similar to a thermostat, supplying a heat flux to a specific region; • FluxBoundary does not affect the simulation but measures flow rates through a given plane.
We refer the reader to the documentation (docs.mercurydpm .org) for further reading.

Contact models
Contact models are used to determine the forces acting between particle pairs. Many different contact forces have been described in literature, which can roughly be classified into three categories: elastic, plastic and dissipative forces f n ij that act in the normal direction to the contact area, n ij ; tangential forces f t ij and torques τ ij due to sliding, rolling and torsion friction; and adhesive normal forces f a ij that may act between nearby particles even if they are not in contact. Which contact model best describes the real contact behaviour depends on the material type and particle size, and on ambient effects such as temperature and moisture. In most cases, a combination of these forces needs to be taken into account, i.e. the total contact force is given as All contact models in MercuryDPM are defined by combining a normal, frictional, and adhesive contact model. The normal, frictional and adhesive contact models currently available in MercuryDPM are summarised in Table 1. The name of a contact model is obtained by concatenating the names of the normal, frictional, and adhesive contact model and adding the word Species. For example, particles of type LinearViscoelas-ticFrictionLiquidBridgeWilletSpecies interact with a linear spring-dashpot normal force, sliding, torsion and rolling friction, and liquid-bridge adhesion forces. All contact models require a normal force, but frictional and adhesive forces are optional. Thus, LinearViscoelasticSpecies denotes the simple linear spring-dashpot contact model.

Normal force models
For a pair of spherical particles i, j with radii a i , a j and position r i , r j , we first compute the distance vector r ij = r j − r i , then the overlap δ n ij = a 1 + a 2 − |r ij |, the unit normal n ij = r ij /|r ij |, the branch vector r ij = (a i − δ i /2)n ij and the contact point c ij = r i + r ij . The same quantities need to be defined for pairs of nonspherical particles and particle-wall contacts, as these quantities are necessary to define the contact force.
We further define the relative velocity at the contact point, The normal contact force f n ij is nonzero if the two particles are in contact, i.e. if their overlap is positive. Several different normal force models are implemented:

Linear spring-dashpot model
The linear spring-dashpot model, implemented as Linear-ViscoelasticSpecies, is defined as with stiffness (or spring constant) k n > 0 and damping coefficient γ n ≥ 0. The force-displacement relation is shown in Fig. 11 left.
The model is efficient and simple to analyse, with collision time t c and restitution coefficient ϵ n only dependent on particle mass, not relative velocity. It is an appropriate contact models for large particles (>100 µm) or upscaled systems, where one particle represents a conglomerate of particles. For large deformation, plasticity needs to be taken into account, see Section 5.2. For sound and compaction experiments, calibrating the stiffness is important. For flows, stiffness just has to be sufficiently high such that the deformations remain small.

Hertz spring-dashpot model
The Hertz elastic normal force (HertzianViscoelastic-Species) is based on the contact between perfectly elastic, spherical particles [39]. It is based on measurable material parameters (the elastic modulus and Poisson ratios) and accurate for powders and small deformations of larger granules, as they appear e.g. in sound propagation experiments. In most other experiments, however, the choice of stiffness is either of little importance, as other effects (such as dissipation, friction or plasticity) become dominant. The Hertz force model is given by (3), but the stiffness k n is now a variable of the radius of the contact area a c ij , with the effective modulus E eff The dissipation coefficient γ n is chosen to satisfy the modeller's assumption on the restitution coefficient. See [40][41][42] for a review of the different dissipation models. In MercuryDPM, the dissipation coefficient γ n is proportional to √ m eff ij k n , resulting in a constant restitution coefficient [43].
The proportionality of stiffness and contact radius has an important effect on the particle behaviour: It makes particles stiffer under pressure (see Fig. 11 middle), thus the speed of a pressure wave increases if the material is compressed. This is important in sound propagation but can usually be neglected in flow situations, where stiffness has only a minor effect. Furthermore, while Hertz models work well for spherical powders, the stiffness measured in macroscopic particles (> 100 µm) is often relatively constant, due to plasticity, surface roughness and particle shape effects.

Linear elasto-plastic cohesive
To mimic plastic deformation, as observed in experiments [44], Walton and Braun [45] and Walton [46] introduced a so-called 'partially latching spring' model that used different normal spring stiffnesses for loading and unloading, To track the plastic deformation, the maximum overlap δ max ij is stored and used to compute the plastic overlap δ 0 ij and minimumforce overlap δ min ij , This model was extended by Luding [32] to allow for slowly changing stiffness: During loading, the stiffness increases linearly with the maximum overlap until a maximum unloading stiffnesŝ k 2 is reached at a fraction φ of the effective radius, ) .
The force-displacement curve for this model (LinearPlasticViscoelasticSpecies), is shown in Fig. 11 right.

Solid-state sintering
Solid-state sintering is a thermal treatment for bonding particles into a solid structure. Particles are sintered by heating particles beyond the glass temperature, but below the melting point of a material. This process is controlled by transport and diffusion of material along the particle's surface and volume, which leads to a reduction of the particle surface area. Solid-state sintering has three stages [47]. The first stage is neck formation: Matter from the particle is transported from regions of high chemical potential (contact region) to regions of low chemical potential (concave neck regions). In the second stage the diameters of the pores channels shrink until the pore structure changes. In the last stage isolated pores form. All stages are dominated by different transport mechanisms, and there is a strong dependence on temperature and initial particle size in the final stage of the process. The solid-state sinter model in MercuryDPM (SinterSpecies) describes the first stage, introducing a gradual increase in plastic overlap between particles at high temperatures. In [48][49][50], the model is applied to sintered polystyrene particles, and numerical and experimental indentation tests are executed and compared, see Fig. 12.

Partial melting
If particles are heated beyond their melting temperature, they start to melt. The particles first melt at the surface, where the heat is applied, forming a melt layer that increases in thickness until the particles are fully melted. This process can be modelled in MercuryDPM using the MeltableSpecies and was initially developed for additive manufacturing processes. In particular, we consider powder bed fusion (PBF), where objects are produced by spreading successive layers of powdered material and hardening selected parts by partially or fully melting them with a laser. PBF processes are highly sensitive to the powder characteristics; therefore, the process parameters need to be optimised for each material. This is typically done by performing costly experimental trials, so developing a computational tool capable of capturing the stochastic nature of the process will help in reducing the amount of trials and thus lower manufacturing costs. In addition, particle-scale simulations of the spreading process can provide information on the powder layer behaviour and quality that is not accessible by experiments (porosity, particles segregation, etc.). A parametric study of the influence of inter-particle friction on the powder layer quality has been done by [51]. MeltableSpecies is based on the model of [52], which was applied to the formation and cyclic melting of faults during earthquakes. It is assumed that solid particles can melt during heating. On cooling these melt layers solidify, potentially forming permanent bonds between the particles. The model was modified and extended to include thermal conduction, radiation and convection. Further extensions will include phase transformation and laser heat input modelling. Fig. 13 shows a snapshot of particles partial melting, with particles diameter range 40-50 µm, and illustrates the heating and solidification of a new layer of particles.

Sliding friction
Similarly to the normal forces, one can define forces in tangential direction. For this, we define the lateral relative velocity, v l ij = v ij − v n n ij , and the tangential elastic displacement δ l ij , which is set to 0 at the initiation of contact, and incremented after every timestep by the formulã The two-step procedure is necessary to keep the tangential displacement in the tangential plane while the particle pair rotates. Note that for a fixed normal direction n ij , we obtain d ij . An elastic and dissipative lateral force can then be defined as If the lateral force exceeds a certain level, the particle begins to slide. This is modelled by a Coulomb yield criterion, cutting off |f l ij | ≤ µ l f n ij .
The model, shown schematically in 14, agrees well with experimental data.
The above sliding friction model is implemented in the Slid-ingFrictionSpecies, and is intended to be used with a linear normal contact force. For the Hertzian normal force, the MindlinSpecies is more appropriate, which has a variable tangential stiffness k t that depends on the effective shear modulus; see [33] for details.

Rolling and torsion torque
Similar to sliding friction resisting lateral motion, rolling and torsion torques are modelled to resist angular motion. Like sliding friction, these torques are modelled as elastic and dissipative with a yield criterion. For this, we define the rolling and torsion velocity, v ro ij = a eff ij (ω ij × n ij ), v to ij = a eff ij (ω ij · n ij )n ij . Their respective displacements δ ro ij , δ to ij are defined equivalently to (6). We then define elastic-dissipative rolling and torsion torques with a eff ij = |r ij | |r ji | |r ij |+|r ji | the effective length of the branch vectors. If these torques exceed a certain fraction of the normal contact force, the particle begins to roll or torque, respectively. Thus, the elastic displacement is cut off to satisfy Similar to the sliding friction model, this model (Friction Species) is intended to be used with a linear normal contact force. For the Hertzian normal force, the MindlinRollingTor-sionSpecies is more appropriate, see [33] for details.

Linear reversible adhesive force
The simplest adhesion model in MercuryDPM, Reversible AdhesiveSpecies, models a linear elastic-dissipative shortrange force, It is called reversible, because the adhesive force is equal during loading and unloading. This model is mostly used to model dry cohesion, i.e. attractive forces due to close proximity between surfaces, such as van-der-Walls interactions.

Linear irreversible adhesive force
IrreversibleAdhesiveSpecies is an irreversible linear elastic-dissipative short-range force, where the short-range adhesive force is active only during unloading (i.e. after the particles have been in contact) Such models are often used to model wet cohesion, i.e. liquid bridges between particle pairs, which form when the particle pair gets in contact, but persist during unloading until the liquid bridge snaps. Both contact models are sketched in Fig. 15.

Liquid-bridge cohesion
LiquidBridgeWilletSpecies is a nonlinear model to model liquid bridges [36,54], based on the theoretical (and experimentally validated) results of Willet et al. [55]. This forcedisplacement relation has been derived from first principles, based on solving the Young-Laplace equation, resulting in the following force model, shown in Fig. 15.
if −S c < δ ij < 0 and was in contact, 0 else, ij cos θ, liquid volume V b , contact angle θ, and surface tension γ . This model has further been extended to account for liquid migration, implemented in LiquidMigrationWilletSpecies.
The methodology is quite straightforward: Particles and liquid are considered as two different entities in the system. Liquid is either attached to the particles (as a thin liquid film), or to the contacts (as liquid bridges). Liquid is transferred whenever contacts are formed or broken. Thus, when a contact is formed between two particles, the liquid attached to the particles can form a liquid bridge [37]. When a liquid bridge is ruptured, the bridge volume is distributed to neighbouring particles and contacts; total liquid conservation is ensured. The microscopic simulations of liquid migration has been used to validate a continuum scale model that describes the migration of liquid as shear-rate dependent diffusion [56].

Bonds and charges
Several other adhesive force models exist: BondedAdhesive Species allows the user to specify a strong adhesive force f b ij between particles, which bonds the particles together. A bond force can be turned on or off for each individual particle pair, allowing for the simulation of soft, bendable particle clusters: In ChargedBondedAdhesiveSpecies, the bond force is combined with a short-range normal force f c ij simulating particle charges, which can be either repulsive or adhesive, depending whether both particles have the same or opposite charge: It has been used to simulate clay particles as elongated, stringshaped particles (in 2D) or oblate particles (in 3D) with opposite charges at centre and edge of the particles [38].

User-defined species
Of course, the user can also define new contact models. For this, it is easiest to modify the most similar existing model. This is described in detail in the For Developers section of the documentation (http://docs.mercurydpm.org).

Non-spherical particles
As DPM studies become more complex and detailed, many users wish to use non-spherical particles in their simulations. MercuryDPM supports several ways to define non-spherical particle shapes, such as multi-spheres [57], superquadrics [58], agglomeration [35], and bonding [38]. We now show different non-spherical particles, using example applications.
The parameters a, b, c determine the particle size in the x, y, z direction, and n 1 , n 2 determine the roundness of the edges. For example, we get ellipsoids for n 1 = n 2 = 2, a cylinder with rounded edges for n 1 ≫ 2, n 2 = 2, and a cuboid with roundededges when n 1 , n 2 ≫ 2. Contact detection and computation of the overlap is implemented similar to [58]: each particle fits into a bounding sphere of radius r = max (a, b, c); based on that radius, the particles are inserted into the hierarchical grid. Whenever the hierarchical grid finds a potential contact, it is tested whether the bounding spheres of the particles i, j intersect. If that is the case, the contact-point c ij is the defined as the x-value minimising the is the shape-function of the particle i, translated and rotated by From there, the direction of the contact and its overlap can be computed. For implementation details, see [59]. Note that, since the coarse-graining tool MercuryCG does not rely on particle shape, the continuum fields for these quantities can automatically be computed without any changes to the code.
To study the influence of particle elongation on segregation, we construct a rotating cylindrical drum made out of small particles. The drum is filled with mixtures of spheres and prolate ellipsoids of equal volume and equal density. After ten rotations, the mixture is coarse-grained over one and a half a rotation period in order to obtain the concentration of spheres throughout the drum. We confirmed the observation of [60] that for a combination of spheres and prolate ellipsoids with aspect ratio 2, the ellipsoids segregate to the core, while for a combination of spheres and prolate ellipsoids with aspect ratio 4, the ellipsoids segregate to the periphery of the flow; more detailed observations will be presented in a follow-up publication. Fig. 16 shows the segregation profile for both these cases.

Multispheres
For some applications, the geometry of the particles is relevant (e.g., railway ballast), and assuming spheres or ellipsoids may be a very simplistic approximation. For this reason, the concept of ''multispheres'' is being implemented in MercuryDPM. A multisphere is a cluster of spheres (or superquadrics) whose relative positions are fixed and are placed such that the boundary of the cluster approximates the intended geometry. The mass and inertia of the particle is computed such that it matches the particle's geometry. In this way, a multisphere is similar to an agglomerate of elementary particles, but no internal deformation and/or breakage is allowed. The elementary particles (slaves) composing a multisphere can overlap and their radius may vary. Due to the possible overlap of slaves composing a multisphere particle, the inertia of the multisphere cannot be calculated internally by the software; instead, it must be specified by the user.
The steps to define a multisphere particle are: • Create a new particle, e.g. SphericalParticle p; • Define its initial position (centre of gravity): p.setPosition(Vec3D pos); • Define its principal axes: p.setPrincipalDirections(Matrix3D dir); Each column of the 3D matrix represents a principal axis; the software enforces orthogonality internally • Define mass and inertia: p.setMass(double mass); p0.setInertia(MatrixSymmetric3D inertia);  • Add as many slaves as needed to achieve the intended geometry: p.addSlave(Vec3D pos, double radius); The position of slaves is defined relative to the centre of gravity and in terms of principal axes.
Contact forces between slave particles of a multisphere and other bodies (not belonging to the same multisphere) are calculated the same way as for any other particle, but the resulting forces and torques are applied to the multisphere's centre-ofmass. The response of the multisphere is ultimately determined by solving the equations of motion of a rigid body [57] for its (linear and angular) accelerations. Fig. 17 depicts the application of multispheres to simulate a compression test of railway ballast material (in 2D). As can be seen, each particle is composed by small spheres delimiting the desired geometry.

Deformable/breakable clusters (agglomerates)
This new feature of MercuryDPM allows the user to create agglomerates (or clusters) composed of individual elementary particles. Clusters are formed by radial isotropic compression, which causes the cluster particles to adhere to one another, as shown in Fig. 18, but their relative position is not fixed, making the agglomerates deformable and breakable. This feature is useful in simulations where such properties are required, such as particle breakage [61], tableting [62], granulation, simulation of clay particles [63], etc.
The LinearPlasticViscoelasticSpecies [32] allows particles to be in mechanical equilibrium despite having a finite overlap, and a proportional finite tensile force is needed to pull them apart. The latter is what keeps agglomerates together, but also allows them to be deformed and broken when sufficiently strong external forces are exerted. As displayed in Fig. 18, clusters are mechanically stable before (left) and after (right) deformation.
Cluster radius R and mass fraction ζ follow an analytical relation dependent on the number N of elementary particles and their plasticity φ [32]. Since Eq. (5) is continuous, the rearrangement of particles inside of clusters due to an external (e.g. uniaxial compressive) force will be steady and gradual. This behaviour characterises plastic materials such as clay and rubber, and is only suited to model such deformation processes. However, to model brittle breakage, the elastic energy release must be sudden and trigger a fracture propagation along the body. To recover a similar behaviour the cohesive branch of the force in (5) must be discontinuous, and for this modelling purposes is substituted by where δ crit is a critical minimum overlap below which the cohesive bond is broken and ϵ is a new dimensionless material parameter. How this parameter affects the breakage behaviour is depicted in Fig. 19 where clusters are tested during uniaxial compression. Full details of this model are being prepared for a separate publication.

Writing applications
To write a MercuryDPM simulation, the class Mercury3D is used: this class contains the algorithm for time integration, contact detection, etc, and containers to store the elementary objects such as particles, walls, boundaries, contact models, etc. To make a new process simulation, the user creates a source file ('driver code'). In this file, one defines an (empty) object of type Mer-cury3D, and defines all the elementary objects that define the process he wants to simulate. One then calls the member function solve(), which calls setupInitialConditions() and continues the simulation. A typical user code is shown below: The distinction between process parameters (contact law, time step, final time, domain size, etc.) and geometric setup (walls, boundary conditions, particle positions) is due to our parallelisation strategy, which requires process parameters to be set first.

General interfaces
To store elementary objects, such as particles, walls and boundary conditions, MercuryDPM uses a series of handler classes. The ParticleHandler, for example, stores all types of particles (spherical, superquadric, etc.), the WallHandler all types of walls, etc. This is shown in the top left of Fig. 20.
All objects in a handler share a common base class. This ensures that the syntax for all objects and handlers is the same.
For example, BaseParticle contains the common properties of all particles, such as position, orientation, and velocity; the same member function, getObject(int), can used to access an object in the Particle-, Wall-, or BoundaryHandler; and the same function, getID(), is used to access the unique id of any particle, wall, or boundary. This is shown in the bottom left of Fig. 20.
Using the inheritance structure, the user can easily define new classes of elementary objects: For example, to define a sinusoidally-shaped wall, the user creates a new class SineWall, inherited from BaseWall, introduces parameters such as amplitude and oscillation frequency, and defines the member functions, such as the getDistanceAndNormal function.

Code samples
Simple parameters can be defined using set-functions; a typical user code specifies the following process parameters

Common geometries
For most processes, the user has to define the full geometric setup from scratch. However, certain setups are so common that we have predefined them, using classes derived from Mercury3D (see right of Fig. 20). The Chute class, for example, contains a function to create an inclined plane, which can be rough or smooth, and has predefined periodic boundaries that allows the user to quickly setup a periodic chute flow simulation. Similarly, ChuteWithHopper can be used to simulate a chute with an inflow hopper. We recommend users to define their own classes with predefined setups, e.g. for parameter studies where simulations only vary slightly, and thus avoid code duplication. An application of the Chute class is shown in [35].

Restarting
Each Mercury3D class has a write function, which stores the current state of a simulation in a text file; and a read function, which reloads the written state of a simulation. This allows simulations to be restarted. This functionality can be executed via a command-line interface: for example, by calling the executable HourGlass2DDemo, a simulation of two seconds is launched.
One can now restart this simulation and run it for a further two seconds by executing the command HourGlass2DDemo -r HourGlass2DDemo.restart -timeMax 4.

Parallelisation
Although MercuryDPM performs well for DPM simulations, due to advanced contact detecting and clever treatment of the walls, the computational power of a single processor (or thread) is limited. Thus, sequential DPM computations are limited to at most a few million particles and a few minutes of process time. In order to finish the simulation in a reasonable time, for larger computations, parallel processing is required. Currently MercuryDPM uses a domain-decomposition based parallel computing algorithm utilising MPI [64]. The current domain decomposition is simple: the process domain is decomposed into n x -by-n y -by-n z sub-domains of equal size (as specified by the user), and each processor computes the movement of particles in one sub-domain.
To determine the contacts with particles from neighbouring subdomains, a communication zone is established in the vicinity of the sub-domain boundaries, in which the processors communicate via MPI the location of the particles to their neighbours. This parallel computing algorithm can handle complex boundaries such as periodic boundaries, insertion/deletion boundaries and maser boundaries [4].
The performance of the parallel algorithm has been tested for the case of a rotating drum of varying width; a snapshot of the simulations is shown in Fig. 21. The serial algorithm shows a near-linear scaling of computing time with the number of particles (Fig. 22a). Weak scaling of the parallel implementation is shown by measuring the efficiency E = T p /T s , the ratio of computing time for the parallel and serial implementation, for a varying number of cores N c , keeping the number of particles per core constant. On a single node, efficiency decreases slowly with the number of cores, but levels off at around 40% for simulations that use hyperthreading (Fig. 22b). On multiple nodes, hyperthreading can be avoided, and the efficiency remains above 60% (Fig. 22c). Thus, the algorithm performs very efficient for large simulations, if the computational load per core is homogeneous.

Visualisation
There are two programs to visualise MercuryDPM output: xBalls and ParaView. xBalls, written by Stefan Luding, is a simple X11-based viewer that allows the user to quickly check the progress of the simulation. It is automatically installed with Mer-curyDPM; to visualise a simulation such as HourGlass2DDemo with xBalls, one simply needs to execute a script file that is part of the default simulation output, in this case HourGlass2D Demo.xballs. A more detailed three-dimensional visualisation of the walls and particles can be obtained via ParaView. For more information, see the MercuryDPM documentation at https: //docs.mercurydpm.org.

Versioning
MercuryDPM is available for download at http://mercurydpm. org. One can download either the latest release, or the developer's version (''Trunk''). The Trunk is updated as soon as a new feature is complete and is intended for developers only. After six months in the Trunk (where the developers' community will be able to debug the feature), a feature is considered save to use and ready to be merged into the next release.

Self-test suite
Developing new features can have unintended consequences.
For example, introducing a new variable in DPMBase could accidentally break the ability to restart simulations. To avoid breaking existing code by introducing new features, MercuryDPM uses the CTest software: Before any new code is committed to the Trunk or Release, the developer calls the command make fullTest, which (a) checks whether all codes in MercuryDPM compile, and (b) runs a series of self-and unit-tests to validate that no existing feature has been broken. Unit tests are designed to test a certain feature (e.g. whether the restitution coefficient is computed correctly) and return true if the test was successful; these are basic simulations that should run in less than one second. Selftests validate more complex features (e.g. restarting), and checks whether the output files have changed; these are slightly more elaborate simulations that should run in less than 10 s. There are now more than 300 unit-and self-tests in current developer's version of MercuryDPM. To ensure that each feature is tested, new tests have to be committed for each new feature.

Documentation and tutorials
The documentation of MercuryDPM is available at docs.mer curydpm.org. All classes of MercuryDPM are documented here, using the Doxygen suite, which extracts documentation from comments written by the developers in the MercuryDPM source files. In addition, the website contains tutorials, a list of demo codes, and a basic manual that will help new users and new developers to get acquainted with MercuryDPM.
Finally, training materials for both C++11 (and extensions), particle simulations in general, and MercuryDPM is freely available from the MercuryDPM website.

Particle-solid interaction
MercuryDPM can now be coupled with oomph-lib, an objectoriented, open-source finite-element library for the simulation of multi-physics problems [65]. This allowed the development of new features, such as surface coupling with FEM walls. The coupled code can simulate walls that deform in response to the forces exerted on them by the particles, and the particle motion updated by the walls, as shown in Fig. 23.
For surface coupling, we implement wrapper classes for the solid elements (and the governing equations), provided by the oomph-lib library, for efficient, I/O-free access and assignment of the nodal/contact forces and displacements. The oomph-lib geometry is mapped onto triangulated MercuryDPM walls that can interact with the discrete particles. For each particle-wall interaction, the contact forces are added as external loads into the weak form of the linear momentum equation, and the wall positions and velocities in MercuryDPM are updated by the finite element approximation.

Multi-resolution particle-fluid coupling
A second feature enabled by the coupling with oomph-lib is fluid-particle coupling. For under-resolved simulations, we use the Anderson and Jackson formulation [70]. This introduces a voidage field, a measure for the fraction of total particle volume inside an element, to simulate the fluid flow. Coupling forces are be defined that specify the interaction between particles and fluid [71]. A different coupling method is semi-resolved, giving more detailed results than an under resolved method but being much faster than a fully-resolved method. The voidage can be described as continuous function by coarse-graining the particles in space. The fluid elements can then simply evaluate this function at their location. For fully resolved simulation there is a noslip boundary condition for the fluid on the particle surface, and the coupling forces can be computed by integrating the pressure along the particle surface.
At the moment the three methods have been independent implemented; however, a multi-resolution is currently under development, where the code adapts between fully-, semi-and under-resolved situations, allowing the simulation of suspensions with arbitrary large particle size-distributions. This is illustrated in Fig. 24.

Better hybrid openMP-MPI parallelisation
The current MPI parallelisation strategy, see Section 7.5, works well for evenly distributed particle systems. However, for inhomogeneous systems, the workload of the processors is unbalanced, resulting in suboptimal scalability. This weakness of the current implementation will be addressed in the near future: one possibility to enhance load balancing is using a cyclic distribution with respect to particle identities (indices); a second possibility is to implement an adaptive mesh of differently-sized domains. These parallel-processing strategies for the MercuryDPM software package will be done using a combination of OpenMP and MPI for CPU parallel computing and OpenACC/CUDA for GPU computing.

STL/STEP readers for reading in industrial geometries
We have an STL reader for MercuryDPM but this by-passed our very nice complex curved wall support. We are currently working on an STEP reader that keeps all the curvature information, leading to quicker and more accurate simulations.

Calibration via grain learning
Calibrating contact models is a huge problem in granular materials. A wide and varied range of granular characterisation machines exists, and many of these machines produce vendorspecific, none-standard measures that are not comparable with each other. To solve this problem, MercuryDPM is integrating Grain Learning (https://github.com/chyalexcheng/grainLearning/) into its software suite. Grain learning uses Bayesian inference and machine learning to train the contact model to reproduce the supplied characterisation data [72,73]. The algorithm further identifies automatically whether the provided characterisation data is sufficient to uniquely determine the parameters of the contact model, or whether more experimentation is required. A major advantage is that the method is independent of the type of characterisation data; this makes it a very valuable tool and it will remain so until a standardisation of the measures characterising granular materials has been established, i.e., until we have granular properties that mirror surface tension, viscosity, etc. in fluid dynamics.

Release strategy
Originally, MercuryDPM has been released once a year; however, this was becoming less practical due to the large number of contributors, so we have moved to an open-development model, i.e. opening the developer's version to public download. For more information about MercuryDPM please visit http://MercuryDPM. org.