Volume-interacting level set discrete element method: The porosity and angle of repose of aspherical, angular, and concave particles

to highly angular, was generated and investigated. Results show that, when combined, angularity and concavity have a disproportionately large effect on the porosity and angle of repose. They also demonstrate the importance of particle shape and that even convex polyhedra may be insufficient to accurately capture the effect of shape on fundamental powder properties.


G R A P H I C A L A B S T R A C T H I G H L I G H T S
• Volume-interacting level-set DEM simulates highly angular and concave shapes.
• Continuous spectrum of particle shapes by hyperbolisation of the Platonic solids.
• Porosity and angle of repose increase with asphericity, angularity, and concavity.
• When combined, angularity and concavity have a large non-additive effect.

Keywords:
Angle of repose Porosity Concave Aspherical Angularity Roundness A B S T R A C T Despite the many advances in the numerical modelling of granular materials using the discrete element method (DEM), the description of realistic particle shapes, especially concave shapes, remains challenging.Little is known about the effect of concavity on properties such as the powder porosity and angle of repose.In this work, the volume-interacting level set DEM (VLS-DEM) was used to study the effect of asphericity, angularity, and concavity on the porosity and angle of repose.By taking the Platonic solids and applying a hyperbolisation procedure, a continuous spectrum of shapes, ranging from rounded to highly angular, was generated and investigated.Results show that, when combined, angularity and concavity have a disproportionately large effect on the porosity and angle of repose.They also demonstrate the importance of particle shape and that even convex polyhedra may be insufficient to accurately capture the effect of shape on fundamental powder properties.

Introduction
A wide range of industries, from food and pharmaceuticals to recycling and construction, are faced with the challenge of processing granular materials.The characterisation of granular materials is therefore also of major importance, as improvements in processing efficiency can result in large savings in energy and time.One property that is frequently used to assess the flowability of powders is the angle of repose, commonly defined as the steepest angle for which a freestanding surface of a pile of granular material is stable [1].Although there are arguably better and less ambiguous measures of flowability, such as the flow function coefficient (FFC), the angle of repose has the advantage of being relatively easy to visualise and measure.The angle of repose is influenced by many different factors such as the interparticle cohesion, friction, and shape of the individual particles.Similarly, the bulk porosity of a material can be equally important, e.g. to describe the dilation or densification of a powder during flow and compaction, and depends on the same factors.This work focuses on unexplored aspects of the effect particle shape on the porosity and angle of repose.
Early experimental studies on the porosity and angle of repose focused on simple geometries such as spheres and cubes [2,3].Since then, the field has expanded with many studies on various factors affecting the porosity or angle of repose, such as particle cohesion [4][5][6], container size [7], dynamics [8], particle-size distribution [9], water content [10], and even gravitational acceleration [6,11].The expansion of the field has also led to an increase in computational studies and a large public data set of angle of repose measurements [12].Computational studies in particular have become popular to investigate the effect of particle shape, including studies on rounded [13], simple polygonal [14][15][16], ellipsoidal [17], super-ellipsoidal [18], or arbitrarily defined smooth shapes [8].Common shape descriptors like angularity, sphericity, and elongation are often considered [19].In spite of the many studies on the dependence of the porosity and angle of repose on particle shape, these have been constrained by what materials are experimentally available or what particle shapes were permitted by the simulation method.In particularl, studies on well-defined concave or sharply angular shapes appear to be missing.This leaves a gap in the literature, because even small protrusions can cause large differences in the angle of repose [20].
Limitations in computational studies typically follow from the simulation method used.The vast majority of studies were performed using the discrete element method (DEM) that explicitly describes the behaviour of distinct particles.Despite the extensive use of DEM, the possibility to represent various particle shapes remained limited.The original DEM studies only used individual spheres, but the method has since been extended to rigid clusters of spheres [14], potential particles [21], and polyhedra [15,16,22,23].Nonetheless, the range of accessible shapes is still restricted, either due to inherent limitations of the chosen shape-representation method or simply by its computational cost.In particular, concave and highly angular particles remain challenging to simulate, often resulting in the necessity for ad hoc artificial rolling friction parameters and the neglect of particle interlocking effects.
A recently developed shape descriptor, using level sets, has shown great promise in being able to accurately describe any particle shape or geometry [24].Particles are described using a discrete level set or signed distance function of which the zero contour marks the particle's surface.This method was named level set DEM (LS-DEM) and can simulate virtually any shape or model by converting it into a level set.However, the original LS-DEM missed some essential physics for properly simulating arbitrarily shaped particles, which was later added by introducing volume interactions, giving rise to volume-interacting level set DEM (VLS-DEM) [25].
These recent developments opened up the opportunity to study systematically the effects of the asphericity, angularity, and concavity of particles on the porosity and angle of repose.The five Platonic solids provide a natural choice of particle shapes for such a study (Fig. 1).When ordered by the number of sides, the solids provide a set of shapes with a monotonically increasing sphericity (0.671, 0.806, 0.846, 0.910, 0.939) and monotonically decreasing angular defect at each corner (180 • , 120 • , 90 • , 60 • , 36 • ).Although there are many possible definitions for these geometrical properties, the Wadell sphericity and angular defect make intuitive sense when respectively describing the sphericity or asphericity and angularity or roundness of these shapes [26].The Platonic solids thus provide a logical variation of sphericities and angularities while also introducing a minimal amount of asymmetries that can otherwise affect results.
The Platonic solids themselves are convex and would not utilise the full potential of VLS-DEM.However, the solids can be made concave (or even more convex) by applying a transformation on each surface point ⃗ .By using the hyperbolic function ⃗  → ⃗ (⃗  ⋅ ⃗ )  with tunable constant  for this transformation, a continuous range of shapes can be generated that are unambiguously related to the original Platonic solid.More importantly, all of the shapes become a sphere for  = −0.5, providing a common reference point for all of the shapes.For positive , the solids become more concave and, because of a larger angular defect, more angular (top row Fig. 1).For negative  below −0.5, the solids also become more concave but, because of an angular excess, more round (bottom row Fig. 1).Hyperbolisation of the Platonic solids therefore provides an excellent range of shapes on a continuous spectrum that vary characteristically in their sphericity, angularity, and concavity.
This study thereby investigates the Platonic solids and variants thereof that have systematically been made more convex or concave than their original shapes.The solids are modified to various degrees by hyperbolising the shape surface with a continuous tunable constant  and, because the hyperbolisation procedure can be used both to deflate and inflate a shape, the variation of angularity is also considered.VLS-DEM is then used to create particle packings, to compute the porosity, and to deposit piles for computing the angles of repose.

