Dense granular flow at the critical state: maximum entropy and topological disorder

TL;DR

This paper models dense granular flows at the critical state using maximum entropy principles to predict topological distributions, aligning well with simulation data.

Contribution

It introduces a maximum entropy framework for predicting topological measures in dense granular flows, incorporating both topological and geometric considerations.

Findings

Predicted void valence distribution matches DEM simulations.

Coordination number distribution aligns with observed data.

Topological entropy maximization effectively models steady-state granular arrangements.

Abstract

After extensive quasi-static shearing, dense dry granular flows attain a steady-state condition of porosity and deviatoric stress, even as particles are continually rearranged. The paper considers two-dimensional flow and derives the probability distributions of two topological measures of particle arrangement---coordination number and void valence---that maximize topological entropy. By only considering topological dispersion, the method closely predicts the distribution of void valences, as measured in discrete element (DEM) simulations. Distributions of coordination number are also derived by considering packings that are geometrically and kinetically consistent with the particle sizes and friction coefficient. A cross-entropy principle results in a distribution of coordination numbers that closely fits DEM simulations.

Tables

Table 1. Table 1: Duality of the valence and coordination number problems

Void valence	Particle coordination number
particles, contacts, voids	voids, contacts, particles
$\left\{s_i\right\}=\left\{(l_i,\Delta M_i)\right\}$	$\left\{s_i^{\prime }\right\}=\left\{(n_i,\Delta M_i)\right\}$
$l_i\in \{3,4,\dots \}$	$n_i\in \{2,3,\dots ,n_{\text{max}}\}$
$\Delta M_i\in \{1,2,\dots ,l_i-1\}$	$\Delta M_i\in \{1,2,\dots ,n_i-1\}$
$P_{l,\Delta M}^{\text{I}}$	$P_{n,\Delta M}^{\text{I}}$
Eq. (3)	$⟨\Delta M⟩=\overline{n}/2$
Eq. (4)	$⟨n⟩=\overline{n}$

Table 2. Table 2: Distances between DEM results and estimates of the coordination number distribution (Eq.  21 )

Estimate	Distance
Model I	0.268
Model II	0.077
Model III	0.019

Equations

$$\overline{n}=2M/N,\overline{l}=2M/L=2\overline{n}/(\overline{n}-2)$$

$$M=\sum _{i=1}^L\Delta M_i,N=1+\sum _{i=1}^L\Delta N_i=1+\sum _{i=1}^L(\Delta M_i-1)$$

$$⟨\Delta M_i⟩\equiv \frac{1}{L}\sum _{i=1}^L\Delta M_i=\frac{\overline{n}}{\overline{n}-2}$$

$$⟨l_i⟩\equiv \frac{1}{L}\sum _{i=1}^L{l}_i=\overline{l}=\frac{2\overline{n}}{\overline{n}-2}$$

$$⟨\Delta M⟩\equiv \sum _{l=3}^{\infty }\sum _{\Delta M=1}^{l-1}\Delta M{P}_{l,\Delta M}^{\text{I}}=\frac{\overline{n}}{\overline{n}-2}$$

$$⟨l⟩\equiv \sum _{l=3}^{\infty }\sum _{\Delta M=1}^{l-1}l{P}_{l,\Delta M}^{\text{I}}=\frac{2\overline{n}}{\overline{n}-2}$$

$$H_l^{\text{I}}=-\sum _{l=3}^{\infty }\sum _{\Delta M=1}^{l-1}{P}_{l,\Delta M}^{\text{I}}\ln \left(P_{l,\Delta M}^{\text{I}}\right)$$

$$P_{l,\Delta M}^{\text{I}}=\frac{1}{Z_l^{\text{I}}\left({\lambda }_{l,1},{\lambda }_{l,2}\right)}\exp \left(-{\lambda }_{l,1}l-{\lambda }_{l,2}\Delta M\right)$$

$$Z_l^{\text{I}}\left({\lambda }_{l,1},{\lambda }_{l,2}\right)=\sum _{l=3}^{\infty }\sum _{\Delta M=1}^{l-1}\exp \left(-{\lambda }_{l,1}l-{\lambda }_{l,2}\Delta M\right)$$

$$P_l^{\text{I}}=\sum _{\Delta M=1}^{l-1}{P}_{l,\Delta M}^{\text{I}}$$

$${\alpha }_{\mathsf{S}}{P}_{D_0}\prod _{k\in \mathsf{S}}{P}_{D_k}$$

$$Q_n=\sum _{\ell =1}^r{P}_{D_{0,\ell }}\sum _{j=1}^{r^n}{\alpha }_{{\mathsf{S}}_j}\prod _{k\in {\mathsf{S}}_j}{P}_{D_k}$$

$$P_n^{\text{II}}=Q_n/\sum _{n=2,3,\dots }{Q}_n$$

$$H_n^{\text{III}}(\vec{P}:\vec{q})=\sum _{n=2}^{\infty }\sum _{\Delta M=1}^{n-1}{P}_{n,\Delta M}^{\text{III}}\ln \left(\frac{P_{n,\Delta M}^{\text{III}}}{q_n}\right)$$

$$q_n=P_n^{\text{II}}$$

$$\sum _{n=2}^{\infty }\sum _{\Delta M=1}^{n-1}n{P}_{n,\Delta M}^{\text{III}}=\overline{n}$$

