A Voronoi strain-based method for granular materials and continua

In a recent article (Frenning in Comp Part Mech 24:1–4, 2021), we demonstrated that a Delaunay-based strain estimate could be used as a starting point for the development of a particle-based method for continua. In this article, we argue that the Voronoi diagram, dual to the previously used Delaunay tetrahedralization, provides a more natural description of the underlying particulate system. For this reason, a Voronoi-based estimate of the deformation gradient is derived and used to the same effect. Although the gradient vectors cease to be antisymmetric, sums over nearest neighbors vanish, which results in a formulation that not only is linearly complete but also satisfies the patch test irrespective of initial particle placement. Pairwise forces, inferred from the local (nonaffine) deformation of each bond or contact, impart a physical stabilization. Forces are obtained from a discrete Lagrangian, thus ensuring that linear and angular momenta are conserved in the absence of external forces and torques. Methods to enforce different types of boundary conditions are described; these are exact for linear displacements, for constant stresses and for free surfaces. The performance of the method is assessed in a number of numerical tests.


Introduction
Particle scale simulations of the mechanics of granular materials at high relative density remain challenging [27,51]. The main reason for this is that particle contacts cease to be independent of each other at a relative density of about 0.7-0.8 [37,49], invalidating the commonly used discrete element method (DEM) [16] in its standard form.
The most common way to resolve this issue is to resort to the multi-particle finite element method (MPFEM; also referred to as the meshed DEM or the combined finite/discrete element method, FEM/DEM), whereby the discretization of each individual particle into finite elements enables an adequate representation of its mechanical behavior [28,55]. Although the MPFEM has successfully been used to study compression of three-dimensional particle assemblies [23,37], it is not practical for large systems due to its high computational cost [27]. As an alternative, there have been some efforts to develop what is commonly referred to as nonlocal contact models for the interaction between particles. The first such model, due to Harthong et al. [36], was based on curve fitting of MPFEM results for elastoplastic particles and utilized a stiffness that tends toward infinity when the local relative density approaches unity, as appropriate for incompressible particles. The local relative density was in turn determined from a Voronoi tessellation of space. Another nonlocal contact model for elastoplastic particles has been developed and calibrated in our laboratory [44]. This model utilizes particles of a (truncated) spherical shape and builds on the conceptions of Arzt [3], who likened the deformation of individual particles at high relative densities to that of an extrusion process. Particle deformation is compensated for by a corresponding increase in the particle radius, so that the particle volume is held constant. In effect, the model utilizes a contact pressure that is considered to be a function of the current particle volume, estimated from a Voronoi construction. Additional work in this field has been carried out by Gonzalez and Cuitiño [34] and Brodu et al. [10], who have devised nonlocal contact models for elastic particles. Related recent work has focused on nonlocal contact models for use in bonded particle packings (the so-called cohesive or bonded DEM) as a basis for fracture mechanics analysis [11,12]. Finally, an interesting contact model has recently been devised by Giannis et al. who included a nonlocal term based on the estimated average pressure in each particle, as obtained from the interparticle forces [31]. Despite their practical usefulness in simulations of powder compression [53], the existing nonlocal contact models are simplistic and typically do not properly account for tangential forces.
A natural way to formulate more satisfactory models would be to combine the DEM with an intrinsically nonlocal method such a particle-based method for continua (smoothed particle hydrodynamics, SPH, and related methods) [32,46,47]. In order to keep the representation of interacting particles intact, it appears desirable to include nearest-neighbor interactions only. If so-called virtual contacts [5] are included, the contact network is often considered to form a Delaunay tetrahedralization of space [17]. Hence, the Delaunay tetrahedralization and its dual, the Voronoi diagram, both constitute natural starting points for such an analysis. (For simplicity, we in this work confine ourselves to monodisperse granular materials; the appropriate generalizations for polydisperse materials are the regular tetrahedralization and the power diagram [4].) To this end, we have recently shown how the most commonly used strain estimate for granular materials, due to Bagi [5], can be translated into a particle-based method for continua [24]. In effect, a compatible strain tensor [38] can be obtained for each particle as the volume average of the local deformation gradient over the corresponding Delaunay cluster (the union of all Delaunay tetrahedra that contain the particle in question; see also [58] where the term contiguous Voronoi cell is used). The Bagi strain thus is a Delaunaybased strain.
However, the Delaunay-based strain does, as any Delaunaybased quantity, suffer from a potentially serious shortcoming: It may change discontinuously when the reference state is updated (when a particle moves over the circumcircle of one of the tetrahedra) [39,64]. For this reason, a Voronoi-based strain would be preferable, similar to the one proposed by Satake [57] (see [6]).
The specific objective of this work therefore is to devise a particle-based method for continua that is based on a Voronoi strain estimate. Such a method is expected to constitute an essential building block for future macroscopically consistent discrete methods for granular materials at high relative density. It is also hoped that the developed method shall be of interest in its own right.
Our method is related to the Voronoi component of the hybrid Lagrangian Voronoi-SPH scheme developed for fluids by Fernández-Gutiérrez et al. [19,20]. Voronoi cells were recently used to estimate strain in DEM simulations via an SPH approach [52], but the underlying methodology is different from ours. In addition, Voronoi cells have been used in . For clarity, the domain boundary ∂D is indicated by a thick red line and the boundary ∂V a of Voronoi cell a by black lines. The thin gray lines drawn between the particle centers represent the underlying Delaunay tessellation. The black dots indicate the referential placement of boundary points (color online) SPH to refine simulations via particle splitting [14], to calculate the volume attributed to each particle in order to improve the accuracy of the SPH method [29,61] and as a means to obtain an optimal particle arrangement via centroidal Voronoi tesselation as obtained from Lloyd's algorithm [29]. Additional related work includes the Voronoi finite element method due to Ghosh et al. [30] and the virtual element method on arbitrary polyhedral meshes [25]. Finally, the proposed method has some, but certainly not all, features in common with the recently proposed continuum-kinematicsinspired or kinematically exact peridynamics [41,42].

