Towards a macroscopically consistent discrete method for granular materials: Delaunay strain-based formulation

We demonstrate that the Delaunay-based strain definition proposed by Bagi (Mech Mater 22:165–177, 1996) for granular media can be straightforwardly translated into a particle-based numerical method for continua. This method has a number of attractive features, including linear completeness and satisfaction of the patch test, exact conservation of linear and angular momenta in the absence of external forces and torques, and anti-symmetry of the gradient vectors for any two points not both on the boundary of the computational domain. The formulation in effect relies on nodal (particle) interpolation of the deformation gradient and is therefore inherently unstable. Drawing on the analogy with granular media, a pairwise interaction between particles is included to alleviate this issue. The underlying idea is to define a local, non-affine deformation of each bond or contact, and to introduce pairwise forces via a stored-energy functional expressed in terms of the corresponding local displacements. In this manner, a generalisation of the Ganzenmüller (Comput Methods Appl Mech Eng 286:87–106, 2015) hourglass stabilisation procedure to non-central forces is obtained. The performance of the method is demonstrated in a range of problems. This work can be considered a first step towards the development of a macroscopically consistent discrete method for granular materials.


Introduction
Smoothed particle hydrodynamics (SPH) is a meshfree particle-based numerical method that was originally proposed by Lucy [36] and Gingold and Monaghan [21] in a cosmological context more than 40 years ago. The heart of the method is a kernel density estimate, i.e. the density of the material at a certain point is calculated from the number, mass and proximity of nearby particles using a smoothing kernel [34]. The method has a number of attractive features such as the ability to straightforwardly incorporate advection and convection, to naturally trace moving boundaries, to handle large deformation and its comparative ease of implementation [34]. Since its inception, SPH has been applied in a number of fields, including solid mechanics, starting with the work by Libersky and Petschek almost 30 years ago [31,32]. However, it is well-known that the standard SPH formulation suffers from a number of issues that primarily are related to the stability of the method and the consistency of its underlying kernel and particle approximations. Firstly, the tensile instability, that causes particle clumping, as first observed by Swegle [52]. Secondly, spurious zeroenergy modes caused by rank deficiency of the stiffness matrix resulting from the commonly used particle (nodal) integration [13,14]. Thirdly, lack of linear completeness of the particle approximation for arbitrarily located particles. This issue is aggravated at boundaries, since the supporting domain of the kernel density estimate will be truncated at boundaries [35]. Fourthly, high-frequency pressure oscillations in shock-dominated problems, as a result of the collocation of mass at the particles [39,40]. These oscillations are particularly evident in the near-compressibility limit.
In the light of the Lax-Richtmyer equivalence theorem, stating that stability and consistency together are equivalent to convergence, these issues would appear to severely limit the practical usefulness of SPH. Fortunately, this is not so, because, there are a number of more or less efficient methods that can be used to resolve them. Whereas spurious zero-energy modes are a generic feature of SPH, the tensile instability is only present when Eulerian kernels are used [2]. It is thus possible to avoid the tensile instability using total [43,56] or updated [54] Lagrangian formulations. It has also been shown that an artificial stress can counteract the tensile instability for Eulerian kernels, and the resulting algorithm, due to Grey et al. [24], has found widespread use. Drawing on the similarity with the hourglass modes commonly encountered in underintegrated finite elements, Ganzenmuller et al. have devised a promising means to alleviate the spurious zero-energy modes in SPH [16,18]. Procedures to restore consistency of the SPH approximations have been discussed in [35]. An often used correction was proposed in [5], which also addressed variational aspects of the SPH in order to derive formulations that preserve linear and angular momentum exactly in the absence of external forces and torques [4]. The first-order hyperbolic framework for large strain computational solid dynamics proposed by Bonet and co-workers [3,20] is also a noteworthy line of developments. This formulation utilises an extended set of strain measures and results in a mixed form that has been implemented using an SPH approach with very good results in the near and truly incompressible limit [28][29][30].
Despite being a particle-based method, SPH has hitherto primarily been used to model continuous matter. For discontinuous matter, such as granular materials, the discrete element method (DEM) is instead commonly used. The DEM was proposed by Cundall and Strack more than 40 years ago [9] and is, like SPH, a Lagrangian particle-based method. However, contrary to SPH, the interparticle forces and torques are generally determined from the (normal and tangential) overlaps between particles, which represent local particle deformations. Hence, the deformation resulting from a certain contact is considered to be local in the sense that it does not affect the overall shape of the particle or the stiffness of any other contact(s) on the same particle. However, more general non-local contact models are emerging [6,19,23,27].
It is interesting to note that the traditional SPH and the DEM thus are complementary in the sense that SPH considers only nonlocal forces and the DEM only local forces. These methods are also strongly related to peridynamics, developed by Stewart Silling at Sandia National Laboratories, which comes in two flavours: bond based [48] and state based [50]. In bond-based peridynamics, particles interact via forces directed along bonds between particles, i.e. via central forces. However, the assumption of central forces leads to severe restrictions of the class of materials that can be modelled using bond-based peridynamics. In particular, the Poisson ratio is always 1/4 for isotropic materials [51]. For this reason, peridynamic states were introduced, which in effect result in a nonlocal coupling between bonds [50]. This means that arbitrary material models can be adapted for use within state-based peridynamics, as for SPH. In fact, it has been pointed out [17] that the commonly used numerical discretisation of state-based peridynamics [50] is equivalent to a total Lagrangian corrected SPH [5].
Acknowledging the complementary nature of SPH and the DEM, the purpose of this work is to take the first steps towards a macroscopically consistent discrete method for granular materials. Specifically, a particle-based method for continua is developed whose SPH-like component is based on the most widely accepted strain definition for granular media, a Delaunay-based strain developed by Bagi [1] and elaborated upon by others [11,12,26]. This strain definition is based on a decomposition of space into simplices and thus has much in common with the Delaunay Density Estimator Method proposed by Schaap and Weygaert [45], suggested for use in SPH by Pelupessy [41]. The DEM-like component is based on pairwise interactions between particles and produces a generalisation of the Ganzenmuller [16,18] hourglass stabilisation procedure to non-central forces. Aiming primarily for a description of solid materials, the energy corresponding to rotation of individual particles is disregarded, as is common practice in SPH. Although the starting points are distinct, it will be seen in the following that the resulting method has strong similarities with the one obtained by nodally integrated finite elements, in particular nodally averaged tetrahedral elements [10,42].

