Efficient Voronoi volume estimation for DEM simulations of granular materials under confined conditions

When the discrete element method (DEM) is used to simulate confined compression of granular materials, the need arises to estimate the void space surrounding each particle with Voronoi polyhedra. This entails recurring Voronoi tessellation with small changes in the geometry, resulting in a considerable computational overhead. To overcome this limitation, we propose a method with the following features:•A local determination of the polyhedron volume is used, which considerably simplifies implementation of the method.•A linear approximation of the polyhedron volume is utilised, with intermittent exact volume calculations when needed.•The method allows highly accurate volume estimates to be obtained at a considerably reduced computational cost.


Method details
We describe a method to approximate the volume of Voronoi polyhedra in situations where recurring updates are needed for small changes in the geometry, as during simulations of confined compression of granular materials with the discrete element method (DEM) [1] (see section 'Additional information: background'). The method utilises a linear approximation of the volume, with intermittent exact volume calculations when needed. The Voronoi polyhedron is specified in terms of n vectors r k (with k =1, . . ., n), see Fig. 1. Hence, the Voronoi polyhedron is represented by a system of inequalities, where x denotes the spatial coordinates, R is an n Â 3 matrix whose kth row equals r k and f is an n dimensional array with components It is clear from (1) and (2) that the polyhedron, and hence its volume V, is a function of R. Moreover, the vectors r k and their magnitudes kr k k are routinely determined during each time step of a DEM simulation.
The polyhedron volume is known to be continuously differentiable with respect to all parameters occurring in the defining inequality system (1) [2]. Hence the volume can be approximated by the linear form where V 0 is the volume at R 0 , is the gradient at R 0 , DR = R À R 0 is the change in R relative to R 0 and the colon indicates double contraction. The approximation (3) is expected to be valid provided that the change in R is small. The polyhedron volume is calculated according to the recursive projection scheme proposed by Lasserre [3], which is convenient to use when the polyhedron has a halfspace representation as in (1). The details are provided in section 'Volume calculation' in order to define the relevant quantities for the subsequent developments. The gradient determination, described in section 'Gradient calculation', utilises formulae derived by Mü ller et al. [2] and Lasserre [4].

