# Stiffness pathologies in discrete granular systems: bifurcation, neutral equilibrium, and instability in the presence of kinematic constraints

###### Abstract

The paper develops the stiffness relationship between the movements and forces among a system of discrete interacting grains. The approach is similar to that used in structural analysis, but the stiffness matrix of granular material is inherently non-symmetric because of the geometrics of particle interactions and of the frictional behavior of the contacts. Internal geometric constraints are imposed by the particles’ shapes, in particular, by the surface curvatures of the particles at their points of contact. Moreover, the stiffness relationship is incrementally non-linear, and even small assemblies require the analysis of multiple stiffness branches, with each branch region being a pointed convex cone in displacement-space. These aspects of the particle-level stiffness relationship gives rise to three types of micro-scale failure: neutral equilibrium, bifurcation and path instability, and instability of equilibrium. These three pathologies are defined in the context of four types of displacement constraints, which can be readily analyzed with certain generalized inverses. That is, instability and non-uniqueness are investigated in the presence of kinematic constraints. Bifurcation paths can be either stable or unstable, as determined with the Hill–Bažant–Petryk criterion. Examples of simple granular systems of three, sixteen, and sixty four disks are analyzed. With each system, multiple contacts were assumed to be at the friction limit. Even with these small systems, micro-scale failure is expressed in many different forms, with some systems having hundreds of micro-scale failure modes. The examples suggest that micro-scale failure is pervasive within granular materials, with particle arrangements being in a nearly continual state of instability.

Stiffness pathologies in discrete granular systems: bifurcation, neutral equilibrium, and instability in the presence of kinematic constraints

Matthew R. Kuhn{}^{1}^{1}^{1}1Correspondence to: Donald P. Shiley School of Engineering,
University of Portland,
5000 N. Willamette Blvd.,
Portland, OR, 97203, USA.
Email: kuhn@up.edu,
Florent Prunier{}^{2}, and Ali Daouadji{}^{2}

{}^{1}Br. Godfrey Vassallo Prof. of Engrg., Donald P. Shiley School of Engrg., Univ. of Portland, 5000 N. Willamette Blvd., Portland, OR, USA 97231 {}^{2}University of Lyon, INSA-Lyon, GEOMAS, F-69621, France

Received …

KEY WORDS: instability; granular material; bifurcation; neutral stability; generalized inverse

## 1 Introduction

In the mid-2000s, several independent works were published on the nature of internal rigidity, uniqueness, and stability of discrete granular materials [Bagi:2007a, Kuhn:2005b, Nicot:2007b]. When violated, these favorable conditions give rise to failure, weakening, and localized deformation, conditions that we broadly designate as stiffness pathologies. The incremental macro-scale behavior of granular materials is known to be exceedingly complex, and the strength, stiffness, and various forms of failure (diffuse, localized, static, dynamic, etc.) at the continuum, macro-scale have received extensive investigation in the past decades. At the risk of making the nearly impenetrable and confounding behavior of these materials even more so, we return to a study of discrete failure — in its many forms — by addressing the internal particle-scale stiffness and rigidity of granular materials.

The works of Bagi [Bagi:2007a], Nicot and Darve [Nicot:2007b], and Kuhn and Chang [Kuhn:2005b], viewed granular materials at the micro-level, treating a material region, which might appear continuous at a larger scale, as a collection of discrete and notionally rigid granules that interact when touching each other at idealized contact points. The paper takes a similar primitive approach, departing from a continuum framework and treating a granular medium as a discrete system of interacting grains. As one difference between discrete and and continuous media, the discrete topology of a granular medium can be expressed as a multi-graph [Satake:1993b, Bagi:1996a] with a finite (or at most, countable) number of vertices (grains) and edges (contacts between grains), and we adopt this elemental view of a granular material. Because a discrete graph expresses the topology of a finite open set (e.g., a planar graph is homeomorphic with a sphere), a graph has no interior, and without an interior, it has no boundary, no unit normal on the boundary, no volume, and no surface area. Our vocabulary is instead of movement and force and of the relationship between them — stiffness. Movements and forces are associated with particles and contacts: the particle forces are external forces, and the contact forces are internal forces, but no distinction is made between boundary and interior forces, as there is no boundary or interior.

