Recent Advances in MercuryDPM

In this paper we introduce the open-source code MercuryDPM: a code for simulating discrete particles. The paper discusses software and management issues that may be interesting for the developers of other open-source codes. Then we review the new features that have been added since the last publication: an improved Hertz-Mindlin model; a new liquid bridge model of Lian and Seville; a droplet-spray model; better support for re-creating complex, measured particle size distributions; a new implementation of rigid clumps; an implementation of elastic membranes; a wear model for walls; a soft-kill feature and a cloud-deployment interface for AWS.


About MercuryDPM
MercuryDPM was started in 2009 by Anthony Thornton and Thomas Weinhart with the aim of creating a discrete particle method software able to solve complex industrial scenarios. This required several features which now form the backbone of MercuryDPM: (i) A flexible implementation allowing complex wall and boundary conditions, (ii) a neighbourhood detection algorithm capable of dealing with highly polydisperse particle packings [8], and (iii) an analysis tool able to extract the most relevant information from the huge amount of data generated by these simulations [9]. The code also has been coupled to the continuum solver oomph-lib [10] to simulate particle interactions with elastic solids, as well as multiscale coupling to simulate granular materials in a computationally efficient way [11].
The code has been open-source since it was started. We started with a GPL license but very quickly moved to a BSD 3-clause license, as this was felt to be more open to external development and simpler to understand. Since its conception, both its user and developer base has grown, with 48 people so far contributing significantly to the code base. Information about past and current contributors can be found on the team page of the MercuryDPM website, www.mercurydpm.org.
MercuryDPM is a versatile, object-oriented C++ code which (we hope) is easy to understand. It is regularly tested on several Linux distributions, Mac OS and Windows 10. To avoid breaking already existing code, a suite of over 260 self-tests have been developed, testing each feature of the code. Developing new applications in the software is straightforward: The user specifies the particulars of their simulation (initial positions, inflow, outflow, walls, interaction parameters) in a single driver file, which calls the MercuryDPM kernel to execute the simulations. All kernel features are documented, and there are many sample driver codes demonstrating the features.
When the code was first started, we used a self-hosted svn server. However, as the user and developer base has grown this was not maintainable; therefore, in May 2022 we moved to a git repository on bitbucket. This version of the code can be found at https://bitbucket.org/mercurydpm/mercurydpm. Moving to git has accelerated the return of features from developers outside the core team.
Building the code is managed using CMake, the test suite via ctest, for code maintenance (bug reporting/tracking, release planning, etc) we use the Atlassian tools: Jira and Confluence. For visualisation we use both Paraview [12] and an in-house code from Stefan Luding: XBalls, which is now distributed with MercuryDPM.

Release Strategy and Version Number
MercuryDPM has a hybrid development pattern: Firstly, we have an open-development model, where the master branch (which includes the newest features) is openly available; plus most of the active feature branches, so code can be seen and used ahead of final testing and journal publication. We encourage new developers to use branches in the main repository so others can see their features in development, but we do not insist on this. The master branch can be accessed at https://bitbucket.org/mercurydpm/mercurydpm/src/master. Secondly, we have returned to doing stable release versions, approximately once a year. After deployment of a release, it will remain in an Alpha state up until 28 days of the last bug fix. At the alpha stage the release notes are drafted. The following 91 days it will be in Beta state until it is considered stable. From version 1.0.0 onward, MercuryDPM will use the following convention of three numbers, a.b.c: The version number, a, is incremented when an interface changes; the second number, b, when new features are added; and the third number, c, for a bug fix to a released version. This scheme complies to semantic versioning rules and for the future, we aim to automate deployment by introduction of conventional commits. The current release can be accessed at https://bitbucket.org/ mercurydpm/mercurydpm/src/1.0.Alpha. 1

Hertz-Mindlin (Improved) No-slip Contact Model
This section discusses changes made to the Hertz-Mindlin contact model, which was changed for 1.x series. The changes are best demonstrated by the driver codes Hertzian2DSelfTest.cpp and MindlinSelfTest.cpp.
The Hertz-Mindlin contact model uses Hertz theory to determine the normal elastic force between contacting spherical particles and a tangential force model established by Mindlin and Deresiewicz [13]. Additionally, rolling friction can be added according to the model of Luding [14]. All three models have been slightly modified:

The Normal Force Model
The MercuryDPM contact model HertzianViscoelastic applies a Hertzian normal force between two particles, given by where the normal stiffness k n and normal dissipation coefficient γ n are and Here, δ n is the normal overlap, v n rel the relative velocity in normal direction, E * the effective Young's modulus, R * the effective radius and M * the effective mass. β n is the normal damping factor which is related to the restitution coefficient as follows, (4.4) Change: Previously, the user directly specified the effective elastic modulus E * using the command setElasticModulus. This caused confusion because most users expect to specify instead Young's modulus E, and Poisson's ratio, ν, and compute the effective elastic modulus as E * = 1 2 E/(1 − ν). Thus, we replaced the function setElasticModulus with setEffectiveElasticModulus and added a helper function, computeEffectiveElasticModulus, to compute E * from E and ν.

The Tangential Force Model
The MercuryDPM contact model Mindlin applies a simplified version of the Mindlin model of sliding friction, as proposed by Di Renzo and Di Maio [15], where the tangential force is given by where the tangential stiffness k t and tangential dissipation coefficient γ t are determined as and In these equations, δ t is the tangential overlap, v t rel is the relative velocity in tangential direction, G * is the effective shear modulus and β t is the tangential damping coefficient.
Change: The dissipation coefficient in Eq. (4.7) previously read γ t = β t √ 8 M * k t . The factor 8 in this expression was removed to allow for the following more compact notation of the tangential damping factor Note: β t is currently an input parameter that the user must set. In the future, a β t will be set to a default value according to Eq. (4.8), and can be overridden by the user. This will change for the 2.x series.

The Rolling Friction Model
The MercuryDPM contact model MindlinRollingTorsion applies a rolling friction according to Luding's model [14]: When a rolling friction coefficient is specified by the user, a corresponding rolling torque is applied, which is calculated as Sketch of a liquid bridge between two unequal-sized particles of radius R 1 and R 2 . Here, θ is the contact angle, a material property, and s the separation distance between the two particles where F ro is the rolling torsion force. Change: In previous versions, the effective diameter was used as a prefactor in (4.9), rather than the effective radius, for calculating the rolling torque. This model has been updated to correctly satisfy Luding's model [14].

Liquid Migration Model of Lian and Seville
The majority of research on granular media has concentrated on dry granular materials. However, in industry and nature, we often encounter wet granular materials. An example currently being investigated with MercuryDPM is the collision between wet (sand) particles, including moist sediment transport by wind and pneumatic conveying. Wet granular materials are cohesive due to surface tension, which is a significant distinction between dry and wet granular materials.
The MercuryDPM contact models LiquidBridgeWillet and LiquidMigrationWillet already have an implementation of liquid bridges. These are based on work of Willet et al. [16], which was developed by experimental data where forces due to liquid bridges between particles were measured.
We have now added a second model, LiquidMigrationLS, which is based on work of Lian and Seville [17], which includes a lubrication viscous force. The viscous force is assumed to be the dominant force during the contact, induced by the squeezing out and pulling in of the liquid. The model was validated in particle simulations against collision experiments between a wet particle and a flat surface by Zhang and Wu [18], and is applicable to describe the dynamic interaction between wet particles; whereas, the previous model was only strictly correct for static (or slow moving) scenarios. Figure 1 shows the capillary liquid bridge between two particles of different sizes.
This code is available in the master version of MercuryDPM and is best demonstrated by TwoParticleElast icCollisionInteractionWithLiquidMigrationLSSelfTest.cpp in the source directory Driver s/SelfTests/Interactions.
It keeps the liquid migration features of Willet's model, but differs in the force calculation: The capillary force is based on the half-filling angle which depends on the value of the contact angle and liquid bridge volume. The viscous force exists both in normal and tangential directions of potential contact, and is active only between a limited separation distance (smaller than rupture distance) and the rupture distance. Its magnitude depends on the relative velocity of the particles with an interstitial liquid bridge.
Within the LiquidMigrationLS model, one can define the liquid viscosity, based on the liquid type they model through setViscosity(). Similar to LiquidMigrationWillet, one also needs to define the volume of the liquid carried by the particles and the maximum and minimum value of the liquid bridge volume. In Fig. 2, particles connected by liquid bridges are visualised in ParaView.  As stated above, many industrial applications include wet granular materials; often, the moisture is added by a spray. In MercuryDPM, the DropletBoundary has been implemented, which can insert liquid droplets to wet the particles during a simulation. The implementation is simple: Each time step, droplets get inserted according to a function specified by the user. Droplets have a mass and velocity, and move under the influence of gravity. They do not interact with each other, but if a droplet contacts a particle, the particle absorbs the liquid and the droplet is removed. Note, as walls currently cannot absorb liquid, the liquid is lost when the droplet contacts a wall.
The user can define freely how droplets are inserted via the function setGenerateDroplets; one can e.g. set the rate of insertion, the droplet radii, the position of insertion (nozzle position and angle), and the velocity of insertion (spray angle and droplet speed). For example, the user can add a droplet boundary that acts as a spray nozzle. In Fig. 3, we show one of the most common nozzles for spraying droplets, a flat fan spray nozzle.
This code is available in version 1.x of MercuryDPM and the figure is created by NozzleDemo.cpp in the source directory Drivers/Demos/IndustrialMixers.