Basic definitions
We consider a system of discrete particles occupying a connected domain D and let lowercase letters in the beginning of the alphabet (e.g. a, b and c) refer to particles. The reference placement of particle a is X a and its current position is x a . We assume the existence of particles located on the boundaries (boundary particles), such that the domain D can be decomposed into simplices via a triangulation (in 2D) or a tetrahedralisation (in 3D) in which the particles constitute the vertices; in effect, we thus assume the existence of an underlying triangular or tetrahedral mesh. We will adapt our language to the three-dimensional case but will, for simplicity, provide schematic illustrations in two dimensions, as in Fig. 1. The particles are assumed to interact via nearestneighbour interactions only, with bonds corresponding to edges of the underlying mesh. The vectors X ab = X b − X a and x ab = x b − x a point from the centre of particle a to the centre of a neighbouring particle b in the reference and current configurations, respectively. These are often referred to as branch vectors in the granular mechanics literature [1,12].  Fig. 1 a Definition of Delaunay clusters based on an underlying tesselation into simplices (triangles in 2D and tetrahedra in 3D). A Delaunay cluster for an internal particle labelled a is shown in red and a cluster for a boundary particle labelled b in blue. b Illustration of the Delaunay cluster volume V (c) a and the particle volume V a (color online) Since distinct simplices are disjoint, the total volume of the domain D, denoted V tot , can be obtained as where the sum extends over all simplices S α in the system. For each particle a, we introduce the Delaunay cluster C a as the union of all simplices that contain particle a, i.e. [45] where the union thus extends over the set of simplices that have particle a as one of its vertices (Fig. 1a). The Delaunay cluster has been referred to as the contiguous Voronoi cell [45] but this denomination is not used here, since the Delaunay cluster represents a set of adjoint Delaunay simplices rather than a Voronoi cell. Since distinct simplices are disjoint, the volume of the Delaunay cluster is obtained as where the sum has the same extent as the union preceding equation. Since each simplex has n + 1 nodes (where n is the number of spatial dimensions), each simplex contributes to n + 1 Delaunay clusters, implying that where the sum extends over all particles in the system. For each particle, we will therefore define two volumes: Firstly the volume V (c) a = Vol C a and secondly the volume V a = Vol C a /(n + 1). We will refer to V (c) a as the Delaunay cluster volume and to V a as the particle volume. The difference between V (c) a and V a is illustrated in Fig. 1b.