Volume-interacting level-set discrete element method (VLS-DEM)
Each particle in VLS-DEM is defined by its own discrete level set or signed distance function [24].Instead of explicitly defining a surface, a function (⃗ ) is defined that gives the shortest distance  to the surface of the particle at any given position ⃗ .The distance  is negative if the point ⃗  is inside the particle and positive outside the particle.The surface of the particle thus follows from the zero contour, (⃗ ) = 0. To avoid constraining VLS-DEM to shapes for which an analytical expression exists, the function (⃗ ) is discretised on a regular grid of points The level set values {(⃗   )} are then saved, implicitly defining the particle's shape, and the value of  at any non-grid point can be obtained by trilinear interpolation of {(⃗   )}.
Contrary to LS-DEM, VLS-DEM does not use surface nodes to define particle-particle interactions but instead directly integrates the overlap volume between the particles.This avoids multiple problems, such as choosing the number and distribution of nodes across the particle surface, and allows the accurate simulation of concave particle geometries [25].An Octree integration algorithm (OIA) is used to determine the overlap volume   between two particles  and .First, the axisaligned bounding boxes (AABBs) of the particles are intersected to define a potential overlap region.The resulting potential overlap region is a large rectangular voxel, which is then recursively divided into eight equally-sized smaller voxels.By definition, each voxel  has a centre point ⃗   with a maximum radius  to the corners of the voxel.The refinement process is then repeated either until a maximum of five refinement steps has been reached or one of two criteria has been met:   (⃗   ) >  and   (⃗   ) > , meaning that the entire voxel is outside of both particles and it can be discarded, or   (⃗   ) <  and   (⃗   ) < , meaning that the entire voxel is inside of both particles and can be fully added to the overlap volume   .If the refinement finishes with neither of the two previous criteria having been met, the remaining voxels are added to   using a smearing rule defined below.
The overlap volume is thus determined according to: where the sum is over all undiscarded voxels ,   is the volume of voxel ,  max , and () is the smeared Heaviside step function given by: where  is the smearing width, which is set equal to the value of .
The contact point ⃗  of the interaction follows as: Each particle also has a contact area vector ⃗   defined by: with  , =   ( ⃗   ) and () being the smeared Dirac delta-like function: Each particle then has its own contact area and where ⃗   is the linear shear displacement,   is the volumetric normal stiffness,   is the linear tangential stiffness, and  the angle of friction.All simulations were implemented using YADE [27].A standard Verlet scheme is used to carry out the discrete time integration.Damping is effected by the default numerical damping scheme that directly adjusts the total force acting upon a particle (full details in YADE documentation).The damping constant was set to   = 0.15.A more extensive description of VLS-DEM can be found in van der Haven et al. [25].

