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 (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 and f is an n dimensional array with componentsIt is clear from (1), (2) that the polyhedron, and hence its volume V, is a function of R. Moreover, the vectors r and their magnitudes ∥r∥ are routinely determined during each time step of a DEM simulation.

Fig. 1

Specification of the Voronoi polyhedron (or polygon in two dimensions) in terms of the vectors r, which in this case are obtained from an analysis of particle packing.

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 V0 is the volume at R0,

is the gradient at R0, ΔR = R − R0 is the change in R relative to R0 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 aswhere A 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 denote the area of face k when projected onto the plane x = 0, as illustrated in Fig. 2a, we obtain

Fig. 2

Projection of (a) face onto the plane x = 0 (to yield ) and (b) edge onto the plane x = 0 (to yield ).

As in Fig. 2a, the unprojected and projected face k will henceforth be denoted by and , respectively. For Eq. (6) to be useful, we must demand that |r| > 0. As is commonly done [5], we select γ so that |r| is maximised (although γ depends on k, we will write γ rather than γ for notational simplicity). Combining Eqs. (2), (5), (6), the polyhedron volume can be expressed as

The next task is to determine the area .

Eliminating x from the inequalities (1) using the equation r · x = f, one finds that is defined by a systems of inequalities of the formHere, S is an n × 2 matrix, g is an n dimensional array and y is a vector containing the remaining two spatial coordinates. Row j of S may be considered as a two-dimensional vector, denoted by s, analogous to the three-dimensional vector r. If we let s be component (j, m) of S and g be component j of g, we specifically obtain

where p = r/r. In Eq. (9), M = m when m < γ and M = m + 1 when m ≥ γ, to accommodate for the projection. Clearly, s and g need only be determined for j ≠ k, since the inequality on row j in (8) will be trivially satisfied.

The special case when ∥s∥ vanishes deserves some special attention. From the two equalities embodied in Eq. (9) together with the definition of p, it is seen that ∥s∥ can only vanish if r and r are linearly dependent (parallel if p > 0 and antiparallel if p < 0). When r and r 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 (the constraint imposed by r is redundant when g > 0 and the one imposed by r when g < 0). Hence, r and r are parallel and r is redundant when p > 0 and g < 0, in which case 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 denotes the index set for the edges on is the length of edge and

is the perpendicular distance from the origin to the edge. The edge length is calculated via Lasserre's projection scheme. If we let denote the length of edge when projected onto the line x = y = 0, as illustrated in Fig. 2b, we obtain

As in Fig. 2b, the once-projected and twice-projected edge will henceforth be denoted by and , respectively. We select b so that |s| is maximised (although b make take on different values for different faces , we will write b rather than b for notational simplicity). Combining Eqs. (11), (12), (13) we find that

Eliminating y from (8) using the equation s · y = g, one finds that is defined by a system of inequalities of the form

where t and h both are n dimensional arrays and z is a scalar. If we let t and h be component i of t and h, respectively, we specifically obtain

where q = s/s and a = 3 − b indicates the remaining component of y. Clearly, t and h need only be determined when i, j and k are all unequal.

Special considerations are needed when t vanishes. From Eq. (16) and the definition of q, it is seen that t can only vanish if s and s are linearly dependent (parallel if q > 0 and antiparallel if q < 0). When s and s 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 (the constraint imposed by s is redundant when h > 0 and the one imposed by s when h < 0). Hence, s and s are parallel and s is redundant when q > 0 and h < 0, in which case 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/t for which t > 0, denoted by u. Analogously, the greatest lower bound equals the largest value of h/t for which t < 0, denoted by u. When and the edge length becomes . Otherwise, and this combination of j and k needs not be considered further. Moreover, if there is no t > 0 or no t < 0, the polyhedron is unbounded.

Gradient calculation

Using the chain rule, one finds thatwhere 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

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 {α, β, γ} be a permutation of {1, 2, 3}, such that the first projection is along γ (the value of γ depends on k) and the second projection is along β (the value of β depends on j and k).

When l ≠ γ, x indeed is a homogeneous function of degree 1 on . 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 . The case l ≠ γ can be further subdivided into l = α and l = β. When l = α,

where . When l = β, the second projection onto the line x = x = 0 is accomplished via the equation sx + sx = g. Consequently

When l = γ, the first projection onto the plane x = x = 0 is accomplished via the equation rx + rx + rx = f. Consequently

where Theorem 2.4 in [4] has been used (note that rx + rx is a homogeneous function of degree 1 on ). A variable substitution has been accomplished via Eq. (13), again noting that , and we have let

where d = rs − rs.

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 (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 d0 (corresponding to the radius r0 = d0/2) and that the maximal particle–particle overlap was δ, corresponding to a minimal dihedral angle between the vectors r of

The value of δmax/d0 was kept fixed at 0.2, corresponding to a minimal dihedral angle φmax of about 47°. First, n random unit vectors were generated, such that the dihedral angle between any pair was at least φmin. Specifically, the components of n were drawn from independent uniform distributions on the interval [−1,1], and was calculated as n/∥ n ∥ whenever ∥n ∥ >0. Thereafter, r was calculated as , where the magnitude r was drawn from a uniform distribution on the interval [r0 − δ/2, r0]. Keeping the intended application in mind, only polyhedra with volumes not exceeding (3/2)V were included in the subsequent analysis, where is the particle volume. Components of ΔR were drawn from a uniform distribution on the interval [− ρ, + ρ], with ρ/r0 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 ρ, 1000 different variations ΔR were generated, and for each variation, both the exact (V) and approximate (Vapprox) Voronoi volumes were determined. As error estimate, the relative error of the Voronoi volume was calculated as

and the average and maximal values of |ɛ| 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/cm3 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 μs, implying that about 3.8 × 106 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 ΔR exceeded a certain predefined threshold, viz.

where ϵ was selected as 0.1 and 1% (as before, r0 denotes the particle radius). Whenever updates were promted by changes in ΔR, the exact and approximate Voronoi volumes were stored. Based on these, the relative error of the volume change was calculated as

where ΔV = V − V0 is the exact volume change and ΔVapprox = Vapprox − V0 is the approximate volume change since the last update that yielded the volume V0. Data were binned based on the magnitude of ΔV and the average and maximal values of |ɛΔ| were determined for each bin.

Numerical results and discussion

The results obtained from the random polyhedra are summarised in Fig. 3, Fig. 4. The average and maximal relative errors of the Voronoi volumes, determined for five different values of ρ/r0 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/n3/2, where the proportionality constant C depended on ρ/r0 (solid lines in Fig. 3a). Since the approximate Eq. (3) is first-order accurate, one expects the truncation error to be proportional to (ρ/r0)2, and the constant C was indeed found to be proportional to (ρ/r0)2 (inset in Fig. 3a). It is clear that an accurate approximation was obtained as long as ρ/r0 remained small. Specifically, a maximal error of about 0.2% was observed when ρ/r0 did not exceed 1% and a maximal error of about 1% was obtained when ρ/r0 was at most 2%.

Fig. 3

Average (a) and maximal relative errors (b) of the Voronoi volumes for the indicated values of ρ/r0. The solid lines in (a) represent fits of the function C/n3/2 to the numerical data, and the inset displays the proportionality constant C vs. (ρ/r0)2.

Fig. 4

Computation times, relative to those for exact calculations, for approximate and full (i.e., exact plus gradient) calculations for the indicated values of ρ/r0 (the solid lines represent power-law fits to the numerical data).

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 ρ/r0. 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 |ΔV/V0|. 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 |ΔV/V0| decreases and may, for sufficiently small values of |ΔV/V0|, 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

Fig. 5

Average and maximal relative errors of the volume changes as a function of the relative volume change for (a) ϵ = 0.1% and (b) ϵ = 1%. The insets in (a) and (b) show the initial and final particle arrangements.

where and represent ‘flattened’ versions of ΔR and G0, i.e., vectors with 3n components rather than n × 3 matrices, 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 . As demonstrated by the numerical simulations, this condition is generally fulfilled provided that the magnitude of all components of are sufficiently small. Exceptions occur when the positive and negative contributions to balance each other, in which case ΔV/V0 becomes small. According to our experience, the proposed method nonetheless appears to work well in practice with ϵ ∼ 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 [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].

Table 1

Number of gradient updates and computational times for the small-scale DEM simulations. For comparison, the total number of iterations was 3,773,748 and the computational time was 8192 s when exact volume determination was used.

Update limit ϵ	No. of updates	Computational time (s)
0.1%	148,238	2676
1.0%	45,271	2671

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.