$$\sum _{n=2}^{\infty }\sum _{\Delta M=1}^{n-1}\Delta M{P}_{n,\Delta M}^{\text{III}}=\frac{\overline{n}}{2}$$

$$P_{n,\Delta M}^{\text{III}}=\frac{q_n}{Z_n^{\text{III}}\left({\lambda }_{n,1},{\lambda }_{n,2}\right)}\exp \left(-{\lambda }_{n,1}n-{\lambda }_{n,2}\Delta M\right)$$

$$Z_n^{\text{III}}\left({\lambda }_{n,1},{\lambda }_{n,2}\right)=\sum _{n=2}^{\infty }\sum _{\Delta M=1}^{n-1}{q}_n\exp \left(-{\lambda }_{n,1}n-{\lambda }_{n,2}\Delta M\right)$$

$$P_n^{\text{III}}=\sum _{\Delta M=1}^{n-1}{P}_{n,\Delta M}^{\text{III}}$$

$$\sum _{n=2,3,\dots }{P}_n^{(\bullet )}\ln \frac{P_n^{(\bullet )}}{P_n^{\text{DEM}}}$$

Full text

1 Introduction

The critical state principle, a unifying concept in geomechanics, holds that dense confined granular materials attain a steady-state condition of flow after extensive slow shearing Schofield:1968a . A distinctive feature of such flow is the continual micro-scale rearrangement of particles even as bulk characteristics — stress, density, and fabric — remain nearly constant. Micro-scale rearrangements are expressed in three ways: (a) statically (kinetically), as alterations of inter-particle contact forces, (b) geometrically, as changes in the particles’ positions or in the local density, and (c) topologically, as changes in the load-bearing contact network among the particles. This Paper addresses the latter form of granular arrangement as expressed in a two-dimensional (2D) setting, focusing on the local coordination number (number of contacts per particle) and the local void valence (number of particles surrounding a void). The analysis applies to the critical state flow of frictional materials of sufficient density to develop a load-bearing network of contacts during slow (quasi-static) shearing.

The discrete 2D topological arrangement of load-bearing particles is conveniently represented with a planar particle graph, a tessellation in which edges, nodes, and faces represent inter-particle contacts, particles, and voids, respectively Satake:1992a . Fig. 1 shows particle graphs from discrete element (DEM) simulations Cundall:1979a of a two-dimensional assembly of circular disks at two states: an initial state of 676 densely arranged disks and a deformed state in which the assembly width has been slowly reduced by a horizontal strain $\varepsilon_{11}=-20\%$ while maintaining the original, constant mean stress  $p=\frac{1}{2}(\sigma_{11}+\sigma_{22})$.

Such biaxial loading departs from the initial isotropic condition, producing deviatoric stress $\sigma_{11}-\sigma_{22}$ and inducing steady-state (critical state) flow in which continued, plastic deformation progresses at constant volume, mean stress, and deviatoric stress Schofield:1968a . Gross topological parameters, such as the average coordination number and average void valence (i.e., the number of sides of an $l$-polygon void), also remain constant during such flow Thornton:2000a ; Pena:2009a . The micro-scale topology, however, is continually altered, to the extent that a further small deformation of the assembly in Fig. 1(b) of as little as $\Delta\varepsilon_{11}=2\%$ will produce a particle graph that is scarcely recognizable from the one shown. During deformation, an existing polygonal void can split into two smaller voids when fresh contacts are newly established among unloaded, rattler particles within the original void (i.e., particles with zero or one contact) Kruyt:2012a . Likewise, adjacent voids can merge when contacts disengage along shared edges. These transmutations, similar to those in hexagonal cellular and froth structures Weaire:1984a , produce a disordered and constantly changing (and seemingly random) topology.

We view this continual transmutation of particle arrangement as a maximally disordered process, in which void polygons split and merge in a random, disordered manner, subject only to bulk topological constraints. The principle of maximum entropy has been applied to granular materials for several decades. In these studies, disorder is usually expressed as the Shannon entropy Shannon:1948a ; BenNaim:2008a . Brown et al. Brown:2000b conducted experiments on two-dimensional assemblies of spheres, and by applying a back-and-forth shearing, disorder in the local density and coordination number increased with each shearing cycle. The Jaynes Jaynes:1957a maximum entropy (MaxEnt) formalism has typically been used to solve the condition of maximum disorder. This approach has been applied to the local fabric of disk assemblies by categorizing voids into several canonical types Shahinpoor:1983b ; Brown:2000a . Similar maximum entropy approaches have also been applied to the local packing density Moroto:1983a ; Edwards:1989a ; Kumar:2005a ; Yoon:2012a , to contact forces Coppersmith:1996a ; Edwards:2001a ; Chakraborty:2010a , to contact orientations Troadec:2002a , and to contact displacements and bulk elastic moduli Rothenburg:2009a .

The Euler equation, $M=N+L-1$, applies to the bulk topological quantities of a particle graph: $M$ is the number of contacts (edges), $N$ is the number of particles (nodes), and $L$ is the number of void polygons (faces). The equation holds for any connected planar graph, but we restrict attention to the load-bearing subgraph of an assembly, which excludes peninsular and island particles (rattlers), such as those apparent in Fig. 1(c). (Such particles, shown in red, are typically “nudged” along by their neighbors until the surrounding void eventually collapse onto them.) For large assemblies, the average coordination number $\overline{n}$ among load-bearing particles and the average void valence $\overline{l}$ are defined as