Particle shapes: Hyperbolised platonic solids
The five Platonic solids and hyperbolisations thereof, created by applying the transformation given by Eq. ( 8), are used for the particle shapes (Fig. 1).First, the vertices, faces, and face normals of a Platonic solid are computed.For convenience, the vertices are chosen to lie on the unit sphere during the hyperbolisation procedure.The next step is to generate a regular grid that forms the basis of the discrete level set function.All points ⃗  in the regular grid are considered one by one to initialise the level set values ( ⃗ ).Only the distance to the closest face is needed to initialise a point ⃗ .This is the face for which the angle between the face normal f and point normal p = ⃗ ∕| ⃗ | is the smallest.Then, the ray that goes both through the origin and ⃗  is taken, giving a point on the face called ⃗ .The centre-to-face distance  = |⃗ | is then remapped using: where  is the distance from the origin to the centre of the face and  the hyperbolisation constant or hyperbolicity.The result is identical to a remapping of all surface points ⃗  with the transformation ⃗  → ⃗ (⃗  ⋅ ⃗ )  .The initial level set value is then given by ( ⃗ ) = | ⃗ | − .The aforementioned procedure gives the correct contour ( ⃗ ) = 0 but does not produce a correct signed distance function, because hyperbolisation also distorts the distance values.Therefore, the level set close to the interface was reconstructed using a probabilistic reinitialisation method (details to be published in the future as [28]).In short, the absolute distances after hyperbolisation may not be correct, but when sufficiently close to the interface, the values are still in the same order when sorted by magnitude.This information was used to force the distance values of cells adjacent to the interface to obey the probability density function of distance values that follows from a random plane separating two points on a regular grid.This can be justified because, at the resolution used in this study, even a sphere, having the smallest surface area of all 3D shapes, already has  ( 10 4 ) grid points adjacent to the interface.The rest of the level set is then reinitialised using a MATLAB implementation of the Fast Sweeping algorithm [29,30].
The hyperbolisation constants considered in this study are  ∈ (0.5, 0.0, −0.5, −1.0).Regardless of the solid chosen, all shapes become a perfect sphere for  = −0.5, immediately providing a common reference point across all shapes and  values.The final grid resolution of all shapes was around 50 × 50 × 50.All shapes were rescaled such that each of them had the same volume as a sphere with a radius of  ref = 0.003 m, resulting in maximum particle radii in the range  ∈ [0.003, 0.010] m.

Simulation set-up and analysis
A container was made using five infinite walls, four side walls and a floor.The side walls were spaced apart by 10, with  = 2 being twice the maximum particle radius of whichever shape is being simulated, giving a container base measuring 10 ×10.A cloud of 4096 particles was then generated by placing the particles with a random position and orientation in a space of 10 by 10×85 and subsequently deposited under gravity.All particles and walls have stiffness constants   = 5 × 10 7 N m −3 and   = 4 × 10 6 N m −1 .All particles have an angle of friction of  = arctan (0.3) rad whereas all walls have  = arctan (0.1) rad.The particles' density is  = 800 kg m −3 .The acceleration due to gravity is 9.81 m s −1 in the negative z direction.The time step is  = 10 × 10 −6 s, equal to 10 μs.
The deposition simulation is continued until the average kinetic energy of the particles drops below ⟨ kin ⟩ = 1×10 −8 J. Once the deposition has finished, the porosity is calculated using  = 1 − 3 3  ref ∕(4 ) with  being the number of particles and  the volume, excluding all volume and particles within 3 ref of the highest particle or any wall.Due to its large radii, despite having the same volume, the tetrahedron packing with  = 0.5 created a very shallow powder bed, so 8192 particles were simulated instead and the data exclusion margin was taken to be 2 ref instead of 3 ref .
After the deposition, one of the side walls was deleted and the simulation continued until the average kinetic energy dropped below ⟨ kin ⟩ = 1 × 10 −9 J.The packing thus collapses, and the resulting pile is equilibrated such that the angle of repose (AoR) can be measured.
To measure the AoR, the centre of mass of all particles is projected on a side wall.This projection gives a pointillistic 2D image of the pile.The surface of the pile can be identified as the interface between the point-dense region and the point-sparse region of the projection.In such cases, highly accurate estimations of the interface location can be obtained by optimising the log-likelihood function: which is based on a maximum-likelihood estimator that maximises the contrast between two regions given a linear interface [31].The interface is defined by the points ⃗  and ⃗  on the boundaries of the domain, where  is the total number of particles,  1 the number of particles on the sparsely populated side of the interface,  the total area of the domain, and  1 the area on the sparsely populated side of the interface.As before, all particles within 3  of any wall are excluded.The AoR is then computed as: However, the interface is not completely linear.The projected pile is therefore divided into ten equally wide vertical bands to which Eq. ( 9) is applied and the AoR is computed.The final AoR is then given by the average of the five largest angles.