Strain measures
For each particle a, we define a discrete mean deformation gradient F a as the average of the continuum deformation gradient, over a suitably selected volume containing particle a. Here x = ϕ(X, t) is a motion, depending on the reference placement X of a material particle and time t. As noted by Bagi [1], Schaap [45] and He [26], the Delaunay cluster represents a natural choice (at least as long as a total Lagrangian formulation is sought). The details of the derivation are deferred to Sect. 4. Similarly, we introduce a mean deformation gradient for each particle pair ab as where an unweighed average is used for simplicity. The mean deformation gradient F ab maps the referential branch vectors X ab to y ab = F ab X ab . Note that y ab = x ab in general, because F ab only accounts for the affine particle deformation. As a measure of non-affine deformation of the bond between particles a and b, we introduce a local displacement vector u loc ab , defined as A geometric interpretation of the local displacement vector is provided in Fig. 2. It can be compared to the corresponding total displacement vector It will also be convenient to introduce the unit normal vectorŝ so that the (total or local) displacement u ab can be decomposed into normal and tangential parts, and respectively. Fig. 2 a Definition of referential contact vectors X ab /2 and X ba /2 from branch vectors where X a and X b are the referential placements of particles a and b. b Definition of the local displacement vector u loc ab from the current particle coordinates x a and x b and the images F a X ab /2 and F b X ba /2 of the referential contact vectors under the affine transformations described by F a and

Variational total Lagrangian formulation
In order to derive the equations of motion, we here consider a conservative system, for which the discrete Lagrangian L takes the form (see, e.g. [22]) The total kinetic energy T is expressed as where m a and v a = |v a | are the mass and velocity of particle a, respectively. For simplicity, the energy corresponding to rotation of individual particles is thus disregarded, as is common practice in SPH. The potential energy V is considered to be the sum of three terms: 1. The internal (strain) energy V int resulting from the global (affine) particle deformation, characterised by the mean deformation gradient F a . Hence, where U a is the internal energy per unit volume in the reference state and V a is the initial volume attributed to particle a. This term is commonly included in SPH (but not in the DEM). 2. The internal (strain) energy V cnt resulting from local (non-affine) deformation at each contact. Noting that both the local displacement vector u loc ab and its normal and tangential components can be expressed in terms of x ab and y ab [compare Eqs. (7), (10) and (11)], we write Here, U ab represents the pairwise interaction energy and the sum extends over all pairs of nearest neighbours a and b. This term is typically not included in SPH, but the stabilisation of SPH proposed by Ganzenmuller [17] is of this form, albeit restricted to normal forces. It is also related to the stabilisation of state-based peridynamics models proposed by Littlewood [33] and Silling [49].
The DEM does, on the other hand, rely on pairwise contact forces, but inferred from the total rather than the local displacement [c.f. Eqs. (8) and (7)]. 3. The external energy V ext , which for simplicity here is stated for the action of a constant gravitational force g, The Euler-Lagrange equations (for particle a) take the form (see, e.g. [22]) Hence, introducing the canonical momentum together with the internal, contact and external forces, one obtains the equations of motion in the familiar form The external force f ext a is immediately obtained as The internal and contact forces ( f int a and f cnt a , respectively) will be elaborated upon once the mean deformation gradient has been defined.