[TABLE]

which will serve as the constraints on possible particle graphs. The corresponding micro-scale quantities are the number of contacts $n_{j}$ of the “$j$”th load-bearing particle and the valence $l_{i}$ of the “$i$”th $l$-polygon void.

We develop three methods for estimating the probability distributions of $n_{j}$ and $l_{i}$. In the next section, a maximum entropy principle (MaxEnt) is used to maximize the disorder of the two distributions. Although this approach leads to a reasonable approximation of the void valence $l_{i}$ distribution, the estimated distribution of coordination numbers $n_{j}$ misses many features observed in discrete element (DEM) simulations. In Section 3, we invoke certain geometric and kinetic constraints to develop a second distribution of coordination numbers. This distribution is an improvement over the first, and it also resolves the effect of inter-granular friction on the coordination number. In Section 4, the two methods are combined with a minimum cross-entropy (MinXEnt) principle, which yields an improved estimate of the coordination number distribution.

2 Maximum disorder theory (Model I)

We first develop a model, based upon a maximum disorder principle, for estimating the probabilities $P_{l}^{\text{I}}$ and $P_{n}^{\text{I}}$ of encountering a void with valence $l$ and a particle with coordination number $n$. The superscript “I” denotes probabilities derived from this first approach. We begin with a method of constructing a particle graph from a journal $\{s_{i}\}$ of integer pairs $s_{i}=(l_{i},\Delta M_{i})$, a journal that represents a single micro-state of topology. Each journal entry $s_{i}$ corresponds to a single “$i$th” void polygon (face) of valence $l_{i}$ that has been appended to the particle graph by adding $\Delta M_{i}$ new contacts (edges). After describing the construction of a journal, we then establish constraints on such journals. Assuming that all journals that meet these constraints are equiprobable during steady-state flow and that the most likely topology maximizes the disorder of the system, we derive the expected probabilities of pairs $(l,\Delta M)$ and, by extension, the probabilities $P_{l}^{\text{I}}$ of valences $l$. A duality principle leads to the complementary probabilities of coordination numbers, $P_{n}^{\text{I}}$. These expected probabilities, $P_{n}^{\text{I}}$ and $P_{l}^{\text{I}}$, are then compared with those measured in DEM simulations.

The disorder of a particle graph is characterized by the number of ways in which similar, equiprobable graphs can be constructed from sets of polygonal faces. In this regard, we start with a scheme for identifying and “counting” these graphs (i.e., micro-states), each micro-state being a nearly random sequence $\{s_{i}\}$ of pairs $s_{i}=(l_{i},\Delta M_{i})$.

We now describe the meaning of such journals (micro-states), referring to Fig. 2a and the construction of this seven-polygon graph. Starting from a seed polygon and one of its vertices (labeled \raisebox{-.7pt} {\small1}⃝ and “1” in Fig. 2), a planar graph can be constructed by progressively appending faces to the graph. Each new face of valence $l_{i}$ adds $\Delta M_{i}$ edges and $\Delta N_{i}$ vertices ($\Delta N_{i}=\Delta M_{i}-1$). A journal $\left\{s_{i}\right\}$ of integer pairs $s_{i}=(l_{i},\Delta M_{i})$ records this process (Fig. 2a).

Each new face \raisebox{-.7pt} {\smalli}⃝ begins from the terminal vertex “$i-1$” of the previous face. In this manner, new faces, having pairs $(l_{i},\Delta M_{i})$, are spiraled counter-clockwise around the perimeter of the existing graph, with each \raisebox{-.7pt} {\smalli}⃝th face sharing $l_{i}-\Delta M_{i}$ edges with the existing graph. The curved arrows in Fig. 2a show that each added face, beginning from the previous $i-1$ terminus, includes the edge of the existing perimeter that lies in front of (i.e. counter-clockwise relative to) the previous terminus, and the new face also includes the edge immediately clockwise around the terminus node. For example, face \raisebox{-.7pt} {\small3}⃝, rather than the later face \raisebox{-.7pt} {\small7}⃝, is added from terminus 2, even though both faces share this node. Several faces can start from the same node, whereas some nodes are not the start of any face.

A unique graph should result from any given journal, and a unique journal should be associated with any given graph in which the seed has been specified. To achieve these qualities, we require three additional conditions with each added face.

No edge of a new face \raisebox{-.7pt} {\smalli}⃝ should lie behind (i.e. clockwise relative to) the $i-1$ terminus along the existing perimeter: all edges of \raisebox{-.7pt} {\smalli}⃝ should be new edges or should be existing edges that lie in front of the $i-1$ node. With some nodes, a pseudo-face with two edges can be used to skip across an existing perimeter edge. For example, starting at node 3, the face \raisebox{-.7pt} {\small9}⃝ could not be added as the fourth face, since this face would include an edge (labeled 3–8) that lies behind node 3. The pseudo-face \raisebox{-.7pt} {\small4}⃝ skips to the next node, labeled 4, from which face \raisebox{-.7pt} {\small5}⃝ emanates. 2. 2.