Basic definitions
As illustrated in Fig. 1, we consider a set of particles that between them form a material domain D. The material and spatial coordinates of (the center of) particle a are denoted X a and x a , respectively. A standard Vornonoi tesselation is used to subdivide the material domain D into Voronoi cells V a , one for each particle, such that each material point X is assigned to its closest particle center. The faces between particles thus are planar polygons, and two particles a and b are considered to be nearest neighbors if and only if they share a common face. To conveniently state and enforce boundary conditions, we introduce ghost particles as mirror images of particles located at the boundary ∂D of D [15,19] and consider ∂D to be the union of a number of planar polygons.
As a result, the boundary ∂V a of V a can be decomposed into a number of planar faces, such that each face ∂V ab either separates two real particles a and b or a real particle a and a ghost particle b. This procedure is straightforward when the computational domain is bounded by flat surfaces, as in Fig. 1, since the Voronoi faces then reproduce the boundaries exactly, and can be extended to complex geometries [1]. Specifically, for each particle a located at a boundary, we determine a boundary point W ab by closest point projection of X a onto the boundary (indicated by black dots in Fig. 1) and determine the coordinates X b of the ghost particle from the equation W ab = (X a + X b )/2. The volume V a of the Voronoi cell V a is a function of the location of particle a and its nearest neighbors b. In the following, this set of particles (including particle a itself) will collectively be denoted by P a . Similarly, P a denotes the set of (proper) nearest neighbors to a, not including particle a itself. When stating boundary conditions (BCs), it will be convenient to subdivide the set P a of nearest neighbors to a into the set E a of external (ghost) particles and the set I a of internal (real) particles. The reduced set P a is similarly divided into the sets E a and I a , where the prime again is used to indicate that particle a itself is not included. It proves convenient to subdivide the set E a into the sets E a and E a . Specifically, a ghost particle b ∈ E a belongs to the set E a if and only if the current location of the boundary point w ab = (x a + x b )/2 is fully prescribed. Consequently, the set E a contains those ghost particles for which w ab is not fully prescribed. Finally, we let E T a denote traction boundaries, i.e., a ghost particle b ∈ E a belongs to the set E T a if and only if the traction is prescribed on the interface between particles a and b. Note that E T a is a subset of E a , since the displacement is not fully prescribed on traction boundaries (most often not prescribed at all). The vectors X ab = X b − X a and x ab = x b − x a that point from the center of particle a to the center of a neighboring particle b in the reference and current configurations, respectively, are often referred to as branch vectors in the granular mechanics literature [5,17].