Volume calculation
The convex polyhedron defined by the inequality system (1) may be decomposed into n pyramids. Hence its volume can be expressed as where A k is the area of face k. Expression (5) is not useful unless a way be found to calculate the face areas. To this end, Lasserre [3] suggested a straightforward projection scheme. If we let A 0 k denote the area of face k when projected onto the plane x g = 0, as illustrated in Fig. 2a, we obtain As in Fig. 2a, the unprojected and projected face k will henceforth be denoted by F k and F 0 k , respectively. For Eq. (6) to be useful, we must demand that jr kg j > 0. As is commonly done [5], we select g so that jr kg j is maximised (although g depends on k, we will write g rather than g k for notational simplicity). Combining Eqs. (2), (5) and (6), the polyhedron volume can be expressed as The next task is to determine the area A 0 k . Eliminating x g from the inequalities (1) using the equation r k Á x = f k , one finds that F 0 k is defined by a systems of inequalities of the form S k y g k : Here, S k is an n Â 2 matrix, g k is an n dimensional array and y is a vector containing the remaining two spatial coordinates. Row j of S k may be considered as a two-dimensional vector, denoted by s jk , analogous to the three-dimensional vector r k . If we let s jmk be component (j, m) of S k and g jk be component j of g k , we specifically obtain s jmk ¼ r jM À p jk r kM ; where p jk = r jg /r kg . In Eq. (9), M = m when m < g and M = m + 1 when m ! g, to accommodate for the projection. Clearly, s jmk and g jk need only be determined for j 6 ¼ k, since the inequality on row j in (8) will be trivially satisfied. The special case when ks jk k vanishes deserves some special attention. From the two equalities embodied in Eq. (9) together with the definition of p jk , it is seen that ks jk k can only vanish if r j and r k are linearly dependent (parallel if p jk > 0 and antiparallel if p jk < 0). When r j and r k are parallel, the 'outer' constraint is redundant and shall not be included in the subsequent analysis. This constraint is in turn identified by the sign of g jk (the constraint imposed by r j is redundant when g jk > 0 and the one imposed by r k when g jk < 0). Hence, r j and r k are parallel and r k is redundant when p jk > 0 and g jk < 0, in which case F k shall not be considered further.
The convex polygon defined by the inequality system (8) may be decomposed into triangles. Hence its area may be expressed as where S k denotes the index set for the edges on F 0 k ; L 0 jk is the length of edge j 2 S k and is the perpendicular distance from the origin to the edge. The edge length is calculated via Lasserre's projection scheme. If we let L 00 jk denote the length of edge j 2 S k when projected onto the line x b = y b = 0, as illustrated in Fig. 2b, we obtain [ ( F i g . _ 2 ) T D $ F I G ] Fig. 2. Projection of (a) face F k onto the plane x g = 0 (to yield F 0 k ) and (b) edge E 0 jk onto the plane x b = 0 (to yield E 00 jk ).
As in Fig. 2b, the once-projected and twice-projected edge E jk will henceforth be denoted by E 0 jk and E 00 jk , respectively. We select b so that js jbk j is maximised (although b make take on different values for different faces j 2 S k , we will write b rather than b jk for notational simplicity). Combining Eqs. (11)-(13) we find that Eliminating y b from (8) using the equation s jk Á y = g jk , one finds that E 00 jk is defined by a system of inequalities of the form where t jk and h jk both are n dimensional arrays and z is a scalar. If we let t ijk and h ijk be component i of t jk and h jk , respectively, we specifically obtain where q ijk = s ibk /s jbk and a = 3 À b indicates the remaining component of y. Clearly, t ijk and h ijk need only be determined when i, j and k are all unequal. Special considerations are needed when t ijk vanishes. From Eq. (16) and the definition of q ijk , it is seen that t ijk can only vanish if s ik and s jk are linearly dependent (parallel if q ijk > 0 and antiparallel if q ijk < 0). When s ik and s jk are parallel, the 'outer' constraint is redundant and shall not be included in the subsequent analysis. This constraint is in turn identified by the sign of h ijk (the constraint imposed by s ik is redundant when h ijk > 0 and the one imposed by s jk when h ijk < 0). Hence, s ik and s jk are parallel and s jk is redundant when q ijk > 0 and h ijk < 0, in which case E 00 jk shall not be considered further.
It is clear from the inequality system (15) that the least upper bound on z equals the smallest value of h ijk /t ijk for which t ijk > 0, denoted by u mjk . Analogously, the greatest lower bound equals the largest value of h ijk /t ijk for which t ijk < 0, denoted by u njk . When u mjk > u njk ; j 2 S k and the edge length becomes L 00 jk ¼ u mjk À u njk . Otherwise, j = 2 S k and this combination of j and k needs not be considered further. Moreover, if there is no t ijk > 0 or no t ijk < 0, the polyhedron is unbounded.

Gradient calculation
Using the chain rule, one finds that where the first derivative is to be taken with the right-hand-side vector f in (1) fixed and the second derivative with the left-hand-side matrix R fixed. Using Theorem 1 in [2], the derivatives in Eq. (18) can be expressed as @V @r kl where a variable substitution has been accomplished via Eq. (6) and where we have let As demonstrated by Lasserre [4], integrals on faces of convex polyhedra can be expressed in terms of sums of integrals on the face edges, provided that the integrand is a homogeneous function on the face. To be able to use this result, we need to distinguish between different cases. We let {a, b, g} be a permutation of {1, 2, 3}, such that the first projection is along g (the value of g depends on k) and the second projection is along b (the value of b depends on j and k).
When l 6 ¼ g, x l indeed is a homogeneous function of degree 1 on F 0 k . Using Theorem 2.4 in [4], we thus obtain where a variable substitution has been accomplished via Eq. (13) and where we have let noting that dL 00 jk ¼ dx a . The case l 6 ¼ g can be further subdivided into l = a and l = b. When l = a, where M 00 When l = g, the first projection onto the plane x l = x g = 0 is accomplished via the equation where Theorem 2.4 in [4] has been used (note that r ka x a + r kb x b is a homogeneous function of degree 1 on F 0 k ). A variable substitution has been accomplished via Eq. (13), again noting that dL 00 jk ¼ dx a , and we have let where d jk = r ka s jbk À r kb s jak .

Numerical simulations
Numerical simulations were performed to test the accuracy and computational efficiency of the described method. To assess the overall characteristics of the method, a set of random polyhedra was investigated. To test the performance of the method in its intended application, a small-scale DEM simulation was performed.