Packing porosities
Several example particle packings are shown in Fig. 2, where it can be seen that a fixed number of particles of different shapes results in different packing heights and porosities.Table 1 and Fig. 3 show the porosities both as a function of the chosen shape (or number of sides) and the hyperbolicity.For the regular Platonic solids with  = 0, the solids in increasing order of porosity are: hexahedron < icosahedron < dodecahedron < octahedron < tetrahedron.This order is in agreement with earlier computational and experimental studies on loose packings of Platonic solids, giving confidence in the validity of the current method [15,16,23].Other studies observe a different order, such as Zhao et al.where the hexahedron (or cube) seems to be out of place [14].However, the hexahedron is perhaps the shape most sensitive to changes in the deposition method as, unlike the other Fig. 3. Porosities of the packing resulting from the gravity deposition of a cloud of randomly placed particles depending on the shape of the particles.Left: porosity per shape as a function of the hyperbolicity.Right: porosity per hyperbolicity value as a function of the number of faces (in parentheses).

Table 1
The packing porosities  for various particle shapes as a function of the hyperbolicity parameter  after the gravity deposition of a cloud of randomly placed particles.Note that  = −0.5 results in a sphere for all shapes.shapes considered, its minimum porosity is zero.Athanassiadis et al. also observed a different order but also reported high errors that make the actual order uncertain [32].An interesting follow up would be to shake the particles to see how the porosities change as this has been shown to strongly affect the porosities and can change the aforementioned order as well [25,33].The porosity values seem to be slightly higher in general, by about 0.05, than other reports [15,16,23].This difference might disappear if the VLS-DEM accuracy is increased to six integration refinements, as VLS-DEM has a tendency to overestimate porosities for a smaller number of refinements [25].Nonetheless, the porosities of the regular Platonic solids and their trends are in generally in good agreement with the literature, lending confidence to the view that the current methodology should also be sufficiently accurate to assess the porosities of the hyperbolised Platonic solids.

Hyperbolicity
Simply from geometry, one might expect that as the shape of the particles becomes more angular and concave that the porosity of the subsequent packing increases.The chosen Platonic solid determines the angularity, which can be exacerbated by a positive hyperbolicity.As a polygon grows in its number of faces, it will increasingly start to resemble a sphere, becoming less angular.The concavity is determined by the hyperbolicity .For −0.5 ≤  ≤ 0.0, the shapes are convex whereas shapes become concave for  < −0.5 and  > 0.0.Looking at the trends in Fig. 3, it seems that a positive hyperbolicity sharply increases the porosity.Moreover, the overall increase appears to be stronger for angular shapes and weaker for nearly spherical shapes such as the dodecahedron and icosahedron.This suggests that there is substantial interplay between the effects of angularity and concavity.The hexahedron (or cube) again seems to break the generally expected trend with respect to angularity or sphericity.However, it is not surprising that any hyperbolisation increases the porosity of the hexahedron packing as the minimum porosity of zero can no longer be reached.
On the other end, making the particles more concave but also more rounded using a negative hyperbolicity seems to have a much weaker effect on the porosity.The porosity of the tetrahedron drops, implying that for this particular shape its angularity is more important than its convexity.The hexahedron has been discussed already, and the octahedron actually becomes rather cube-like and obtains a lower porosity Fig. 1.For the most spherical shapes, the dodecahedron and icosahedron, it seems that the concavity effect overpowers the reduction in angularity, causing an increase in the porosity.
To conclude, angularity and concavity increase the porosity and their effects can be combined to create a disproportionally strong increase in porosity.Such knowledge may be used to tune the porosity of future man-made or artificial materials.However, more investigations into the interplay of angularity and concavity are needed, especially with regard to the cohesion and friction as these are other main determinants of the porosity [34].