The mean deformation gradient
As already mentioned, we define the discrete mean deformation gradient F a as the average of the continuum deformation gradient (5) over the corresponding Delaunay cluster. To simplify the ensuing developments, we introduce the following notation, illustrated in Fig. 3: First, for a given particle a, we let P a denote the set of all particles that form the vertices of the associated Delaunay cluster C a . Hence, for a not on the boundary, P a corresponds to the set of nearest neighbours of a. For a on the boundary, P a in addition contains (a) and Fig. 3 Definition of the sets of particles a P a and b P aβ and of faces c F a and d F ab described in the text (color online) particle a itself. Occasionally, it will be convenient to refer to the proper set of nearest neighbours of a certain particle a, which we denote as P a . For particle a on the boundary, P a thus contains all particles in P a except particle a itself. For particle a not on the boundary, P a and P a coincide. Second, for a given face β of a Delaunay cluster C a , we let P aβ denote the set of particles that are located at the vertices of the face. Third, for a given particle a, we let F a denote the set of all faces of the associated Delaunay cluster C a . Fourth, for a given particle b of a Delaunay cluster C a , we let F ab denote the set of faces of C a that have particle b as one of its vertices. With this notation in hand, we proceed with the derivation of the mean deformation gradient. Using Eq. (5) and the divergence theorem, we obtain whereN is the referential unit outward normal to the surface ∂C a of C a . We let ∂C aβ denote the planar faces that between them form ∂C a , so that we can write SinceN is constant in each of the integrals in Eq. (23), it suffices to integrate ϕ over the faces ∂C aβ . Assuming that ϕ(X, t) varies linearly between the values ϕ(X b , t) at the vertices, we obtain where A aβ is the area of the face ∂C aβ . Using Eq. (24) and introducing the referential area vector A aβ = A aβN aβ , Eq. (23) takes the form Substituting Eq. (24) in Eq. (25) and interchanging the order of summations, we obtain a β∈F a b∈P aβ a b∈P a β∈F ab Introduction of the referential vector B ab = 1 n(n + 1) β∈F ab A aβ (27) yields the final result where we have used the fact that V a = V (c) a /(n +1). The thus defined vector B ab has been referred to as the complementary area vector by Bagi [1]. However, we will henceforth refer to B ab as the gradient vector, because of its strong similarity to the gradient vectors commonly employed in SPH (see below).
The gradient vectors exhibit a number of useful properties, as we next demonstrate. First, summation over the second index yields b∈P a since the sum of outwards directed area vectors over any closed surface vanishes according to the divergence theorem. Similarly, provided that not both particle a and b are on the boundary of the domain D, because the parenthesis contains a sum of outwards directed area vectors over a closed surface, which vanishes. Provided that not both particle a and b are on the boundary of the domain D, we thus from Eqs. (29) and (30) conclude that the gradient vectors are antisymmetric and that the sum over either index vanishes, i.e. that B ba = −B ab and that Equation (29) enables us to write Eq. (28) in the more symmetric form where x ab = x b −x a is the current branch vector for particles a and b. Since the deformation gradient reduces to the identity tensor I when the reference coordinates are substituted for the current coordinates in either of Eqs. (28) or (32), we immediately obtain the following duality relation where X ab = X b − X a is the referential branch vector.
Noting that the gradient vectors are referential and hence time independent, one immediately obtainṡ where the superposed dot denotes the material time derivative and where v b is the velocity of particle b. Hence, the velocity gradient l a and the rate of deformation tensor d a (see, e.g. [25]) can be determined as and respectively. The rate of deformation tensor can be used, e.g. to define viscous stresses.

Preliminaries
For future reference, we notice that where δ ca is the Kronecker delta. Using Eqs. (6) and (28) and the antisymmetry of X bc (i.e. that X bc = −X cb ), we may expand y bc = F bc X bc as where Since H abc is referential and thus independent of x a , this yields In order to evaluate the derivative of the mean deformation gradient, we resort to an orthonormal Cartesian coordinate system, with components labelled by i, j and k. Since B bc is referential, we find using Eqs. (28) that where the properties of the Kronecker delta have been used.

Internal force
From Eq. (19) and the chain rule, the internal force is obtained as where the colon indicates double contraction and where is the first Piola-Kirchhoff stress tensor for particle b. Using Eq. (41) in Eq. (42), expressed in component form, we obtain Hence, in direct notation, where we have explicitly indicated that the sum is to be taken over the set of nearest neighbours to particle a (including particle a itself if the particle is located on the boundary). One can note that the internal forces vanish for a constant stress field provided that the sum of the gradient vectors over the first index vanishes. Hence, fulfilment of Eq. (31) for particles not on the boundary is required for satisfaction of the patch test. This is nothing else than the integration constraint derived by Chen et al. [7] (see also [42]).
Each term can be made antisymmetric as follows: since the sum of the added terms vanishes as a result of Eq. (29). Hence, we may interpret as the contribution to the internal force on particle a due to contact with particle b.