Another difference between discrete and continuous systems arises in the choice of a reference configuration of stable behavior. With continuous systems having simple boundaries, one can usually distinguish a “fundamental deformation” with which a buckled or bifurcated deformation can be compared (for example, the fundamental deformation of a persistently straight column or of a region that deforms in a fundamental, uniform mode without shear bands). Sliding between discrete particles at their contacts, a primary mechanism of deformation and failure, usually occurs in only a subset of the contacts, and predetermination of this sliding subset is usually quite difficult. As such, the fundamental deformation of a discrete granular system can rarely be presaged.

The paper extends the current understanding of discrete systems, by presenting a thorough accounting of four different types of geometric effects in granular media, for which three of the effects behave as internal follower forces. These geometric effects are a results of the internal geometric constraints that are imposed by the particles’ shapes, in particular, by the surface curvatures of the particles at their points of contact. The paper also describes four types of external displacement constraint, using the theory of generalized inverses to address three of the types. We also extend an accepted definition of stability by incorporating the internal geometric effects within granular systems and give a comprehensive accounting of three types of pathologic behaviors in discrete systems: neutral equilibrium, bifurcation and path instability, and instability of equilibrium, with Bažant–Petryk theory applied to the problem of path instability. This accounting is detailed for each of the four types of displacement constraints. An analysis of potential pathologies is particularly vexing with granular materials, as behavior is incrementally non-linear, and we present a systematic means of determining the consistency of a possible pathology with the assumed contact-level stiffness conditions. We also present examples of simple granular assemblies and show how the various pathologies arise in these systems.

The plan of the paper is as follows. In the next section, we develop essential elements for characterizing the stiffness of discrete systems. These elements include the derivation of the stiffness matrix of a granular assembly, including those stiffness components that depend upon the current forces and their geometric alteration (Section LABEL:sec:equilibrium). A thermodynamic framework for analyzing the stability of a granular system is then developed in Section LABEL:sec:stability. In Section LABEL:sec:constraints, we consider possible constraints placed upon a granular system by its surroundings, as these limitations will restrict the available modes of deformation and instability. Four categories of displacement constraints are developed in this section. To provide a specific framework for these principles, Section LABEL:sec:contacts describes a standard two-branch frictional model for the contact interaction between particles, and Section LABEL:sec:branches places this model in the context of an entire assembly’s stiffness matrix.

After establishing these principles, we define various stiffness pathologies in Section LABEL:sec:Pathologies (controllability, bifurcation, and instability). Examples of several simple granular systems are then analyzed in Section LABEL:sec:Examples. These examples are engaged using the methods and language of linear algebra, which we hope will clarify distinctions among the different pathologies.

With occasional exceptions, we use vector and matrix, rather than index, notation for most objects and operations. Vectors and matrices are enclosed in brackets when they contain information for an entire assembly; but brackets are usually excluded when the vector or matrix is referenced to a single contact. Vectors are written with bold lower case letters; matrices are with bold upper case letters; and scalars are normal lower case letters. Inner, outer, and dyad products are denoted as follows: inner products \mathbf{x}\cdot\mathbf{y}=x_{i}y_{i} and \mathbf{A}\cdot\mathbf{x}=A_{ij}x_{j}; outer product \mathbf{x}\times\mathbf{y}=e_{ijk}x_{j}y_{k}; and dyad product \mathbf{x}\otimes\mathbf{y}=x_{i}y_{j}.

## 2 Stiffness framework for discrete systems