When considering the ordered list of all nodes of a new face \raisebox{-.7pt} {\smalli}⃝, beginning with the previous $i-1$ terminus and moving clockwise around the face, these nodes should consist of existing (perimeter) nodes followed by added nodes. A counter-example is shown in Fig. 2b, where the six nodes (solid dots) of an improper face \raisebox{-.7pt} {\smalli}⃝ (crossed) have the pattern T-T-F-T-T-T ($\text{T}=\text{ }$existing, $\text{F}=\text{ }$new). If this face in Fig. 2b was added to the graph, the interior triangular face would be stranded and could not later be added to the journal. A proper pseudo-face \raisebox{-.7pt} {\smalli}⃝ (not crossed) must be used to skip to its terminus, after which the triangular face $i+1$ would be added. The improper face \raisebox{-.7pt} {\smalli}⃝ (crossed) could be added later, after the process has spiraled fully around the graph once again. 3. 3.

When considering the ordered list of nodes in the existing perimeter, beginning with the previous terminus $i-1$ and moving counter-clockwise around the perimeter, these nodes should consist of nodes that are part of the new face \raisebox{-.7pt} {\smalli}⃝ followed by nodes that are not in \raisebox{-.7pt} {\smalli}⃝. A counter-example is shown in Fig. 2c, with an improper face \raisebox{-.7pt} {\smalli}⃝ (crossed). Beginning with the node $i-1$, we have the pattern T-T-F-F-T-T-T-T-$\ldots$ along the existing perimeter ($\text{T}=\text{ }$included in the nodes of \raisebox{-.7pt} {\smalli}⃝, $\text{F}=\text{ }$not included in the nodes of \raisebox{-.7pt} {\smalli}⃝), where the two bottom-most nodes are not part of the improper face \raisebox{-.7pt} {\smalli}⃝ (crossed). To avoid stranding the bottom quadrilateral face, a pseudo-face \raisebox{-.7pt} {\smalli}⃝ skips ahead so that the proper bottom face $i+1$ (circled) is added to the graph.

The journal of a given graph can be efficiently constructed when the graph is represented with the data structure of a doubly-connected edge list (DCEL) Preparata:1985a , which permits rapid counter-clockwise traversals of the faces and nodes.

Every journal $\left\{s_{i}\right\}$ corresponds to a unique planar graph (a topological micro-state), although the same graph can be constructed from different journals, depending upon the chosen seed. The total number of faces, $L$, in the load-bearing graph is equal to the numbers of pairs in the journal, excluding those of valence 2. The set of all possible journals of length $L$ comprises a topological configuration space. We require, however, that a journal $\left\{s_{i}\right\}$ be consistent with the known coordination number $\overline{n}$ and valence $\overline{l}$. The total number of contacts (edges) and particles (vertices) are

[TABLE]

excluding faces of valence 1 or 2. Substituting Eqs. (11) and (21) into the Euler equation, dividing by $L$, and neglecting the insignificant $1/L$ yields a constraint on the journal $\left\{s_{i}\right\}$,

[TABLE]

The journal must also be consistent with the average valence $\overline{l}$ [Eq. (12)]:

[TABLE]

Elements in a journal $\left\{s_{i}\right\}$ correspond to the faces (voids) of a particle graph. The journals of small assemblies can be affected by the choice of seed, a matter addressed below. In the limit of large assemblies, we can assign discrete probabilities $P_{l,\Delta M}^{\text{I}}$ to the likelihood of each “$l$ species” and its “$\Delta M$ sub-species” among the $L$ faces. The random variables $l$ and $\Delta M$ are assumed to be independent. Eq. (3) implies a constraint on these probabilities:

[TABLE]

The outer sum explicitly excludes valences 1 and 2, as these values apply only to non-load-bearing particles, such as peninsular particles (rattlers). The inner sum ignores pairs with $\Delta M=l$, as this situation arises only with the initial, seed face. A second constraint results from the average valence $\overline{l}$ of an assembly (Eq. (4)):

[TABLE]

The full mechanical description of an assembly’s incremental evolution requires solution of an $N$-body problem, involving $3N$ equilibrium equations, while applying $2M$ contact elastic-frictional rules and any applicable boundary constraints Kuhn:2005b ; Agnolin:2007c . Although this non-linear problem is solvable, its solution also requires full information about the initial particle positions and contact forces. Lacking such information, we simply estimate the probabilities $P_{l,\Delta M}^{\text{I}}$ that describe the most likely topological condition. The most likely probability set — the set that optimally respects the missing information and, hence, corresponds to the greatest number of similar journals within the configuration space — is the one which maximizes the topological disorder (Shannon entropy) Jaynes:1957a ; BenNaim:2008a ,

[TABLE]

such that $\partial H_{l}^{\text{I}}/\partial P_{l,\Delta M}^{\text{I}}=0$, while satisfying the moment constraints of Eqs. (5) and (6). By applying the Jaynes formalism Jaynes:1957a in maximizing $H_{l}^{\text{I}}$ (i.e., using the maximum entropy “MaxEnt” principle), the probabilities $P_{l,\Delta M}^{\text{I}}$ and partition function $Z_{l}^{\text{I}}$ are

[TABLE]

with two Lagrange multipliers, $\lambda_{l,1}$ and $\lambda_{l,2}$, that satisfy Eqs. (5) and (6). The probability $P_{l}^{\text{I}}$ of encountering a void of valence $l$ is the exponential marginal probability

[TABLE]

