A contact dynamics approach to the Granular Element Method

Highlights

• We develop a new contact dynamics approach to the Granular Element Method (CD–GEM).

• CD–GEM combines shape flexibility using NURBS with implicit time discretization.

• We develop a novel separating hyperplane algorithm to simplify contact detection.

• CD–GEM performs robustly with rigid particles in the quasi-static and flow regimes.

• CD–GEM significantly simplifies current CD approaches while maintaining performance.

Abstract

We present a contact dynamics (CD) approach to the Granular Element Method (GEM) Andrade et al. (2012) [1], abbreviated here as CD–GEM. By combining particle shape flexibility through Non-Uniform Rational Basis Splines, properties of implicit time-integration discretization (e.g., larger time steps) and non-penetrating constraints, as well as a reduction to a static formulation in the limit of an infinite time step, CD–GEM targets system properties and deformation regimes in which the classical discrete element method either performs poorly or simply fails; namely, in granular systems comprising of rigid or highly stiff angular particles and subjected to quasi-static or intense dynamic flow conditions. The integration of CD and GEM is made possible while significantly simplifying implementation and maintaining comparable performance with existing CD approaches.

Introduction

The objective of this paper is to develop a contact dynamics (CD) approach to the Granular Element Method (GEM) [1], abbreviated here as CD–GEM. CD–GEM targets system properties and deformation regimes in which the classical discrete element method (DEM) either performs poorly or simply fails; namely, in granular systems comprising of rigid or highly stiff angular particles and subjected to quasi-static or intense dynamic flow conditions. Within the context of such applications, we describe how CD–GEM offers a better solution in terms of particle morphology or shape representation and ease of implementation, while maintaining comparable performance with existing CD approaches. To motivate the development of our approach, we first refer the reader to Table 1 for a brief summary of the key features of and differences between CD and the classical DEM by Cundall and Strack [2]. In the following, we highlight the difficulties associated with CD and DEM followed by a description on how we eliminate them through CD–GEM.

The so-called Non-Smooth CD, originally developed by Moreau [6], [7], [8], [9], is an alternative discrete approach to the DEM. The most prominent feature of CD, in contrast to that of classical DEM, is that the particles are considered perfectly rigid and the contact forces are determined as those that prevent interparticle penetration and at the same time satisfy the frictional stick–slip constraints. In their simplest forms, these contact laws are embodied in the so-called Signorini unilateral contact condition and classical Coulomb law, as shown in Fig. 1(a) and (b), respectively. Commensurate with these physical enhancements, however, is the need for both contact and constraint forces to be solved simultaneously or implicitly since the problem is nonlinear. The need for an implicit solution procedure till today remains the primary reason why CD is deemed much more complicated to implement than DEM. This has thwarted the wide adoption of CD despite the favorable performance that has been shown through a number of studies [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20].

While there is wide applicability of DEM, its application has gone beyond its restriction as a tool that is strictly applicable only to materials with finite elasticity. For example, DEM is widely used as a tool to study real granular materials that are almost rigid or highly stiff in nature. Here, finite elasticity means that the contact interaction is essentially modeled using springs. Under explicit time integration algorithms that are typically used in DEM, the stable time step is restricted by the critical time step, which scales with the inverse of the contact spring-particle mass frequency. This results in infinitesimally small time steps if material parameters corresponding to highly stiff particles (e.g., rocks, sand, steel) are used. Although explicit integration algorithms can be easily parallelized, the runtime for stiff systems remains computationally prohibitive. One modeling technique commonly employed in practice to overcome this restriction is to simply reduce the contact stiffness, usually by two to four orders of magnitude, to the extent that particle kinematics obtained from simulations are still somewhat representative of the overall response of the actual system of interest. If quasi-static behavior is assumed to hold, usually used in combination with stiffness tuning is mass scaling in which the particle masses are adjusted (usually increased) such that the combined spring-particle system frequency is lowered, increasing the time step size. In practice, model calibration by means of mass scaling and/or stiffness tuning is a delicate and cumbersome process. Another problem that is associated with the presence of contact springs and the dynamic nature of DEM is the introduction of unwanted oscillations or noise, with frequencies that increase with spring stiffness. This requires additional calibration of the global and/or local damping parameters. Moreover, under certain loading conditions (for e.g., strain-controlled and dynamic), either particle kinematics or contact forces obtained under such calibration procedures can be highly inaccurate [21].

