CLUMP: A Code Library to generate Universal Multi-sphere Particles

Particle shape plays a key role in the mechanical and rheological behaviour of particulate and granular materials. The simulation of particulate assemblies typically entails the use of Molecular Dynamics, where spheres are the predominant particle shape, and the Discrete Element Method (DEM). Clumps and clusters of spheres have been used to simulate non-spherical particles, primarily due to the simplicity of contact detection among spheres and their ability to approximate practically any irregular geometry. Various approaches have been proposed in the literature to generate such clumps or clusters, while open-source numerical codes applying these are scanty. The CLUMP code, proposed in this paper, provides a unified framework, where a particle morphology can be approximated using different clump-generation approaches from the literature. This framework allows comparing the representations of the particle generated by the different approaches both quantitatively and qualitatively, providing the user with the tools to decide which approach is more appropriate for their application. Also, one novel generation technique is proposed. Outputs are provided in formats used by some of the most popular DEM codes. Moreover, the resulting clumps can be transformed into surface meshes, allowing for easy characterisation of their morphology. Finally, the effect of clump-generation techniques on the mechanical behaviour of granular assemblies is investigated via triaxial compression tests.


Motivation and significance
Particle morphology is believed to be the controlling factor in spreading of powder bed in additive manufacturing [1], wave velocity anisotropy [2], railway track sub-structure design [3] and pharmaceutical tableting [4], just to name a few. Modelling real particles as clumps of overlapping spheres or clusters of non-overlapping spheres is common practice in both academic and industrial studies, when particle morphology is of interest. Clustering spheres to represent a single particle offers versatility within a simulation, as the fine details of particle morphology such as local roundness, roughness or concavities can be introduced. Moreover, the computational cost is affordable due to the simplicity of performing contact detection among spheres, given the total number of spheres considered. An abundance of methods has been proposed in the literature to perform such  a feat, including Favier et al. [5], Matsushima and Saomoto [6], Bradshaw and O'Sullivan [7], Price et al. [8], Wang et al. [9], Garcia et al. [10], Ferellec and McDowell [11], Taghavi [12], Gao et al. [13], Li et al. [14], Zheng and Hryciw [15], Haeri [16], and Katagiri [17]. These methods have been employed in multiple research fields, namely physics [18], pharmaceutical engineering [19], chemical engineering [20], agriculture [21], civil engineering [9,22], mining engineering [23] and geosciences [24]. However, only a limited number of implementations are available open access, such as Haeri [16]: https://github.com/sihaeri/DEM-ClumpedSphere and Bradshaw and O'Sullivan [7]: https://github. com/mlund/spheretree.
In a typical numerical study employing clumps, the researcher employs an available method to approximate the real particle morphology. Most commercial DEM software offer one specific approach, developed in-house, to assist users in defining such clumps. For example, PFC3D (Itasca) provides the Taghavi [12] method and EDEM (Altair Engineering) offers an in-house method in their most recent version. Open-source DEM codes (such as YADE, LAMMPS, etc.) often do not provide any methods to generate clumps, although YADE offers functions to fill a target geometry with non-overlapping spheres. Note that the quality of the particles generated by the available approaches varies with the type of application/problem investigated and with particle morphologies. At present, a quantitative analysis of the effect of clump generation methods is missing from the literature. Also, no or very little information is provided in the publications that introduce a new method for generating a clump about the performance of the method in preserving morphological characteristics. To this end, CLUMP provides a unified framework for the comparison of different multi-sphere particle generation methods and it facilitates morphology characterisation of these particles. The code is easy to use and it can provide an automated workflow to generate thousands of particles within reasonable time runs.

Software description
The code comprises three main modules: (1) Functions to generate clumps, (2) Functions to export clumps and (3) Transformation of clumps to surface meshes, allowing for their morphological characterisation and comparison with the original particle. The first module is the main core of the software, while the second and third ones have an auxiliary role, functioning as utility scripts.

Software architecture
This section details the main structure of the code, which is also visualised in Fig. 1. Module I: Generate clump. CLUMP supports three approaches to generate clumps, namely Favier et al. [5], the first multi-sphere approach proposed in the literature; Ferellec and McDowell [11], a widely used approach; and a newly proposed approach based on the Euclidean distance transform of 3D images. Fig. 2 shows the main objects of the clump generation module.
Module II: Export clump. The generated clumps represented by the centroids and radii of multi-spheres can be exported in various formats (*.py, *.csv), compatible with some of the most mainstream commercial and open-source DEM codes, including EDEM, PFC3D, YADE, and LAMMPS.

Module III: Characterise clump.
A surface extraction module is developed, transforming each clump into a surface mesh, which in turn can be imported in SHAPE [25], to perform shape characterisation.
CLUMP relies on several external functions available within the Matlab FEX community. In particular, the outsourced tasks include the handling of stereolithography files [26], mesh manipulations [27], and calculation of rigid body parameters [28].