Dealing with True Particle Size
The InsertionBoundary in MercuryDPM was reworked for version 1.x. It is now capable of inserting particle mixtures composed of a variety of different materials and particle size distributions (PSD) into a single insertion volume. Moreover, the new interface accepts different types of PSD functions, based on volume, length, surface area or number of particles. This will allow for a more accurate representation of particle mixtures (see Fig. 4), which is crucial for several applications and phenomena, such as segregation and mixing. Figure 4 is created by a Fig. 4 Cumulative PSD of inserted particles, compared to the input PSD script comparing a cumulative number PSD to the same one inserted by the driver code PSDSelfTest.cpp. The script is available in the master.
Furthermore, polydisperse particle packings used to be insertable only into cubic volumes. However, several applications require the insertion of particles into more complex geometries. Therefore, the insertion routine was reworked to allow the geometry of the fill volume and the properties of inserted particles to be set independently. It is now straightforward to create complex insertion regions. For a new insertion boundary object, only the geometry has to be defined.

Dealing with True Particle Shape-Rigid Clumps in MercuryDPM
Rigid clumps of spherical particles are an important tool to analyse the behaviour of granular materials consisting of particles of irregular shapes with the discrete element method. By rigid clump (just clump, or multiparticle) we imply the aggregate of N rigid (spherical) particles of a given density, that are rigidly linked to each other at given relative translational and rotational positions. Note, currently the code has only been demonstrated for spherical particles but in principle you should be able to form rigid clumps from more general super-quadric particles. The constituent particles of a clump will be referred to as pebbles. The pebbles may (or may not) have overlaps, introducing volumes within a clump that belong to more than one pebble. The contact detection algorithm treats contacts and corresponding forces/moments, as well as the forces arising from the force fields at the pebble level, while the dynamics computation algorithm treats a clump as a single rigid body, that is accelerated by the resultant force/moment, properly summed up over the pebbles.
The rigid clump function in MercuryDPM is implemented as a multilevel structure. A third-party library CLUMP [19] is used to generate positions and radii of pebbles that describe the given nonspherical shape (Fig. 5A). The CLUMP tool provides pebble data, which, along with the optionally provided initial STL-format shape of the clump, constitute an input of the MClump pre-processing tool (part of MercuryDPM). It centres and rotates the clump, aligning its principal axes with global Cartesian axes, and computes clump's tensor of inertia using the prescribed algorithm (summation over (non-overlapping) pebbles, summation over voxels, summation over tetrahedrons using STL representation (Fig. 5B)). A special header library for the driver file introduces necessary modifications of MercuryDPM virtual members, enabling clump dynamics. The driver file loads the list of clump instances generated by MClump, and, using them, generates necessary distributions of nonspherical particles and computes their dynamical evolution. Figure 5B highlights an example of using nonspherical particles in MercuryDPM (driver file Drivers/Multi Particle/BulkTs.cpp). The rigid clump feature is currently available in the master (developer's version), and will be included in future releases of MercuryDPM and is best demonstrated by the codes in the directory Drivers/MultiParticle.  5 (A) STL model of a complex particle shape, the corresponding rigid clump, generated by [19] and the voxel model. (B) Sample simulation -100 nonspherical particles in a gravity field, within a box with fully elastic walls

