Minkowski-Voronoi diagrams as a method to generate random packings of spheropolygons for the simulation of soils

Minkowski operators (dilation and erosion of sets in vector spaces) have been extensively used in computer graphics, image processing to analyze the structure of materials, and more recently in molecular dynamics. Here, we apply those mathematical concepts to extend the discrete element method to simulate granular materials with complex-shaped particles. The Voronoi-Minkowski diagrams are introduced to generate random packings of complex-shaped particles with tunable particle roundness. Contact forces and potentials are calculated in terms of distances instead of overlaps. By using the Verlet method to detect neighborhood, we achieve CPU times that grow linearly with the body’s number of sides. Simulations of dissipative granular materials under shear demonstrate that the method maintains conservation of energy in accord with the first law of thermodynamics. A series of simulations for biaxial test, shear band formation, hysteretic behavior, and ratcheting show that the model can reproduce the main features of real granular-soil behavior.

Figure 1

Dilation of a square by a disk.

Figure 2

Erosion of a square by a disk.

Figure 3

Opening of a square by a disk.

Figure 4

A typical Voronoi construction in two dimensions. Exactly one Voronoi site (dot) is randomly generated within each square of a plane grid. Each Voronoi polygon is the set of points closer to a given site than to any other site.

Figure 5

For each vertex $\vec{x}$ of the Voronoi cell (see Fig. 4), one finds a vertex candidate ${\vec{x}}_e$ at a distance $d$ of the cell edges intersecting at $\vec{x}$. The vectors ${\vec{E}}_1$ and ${\vec{E}}_2$ are two adjacent edges of the Voronoi polygon.

Figure 6

The initial Voronoi cell $A$ (left) is eroded by a disk element $B_d$ of radius $d$. The eroded polygon $A\ominus B_d$ (center) is then dilated by the disk element $B_R$, producing the desired spheropolygon $\left(A\ominus B_d\right)\oplus B_R$ (right). Note the pathological case of points 1 and 2; since they are closer than $d$ to another cell side, they are substituted by point 3.

Figure 7

An array of Minkowski-Voronoi spheropolygons constructed by the method described in the text. To allow void spaces between particles, the disk used for dilation is slightly smaller than the one used for erosion.

Figure 8

Each spheropolygon is divided into regions to compute its mass properties.

Figure 9

Snapshot of shear-cell simulation with regular spheropentagons (above) and Minkowski-Voronoi spheropolygons (below).

Figure 10

Relative error in the energy balance versus time step (in units of seconds) in simulation of shear cells with Minkowski-Voronoi spheropolygons.

Figure 11

(Color online) Inverse of the Cundall number (in units of seconds) versus number of vertices of the regular spheropolygons employed for the simulation of the shear band (log-log scale). When the Verlet lists are almost switched off $\left(\alpha =0.9\right)$, one observes a power-law dependence with exponent 1.56(23). Switching on the optimization scheme by taking $\alpha =0.3$ reduces this exponent to a linear relationship [exponent 0.98(31)].

Figure 12

(Color online) Normalized CPU time versus Verlet distance $\alpha$ (in centimeters) for the simulation of a shear cells with Minkowski-Voronoi spheropolygons. The model of Eq. (23) fits this dependence quite well, with parameters $\beta =0.704$, $\gamma =0.3573$, $A=4.34\times {10}^{-7}$, and $B=9.91\times {10}^{-5}$.

Figure 13

(a) Deviatoric stress $q$ and (b) volumetric strain ${ϵ}_v$ as functions of the axial strain ${ϵ}_a$ for both a dense (continuous line) and a loose (dashed line) packing of 1200 spheropolygons under a biaxial test (horizontal confining pressure ${\sigma }_a=0.256K_n$ and vertical constant strain rate ${\dot{ϵ}}_l=0.005/\text{s}$).

Figure 14

Deviatoric stress at failure as a function of the confining pressure. A nonlinear power dependence, with exponent 0.85(2), is obtained, in contrast with the linear relationship established by the Mohr-Coulomb failure criterion.

Figure 15

(Color online) (a) (Color online) The shear bands are identified by plotting the cumulated displacement vector for each spheropolygon. (b) The measured shear band orientations (diamonds) for several confining pressures, compared with those predicted by the Mohr-Coulomb (circles) and the Roscoe (triangles) criteria.

Figure 16

(Color online) Mohr-Coulomb failure line for Voronoi spheropolygon packings with three different erosion parameters: $r=0.5$, 1, and 1.5 mm. Each failure line is the envelope of the Mohr circles at failure (curves in the figure). As erosion is larger the strength of the packing decreases.

Figure 17

(Color online) (a) The hysteresis loops for different oscillation amplitudes in the stress-strain space. The stress difference is given in units of $p_0$. (b) Hysteresis area versus amplitude.

Figure 18

(Color online) Damping coefficient $D$ versus tensile modulus $E={\sigma }_m/{ϵ}_m$.

Figure 19

(a) The relation between the axial stress and axial strain and (b) the accumulated strain at each cycle showing granular ratcheting.

Figure 20

Accumulated strain per cycle versus the oscillation amplitude.

Figure 21

The origin of numerical ratcheting in a spheropolygon. Left: a closed ABCD loop in the trajectory of the contact leads to a spurious change in potential energy, as shown by McNamara et al. [32]. The reason is that the energy stored in AB is not the same as the energy released in CD. This spurious energy is absent in vertex-edge contacts (right), because AB and CD have the same length. Therefore, the spheropolygons only partially solve the problem of numerical ratcheting.

Figure 22

(Color online) Permanent accumulation of plastic deformation versus the number of cycles for three different particle roundness. The roundness is defined by the erosion parameters, which are taken as $r=0.5$, 1, and 1.5 mm.