Global particle deformation: the mean deformation gradient
For each particle a, we define a discrete mean deformation gradient F a as an average of the continuum deformation gradient, F = ∂ϕ/∂ X = Grad ϕ, over the corresponding Voronoi cell V a . Here, x = ϕ(X, t) is a motion that expresses the spatial coordinates x as a function of the material coordinates X and time t. From the definition of the deformation gradient and the divergence theorem, one obtains whereN is the material outward unit normal to the surface ∂V a of V a . The final equality follows from a decomposition of the surface integral into integrals over the planar Voronoi faces ∂V ab with material outward unit normalsN ab . To proceed further, we assume that ϕ(X, t) varies linearly either on each face or in the whole Delaunay cluster [24] for particle a. The Delaunay cluster represents the union of all Delaunay tetrahedra that contain a certain particle. As such, it could be interpreted as the convex hull of the set of nearest neighbors of the particle. Assume first that ϕ(X, t) varies linearly on each face. Introducing the material centroid C ab of face ∂V ab with area A ab as we can write ϕ = c ab + G ab (X − C ab ), where c ab = ϕ(C ab , t) and G ab is the constant gradient. For the assumed linear variation of ϕ(X, t) on each face, Eq. (1) takes the form where A ab = A abN ab is the material area vector. However, this equation is not useful unless a practical way to determine c ab = ϕ(C ab , t) is found. For this reason, we instead follow Springel [64] and assume that ϕ(X, t) varies linearly in the whole Delaunay cluster, so that ϕ = x a +G a (X−X a ), where x a = ϕ(X a , t) and G a is the constant gradient. Hence, the integrand in Eq. (1) takes the form It is realized that only the term [G a (X − X a )]⊗N produces a nonzero result in Eq. (1), because a straightforward application of the divergence theorem demonstrates that the surface integral of the first term can be transformed to the volume integral of Grad x a = 0. To proceed further, we utilize to following tensorial analogue of the vectorial triple product expansion ('BAC-CAB' rule; see Appendix A): Only the term (G aN ) ⊗ (X − X a ) produces a nonzero result in Eq. (1), because another application of the divergence theorem shows that the surface integral of the second term corresponds to the cross product of G a and the volume integral of Curl(X − X a ) = 0. It can be noted that G aN ab is the directional derivative of ϕ in the directionN ab , which for a linearly varying function can be expressed as x ab /L ab , where x ab = x b − x a is the spatial branch vector and L ab = |X ab |.
Since each directional derivative G aN ab is independent of X, Eq. (1) becomes where we have introduced the material vector henceforth referred to as the gradient vector.
Since the deformation gradient reduces to the identity tensor I when the material coordinates are substituted for the spatial coordinates in Eq. (6), we obtain the following duality relation so that Eq. (6) can be restated as where the sum now includes particle a itself.

The gradient vectors
As demonstrated by Flekkøy et al. [21] and Serrano and Español [60], the gradient vectors B ab exhibit interesting relationships with derivatives of the Voronoi cell volume V a with respect to its defining points b ∈ P a (see Appendix B). It is also instructive to decompose the gradient vectors B ab into components that are normal and tangential to the face ∂V ab , denoted by B n ab and B t ab , respectively. From definition (7) it is realized that the normal and tangential components behave differently upon interchange of the indices a and b: As illustrated in Fig. 2, the normal component is antisymmetric (i.e., B n ba = −B n ab ) and the tangential component is symmetric (i.e., B t ba = B t ab ). Moreover, it is realized that is the material area vector. This result has two important consequences. First, the material area vector can be expressed as a difference between oppositely directed gradient vectors: Second, summing Eq. (11) over b ∈ P a , noting that the sum of outwards directed area vectors A ab over a closed surface vanishes as a consequence of the divergence theorem, it is seen that Adding B aa to each member of Eq. (12), we obtain according to definition (9) a vanishing sum, i.e., This result carries significant practical implications. Since the deformation gradient shall be invariant under translations, the substitution of (10) shall have no effect. Clearly, this is true only if the sum of B ab over b vanishes. As demonstrated below, the patch test will not be fulfilled in general unless the sum of B ba over b vanishes.