The probabilities of coordination numbers $n$ are developed with the dual of the particle graph (the void graph), in which particles and voids assume the complementary roles of faces and vertices Satake:1992a . The alternative journal $\{s_{j}^{\prime}\}=\{(n_{j},\Delta M_{j})\}$ is constructed by adding particles (now faces in the dual, void graph) of coordination number $n_{j}$ around a seed void (now a vertex). In this manner, we can develop complementary statistics on $n$ (see Table 1).

This approach differs in three respects from the previous analysis of void valence $l$. First, the minimum coordination number is 2, with $n_{i}\in\{2,3,\ldots\}$, since a particle can have two contacts, but particles with fewer than two contacts cannot serve within a load-bearing contact network. Second, among the faces (particles) with $n_{i}=2$, some are truly particles with two contacts; others can be pseudo-faces. A list must be maintained to distinguish between these two cases. Finally, the maximum coordination number $n_{\text{max}}$ is limited by geometric exclusion and will depend upon the range of particle sizes, a matter more fully developed in Section 3.

2.1 Simulation methods and results — Method I

Statistics of valence and coordination number were measured in assemblies of 676 bi-disperse disks that were sheared in biaxial compression (Fig. 1). The two disk varieties have ratios of 1.5:1 in size, 1:2.25 in number, and 1:1 in cumulative area (that is, 468 particles of size $1.0$ and 208 of size $1.5$). The assemblies were small enough to prevent gross non-homogeneity in the form shear bands, yet large enough to capture the average, bulk material behavior. To develop more robust statistics, 168 different assemblies were created by compacting random sparse mixtures of the two disk sizes into dense isotropic packings within periodic boundaries until the average contact indentation was $0.0002$ times the average radius. Linear contact stiffnesses were applied between particles with equal tangential and normal coefficients ($k^{\text{t}}=k^{\text{n}}$), and the friction coefficient $\mu=0.50$ was enforced during the pair-wise particle interactions daCruz:2005a . Using the discrete element (DEM) algorithm, the initially square assemblies were horizontally compressed in increments $\Delta\epsilon_{11}=1\times 10^{-6}$ while maintaining a constant mean stress of $2\times 10^{-4}k^{\text{n}}$. Stress and volumetric behavior are shown in Fig. 3. A large initial stiffness cause the deviatoric stress to rise quickly from zero to a peak stress at strain $-\varepsilon_{11}=2\%$. The critical state condition was attained at compressive strains $-\varepsilon_{11}$ of 16–18%. During subsequent steady-state deformation, the particle graph of each assembly was interrogated at five strains between $-\varepsilon_{11}=16\%$ and $25\%$ (Fig. 3).

Applying the ergodicity principle, micro-sate statistics were averaged across the five strains and the 168 assemblies, involving 840 particle graphs containing over 300,000 void faces.

At the critical state, the average valence $\overline{l}=5.24$, and the average coordination number $\overline{n}=3.23$ among load-bearing particles (i.e., excluding rattlers). About 11% of particles were non-load-bearing rattlers. The probabilities $P_{l,\Delta M}^{\text{I}}$ are nearly the same across all assemblies and at all strains, although averages $\overline{l}$ and $\overline{n}$ will depend upon the particular material properties (for example, a larger coefficient $\mu$ will increase $\overline{l}$ Kruyt:2013a ).

Predictions of topological statistics are based upon Eqs. (5)–(10) and Table 1, using the measured averages $\overline{l}$ and $\overline{n}$. An upper bound of 14 was applied to $l$, as larger voids were not observed. Micro-state (journal) statistics are shown in Fig. 4 for two void valences, $l=5$ and 7. This figure was produced by reconstructing the journals of 840 particle graphs.

Eqs. (5)–(9) predict varying probabilities $P_{l,\Delta M}^{\text{I}}$ across valences $l$, but they predict a uniform probability for each $\Delta M$ within a single valence $l$. The latter prediction is not realized in the simulations (Fig. 4), which exhibit a variation of $P_{l,\Delta M}^{\text{I}}$ with $\Delta M$. Additions $\Delta M$ of one and $l-1$ are more frequent than those with $\Delta M$ near $l/2$. This result suggests that a somewhat greater order is realized in the simulations than predicted by the theory (i.e., the assumption that $l$ and $\Delta M$ are independent is not supported by the data). The greater topological order is reflected in a measured entropy $H_{l}^{\text{I}}$ of 3.191 compared with a prediction of 3.306 (Eq. 7).

The box-plot within Fig. 4 also shows the effect of the choice of a seed on the statistics of voids having five edges ($l=5$). A single graph was chosen among the 840 graphs, and fifty seeds were applied. Although the choice of seed does affect the distribution of $\Delta M$ values, the effect (expressed as a standard deviation) is less than 9% of the mean probabilities, and the resulting entropy $H_{l}^{\text{I}}$ is affected by less than $0.5\%$.

Fig. 5a shows the predicted and actual probabilities $P_{l}^{\text{I}}$ of void valence, based upon the average $\overline{l}$ (Eq. 10). The valence distributions are in general agreement, although the theory predicts a greater frequency of triangular voids and of voids having valence greater than 9, whereas the theory slightly under-predicts the frequencies of valences 4 through 8. Experiments with assemblies having different friction coefficients $\mu=0.1$–$0.7$ show similar features as Fig. 5a, although $\overline{l}$ increases with $\mu$, and the histograms, both measured and predicted, broaden to the right.