Recent CD and DEM approaches include techniques to represent complex particle morphology or shape. In this aspect, recent trends show a clear dichotomy between the choice of shape representation technique. In DEM, the disc- or sphere-clustering technique (see for e.g., [22], [23]) appears to have become fairly widespread because of its simplicity. The clustering technique approximates the particle shape by overlapping discs or spheres and the same disc-disc or sphere-sphere contact algorithm is reused over all potentially contacting pairs. In CD, the polyhedra approach has been become quite popular [24]. A key component in the CD formulation is the signed separation or gap, which is used in the determination of constraint forces to prevent particle interpenetration. Algorithms for the determination of signed separation operate under the assumption of the existence of a hyperplane separating two disjoint convex sets [25]. One such algorithm is the common plane approach [26] used in early DEM for convex polyhedra. The algorithm works by finding an initial separating plane that is approximately equidistant from each polyhedron of a potentially contacting pair. The orientation and position of the separating plane are then updated in subsequent time steps using an iterative technique so that the equidistant relationship is maintained [27]. More efficient and robust solution algorithms, however, have since been identified (see for e.g., [28], [29], [30]). While the polyhedra approach is considered as a more accurate shape representation technique than the clustering approach, its contact implementation remains complicated due to the need to enumerate all the various combinations of contact entities (node, edge, surface), as well as the complexities associated with the algorithms to estimate the penetration depth. As such, the simpler disc/sphere-clustering is favored over the more accurate polyhedra-based approach.

In the design of CD–GEM, we combine and refine two important developments that allow us to eliminate all the above difficulties:

• We simplify the formulation and implementation of CD significantly by generalizing a variational CD formulation recently developed for discs and spheres [3], [31], [32]. This particular formulation, which is employed in this paper, is appealing because it provides a way for CD to be easily implemented and solved using off-the-shelf mathematical programming solvers. The most prominent advantage of this formulation is its automatic inclusion of the quasi-static limit, enabling quasi-static modeling without the need for adjusting damping parameters or time step.

• We remove the complexities associated with polyhedra-based contact detection algorithms by adopting Non-Uniform Rational Basis Splines (NURBS) to describe arbitrary particle geometries, as used in GEM [1]. In addition to providing a significant enhancement in the representation of particle morphology (namely, sphericity and angularity [33]), NURBS also facilitates contact calculations through its smooth boundaries. The original contact algorithm as described in [1], however, is based on an explicit calculation of intersection between two arbitrary NURBS curves, which does not provide a signed separation and cannot be used in a CD setting. To overcome this problem, we develop a supporting separating hyperplane (SSH) contact algorithm for particles described using NURBS, which provides the required signed separation or gap required in the CD formulation. The proposed algorithm determines, in a consistent manner and simultaneously, the signed separation and contact normal. To further simplify its implementation, we cast the contact algorithm as a constrained minimization problem, which can then be solved using standard constrained optimization techniques. With the SSH contact algorithm in place, the integration of GEM into the CD formulation is shown to be simple and straightforward.

This paper describes the details on how each of the above items are implemented and is structured as follows. In Section 2, we present the CD formulation for frictionless particles with only translational degrees-of-freedom followed by an extension to frictional particles including rotational degrees-of-freedom in Section 3. The formulation will show that the signed separation function provides the way through which GEM can be incorporated into CD. In Section 4, we describe the implementation of the proposed SSH contact procedure for particles described using NURBS. Finally, in Section 5, we present numerical examples to demonstrate the capabilities of CD–GEM before conclusions are drawn in Section 6. We limit our discussion to two-dimensional particles of strictly convex shapes in which the SSH theorem holds, and the contact algorithm can be concisely described and implemented.

Remark 1

An extended variational version of CD, that includes (as a regularization technique) particle elasticity, is presented in [3]. We focus our attention here on the case of perfectly rigid particles, but we note that the contact algorithm described in this paper can be directly incorporated in the extended case without any change.

Remark 2

There is another method, also called the Granular Element Method, which was proposed by Kishino [34]. The method, however, applies only to disks and spheres, and its solution procedure is different than the procedures described in this paper, and in [1], [31].