Bagi [Bagi:2007a] and Kuhn and Chang [Kuhn:2005b], approached the discrete mechanics of granular materials from a structural mechanics perspective. These concurrent works developed a stiffness relationship between movement and force in matrix form, so that established concepts of stability, bifurcation, softening and controllability, already used in structural analysis, could also be applied to granular assemblies (see, for example, Bazant [Bazant:1991a]). Both works were preceded by others that developed stiffness matrices for granular systems [Serrano:1973a, Kishino:1988a, McNamara:2006a], but these earlier works neglected second-order geometric changes, which are essential to the developments that follow. A general stiffness relationship was also developed by Agnolin and Roux [Agnolin:2007a, Agnolin:2007c], who had considered geometric effects in the appendices of these works. We will use the notation of Kuhn and Chang [Kuhn:2005b] as it makes useful distinctions between objective and non-objective quantities, which become relevant when distinguishing various stiffness pathologies.

### 2.1 Notation

The shapes, positions, orientations, contact forces, and loading history of the N particles in a granular assembly are assumed known at time t, and we develop the conditions for equilibrium in a deformed state at t+dt. Each particle is assumed a hard body that interacts with neighboring particles at their shared compliant contacts. Movement and deformation are assumed slow and quasi-static, so that we can neglect viscous or gyroscopic forces. As shown in Fig. 1, a particle p touches particle q at a contact c, with q being a member of the set \mathcal{C}_{p}^{q} of p’s neighboring particles, q\in\mathcal{C}_{p}^{q}, and with c being a member of the set \mathcal{C}_{p}^{c} of p’s contacts, c\in\mathcal{C}_{p}^{c}.

In this work, an assembly’s contact topology is represented as a directed multi-graph (a multi-digraph). The graph is “multi”, because a pair of particles can share multiple contacts, as can occur with non-convex particles. The graph is “directed”, because we treat the pq contact (from p to q) as distinct from the contact qp (from q to p). This approach allows a more straightforward derivation and is consistent with the non-symmetry of the stiffness matrix, as we will see that the pq stiffness can differ from the qp stiffness. The system can also include isolated “rattler” particles that are not in contact with other particles. With this situation, we can adopt two different approaches. We can include these rattlers as isolated nodes of a complete topology, which will lead to numerous zero-stiffness modes of neutral equilibrium, or we can consider only the load-bearing network of contacting (non-rattler) particles and then update the topology whenever rattlers are freshly released from or incorporated into the network. For reasons given in Section LABEL:sec:equilibrium, we take the latter approach and consider only the current load-bearing network of contacts and particles.

Vector \mathbf{u}^{p} is the location of a material reference point \chi^{p} attached to p; \boldsymbol{\theta}^{p} are the orientation cosines of p; and \mathbf{b}^{p} and \mathbf{w}^{p} are the external force and moment that act upon p at the point \chi^{p}.

The entire particle system has M contacts. The single contact vector \mathbf{r}^{c,p}, is from \chi^{p} to its contact c; \mathbf{n}^{c} is the outward unit vector normal to the surface of p at contact c; and \mathbf{f}^{c} and \mathbf{m}^{c} are the contact force and moment exerted upon p at c.

Although the two particles, p and q, can share multiple contacts c, we will use the more convenient notations \mathbf{r}^{pq,p}, \mathbf{n}^{pq}, \mathbf{f}^{pq}, and \mathbf{m}^{pq} with the understanding that the pair pq can represent one of several contacts c between p and q. With this notation, force \mathbf{f}^{pq} and moment \mathbf{m}^{pq} act upon particle p, and \mathbf{n}^{pq} is directed outward from p. We must distinguish, however the two “\mathbf{r}^{pq}” contact vectors for contact pq: vector \mathbf{r}^{pq,p} is directed from p to the contact; whereas vector \mathbf{r}^{pq,q} is directed from q to the contact.

We gather the particle positions and orientations of all N particles into a stacked column vector [\mathbf{u}/\boldsymbol{\theta}], the external forces and moments into the stacked complementary vector [\mathbf{b}/\mathbf{w}], and the internal contact forces and moments into the stacked vector [\mathbf{f}/\mathbf{m}]:

\left.\begin{array}[]{r}\mathbf{u}^{p}\\ \boldsymbol{\theta}^{p}\cr\@@LTX@noalign{ }\omit\\ \end{array} |