A similar analysis of coordination number is based upon Table 1 and the observed average $\overline{n}$, with an upper bound of 7 on $n$, since no more than seven disks of radius 1.0 can touch a central disk of radius 1.5 (Section 3). The comparison of theory and simulations (Fig. 5b) is less supporting than that of the valences. Coordination numbers of 2, 5, 6 and 7 are much less frequent than the predictions; coordination numbers of 3 and 4 are under-predicted in their frequency.

2.2 Discussion of Model I results

The theory captures features of the void valence probabilities, but gives a rather poorer prediction of the distribution of coordination numbers. Discrepancies result from a theory that addresses topological disorder alone, but is uninformed by the geometric or kinetic aspects of granular flow. That is, the theory is unburdened by the reality that each vertex or face is a real particle or void, having geometric character (size and shape) and an obligation to interact with other vertices or faces while respecting kinematic (geometric), equilibrium, and boundary constraints. These constraints are manifested in several ways. Large coordination numbers $n$ are discouraged by the geometric impossibility (or unlikelihood) of fitting numerous particles around a central particle; whereas, an $n=2$ is discouraged by a frictional limit that can cause a central particle to “squirt” from its two neighbors. These two constraints are examined in the following section.

In regard to the void valence distribution (Fig. 5a), the occurrence of high-valence voids is discouraged by their inherent instability Hunt:2010a , since large voids are short lived and readily collapse onto their inner rattler particles; whereas, small voids ($l=3$) are discouraged by the geometric necessity that sets of three particles be coordinated with small interior angles (less than $\pi/2$, as evident in Fig. 1b) while avoiding interior rattlers.

One also notes the spatial patterning of void valences in Fig. 1b: voids of valence 3 and 4 are typically clustered to form ladder-like chains, and voids of valence 6 and greater are also clustered near each other. Examples of such chains and clusters are shaded orange and green in Fig. 1b (these patterns are more apparent when the periodic figure is multiply tiled). Stress transmission and deformation are also know to be spatially organized into “force chains” and “micro-bands” Radjai:1996a ; Kuhn:1999a ; Azema:2007a , and these patterns affect the statistics of force and motion Radjai:1998a in a manner that is not yet fully understood Chakraborty:2010a . The spatial ordering of void valences is certainly beyond Eqs. (5)–(10).

With its limitations, the model is not fully predictive, but even this primitive view of granular flow as a solely topological process yields a reasonable prediction of the distribution of void valences. An entropy-based theory can, of course, be improved by supplying additional information (pre-knowledge) that further restricts the configuration space, and such information, in the form of geometric and equilibrium biases, is applied to the distribution of coordination numbers in the next two sections.

3 Geometric and kinetic constraints (Model II)

Recognizing that geometric limitations can affect particle packings, we now consider the ease (and, by extension, the likelihood) with which $n$ outer particles of sizes $\mathsf{S}=\{D_{1},D_{2},\ldots,D_{n}\}$ can be packed around an inner particle of size $D_{0}$, such that all particles of the ordered set $\mathsf{S}$ touch the inner particle (Fig. 6).

Unlike the dense packing algorithm in Hihinashvili:2012a , we permit loose packings in which gaps occur among the outer (shielding) particles. We evaluate the angle $\theta_{j}$ between the center of particle $j$ in $\mathsf{S}$ and its neighbor $j+1$, when the two particles, $D_{j}$ and $D_{j+1}$, are temporarily placed together and in contact with $D_{0}$ ($\theta_{n}$ is the final angle between $D_{n}$ and $D_{1}$). In some cases, it is not possible to pack the $n$ particles around $D_{0}$: if the sum $\sum_{j=1}^{n}\theta_{j}$ exceeds $2\pi$, we assign a probability of zero for this combination of particles. In other cases, the $n$ particles will fit with excess space. We contend that the likelihood of these arrangements is proportional to the excess angular space (gap) $\alpha_{\mathsf{S}}\equiv 2\pi-\sum_{j=1}^{n}\theta_{j}$, as in Fig. 6. The probability of encountering a particular arrangement $\mathsf{S}=\{D_{1},D_{2},\ldots,D_{n}\}$ around a central particle $D_{0}$ will depend on their gap $\alpha_{\mathsf{S}}$ and on the distribution of particle sizes in the entire assembly (i.e., the probabilities of the individual sizes $D_{0},D_{1},D_{2},\ldots$), as described below.

Besides this geometric influence, a kinetic constraint applies in the case of only two outer particles, $n=2$, due to the limiting friction coefficient $\mu$ (Fig. 7).

Unless the exterior angle between the two outer particles is less than $2\tan^{-1}\mu$, equilibrium can not be maintained, and the inner particle will squirt from its two neighbors. As a result, an angular range (gap) $\alpha_{\mathsf{S}}=4\tan^{-1}\mu$ is available for the arrangement of two particles around an inner particle. Within this available angular range, the mobilized friction can be as small as zero (when the particles are diametrically opposed) and as large as $\mu$ (when the particles are in the limiting arrangements shown in Fig. 7), an assumption that is consistent with experimental results Majmudar:2005a .