Contact force
From Eq. (19) and the chain rule, the contact force is obtained as Inserting the expressions for ∂ x bc /∂ x a and ∂ y bc /∂ x a provided by Eqs. (37) and (40), respectively, one obtains Noting that the derivative ∂U bc /∂ x bc and ∂U bc /∂ y bc and the pre-factors δ ca −δ ba and H bca − H cba all change sign when b and c are interchanged, the sum over all particle pairs b < c is converted to a double sum over nearest neighbours b and c. Moreover, exploiting the properties of the Kronecker delta, the final result is obtained as whereH bca = H bca /2 has been introduced for convenience.
Here we have written P a and P b to indicate the proper set of nearest neighbours of a and b, respectively.

Material model
We assume a neo-Hookean material model with strain energy where C b = F T b F b is the mean right Cauchy-Green deformation tensor for particle b and J b = det F b = (det C b ) 1/2 . Using standard solid mechanics results (see, e.g. [25]) one obtains the second Piola-Kirchhoff stress tensor S b as (52) where C −T b is the transpose of the inverse right Cauchy-Green deformation tensor and I is the second-order identity tensor. Moreover, the first Piola-Kirchhoff stress tensor P b is obtained as The first Piola-Kirchhoff stress is used to calculate the internal force via Eq. (45). The Lamé parameters may be expressed as where E is Young's modulus and ν Poisson's ratio. We also note that the longitudinal (primary) and the transverse (secondary) wave speeds are obtained as respectively.

Contact model
We will here assume a linearly elastic contact model with normal and tangential stiffness K n ab = K n ba and K t ab = K t ba for contact between the nearest neighbours a and b. Hence the contact energy is written Using Eqs. (10) and (11), we thus obtain and ∂U bc ∂ y bc = − K n bc u loc,n bc + K t bc u loc,t bc .
In the special case that the normal and tangential contact stiffness coincide, i.e. when K n ab = K t ab = K ab , the above equations reduce to In this work, the normal and tangential contact stiffness are expressed as where E is Young's modulus and G = μ is the shear modulus. Two geometric quantities are included, namely A bc = 4×(|B bc |+|B cb |)/2, which is an estimate of the area of the interface between particle b and c [cf. Eq. (27)] and L bc = |X bc |, which is the initial distance between particle b and c. Finally, ξ n and ξ t are non-dimensional parameters.

Artificial viscosity
An artificial viscosity can be used to reduce oscillations in shock-dominated problems [39,40]. We here use a viscous (Cauchy) stress of the form [37] where c 1 and c 2 are two non-dimensional constants, c p the longitudinal wave speed and h is a characteristic size, here taken as the minimum length of the referential branch vectors. The viscous stress is incorporated as an addition to the elastic first Piola-Kirchhoff stress, obtained as (see, e.g. [25]) This type of dissipation can also be added constitutively, i.e. be considered to be part of the material rather than algorithmic behaviour [44]. When artificial viscosity is included, the sum of the elastic and viscous first Piola-Kirchhoff stress, Eqs. (53) and (62), is used to calculate the internal force via Eq. (45).

Implementation details
Based on the initial particle arrangement, a Delaunay tetrahedralisation was obtained using Tetgen, developed and maintained by Si [47]. As is commonly done in SPH [38], a standard leapfrog scheme was used to integrate the equations of motion in time. When damping was applied, the velocity was calculated as described in [15].

Numerical tests and examples
As described in Sect. 2.1, it is assumed that the computational domain can be decomposed into simplices with particles as vertices but the particle placement is otherwise arbitrary. In practice, however, it is preferable to avoid very short inter-particle distances, since the stable time step is proportional to the minimum particle separation when a Courant-Friedrichs-Lewy condition is used [8]. In all numerical tests and examples, the initial particle positions were based on a regular cubic lattice. In most of the examples, randomness was introduced by perturbing the particles in all directions, while respecting the location of faces and edges, by an amount not exceeding 25% of the regular particle spacing in each direction. The constitutive equations summarised in Sect. 6 were used throughout.