Angles of repose
Examples of the piles used to compute the angle of repose can be seen in Fig. 4. Table 2 and Fig. 5 show the angle of repose both as a function of the chosen shape (or number of sides) and the hyperbolicity.Again, first looking at the regular Platonic solids with  = 0, the solids in increasing order of the angle of repose are: icosahedron < tetrahedron < dodecahedron < hexahedron < octahedron with a range of 5.1 • between the lowest and highest angles of repose.
Studies on the angle of repose of Platonic solids are much rarer than studies on the porosity, making a systematic comparison with literature results more difficult, but two suitable studies were found.Landauer et al. reported an order of hexahedron < octahedron < tetrahedron [16].The tetrahedron falls out of place but the order of the other two is preserved.The reported range of angles spanned more than 10 degrees under static conditions but reduced to about 6 degrees in a slowly (0.15 Hz) rotating cylinder.Zhou et al. reported icosahedron < dodecahedron < tetrahedron < octahedron < hexahedron with a range spanning just under 6 degrees [14].Compared to our results, the dodecahedron and tetrahedron as well as the octahedron and hexahedron seem to have swapped places.Moreover, Landuaer and Zhou contradict each other, reporting completely reversed trends.It thus remains difficult to draw any definite conclusions as the differences in friction coefficients, depositions methods, and boundary conditions are obfuscating the underlying trend.Nonetheless, by maintaining identical

Table 2
Angle of repose in degrees ( • ) of the pile following the collapse of a particle packing depending on the shape of the particles.Note that  = −0.5 results in a sphere for all shapes.conditions amongst all simulations within the current study, it should still be possible to use the regular Platonic solids as a standard for comparison with their hyperbolised versions.It appears that the angle of repose is generally minimal for the convex hyperbolicities, −0.5 ≤  ≤ 0.0, and sharply increases outside this range.A possible explanation may be that hyperbolisation increases the effective rolling friction of the particles when moving further away from the convex range, thereby increasing the angle of repose, whilst also creating the potential for interlocking of the concave particles.Regardless of the explanation, this implies that even convex polyhedra might need a rolling friction parameter to accurately capture the angle of repose of otherwise concave shapes.It is also noticeable that, contrary to the porosities, the increases in the angles of repose are similar for  = −1.0 and  = 0.5.Seeing that these trends differ between the porosity and angle of repose further indicates that, even with a rolling friction parameter, convex polyhedra may not be capable of capturing both material behaviours at once.

Conclusions
The porosity or solid fraction of granular materials has previously been studied as a function of friction, cohesion, deposition method, and particle geometry.However, a large part of the parameter space has always been unreachable due to the inability to simulate concave particle geometries.By using the volume-interacting level discrete element method (VLS-DEM), this study provides the first systematic investigation into the angle of repose of concave particles, showing that not only the angularity or sphericity but also the concavity of particles strongly affects the angle of repose.Angularity and concavity further seem to increase each other's influence on the porosity and angle of repose, generally leading to a larger than expected increase in both properties.Finally, the trends in the porosity and angle of repose with respect to angularity and concavity are not the same, indicating that even convex polyhedra with a rolling friction parameter may not be able to capture all the effects of particle geometry.This highlights the need for simulation techniques such as VLS-DEM that are capable of representing all particle geometries when considering realistic problems.

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.

Fig. 1 .
Fig. 1.The five Platonic solids and the hyperbolised variants thereof considered in this study.The columns from left to right are: tetrahedron, hexahedron (cube), octahedron, dodecahedron, and icosahedron.The original Platonic solids in the middle row are convex, whereas the hyperbolised variants in the top and bottom rows are concave.

Fig. 5 .
Fig.5.Angle of repose of the pile resulting from the collapse of a particle packing depending on the shape of the particles.Left: angle of repose per shape as a function of the hyperbolicity, right: angle of repose per hyperbolicity value as a function of the number of faces (in parentheses).