Random polyhedra
Each random polyhedron was specified in terms of n vectors r k (cf. 'Method details' section). The geometry was selected so that it would be relevant for volume estimation in conjunction with simulations of particulate systems with the DEM. It was assumed that all particles had the same diameter d 0 (corresponding to the radius r 0 = d 0 /2) and that the maximal particle-particle overlap was d max , corresponding to a minimal dihedral angle between the vectors r k of ' min ¼ arccos 1 À The value of d max /d 0 was kept fixed at 0.2, corresponding to a minimal dihedral angle w max of about 478. First, n random unit vectorsn k were generated, such that the dihedral angle between any pair was at least w min . Specifically, the components of n k were drawn from independent uniform distributions on the interval [À1,1], andn k was calculated as n k /k n k k whenever kn k k >0. Thereafter, r k was calculated as r knk , where the magnitude r k was drawn from a uniform distribution on the interval [r 0 À d max /2, r 0 ]. Keeping the intended application in mind, only polyhedra with volumes not exceeding (3/2)V sph were included in the subsequent analysis, where V sph ¼ pd 3 0 =6 is the particle volume. Components of DR were drawn from a uniform distribution on the interval [À r, +r], with r/r 0 ranging from 1 to 5%.
For each value of n between 6 and 14, 1000 random polyhedra were generated as described above.
For each random polyhedron and value of r, 1000 different variations DR were generated, and for each variation, both the exact (V) and approximate (V approx ) Voronoi volumes were determined. As error estimate, the relative error of the Voronoi volume was calculated as and the average and maximal values of je V j were determined. The computation times needed for full (i.e., exact plus gradient), exact and approximate volume determinations were recorded. A special routine implemented according to the description in section 'Volume calculation' were used for the exact volume determinations.

DEM simulation
Uniaxial compression of 100 spherical particles (diameter 1.0 mm) in a cylindrical die (diameter 5.0 mm) was simulated with the DEM [1]. The initial particle bed had a height of about 5.25 mm and the upper punch was lowered at a rate of 5 mm/s for a period of 0.52 s until a nonporous compact was formed. For illustrative purposes, a standard contact model of the linear spring-dashpot type was used as described in [6]. The particle density was 1.45 g/cm 3 and the normal and tangential stiffness were 100 N/mm. The sliding and rolling friction coefficients were 0.5 and 0.001 for contact between particles whereas five times smaller values were used for contact between particles and confining surfaces. The normal and tangential damping coefficients were chosen so that the fractional damping was 0.3. The time step was 0.138 ms, implying that about 3.8 Â 10 6 time steps were needed to complete the simulation.
During the course of the simulation, Voronoi volumes were determined as described in section 'Method details'. Exact volume determinations and gradient updates were made when contacts were formed or broken and when the magnitude of any component of DR exceeded a certain predefined threshold, viz.
where e was selected as 0.1 and 1% (as before, r 0 denotes the particle radius). Whenever updates were promted by changes in DR, the exact and approximate Voronoi volumes were stored. Based on these, the relative error of the volume change was calculated as where DV= V À V 0 is the exact volume change and DV approx = V approx À V 0 is the approximate volume change since the last update that yielded the volume V 0 . Data were binned based on the magnitude of DV and the average and maximal values of je DV j were determined for each bin.