Software functionalities
CLUMP supports the following clump-generation approaches: 2.2.1. Favier et al. [5] This method only applies to symmetrical particles which are constructed of spheres whose centres are located on the particle axis of symmetry. Spheres may overlap and may vary in diameter along the length of the axis of symmetry. The surface of a particle is approximated by inscribing spheres such that the surface of each sphere is tangent to the surface of the particle at the point of contact. The position of each element sphere is fixed relative to the other elements within a particle.   [11] The clump generation methodology of Ferellec and McDowell [11] selects random points of the mesh of the particle surface and generates tangent spheres of increasing radii, until they intersect a different vertex of the particle surface. This methodology was designed to eliminate artificial roughness of flat faces, through the consideration of large numbers of spheres per particle. Due to this, the approach of Ferellec and McDowell will not work as well if we want to generate 2 or 3 spheres per clump, as the selection of the initial points is random. This happens because the methodology lacks a deterministic criterion of which vertex to choose every time to generate spheres. In our implementation of this method, a seed parameter is introduced in the selection of vertices to achieve reproducible clumps. Though, this seed is not guaranteed to lead in any way to the best-possible clump, for the given number of spheres.

Ferellec and McDowell
In more detail, the volume of the real particle is filled optimally with spheres of different sizes: from a point chosen randomly on the surface, a sphere is grown internally along the normal vector at that point to the maximum extent possible inside the boundary of the particle. In other words, the expansion of the sphere continues until it reaches another point on the surface of the particle. Then, the process is repeated for other points on the surface of the particle, which must be farther by a distance d min to any existing sphere. The number of spheres inside the clump is directly related to the number of points on the surface, by monitoring a percentage p max , which if met, terminates the generation process. If the density of points on the surface of the particle is high enough, then the process can theoretically lead to a clump composed of thousands of spheres, although Ferellec and McDowell [11] set criteria that can reduce local crowding of spheres, using a minimal radius r min , or the aforementioned d min and p max variables. If a large number of spheres per clump is computationally acceptable, this approach can lead to smaller artificial roughness of flat faces, compared to other methods, as Ferellec and McDowell report [11]. This approach was designed to be applied on surface meshes, but is also applicable directly on point clouds, since techniques do exist to estimate the normal vector of each member point of a cloud.

Euclidean distance transform
A new approach is presented in this paper to generate multisphere particles, based on the concept of the Euclidean distance transform of 3D images. Initially, the particle morphology is transformed from a surface mesh to a voxelated representation of a minimum dimension div, where the value of all voxels belonging to the particle is set equal to zero. Then, the inscribed sphere of the particle is found as the maximum value of the Euclidean distance transform of the voxelated image. The voxels corresponding to the inscribed sphere are then set equal to one and the Euclidean distance map is calculated for the new voxelated image. This process is repeated until a user-defined number of spheres is found (N) or until the user-defined minimum radius (r min ) has been reached, as the spheres are generated in decreasing sizes. This methodology can also generate overlapping spheres, if only a percentage of the voxels making each new sphere is set equal to one, instead of all of them, expressed by a variable overlap, which takes values within the range [0, 1). Fig. 3 shows the algorithms behind the Ferellec and McDowell [11] and Euclidean distance transform implementations. The output of the above approaches is a set of spheres, represented by their centroids and radii, as needed in commercial and opensource codes. More information on these approaches is provided as supplementary material.

Particle surface extraction
One of the main objectives of CLUMP is to allow for a direct comparison among the morphologies of multi-sphere particles generated using different approaches. To achieve a quantitative comparison, particle shape-characterisation tools can be used, such as SHAPE [25]. To this end, the surface of each clump has to be extracted in the form of a surface mesh, allowing for compatibility between the clump generation and the shape characterisation codes.
Extracting the surface of a multi-sphere particle is not a trivial task. The problem can be tackled using many different approaches. In this implementation, the surface of each clump is extracted following the steps below: 1. Contact detection is performed between all possible pairs of the spheres making the clump, identifying the pairs with geometric overlap. 2. A point cloud is generated on the surface of each sphere. 3. For each pair of overlapping spheres, the points within the overlap region are deleted (Fig. 4). 4. A circle is found as the intersection of each pair of overlapping spheres and points are generated along each circle, which are added to the point cloud.
5. The point cloud is then tessellated using a surface reconstruction technique, making the surface of the clump.
It should be noted that generating points on the surface of each sphere introduces mesh-dependency in the calculation of local shape parameters, like roundness and roughness. Therefore a dense mesh should be chosen when these shape aspects are of interest. The local morphology can also be affected by the technique chosen to perform surface reconstruction, which in this case is an implementation of the Crust method [29] provided by Giaccari [30]. The points on the circles calculated as the overlap of each pair of spheres are added to the point cloud, aiming to minimise the numerical noise of the clump surface near the region of an overlap and provide a smooth transition from the surface of a sphere to its adjacent ones. Fig. 4 shows the surface of a clump made of two spheres, where the points inside the overlap region (removed from the point cloud) and on the intersection circle (added to the point cloud) can be observed.

Sample code snippets analysis
This section highlights two implementation details of CLUMP. First, Fig. 5(a) demonstrates how the newly introduced approach based on the Euclidean distance transform allows the generation of overlapping spheres, while Fig. 5(b) shows a typical loop of how a tangent sphere is grown in the methodology of Ferellec and McDowell [11].