Boundary treatment
The location of ghosts needs to be updated when the material deforms or undergoes rigid-body motion. This is straightforward when the location of boundary points is prescribed but less so on free boundaries. The location of free (coordinates of) boundary points, or equivalently of free (components of) branch vectors x ab , can be determined from a least-squares minimization [15,19]. Specifically, the location is here determined from minimization of the functional b∈P a |x ab − F a X ab | 2 over the space of branch vectors x ab subject to the constraints resulting from known locations of internal points or imposed displacement BCs. In order to enforce (total, symmetric and antisymmetric) displacement BCs, we introduce for each particle a located at the boundary of the domain D and each ghost particle b ∈ E a a projection D ab onto the constrained degrees of freedom and for now assume that the remaining degrees of freedom are free. Specifically, whereN ab is a unit normal to the face (for simplicity assumed to be time independent, thus coinciding with the material normal). Likewise, we let D ab = 0 on free boundaries (since there are no constraints) and D ab = I for b ∈ I a (since the corresponding branch vector then is known). With the aid of the projection D ab , we may specify constraints as D ab x ab = D abxab , where the relevant components ofx ab are considered known. Specifically, for displacement BCs, x ab = 2(W ab +ū ab − x a ), whereū ab the prescribed displacement of the boundary point.
Introducing Lagrange multipliers λ ac , the constrained minimization stated above is transformed into an unconstrained minimization of Using Eq. (6), we may expand F a X ac as where H acd = X ac · B ad /V a . The stationary values of x ab are obtained by letting ∂L a /∂ x ab = 0. A straightforward but fairly lengthy calculation shows that where Letting ∂L a /∂λ ab = 0 reproduces the constraint, i.e., Since D ab is a projection and hence symmetric and idempotent, we find by multiplication of Eq. (17) by D ab from the left that which when substituted back in Eq. (17) produces Clearly, Eq. (21) is trivially satisfied when the branch vector x ab is known, i.e., when D ab = I, which is the case for interior particles (b ∈ I a ) or for total displacement bound- , at least one component remains after projection, implying that the sum must vanish. We may form a coefficient matrix K a from the elements K abc with b ∈ E a and c ∈ P a and extract a symmetric matrixK a from the elements K abc with b, c ∈ E a . Provided that this matrix is nonsingular, we may for each b ∈ E a solve for x ab . Letting x ab denote the obtained solution, we may write where the coefficients K abc are obtained from the matrix It is sometimes convenient to separate the sum in the above equation into sums over internal (real) and external (ghost) particles and to write where we have let K aba = − c∈I a K abc . Pre-multiplying x ab , as obtained from Eq. (22) or (23), with I − D ab , we obtain the free components of x ab . The constrained components are obtained from Eq. (19), and hence, we finally obtain where the required components ofx ab are obtained from the displacement BCs.

A possible static condensation
We may separate the sum in Eq. (10) into sums over internal and external particles and use Eq. (24) in the latter. For convenience introducinḡ Eq. (10) can be restated as It can be noted that Eq. (26) reduces to Eq. (10) for internal particles (i.e., particles not located at the boundary). For boundary particles that are not subject to any displacement BCs, D ab = 0 in Eq. (24) and the second sum in Eq. (23) vanishes. Inserting the resulting expression in Eq. (26), we obtain wherẽ Equation (27) represents a static condensation that can be used to evaluate the mean deformation gradient for boundary particles not subject to any essential boundary conditions.

Local particle deformation: local displacement
As in our previous work [24], we introduce a local displacement vector u loc ab as a measure of nonaffine deformation of the bond between particles a and b. For contact between two (real) particles (i.e., for b ∈ I a ), this vector is defined as where y ab = F ab X ab represents the image of the material branch vector X ab under the affine transformation represented by the mean deformation gradient F ab = (F a + F b )/2. In order to avoid complications resulting from a nonzero local displacement on free boundaries, we include the projection D ab in the definition for contact between a real particle a and a ghost b (i.e., for b ∈ E a ), i.e., where y ab = F a X ab . We can summarize Eqs. (29) and (30) by writing u loc

Variational total Lagrangian formulation
As in our previous work [24], we derive the equations of motion in a total Lagrangian framework. The starting point is a discrete Lagrangian L = T − V, where T is the total kinetic energy of the particle system and V is its potential energy (see, e.g., [33]). Specifically, where m a and v a = |v a | are the mass and velocity of particle a, respectively. The potential energy V is expressed as where V int is the internal strain energy resulting from the global (affine) particle deformation, V cnt is the strain energy resulting from local (nonaffine) deformation at each contact and V ext is the energy imparted by agents external to the system, such as gravity. Specifically, where U a is the internal energy per volume unit in the reference state, considered to be a function of the mean deformation gradient F a . Similarly, where U ab represents a pairwise interaction energy and the sum extends over all pairs of nearest neighbors a and b. For contact between internal (real) particles, U ab is considered to be a function of the current branch vector x ab (=x ab ) and the vector y ab (=ỹ ab ), defined in Sect. 2.6, in order to facilitate a discrimination between the normal (x ab ) and tangential directions. In order to obtain a unified treatment, U ab is expressed a function ofx ab = D ab x ab andỹ ab = D ab y ab for contact between an internal (real) and an external (ghost) particle. This is possible because the normal and tangential directions are in this case defined by the normal to the boundary and not by the branch vector x ab . Considering a body force resulting from gravity (gravitational constant g) and nominal tractionsT ab acting on the boundary between the internal (real) particle a and the external (ghost) particle b, the energy due to external forces can be expressed as The Euler-Lagrange equations (see, e.g., [33]) result in the equations of motion where are referred to as the internal, contact and external forces, respectively. We remark that some care is needed to ensure a proper distinction between internal and external forces when ghost particles are used, as elaborated further upon below.