Membranes
Interactions between membranes and granular particles may occur in systems like triaxial tests but also in applications from the soft robotics community, such as the granular gripper.
One possibility to represent membranes within granular simulations is a mass-spring system. Due to its particlebased nature it is easily integrated in discrete particle simulations. Depending on the chosen mesh and the nature of the spring, it is possible to approximate all material properties [20].
MercuryDPM version 1.x now supports the use of mass-spring systems containing distance springs for the in-plane dynamics and a bending force penalty originating from cloth simulations [21]. If the connectivity of the masses is given by a triangular mesh with a hexagonal unit cell, the spring constant of the distance springs may be , where t is the thickness and E the elastic modulus of the membrane [22,23]. Note that the simplicity of this setup restricts us to a Poisson's ratio of ν = 1/3.
In MercuryDPM the mass-spring system is constructed by inserting particles in place of the masses and connecting these according to a given mesh connectivity, which may be specified via a STL file or given directly in the code. As the connectivity has to be specified, the implementation does not rely on any contact detection algorithm. Although the code supports the usage of the vertex particles to detect contacts between the membrane and the granulate, as it was originally introduced by de Bono et al. [24] for triaxial test simulations, it is advised to deactivate these by a proper choice of species. In that case, the contacts are instead calculated using additional triangular wall elements  inserted between the vertex particles. The position and dynamics of these wall elements are purely defined by the vertex particles. Quantities such as local velocities and the calculated contact forces are transferred between the triangles and the vertex particles using barycentric interpolation. Figure 6 shows an example of the triaxial shear test and is created by a private driver code. Figure 7 shows the granular gripper application from the published paper [25]. The membrane is best demonstrated by the code MembraneDemo.cpp in the directory Drivers/Demos/Membrane.

Modelling Surface Wear
Often in industry, granular flows lead to the wearing of surfaces over time due to abrasion, causing apparatus to misperform or even fail. One such apparatus is a vibrating screen, used to grade materials in various industrial applications including steel making. The left panel in Fig. 8 shows a typical vibrating screen from the front. Grains are dropped at the back of the screen and slide forward, falling through the increasingly widening gaps. As the material slides across, the surface is worn down, increasing the size of the particles that are able to fit through the screen. This feature will be available in future releases.
In MercuryDPM, the Reye-Archard-Khrushchov wear model has been implemented [26], which states that the volume of material removed is proportional to the work done by the friction forces. The proportional is normally set to higher values than for the real material in order to accelerate the process and reduce the simulation time. Currently, only for triangulated walls, the wear model can be used to update the surface and hence, predict the long term effect of particle abrasion on walls. The process is summarised as follows: 1. Compute the (local) frictional forces on each wall segment. 2. For each triangle compute the (local) volume to be removed via the Reye-Archard-Khrushchov model. 3. Iteratively move the mesh vertices to match each local volume of removed material.
In Fig. 8 you can see the shape of a vibrating screen before (left) and after (right) a granular material has flowed over it.

Soft Kill
MercuryDPM now has an internal signal handler, which allows us to control what happens when external signals are received by the code. This enables the code to be run on resources with limited time available, as it will allow the code to be stopped (and restarted) automatically based on external triggers. Thus, code can be deployed on cheap, unused, temporary computing resources like AWS spot instances.
The signal handler treats various interrupt signals that are sent to the processors. A common example is SIGINT, which can be sent by the user to a running program by pressing CTRL+C. This interrupt is caught by the signal handler, triggering a soft-kill routine: instead of stopping the code immediately, this feature allows MercuryDPM simulations to continue until all operations necessary to write the restart file are finished. This way the output data (e.g. *.restart, *.fstat, *.data) will be consistent (same timestep) and complete (no partly written files). The simulation can thus be restarted safely at a later time.
Currently this feature catches the following signals, and triggers the actions described: • SIGINT (signal 2): This signal is sent by CTRL+ C. Soft kill continues the simulation until the current time step is completed. Then it forces MercuryDPM to write to all output files, in particular the restart file. • SIGTERM (signal 15): This is generated in linux by kill. It triggers the same action as SIGINT.
• SIGKILL (signal 9): This is generated in linux by kill -9. It terminates MercuryDPM immediately. In other words, this signal stops the soft-kill feature.
Note: Using tools like htop and F9 in Linux, the user can send specific signals or kill commands to MercuryDPM. Also, one can use $kill -SIGNAL PID, where SIGNAL is the signal number listed above and PID is the process id of the MercuryDPM executable.

Mercury Cloud
MercuryDPM also has its own official spin-off company, MercuryLab, which aims to facilitate access to the advanced features of MercuryDPM for industry and academia. MercuryLab often develops new features for industrial clients, One of the developments of MercuryLab is its AWS-deployed cloud interface for the code, MercuryCloud (https:// cloud.mercurylab.org/). This interface allows scalable, GUI-based, on-demand access to (certain) features of the code via an internally developed web interface. The architecture of the cloud is shown in Fig. 9 and is available from any device with internet access and a web browser.
MercuryCloud is still in development, thus its capabilities are limited, but growing. Academic users are able to access the facilities at a considerable discount. Utilising the soft-stop feature described in the previous section, we hope to make use of cheap AWS spot instances, which are often at least 60% cheaper.