Illustrative examples
This example demonstrates the mechanical behaviour of assemblies of rice grains and soil grains under triaxial compression. The morphological characteristics of clumped spheres used in this illustrative example are presented in Table 1. Hereinafter, EU, FM and FA refer to grains represented by 25 spheres using the Euclidean transform, Ferellec and McDowell [11] and Favier et al. [5] methods, respectively. The samples consist of (1) rice grains using all three approaches, and (2) sand grains using the Euclidean transform and Ferellec and McDowell [11] approaches. The major axis length of the rice grain and silica sand was approximately 18 mm and 6 mm, respectively. The particle shear modulus and particle Poisson's ratio were considered as 20 MPa and 0.1 for rice grains and 29.1 GPa and 0.23 for sand grains, respectively. The particle density was 1470 kg/m 3 for the rice grains and 2650 kg/m 3 for the sand grains; these values were increased by 1000 times to save computational time.
The LAMMPS molecular dynamics code [31,32] is employed with a modified servocontrol, as used by Hanley et al. [33] and Morimoto et al. [34]. Each sample is composed of 5000 clumped grains, which are generated randomly without initial contacts with other particles, under no gravity conditions. The sample is subjected to an initial isotropic compression by moving the boundaries with a controlled velocity, to achieve an isotropic stress of 50 kPa. A simplified Hertz-Mindlin contact model with Coulomb friction is adopted. Dense, medium and loose samples are prepared using inter-particle friction coefficients (µ) of 0.0001, ≈0.1 and 0.35, respectively, during the isotropic compression, giving relative density (Dr) values of 100%, 50% and 0%, following the approach of Morimoto et al. [34]. In the subsequent triaxial compression, the µ value was set to 0.35 for all the cases as a material constant. In the present simulations, no damping was used.
The cubical samples composed of rice grains in Fig. 6(a) were subjected to quasi-static axial compression in the Z -axis while keeping the lateral stresses in the X-and Y-axes at 50 kPa. The overall trend of stress-strain relationships is similar for each density level (Fig. 6(b)); however, the FA Dr=100% case underestimates the peak and residual stresses, while the FA cases overestimate the initial stiffness, i.e. the initial increase in deviator stress (Fig. 6(c)). A clear difference between FA and the other approaches is evident in the variation of mean coordination number (Fig. 6(d)). A representative sample composed of sand grains is illustrated in Fig. 7(a). Similar variations of stressstrain responses are observed between the EU and FM approaches ( Fig. 7(b)), although measurable differences in the variation of void ratio are evident (Fig. 7(c)).
For the illustrative examples presented here, the EU and FM approaches gave comparable results, while the results of the FA approach differed. Referring to Table 1, this discrepancy can be attributed to the difference in the shape parameters of the generated multi-sphere particles. EU and FM provide similar shape parameters, indicating that quantifying shape parameters of the generated clumps is recommended as an effective measure to ensure the quality of modelling non-spherical grains. On the other hand, FA cannot represent particles with flatness, as the axial symmetric clumps it generates can only result in compact or elongated particle shapes, a fact that shows a distinct effect on the macro-mechanical behaviour of the simulated triaxial tests on rice. Interesting to note that apart from flatness, the rest of the studied morphological parameters between FA, and EU and FM are similar.

Impact
CLUMP enables engineers, researchers and scientists from a variety of disciplines and industries to examine the dependency of their application to particle shape effects, through the generation of irregular particle shapes. Recent advances in granular physics demonstrate a clear need for consideration of particle shape in numerical modelling. In engineering practice, the clumped-sphere approach has been adopted in additive manufacturing, railway engineering, pharmaceuticals, etc. However, there is currently little evidence regarding how many spheres need to be generated. This is, mainly, due to the limited accessibility of open-source implementations of different clump-generation approaches. To this end, CLUMP offers two widely-used approaches to generate clumps and proposes a new one, aiming to provide an extendable library of available scripts. The implementations are computationally efficient, i.e. it takes few minutes to generate hundreds of spheres per grain. CLUMP provides an option to extract the surfaces of the generated multi-sphere particles, allowing for the direct calculations of their shape descriptors, using available, open-source shapecharacterisation codes. Comparing the morphological characteristics of the generated clumps to those of the original particle morphology, one can decide the number of spheres and method of application based on a quantitative and rigorous particle-generation approach.

Conclusions
We present CLUMP, a code library to approximate threedimensional particle morphologies using different multi-sphere generation methods. The software allows exporting the results in formats compatible with several commercial and open-source DEM codes. CLUMP enables to assess the suitability of a clumpgeneration method for different applications and characterise the degree of manipulation of the original morphology in a quantitative manner.
The examples illustrated provide insight into the effect of the clump generation methods on the mechanical behaviour of two granular assemblies subjected to triaxial compression. We found that the adoption of particles with near zero flatness leads to the underestimation of the peak and residual stresses, while the adoption of larger contacting spheres leads to the overestimation of initial stiffness.
We hope code developers will find the framework provided here useful so that other existing clump-generating approaches and future ones will be implemented too.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.