External force
Combining Eqs. (38c) and (35), the external force is immediately obtained as Moreover, for displacement BCs, reaction forces yield an additional contribution to the external force, as elaborated further upon below.

Internal force
Using Eqs. (38a) and (33), together with the chain rule, the internal force is obtained as where the colon indicates double contraction and P b = ∂U b /∂ F b is the first Piola-Kirchhoff stress tensor for particle b. If particle a is located in the interior of the domain (i.e., if none of the particles b ∈ I a is located at the boundary), there is no need to explicitly take boundary effects into account. Hence, differentiation of expression (10) (most easily done in Cartesian components; see [24]) and substitution of the result in Eq. (40) produces the internal force in the form For interior particles, it follows from Eq. (41) that the force resulting from a constant stress field vanishes provided that the sum of B ba over b vanishes. This results is a special instance of the integration constraint derived by Chen et al. [13,56]. That the integration constraint indeed is fulfilled follows from Eq. (13). Using definition (9), noting that P a = I a for particles in the internal of the domain, the internal force may alternatively be expressed as where f int ab = P a B ab − P b B ba represents the force on particle a caused by contact with particle b. When written in this form, the expression for the internal force embodied in Eq. (42) is valid also for particles located at boundaries. In the special case that P b = P a , it follows from Eq. (11) that the internal force can be expressed as f int ab = P a A ab = A ab T ab , where T ab = P aN ab represents the nominal traction on the interface between particles a and b. Likewise, for a total displacement BC, the reaction force on particle a caused by contact with ghost particle b is obtained as P a A ab . The reaction force represents a contribution that is to be added to expression (39). Hence, for a state of constant stress ( P b = P a , without body forces), the total force f int a + f ext a vanishes, as it should. When only some degrees of freedom are prescribed, the projection operator D ab shall be included to ensure that no spurious forces are picked up from free boundaries, i.e., the reaction force takes the form D ab P a A ab . The formulation described this far is fully satisfactory for static and quasi-static problems. In particular, it can be noted that identical expressions are obtained for forces due to applied tractions and for reaction forces (compare f ext ab = A abT ab with f ext ab = A ab T ab where T ab = P aN ab ). However, an alternative expression for the internal force is preferable in highly dynamic situations. For particles located at free boundaries, the internal force provided by Eq. (41) is not in general derivable from a rotationally invariant potential energy, and as a consequence, the total angular momentum will not be conserved when this expression is used. To alle-viate this issue, one can base the derivation on Eq. (26) rather than on Eq. (10). Moreover, (components of) branch vectors corresponding to prescribed displacements [i.e., D abxab , c.f. Eq (24)] are treated as constants to ensure that no spurious external forces are obtained. Hence, differentiation of Eq. (26) and substitution of the derivative into Eq. (40) produces Clearly, Eq. (43) reduces to Eq. (41) in the interior of the domain, as is should. It is also instructive to consider the special cases that all degrees of freedom are free or constrained.
which is identical in form to expression (41)

Contact forces
The contact forces are obtained from Eqs. (38b) and (34) together with the chain rule. The derivation parallels the one outlined in [24] but is more technical and is therefore deferred to Appendix C. The final result can be expressed as where M bca is a symmetric second-order tensor defined in Appendix C and where we have let (47) in order to be able to combine the terms obtained for internal (real) and external (ghost) particles. Alternatively, one can letM bca = M bca E bc /2 to simplify the appearance of the second sum somewhat. If b is located in the interior of the domain (i.e., not at the boundary),M bca reduces toH bca I withH bca ≡ H bca /2 so that the matrix-vector multiplication can be avoided.

Constitutive equations
We use the same constitutive equations as in our previous work [24], to which the reader is referred for further details.
The mean (first Piola-Kirchhoff) stress was inferred from a neo-Hookean material model with strain energy of the form where C a = F T a F a is the mean right Cauchy-Green deformation tensor for particle a, J a = det F a is the determinant of F a and tr C a is the trace of C a . The Lamé parameters λ and μ were calculated from Young's modulus E and Poisson's ratio ν in a standard manner. The contact forces were derived from an elastic contact model with contact energy The normal (K n ab ) and tangential (K t ab ) contact stiffness were expressed as where E is Young's modulus, G = μ is the shear modulus, A ab is the initial area of the interface and L ab = |X ab | is the initial distance between particle a and b. Finally, ξ n and ξ t are nondimensional parameters. Artificial viscosity was included in some simulations, in which case a viscous stress was incorporated as an addition to the elastic first Piola-Kirchhoff stress. This addition has the form P visc a = J a σ visc a F −T a (see, e.g., [35]), where σ visc a is the viscous Cauchy stress for particle a. Assuming a homogenous material, the latter was calculated as [48] σ visc where d a is the rate of deformation tensor for particle a, ρ is the material density, c p = √ (λ + 2μ)/ρ is the longitudinal wave speed and h is a characteristic size, here taken as the minimum length of the material branch vectors. Finally, c 1 and c 2 are two nondimensional constants.