Patch test
The patch test in the forms described by Taylor et al. [53] was utilised to assess the basic characteristics of the proposed method. As illustrated in Fig. 4, an irregular arrangement of 4 × 4 × 4 particles was used to represent a cubic specimen with a 10 mm side. For clarity, the particles are in this example displayed as red spheres and the bonds or contacts between   Normalised displacement error (in the L 2 norm) versus mean particle spacing for the swinging-cube test them as black solid lines. The material parameters are E = 7.0GPa, ν = 0.25 and ρ = 1.0 × 10 3 kg/m 3 and the contact parameters as set to unity (i.e. ξ n = ξ t = 1). Following [55], the displacement field is postulated as where u 0 = 5 × 10 −4 . In test A, the displacement is prescribed for all particles and force equilibrium is tested for internal particles. In test B, the displacement is prescribed for all boundary particles, and the displacement error is determined for internal particles. In test C, a minimal number of essential boundary conditions needed to avoid rigid body translation or rotation are enforced and the traction corresponding to the displacement field (63) is prescribed on all boundaries. Specifically, all displacement components were prescribed on one corner particle and the vertical displacement was prescribed on two additional corner particles. Since an explicit solver was used, the displacements and/or loads were applied gradually during 5 ms for tests B and C and damping was used to obtain a static solution. The results obtained are summarised in

Swinging cube
The purpose of this test is to assess the convergence of the proposed method in a dynamic three-dimensional setting. A cube with side-length L = 1 m is considered, occupying the referential domain 0 ≤ X 1 , X 2 , X 3 ≤ 1 m. Zero normal dis- placement is enforced for X 1 = X 2 = X 3 = 0 (symmetric boundary conditions) whereas zero tangential displacement is enforced for X 1 = X 2 = X 3 = 1 (anti-symmetric boundary conditions). The following analytic solution for the displacement field u is valid [3,46] u = cos(ωt) provided that the magnitude of the oscillation is sufficiently small and that the amplitude coefficients fulfill the auxiliary condition U 1 + U 2 + U 3 = 0. Here, ω = √ 3π c s /(2L) and k = π/(2L). Specifically, we let U 1 = U 2 = 5.0 × 10 −4 m and U 3 = −1.0 × 10 −3 m and prescribe the initial displacement (and velocity) according to Eq. (64). As previously done [3], the error of the numerical solution after 2 × 10 −3 s was benchmarked against the analytical result. The material parameters are E = 17.0 MPa, ν = 0.3 and ρ = 1.1 × 10 3 kg/m 3 and the contact parameters are set to unity (ξ n = ξ t = 1). Figure 5 displays the normalised displacement error, measured in the L 2 norm, as a function of the mean particle spacing. The displacement is seen to exhibit a quadratic convergence towards the correct solution. This result is anticipated for a linearly consistent method, since the truncation error then will be of second order.

Bending and rotation of free beam
The purpose of this test is to assess the momentumpreservation properties of the proposed method. A beam of length 6 m, with a square cross section (side length 1 m) is considered, occupying the referential domain − 0.5 m ≤ X 1 , X 2 ≤ 0.5 m and 0 ≤ X 3 ≤ 6 m. The beam is free, i.e. no tractions or displacements are prescribed. The initial velocity is where set to unity (ξ n = ξ t = 1). The nonzero linear and angular momentum components (with respect to the origin) can readily be calculated analytically as p 1 = 8.8 × 10 4 kg m/s and 2 = 2.97 × 10 5 kg m 2 /s. The location and shape of the beam at certain instants of time are displayed in Fig. 6. Since the Young's modulus is relatively small, the beam is seen to undergo a pronounced bending in addition to translational and rotational motion. In Fig. 7a, the linear and angular momentum components are shown as a function of time. As a result of collocation of mass, the numerically computed momentum components differed slightly from their analytical counterparts (e.g. p 1 ≈ 8.802 × 10 3 kg m/s and 2 ≈ 2.972 × 10 5 kg m 2 /s initially), but this difference is too small to be noticeable in the figure. In Fig. 7b, the magnitude of the change in each momentum component from its initial value is displayed on a logarithmic scale. The magnitude of the change is 1 × 10 −9 kg m 2 /s for 2 , 1 × 10 −10 kg m/s for p 1 and one to two orders of magnitude smaller for the remaining components. Hence, both the linear and angular momentum can for all practical purposes be considered to be conserved. This result is anticipated, since the equations of motion have been derived from a discrete Lagrangian.