Numerical results and discussion
The results obtained from the random polyhedra are summarised in Figs. 3 and 4. The average and maximal relative errors of the Voronoi volumes, determined for five different values of r/ r 0 between 1 and 5%, are displayed as a function of the number of faces, n, in Fig. 3a and b. In general terms, both the average and maximal errors decreased with increasing n and the maximal error is seen to be about one order of magnitude larger than the average error. It was empirically found that the average error could be described by a function of the form C/n 3/2 , where the proportionality constant C depended on r/r 0 (solid lines in Fig. 3a). Since the approximate Eq. (3) is first-order accurate, one expects the truncation error to be proportional to (r/r 0 ) 2 , and the constant C was indeed found to be proportional to (r/r 0 ) 2 (inset in Fig. 3a). It is clear that an accurate approximation was obtained as long as r/r 0 remained small. Specifically, a maximal error of about 0.2% was observed when r/r 0 did not exceed 1% and a maximal error of about 1% was obtained when r/r 0 was at most 2%.
[ ( F i g . _ 3 ) T D $ F I G ] Computation times for approximate and full (i.e., exact plus gradient) calculations are displayed as a function of n in Fig. 4. To enable a clearer view, all values have been normalised by the time required for the corresponding exact volume calculation (depending on n, this value ranged from a few seconds to about 1 min on the particular hardware used). As seen, consistent computation times were observed for all five values of r/r 0 . It is evident that the proposed approximation considerably reduced computation times, especially for large values of n, where a speedup around 450 was observed. The computational overhead introduced by gradient calculation ranged from about 50% for n = 6 to about 10% for n = 14.
The results from the small-scale DEM simulations are summarised in Fig. 5, which displays the average and maximal relative errors of the volume changes as a function of the magnitude of the volume change, expressed as jDV/V 0 j. The insets in Fig. 5a and b show the initial and final particle arrangements. As can be clearly seen, the magnitude of the relative errors of the volume changes (Fig. 5) are considerably larger than the magnitude of the relative errors of the Voronoi volumes (Fig. 3). Moreover, the relative errors increase when jDV/V 0 j decreases and may, for sufficiently small values of jDV/V 0 j, exceed unity. This behaviour is not unexpected but is rather a consequence of the nature of the approximation. Using matrix notation and expanding the volume to second order, one obtains where DR and G 0 represent 'flattened' versions of DR and G 0 , i.e., vectors with 3n components rather than n Â 3 matrices, H 0 is the corresponding Hessian matrix and the superscript 'T' indicates the matrix transpose. Hence, the relative error of the volume change can be approximated as The truncation error will be of the first order but will nevertheless be small provided that jG generally fulfilled provided that the magnitude of all components of DR are sufficiently small.
Exceptions occur when the positive and negative contributions to G T 0 DR balance each other, in which case DV/V 0 becomes small. According to our experience, the proposed method nonetheless appears to work well in practice with e $ 0.1%, likely because relative changes of the order of 10 À4 in the Voronoi volume are insignificant.
When the described method is used in DEM simulations, the overall savings in computation times will depend on how frequently the gradient needs to be updated on average. However, given the relatively modest overhead introduced by gradient calculation, significant speedups will result also for relatively frequent updates. This is substantiated by Table 1, which shows the total number of gradient updates and computational times for the small-scale DEM simulations. It might be possible to optimise the procedure used for the exact volume determination even further [ ( F i g . _ 5 ) T D $ F I G ]   [7]. However, judging from the data presented in Table 1, such an optimisation would have a modest impact on the overall computational times. Moreover, the proposed local procedure circumvents the costly determination of Voronoi diagrams (or the dual Delaunay tetrahedralization) inherent in alternative methods [8].

Additional information: background
The discrete element method (DEM), developed by Cundall and Strack [1], has established itself as the de facto standard technique for micromechanical simulations of granular systems in diverse application areas (e.g. [9,10]). The method is fairly efficient, at least in relative terms, owing to the generally employed simplified contact models. Assuming mechanical independence of contacts, interparticle forces are defined in terms of certain functions of the particle overlap (e.g. [11,12]). Although the majority of applications of the DEM have been related to dynamical processes with limited particle deformation, the method is increasingly used to study the behaviour of granular materials under confined conditions, as during manufacturing of tablets or machine parts by compression/compaction (see [9] and references therein).
Compression/compaction can often be considered as a macroscopically quasistatic process [13], characterised by extensive particle deformation. As a result, the assumption of contact independence will not be valid at high relative densities (exceeding about 0.85-0.90 for monodisperse spherical particles [14,15]). The combined finite/discrete element method (FE/DE method) [16,17] sometimes also referred to as the multiparticle finite element method (MPFEM) [18] or the meshed discrete element method (MDEM) [19] -has been proposed to overcome this limitation. When the FE/DE method is used, each particle is meshed by finite elements that enable a superior representation of particle deformation, but also result in a significantly higher computational cost. The FE/DE method is highly valuable in the detailed study of systems comprising a few particles, but is unpractical for large-scale simulations due to its prohibitive computational cost. Hence simplified models for the interaction between particles under confined conditions, suitable for implementation in the DEM, are needed.
The most important ingredient in such models appears to be the constraint imposed by plastic incompressibility [17,19,20]. Plastic deformation can only proceed as long as there is a void space that can accommodate the displaced material. This in turn necessitates that a way be found to estimate the volume of the available void space. The most natural and promising approach seems to be to use volume estimates based on Voronoi constructions, as originally proposed by Arzt [14] and subsequently used in the DEM by Donzé et al. [19]. Since the methodology is well developed and open-source software exists, such as CGAL [21] and Voroþþ [22], the main challenge is to retain the computational efficiency of the DEM despite the overhead introduced by Voronoi volume determination. Voronoi volumes can readily be determined from half-space representations using the mentioned codes, but other representations are used internally.