Implementation details
Delaunay tetrahedralization and Voronoi tessellation were performed using TetGen (version 1.6) [62]. A standard leapfrog scheme was used to integrate the equations of motion in time, as is commonly done in SPH [50], and the velocity was calculated as described in [22] when damping was applied.

Patch test
A patch test was used to assess the basic characteristics of the method. An irregular arrangement of 3 × 3 × 3 particles was considered. The particles constituted a cube with side length 1 m, occupying the referential domain 0 ≤ X 1 , X 2 , X 3 ≤ 1 m. For illustrative purposes, the particles are displayed as gray spheres in Fig. 3 and the Voronoi cells formed from them and their ghost particles (not shown) are drawn in blue where α = 1 × 10 −4 . The material parameters were E = 1.0 MPa, ν = 0.3 and ρ = 1.0 × 10 3 kg/m 3 and the contact parameters were set to unity (i.e., ξ n = ξ t = 1). Following Taylor et al. [66], three different variants of the test were performed. First, the displacement was prescribed for all particles and force equilibrium was tested for each particle (test A). Second, the displacement was prescribed on all boundaries, and the displacement error was determined for each particle (test B). Third, the displacement was prescribed on the bottom face (i.e., for X 3 = 0) and the traction corresponding to the displacement field (52) was prescribed on all other boundaries (test C). When equilibrium solutions were sought (tests B and C), the displacements and/or tractions were gradually applied during 50 s and damping was used. The obtained results are summarized in Table 1, in which the maximal errors (in the L 2 norm) in the force on or displacement of each particle are provided.
It is evident that all versions of the patch test are passed. (The error in test A is somewhat larger than those in tests B and C, since a minute error in displacement is translated into a larger error in force unless the material is very soft.) Satisfaction of test A demonstrates, firstly, that the expressions for the mean deformation gradient embodied in Eqs. (6) and (10) are exact for linear displacements and, secondly, that the Chen integration constraint [13] indeed is fulfilled. In essence, fulfillment of these requirements follows from the duality relation (8) and the fact that the sum of the gradient vectors over the first and second index vanishes, Eq. (13). Satisfaction of test B in addition shows that the stiffness is nonsingular, i.e., that no spurious zero-energy modes exist [65]. Finally, satisfaction of test C demonstrates that constant tractions can be prescribed exactly, as anticipated from the discussion following Eq. (42).

Swinging cube test
In this example, a cube with side length L = 1 m was considered, occupying the referential domain 0 ≤ X 1 , X 2 , X 3 ≤ 1 m. The normal motion was constrained on the boundaries at X 1 = X 2 = X 3 = 0 (symmetric BCs), whereas the tangential motion was constrained for X 1 = X 2 = X 3 = 1 m (antisymmetric BCs). The initial displacement field u 0 was prescribed according to the following analytic solution [8,59] which is valid in the linear regime for volume-preserving motions (U 1 + U 2 + U 3 = 0). Specifically, we let U 1 = No damping was used. This example is typically used to test the convergence of the solution [8,59] and it was here also used to assess the boundary treatment (cf. Sect. 2.4). The four slightly irregular configurations displayed in Fig. 4 were considered, in which particles are colored according to the magnitude or their initial shear (von Mises) stress. The error of the numerical solution after 5 × 10 −3 s, which corresponds to about 1/6 of the period of oscillation, was benchmarked against the analytical result. The results are summarized in Fig. 5, which displays the normalized displacement error, measured in the L 2 norm, as a function of the mean particle spacing. It is evident that the displacement exhibits a quadratic convergence toward the correct solution, as expected for a second-order accurate method. This result is anticipated in the present context, since our Voronoi gradient estimate [Eqs. (6) or (10)] provides the material gradient. It deserves to be mentioned, however, that convergence issues have been identified for Voronoi gradient estimates on moving meshes [54] unless proper account is made for the change in geometry that occurs during each time step. In fact, the original spatial Voronoi gradient estimate derived by Springel [64], upon which our work is based, was abandoned in favor of a least-squares gradient estimate in the public release of the cosmological code AREPO [67].