A poly-disperse assembly with $r$ different particle sizes admits $r^{n}$ possible ordered sets $\mathsf{S}$ of cardinality $n$. We assume knowledge of the distribution of the $r$ sizes: size $D_{k}$ comprises fraction $P_{D_{k}}$ of an entire assembly, such that $\sum_{k=1}^{r}P_{D_{k}}=1$. We contend the likelihood of encountering a particular ordered combination $\mathsf{S}$ (of $n$ particles around an inner particle of size $D_{0}$) among all such combinations of cardinality $n$ is proportional to three factors: to the gap $\alpha_{\mathsf{S}}$ among the $n$ particles in $\mathsf{S}$, to the probability $P_{D_{0}}$ of the inner size $D_{0}$, and to the product of the $n$ probabilities of the outer sizes comprising $\mathsf{S}$:

[TABLE]

If so, the probability of encountering a particle with $n$ neighbors in the entire assembly is proportional to the sum $Q_{n}$ of all such $r^{n}$ combinations of $\mathsf{S}$ and all $r$ sizes $D_{0}$:

[TABLE]

The sum of the probabilities for all coordination numbers, $n=2,3,\ldots$, must equal one, so that the estimated probability of coordination number $n$ is

[TABLE]

This equation gives the probabilities of coordination numbers $n$ for the second “II” model. Although no entropy principle is explicitly applied, we assume that the probabilities in Eq. (13) are not biased by factors other than geometric packing constraints and the kinetic, frictional constraint. Within these constraints, maximum entropy is attained for each set of sizes $\mathsf{S}$ by a uniform distribution of gaps $\alpha$ within the range $|\alpha|\leq\alpha_{\mathsf{S}}$.

Anisotropy of particle arrangements is a dominant feature of granular flow at the critical state Radjai:2012a ; Kuhn:2010a . The model does not explicitly incorporate macro-scale geometric anisotropy, but the model does permit anisotropy at the micro-scale: for example, the particles $D_{1}$ and $D_{2}$ in Fig. 7 are vertically aligned, and during vertical compression this arrangement would be more likely than a horizontal alignment.

3.1 Simulation results and Model II

With a bi-disperse assembly ($r=2$), the gaps and probabilities of all $r^{n}=2^{n}$ combinations of sizes can be computed in reasonable time. The two particle sizes sizes, $D_{a}=1$ and $D_{b}=1.5$, were in proportions $P_{D_{a}}=2.25/3.25$ and $P_{D_{b}}=1/3.25$, such that both species comprise the same total area. Figure 8 shows the estimated probabilities $P_{n}^{\text{II}}$ for assemblies with $\mu=0.50$.

This figure compares the estimated distribution $P_{n}^{\text{II}}$ with the DEM simulations of dense granular flow at the critical state. These DEM results are also shown in Fig. 5b. The model II yields a better representation of the coordination number probabilities than model I, capturing the general trend of the simulation results and yielding a smaller fraction having $n=2$. The model, however, over-estimates the average coordination number: the average $\overline{n}$ of the simulations is $3.23$; whereas, the estimated average is $3.40$.

The arguments that were used in developing this method would suggest that the estimated probability of a particle having only two outer contacts is affected, in part, by the friction coefficient $\mu$: larger coefficients permit a greater range of orientations of two outer particles (see Fig. 7) and should increase the likelihood of coordination number 2. The average coordination number should, therefore, be smaller for assemblies having a larger friction coefficient. We conducted additional simulations with five different coefficients $\mu$. These coefficients lead to different topological arrangements, as expressed in different average coordination numbers $\overline{n}$ and in different distributions of $n$. The average coordination numbers are shown in Fig. 9 and are compared with the estimates of Eqs. (12)–(13).

The decrease in $\overline{n}$ with increasing $\mu$ is consistent with results in Shaebani:2012a . The estimates of the average coordination number $\overline{n}$ follow the trend of the simulation data, but $\overline{n}$ is over-predicted for all but the lowest friction coefficients, $\mu\approx 0.10$.

Other simulations of moderately poly-disperse materials Shaebani:2012a have shown that the average deficit angle $\alpha_{\mathsf{S}}$ is about the same for the smaller and larger particles in an assembly. These results were reported for the initial, compacted state, prior to deviatoric loading. At the critical state for the bi-disperse simulations, we measured average deficits of 162*∘* and 172*∘* for the small and large particles. Eqs. (11)–(13) include all possible clustered arrangements $\mathsf{S}$, with the probability of each arrangement being naively weighted by its deficits $\alpha_{\mathsf{S}}$. An estimate of the coordination number distribution can be improved by including information (pre-knowledge) of the average deficit for each central particle size. In the next section, we instead apply information of the average coordination number $\overline{n}$ to improve the estimated distribution of $n$.

4 Disorder with geometric and kinetic biases (Model III)

In this section, we develop an amalgam of the previous two approaches by applying Kullback’s minimum cross-entropy (MinXEnt) principle Kullback:1951a ; Kapur:1992a . To the author’s knowledge, this principle has not yet been applied to granular materials. As with Jaynes’ maximum entropy (MaxEnt, model I) approach in Section 2, the average $\overline{n}$ is applied as a rigid moment constraint on the probability distribution $P_{n}$, but we also profess certain inclinations of the probabilities in the form of “a priori estimates” $q_{n}$. The directed distance $H^{\text{III}}_{n}(\vec{P}:\vec{q})$ between the two distributions — from $\vec{q}$ to $\vec{P}$ — is the quasimetric