Large-amplitude vibrations of cantilever beam
In this example, a cantilever beam of length 6 m with a square cross section (side length 1 m) is considered, occupying the referential domain − 0.5 m ≤ X 2 , X 3 ≤ 0.5 m and 0 ≤ X 1 ≤ 6 m. The beam is clamped at X 1 = 0 (i.e. zero displacement is prescribed) and all other boundaries are free. The initial velocity of the beam increases with increasing distance from the clamped end, according to with v 0 = 10 m/s, so that the beam undergoes largeamplitude oscillations in the X 1 − X 3 plane. The material parameters are E = 17.0 MPa, ν = 0.3 and ρ = 1.1 × 10 3 kg/m 3 and the contact parameters are set to unity (ξ n = ξ t = 1). Figure 8 displays the coordinates of the top right corner of the free end as a function of time for four different particle densities (2/m, 4/m, 8/m and 16/m). The period of oscillation is somewhat too small for the lowest particle density, indicating a too stiff response, but the results are similar for the remaining particle densities. The deformed shape of the beam after 0.5 s, which corresponds the maximal deflection (see Fig. 8), is provided in Fig. 9. The results are virtually identical for particle densities ≥ 4/m but the deflection is somewhat smaller for the lowest particle density (2/m) and the stress is less uniform in this case. Punch displacement (mm) Vertical force (N) 3/mm 6/mm 9/mm 12/mm Fig. 11 Vertical force vs. punch displacement for four different particle densities

Punch problem
In this example, the central part of a rectangular block with dimensions 2 × 2 × 1 mm (length × width × height) is compressed by a square punch with 2/3 mm side length. Compression is done under displacement control at a rate of 0.1 mm/s for 5 s, resulting in a maximal punch displacement of 0.5 mm. The bottom face of the block is kept fixed. The material parameters are E = 240 MPa, ν = 0.4 and ρ = 1.0 × 10 3 kg/m 3 and the contact parameters are set to unity (ξ n = ξ t = 1). No damping was used. Considering the symmetry of the problem, only one quarter is considered, occupying the referential domain 0 ≤ X 1 , X 2 , X 3 ≤ 1 mm. Zero normal displacement is enforced for X 1 = X 2 = 0 (symmetric boundary conditions). A regular particle arrangement was used in order to be able to clearly see any tendency for hourglassing. The deformed shape at maximal punch displacement for four different particle densities (3/mm, 6/mm, 9/mm and 12/mm) are displayed in Fig. 10. The overall shape of the block is similar in all cases, but the significant deformation in the vicinity of the punch boundary cannot be resolved unless a relatively high particle density is used. No tendency for hourglassing can be seen, indicating that the contact forces provided sufficient stabilisation in this case. The total vertical (punch) force is displayed as a function of punch displacement in Fig. 11. Consistent with the qualitative observations made above, the results obtained indicate that a relatively large particle density (at least 10/mm) is required to obtain an adequate solution.

Conclusions
The most widely accepted Delaunay-based strain definition for granular media has been translated into a particle-based method for simulation of the mechanics of continua. The structure of the method resembles the one obtained from smoothed particle hydrodynamics or peridynamics. The determination of the deformation gradient is based on an underlying tetrahedral mesh, resulting in a formulation that is linearly consistent and satisfies the patch test. Only nearestneighbour interactions are considered, which considerably simplifies the implementation of boundaries. The interactions between particles are assumed to have two components: The first component originates from the mean stress that results from the affine (global) particle deformation (as per the standard smoothed-particle hydrodynamics). The second part emerges from non-affine (local) deformation of the contact between any two particles (in a similar manner as for the discrete element method). The particle interactions are derived from a discrete Lagrangian, thus ensuring that linear and angular momentum are conserved in the absence of external forces or torques. The performance of the method has been demonstrated in a number of numerical examples. This work can be considered a first step towards the development of a macroscopically consistent discrete method for granular materials.