L-shaped block
In this test, a free L-shaped block with dimensions provided in Fig. 6a was subjected to oppositely directed transient nominal tractions on two of its faces (shown in gray in the figure).
where p(t) = max(2.5 − |kt − 2.5| , 0) is a triangular impulse function (the rate constant k = 1 s −1 has been included for dimensional consistency). The test was originally proposed by Simo [63] and has subsequently been used by many others ( [2] and references therein). Its main purpose is to assess momentum preservation for t > 5 s when no external forces are applied and, for this reason, the test was continued for an additional 50 s. The material parameters were E = 50.0 kPa, ν = 0.3 and ρ = 1.0 × 10 3 kg/m 3 and the contact parameters were set to unity (ξ n = ξ t = 1). No damping was used. As shown in Fig. 7, the block starts to rotate and undergoes significant deformations as a result of the applied forces. The linear and angular momenta (about the origin) are provided in Fig. 8a. As expected, the components of linear momentum remain negligible throughout the test and the components of angular momentum increase in magnitude during the first 5 s and thereafter stay constant. This claim is substantiated by the data presented in Fig. 8b, in which the change in linear momentum components from the start of the test and the change in angular momentum components from their value at the end of loading are displayed. Clearly, all momentum components can be considered to be conserved. This result is expected when the internal forces are calculated with Eq. (43) [but not with Eq. (41)], since the internal force is then derivable from a rotationally invariant potential energy.

Punch test
In this test, a rectangular specimen (width 1 m, height 0.5 m and thickness 0.1 m) was considered. As illustrated in Fig. 9, a rigid punch indented the middle third of the top face (shaded in the figure). The velocity of the punch was 0.1 m/s and it was assumed that friction was large enough to prevent any slip at the punch-specimen interface. The remainder of the top face was free. Symmetric BCs were enforced on all other boundaries (i.e., zero normal displacement). This test has previously been used to assess the performance of different formulations in the near incompressible limit (for ν → 1/2) [18,45]. However, our objective here rather is to assess the performance of the proposed method in a compressiondominated problem, and therefore, a compressible material is considered. Specifically, the material parameters are E = 1.0 MPa, ν = 0.45 and ρ = 1.0 × 10 3 kg/m 3 and the contact parameters are set to unity (ξ n = ξ t = 1). Viscous damping was included in this example (c 1 = c 2 = 0.1). Following [45], a discretization using 21 × 11 × 3 particles was used. A regular particle arrangement was used in order to facilitate the identification of potential instabilities, which would be seen as irregular particle displacements and/or as hourglass patterns [26]. However, no sign of instability can be observed in Fig. 10, indicating that the contact forces, obtained from Eq. (46), imparted sufficient stabilization to the formulation. It is commonly observed that nodally integrated meshless methods exhibit severe pressure oscillations for highly constrained problems [18,45]. However, no pressure oscillations are evident for this relatively compressible material, although a strong pressure gradient exists in the vicinity of the punch boundary. This is not surprising, since an analysis based on linear elasticity reveals that the pressure not only becomes infinite at this point but also exhibits an oscillatory behavior indicative of the breakdown of the linear theory [43, p. 49]. Note, however, that the described displacement-based method is not optimal in the near incompressible limit (when ν → 1/2).

Conclusions
A Voronoi-based estimate of the mean deformation gradient has been derived and translated into a particle-based method for continua. The gradient estimate is expressed in terms of material gradient vectors, defined for each particle and its nearest neighbors, and has a structure that resembles strain estimates used in smoothed particle hydrodynamics and peridynamics. For each particle, the gradient vectors can be computed from the geometry of its material Voronoi cell (particle placement, face centroids and normals and cell volume). They can also be expressed in terms of the material gradient of the Voronoi cell volume. A decisive advantage of this Voronoi-based strain estimate over Delaunay-based estimates is that it is guaranteed to evolve smoothly when the reference state is updated. It is therefore expected to provide an important corner stone in future developments geared toward macroscopically consistent discrete method for granular materials. Moreover, the sum of the gradient vectors over either index vanishes, a fact that results in a formulation that not only is linearly complete but also satisfies the patch test when forces are derived from a discrete Lagrangian. Pairwise forces between particles, inferred from the local (nonaffine) deformation of each bond or contact, provide a physical stabilization. Procedures to enforce different types of boundary conditions in a consistent manner are provided. The presented numerical tests demonstrate that the formulation satisfies different versions of the patch test, exhibits a second-order convergence toward the correct solution and conserves linear and angular momenta in the absence of external forces and torques.

Conflict of interest
The author declares no conflict of interest.