[TABLE]

In our case, the probabilities $P_{n}^{\text{II}}$ derived in the previous section from geometric and kinetic considerations are applied as the a priori estimates $q_{n}$ of coordination numbers $n$:

[TABLE]

The MinXEnt principle calls for minimizing the distance $H^{\text{III}}_{n}$ subject to the applicable moment constraints, which will lead to a new set of probabilities $P_{n,\Delta M}^{\text{III}}$. As before, the average coordination number $\overline{n}=\langle n\rangle$ and contact additions $\overline{\Delta M}=\langle\Delta M\rangle=\overline{n}/2$ are the moment constraints that apply to the problem of coordination numbers (see Section 2 and Table 1):

[TABLE]

Minimizing (14) with respect to the $P_{n,\Delta M}^{\text{III}}$ subject to Eqs. (16) and (17) leads to the probabilities of the third “III” model:

[TABLE]

with Lagrange multipliers, $\lambda_{n,1}$ and $\lambda_{n,2}$, that satisfy Eqs. (16) and (17).

The estimated probability $P_{n}^{\text{III}}$ of encountering coordination number $n$ is the marginal probability

[TABLE]

This approach incorporates more information than either of the previous methods: it combines the rigid moment constraints of the average coordination number $\overline{n}$ and the average added contacts $\overline{\Delta M}$, along with the geometric and kinetic inclinations $q_{n}$, which are intentional biases that arise from the friction coefficient and distribution of particle sizes. If the inclinations $q_{n}$ were to yield, by themselves, the same averages $\overline{n}$ and $\overline{\Delta M}$ (that is, if $\sum n\,q_{n}=\overline{n}$), then the moment constraints of Eqs. (16) and (17) bring no further information beyond that provided by the $q_{n}$. In this case, the distribution $P_{n}^{\text{III}}$ will coincide with the $P_{n}^{\text{II}}$ of Section 3. On the other hand, if the geometric and kinetic inclinations do not bias the results (that is, if $q_{2}=q_{3}=\ldots$), then the inclinations of Eq. (15) bring no further information beyond that of Eqs. (16) and (17), and the distribution $P_{n}^{\text{III}}$ will equal the $P_{n}^{\text{I}}$ of Section 2. (We note a subtle inclusion of such inclinations in that section, when we disallowed coordination numbers greater than seven. In a sense, we had assumed $q_{2}=q_{3}=\ldots=q_{7}=1/6$, $q_{8}=q_{9}=\ldots=0$.)

4.1 Simulation results and Model III

The distribution of coordination numbers predicted by model III (Fig. 10) is fairly close to the DEM results and is better than those predicted with models I and II (Figs. 5b and 8).

As with model I, which predicts a uniform distribution across $\Delta M$ for each $l$ (Fig. 4), model III predicts a uniform distribution across $\Delta M$ for each $n$.

All three estimates of coordination number distribution are summarized in Table 2, which compares them with the DEM results. In a recursion of the directed distance of Eq. (14), we use the following sum as the distance of an estimated distribution from the measured DEM distribution $P_{n}^{\text{DEM}}$,

[TABLE]

where $P_{n}^{(\bullet)}$ is a particular estimate ($P_{n}^{\text{I}}$, $P_{n}^{\text{II}}$, or $P_{n}^{\text{III}}$).

In regard to the distribution of coordination numbers, the model II, which respects only geometric and kinetic information, yields a better estimate than a model I, which disregards such information and only considers topological dispersion. Model III, which combines these factors, gives the best estimate.

5 Conclusion

When a granular material is slowly sheared from an initial rested condition, the critical state is eventually attained at large shear strains, and the bulk characteristics of this state — density, fabric, and strength — are insensitive to the initial particle arrangement. Because of this resilience, the critical state is used as a reference state against which other conditions are compared (e.g., the jamming threshold in powder flows or with the so-called state parameter in geomechanics Been:1991a ). The critical state seems to be characterized by maximum disorder, as simple disorder models are shown to be modestly successful in predicting the distributions of void valence and coordination number. Such disorder is moderated by a tendency toward spatial patterning of the topological arrangement (Section 2.2) and by a bias toward certain $(l_{i},\Delta M_{i})$ pairs in the particle graph (Fig. 4). We should also note that mono-disperse assemblies have a tendency to crystallize into hexagonal arrangements of particles, although this was not observed in the bi-disperse simulations of this study.

Notwithstanding these tendencies for greater order, the view of the critical state as a maximally disordered condition could also be extended to other micro-scale characteristics, such as contact forces, inter-particle motions, and contact orientation, and these matters should be the focus of future study. The method should also be extended to three-dimensional assemblies, although this poses new difficulties: in particular, a journaling scheme for constructing a given granular topology by progressively adding particles, contacts, and polygonal faces around an existing arrangement to form new volume cells. Finally, greater understanding should be sought for the subtle tendency of greater order among the occurrences of $\Delta M$ for each $n$ and $l$ (as in Fig. 4) and for the effects of macro-scale anisotropy on the topological entropy. These and other extensions of maximum disorder models will lead to a richer understanding of critical state flow.

Acknowledgements.

This work is dedicated to the memory of Prof. Colin B. Brown (1929–2013), who made significant contributions to the understanding of granular entropy.