Availability of data and material Not applicable.
Code availability Available from the author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Tensorial cross products
Although generalizations of the vectorial cross product to tensors have been discussed in textbooks [7,40] and in recent scientific literature [8,9], they may be less familiar to the reader. For our purposes, it is sufficient to note that the cross product between a second-order tensor T and a vector v produces a second-order tensor, denoted T × v, defined so that the identity holds for all vectors a. This definition is consistent with (and in fact follows from) the usual definition of the cross product between vectors, as can be seen by expanding T , v and a in an orthonormal Cartesian bases (i.e., T = T i jêi ⊗ê j , v = v kêk and a = a ê whereê i , etc., denote basis vectors). To prove Eq. (5) in the main text, we let the second-order tensor G a × [N × (X − X a )] operate on an arbitrary vector V , to obtain where the first equality follows from the definition (54), the second from the triple product expansion ('BAC-CAB' rule) for vectors and the third from the definition of the tensor product. Since the vector V is arbitrary, Eq. (5) follows.

B Gradient vectors as derivatives of Voronoi cell volumes
Using a smoothed characteristic function of the Voronoi cells, Flekkøy et al. [21] and Serrano and Español [60] have provided elegant derivations of a number of Voronoi cell properties. Of particular relevance to this work are the derivatives of the Voronoi cell volume V a with respect to the coordinates of its defining points X b , with b ∈ P a . These results can also be straightforwardly derived as follows. For b = a (i.e., b ∈ P a ), it proves convenient to introduce a local coordinate system with origin at (X a + X b )/2, one axis (ζ ) directed along the material branch vector X ab = X b − X a and the other two mutually perpendicular axes (ξ andη) in the plane of the face ∂V ab , see Fig. 11. According to the definition of Voronoi faces, a change of X b in the normal direction ζ by a small amount Δζ produces a translation of the face ∂V ab in the same direction by an amount Δζ /2 (Fig. 11a). Edge effects are of second order, and the change in volume of the Voronoi cell is obtained as ΔV = A ab Δζ /2. The directional derivative of V a with respect to X b in direction ζ , denoted by Dζ V a (X b ), thus becomes A change of X b in the tangential directionξ by a small amount Δξ produces a rotation by an angle Δξ/L ab of the face ∂V ab around an axis that is parallel to theη axis and passes through particle a (Fig. 11b). To the first order, the change in volume of the Voronoi cell can therefore be calcu- where C ξ is the ξ coordinate of the centroid of the face ∂V ab . Hence, An analogous calculation shows that Combination of Eqs. (56), (58) and (59) produces Comparing this result with the definition (7) of the gradient vectors, one arrives at the interesting conclusion that In addition where the sum does not include particle a itself.

C The contact force
The contact forces are obtained by differentiation of the pairwise interaction energy (34) according to Eq. (38b). It proves convenient to decompose this energy into contributions resulting from contact between internal (real) particles and between an internal and an external (ghost) particle. In this manner, one obtains where the factor 1/2 is included because each pair is otherwise counted twice and where it is understood that a = b in the first sum. A straightforward application of the chain rule shows that The derivatives ∂U ab /∂x ab and ∂U ab /∂ỹ ab are obtained from the material constitution; cf. Sect. 4. Moreover, the derivatives ∂x ab /∂ x ab and ∂ỹ ab /∂ y ab equal I if b ∈ I a and D ab if b ∈ E a ; cf. Sect. 2.6. The task thus boils down to determining the derivatives ∂ x ab /∂ x c and ∂ y ab /∂ y c . If b ∈ I a , straightforward differentiation of the branch vector where δ ab is the Kronecker delta. The result is less immediate if b ∈ E a , in which case we introduce the shorthand notation M abc = ∂ x ab /∂ x c . When all displacement components are prescribed (i.e., for b ∈ E a ), D ab = I in Eq. (24). Hence, x ab =x ab and the derivative becomes Otherwise (i.e., for b ∈ E a ), we find by differentiation of Eq. (24) that whereK abc is defined as In order to determine the derivative of y ab , defined in Sect. 2.6, we introduce the vector z ab = F a · X ab . Since the material branch vector X ab is antisymmetric (X ba = −X ab ), we find that y ab = z ab − z ba if b ∈ I a and that y ab = z ab if b ∈ E a . Using Eq. (26), z ab can be expanded as whereH abd = X ab ·B ad /V a , H abd is defined in conjunction with Eq. (16) and where x ad is given by Eq. (24). Differentiation of the above equation with the help of Eqs. (66) and (67) With these results in hand, the contact force is obtained from the chain rule, cf. Eq. (46) in the main text. The tensor M abc does not appear in the final expression, since it enters in the form of the product D ab M abc = −2δ ac D ab as a result of the orthogonality of D ab and I − D ab .