Maximum disorder model for dense steady-state flowof granular materials

Maximum disorder model for dense steady-state flow
of granular materials

Matthew R. Kuhn Donald P. Shiley School of Engineering, University of Portland, 5000 N. Willamette Blvd., Portland, Oregon, USA 97203,, (503)943-7361

A flow model is developed for dense shear-driven granular flow. As described in the geomechanics literature, a critical state condition is reached after sufficient shearing beyond an initial static packing. During further shearing at the critical state, the stress, fabric, and density remain nearly constant, even as particles are being continually rearranged. The paper proposes a predictive framework for critical state flow, viewing it as a condition of maximum disorder at the micro-scale. The flow model is constructed in a two-dimensional setting from the probability density of the motions, forces, and orientations of inter-particle contacts. Constraints are applied to this probability density: constant mean stress, constant volume, consistency of the contact dissipation rate with the stress work, and the fraction of sliding contacts. The differential form of Shannon entropy, a measure of disorder, is applied to the density, and the Jaynes formalism is used to find the density of maximum disorder in the underlying phase space. The resulting distributions of contact force, movement, and orientation are compared with two-dimensional DEM simulations of biaxial compression. The model favorably predicts anisotropies of the contact orientations, contact forces, contact movements, and the orientations of those contacts undergoing slip. The model also predicts the relationships between contact force magnitude and contact motion. The model is an alternative to affine-field descriptions of granular flow.

granular matter, entropy, critical state, fabric, anisotropy

1 Introduction

The critical state concept in geomechanics holds that dense granular materials, when loaded beyond an initial static packing, eventually attain a steady state condition of constant density, fabric, and stress (Schofield and Wroth, 1968). This condition is often associated with shear-driven flow and failure: granular avalanches, landslides, tectonic faults, and failures of foundation systems and embankments. As such, the critical state has received intense interest from geologists, engineers, and physicists, who have devoted great effort in understanding the state’s underlying mechanics. Density, fabric, and deviatoric stress at the critical state are known to depend upon the particles’ shapes and contact properties as well as on the mean stress and intermediate principal stress (Zhao and Guo, 2013). Even so, the eventual bulk characteristics for a given assembly are insensitive to the initial particle arrangement and to the stress path that ends in the critical state: for example, materials that are initially either loose or dense eventually arrive at the same density condition after sufficient shearing. This convergent characteristic resembles that of thermal systems that approach an equilibrium condition with sufficient passage of time.

Another pervasive feature of critical state flow is the continual and intense activity of grains at the micro-scale, yet this local tumult produces a monotony in the bulk fabric, stress, and density. Micro-scale activity occurs in three ways: (a) statically, as alterations of the inter-particle contact forces, (b) geometrically, as changes in the particles’ configuration and local density, and (c) topologically, as changes in the load-bearing contact network among the particles. In these three respects, we view the critical state as the condition of maximum disorder that emerges during sustained shearing. Early work by Brown et al. (2000) investigated disorder in the local density, and a recent paper by the author explored topological disorder at the critical state (Kuhn, 2014). The current paper addresses statical disorder as expressed in a two-dimensional (2D) setting, focusing on anisotropies and distributions of the contacts’ orientations, forces, and movements. The analysis applies to the critical state flow of dry unbonded frictional materials of sufficient density to develop a load-bearing (persistently jammed) network of contacts during slow (quasi-static, non-collisional) shearing. Particles are assumed durable (non-breaking) and nearly rigid, such that deformations of the particles are small, even in the vicinity of their contacts.

During flow, granular materials have an internal organization of movement and force, an organization that is in some respects pronounced but in others subtle. We briefly review these characteristics of the critical state, as observed in laboratory experiments, numerical simulations, or both.

  1. The motions of individual particles do not conform to an affine, mean deformation field, and fluctuations from the mean field are large and seemingly erratic (Kuhn, 2003). The rates at which contacting particles approach or withdraw from each other (i.e. contact movements in their normal directions) are generally much smaller than those corresponding to an affine field. In contrast, the transverse, tangential movements between contacting pairs are much larger than those of affine deformation (Kuhn, 2003). As a result, bulk deformation is almost entirely attributed to the tangential movements of particles (Kuhn and Bagi, 2004).

  2. Strength, expressed as a ratio of the principal stresses, is insensitive to the contacts’ elastic stiffness and to the mean stress, such that simulations of either soft or hard particles exhibit similar strengths at the critical state (Härtl and Ooi, 2008; Kruyt and Rothenburg, 2014). Fabric measures at the critical state (fraction of sliding contacts, contact anisotropy, etc.) are also insensitive to contact stiffness (da Cruz et al., 2005).

  3. Particle rotations are large when compared with the bulk deformation rate (Kuhn and Bagi, 2004). In particular, the rolling motions between particles are much larger than the sliding movements (Kuhn, 2004b).

  4. The contact network is relatively sparse, in that the number of contacts within the load-bearing contact network is sufficient to produce a static indeterminacy (hyperstaticity) but with only a modest excess of contacts (Thornton, 2000). The excess very nearly corresponds to the number of contacts that are sliding (Kruyt and Antony, 2007). This modest indeterminacy is consistent with observations of intermittent, sudden reductions of stress, which result from periodic collapse events that are occasioned by fresh slip events or loss of contacts (Peña et al., 2008). This condition of marginal hyperstaticity is referred to as “jammed” within the granular physics community.

  1. During critical state flow, the contact fabric is anisotropic, with the normals of the contacts oriented predominantly in the direction of the major principal compressive stress (Rothenburg and Bathurst, 1989).

  2. The normal contact forces are larger among those contacts oriented in the direction of the major principal compressive stress; whereas, the averaged tangential forces are larger for contact surfaces oblique to the principal stress directions (Rothenburg and Bathurst, 1989; Majmudar and Bhehringer, 2005).

  3. Deviatoric stress is primarily borne by the normal contact forces between particles; whereas, the tangential contact forces make a much smaller contribution to the deviatoric stress (Thornton, 2000).

  4. When considering only the normal forces, deviatoric stress is primarily carried by those contacts with forces that are larger than the mean force (strong contacts), whereas the remaining (weak) contacts contribute far less to the deviatoric stress (Radjai et al., 1998; Kruyt and Antony, 2007).

  5. Anisotropy of the contact network is also largely attributed to strong contacts, which are predominantly oriented in the direction of the major principal stress (Radjai et al., 1998).

  6. Many contacts slide in the “wrong direction” with respect to the direction that corresponds to an affine deformation (Kuhn, 2003).

  7. Compared with other orientations, contacts that are oriented in the direction of extension have a greater average magnitude of slip, but the mean slip velocity is largest among contacts that are oriented obliquely to the directions of compression and extension (Kuhn and Bagi, 2004).

  8. Frictional sliding is more common among those contacts with a smaller-than-mean normal force (Radjai et al., 1998).

  9. The more mobile contacts — those with large sliding movements — tend to be those that bear a smaller-than-average normal force (i.e., weak contacts) (Kruyt and Antony, 2007).

  10. Deformation, when measured at the meso-scale of particle clusters, is related to contact orientation: contacts with branch vectors that are more aligned with that of bulk compression tend to produce local dilation; whereas, contacts that are more aligned with the direction of bulk extension tend to produce local compression. These trends have been determined by studying the elongations of voids that are surrounded by rings of particles and their branch vectors (Nguyen et al., 2009).

  11. The probability density of the normal contact forces usually decreases exponentially for forces that are greater than the mean (Majmudar and Bhehringer, 2005). With forces less than the mean, however, the density is more uniform than exponential.

  12. Strength at the critical state increases with an increasing inter-granular friction coefficient, but the relation is non-linear, and little strength gain occurs when the friction coefficient increases beyond 0.3 (Thornton, 2000; Härtl and Ooi, 2008). Fabric anisotropy also increases with an increasing friction coefficient (da Cruz et al., 2005; Kruyt and Rothenburg, 2014).

  1. Relatively few contacts attain the frictional limit at any particular moment during flow, and sliding typically occurs among only 8%–20% of the contacts in assemblies of disks and spheres (Radjai et al., 1998; Thornton, 2000). The fraction of sliding contacts is reduced when the friction coefficient is increased (Thornton, 2000).

  2. For contacts that are not sliding, the probability density of the mobilized friction is greatest for the nearly neutral condition of zero tangential force, and the density decreases with increasing mobilized friction (Majmudar and Bhehringer, 2005).

  1. Internal force, movement, and deformation are highly heterogeneous and are spatially organized across scales of ten or more particles (Kuhn, 2003). Contact forces are patterned in “force chains” of highly loaded particles that are roughly aligned with the major principal compressive stress (Majmudar and Bhehringer, 2005). Deformation is localized into obliquely oriented microbands of thickness 1–3 particle diameters (e.g., Kuhn, 1999) and into shear bands with a thickness of 8–20 diameters (Desrues and Viggiani, 2004). Local stiffness (Tordesillas et al., 2011) and local dilation and contraction are often clustered (Kuhn, 1999), and chains of rapidly rotating particles are organized obliquely to the major principal stress axes (Kuhn, 1999). The non-affine particle movements in 2D assemblies are coordinated in large circulating vorticity cells that encompass several dozens of particles (Williams and Nabha, 1997). In 2D materials, the topology of the contact network is also spatially organized: voids surrounded by 3-5 particles form elongated chains, and larger voids are often organized as clusters (Kuhn, 2014).

These and other characteristics, as determined from new DEM results, will be discussed and compared in relation to a proposed model of critical state flow. Whereas the results listed above were gained from descriptive, observational studies, the current work is largely predictive. The above characteristics are arranged as follows: “A” characteristics serve as the basis in developing the model, “B” characteristics are favorably predicted by the model, “C” characteristics are not well predicted or must be forced into agreement with an ad hoc intervention, and “D” characteristics are inaccessible, due to the nature of the model.

The proposed model focuses on the movements, forces, and orientations of inter-particle contacts during shear-driven flow. The model is premised on the conjecture that the critical state is a condition of maximum micro-scale disorder among these contact characteristics. This maximally disorder state is extracted, and its resulting conditions are shown to exhibit many of the behaviors listed above. With some reluctance, we use the term “entropy” interchangeably with disorder, with entropy connoting a statistical measure of micro-scale unpredictability (i.e., as in Shannon, 1948) rather than an extensive thermodynamic quantity having an intensive temperature-like dual (in contrast with the compactivity or angoricity settings, e.g. Blumenfeld and Edwards, 2009).

Principles of maximum disorder have been applied to granular materials for several decades. Brown et al. (2000) conducted experiments on two-dimensional assemblies of spheres, and by applying a back-and-forth shearing, found that disorder in the local density and coordination number increased with each shearing cycle. The Jaynes maximum entropy (MaxEnt) formalism has typically been used to predict the condition of maximum disorder (Jaynes, 1957). This approach has been applied to the local fabric of disk assemblies by categorizing voids into several canonical types (Brown, 2000). Similar maximum entropy approaches have also been applied to local packing density (Yoon and Giménez, 2012), to the contact forces (Chakraborty, 2010), to the contact orientations (Troadec et al., 2002), and to the contact displacements and bulk elastic moduli (Rothenburg and Kruyt, 2009). Among these aspects, the distribution of contact force has recently received the greatest attention. Early theories derived probability densities of force distribution by addressing the bulk mean stress but without respecting the local force equilibrium (Goddard, 2004). More recent approaches to the distribution of contact force have enforced the local equilibrium of particles and examined the statically admissible space of contact forces within a granular system (Chakraborty, 2010). These disorder theories address frictionless, non-flowing assemblies and do not address the anisotropies of force, fabric, and movement that accompany flow. The paper not only admits flow but relies upon the shearing deformation to drive the contact movements and the contact forces that generate the bulk anisotropies of fabric and force.

Our approach also differs from mean field theory, which is currently the dominant paradigm for relating the micro- and macro-scale behaviors of dense materials, primarily at small strains. Mean-field theories include the upscaling methods of (Jenkins and Strack, 1993), in which the small-strain bulk stiffness is estimating by assuming that the particle motions or contact forces around a central particle conform to an affine, homogeneous kinematic (motion) or static (stress) field. We will not impose such affine restrictions, as doing so would imply a strong order in the critical state condition — an order that simply does not exist (Kuhn, 2003). Most mean-field approaches are also inherently deterministic insofar as they result in the unique response of a particle for a given description of its neighborhood. Rothenburg and Kruyt (2009) extended such methods by adopting probability densities for both the surroundings and the response. The current work also adopts a probabilistic framework, but the approach is intended for granular flow in which the micro-scale interactions are dominated by tangential contact motions.

The plan of the paper is as follows. We build a micro-scale flow model by first identifying an essential set of contact quantities, forming a rather comprehensive phase space of motion, force, and arrangement (Sec. 2.1). We regard these quantities as random variables with a bulk probability density distribution (Sec. 2.2) and apply constraints, both rational and empirical, to this distribution (Secs. 2.32.7). By themselves, these constraints do not describe a unique distribution. We assume that all such micro-states satisfying the constraints are equiprobable during critical state flow and that the most likely macro-state (i.e., density distribution) is the one that encompasses the greatest breadth of micro-states, maximizing the system’s disorder. A maximum entropy condition is applied to the contact attributes to arrive at a predicted distribution (Sec. 2.8). Section 3 describes the model’s solution and compares it with past observations and with the results of new DEM simulations.

Although the principles presented in these sections are fundamental, the set of constraints sparse, and the results promising, their evaluation is computationally taxing. Complete computational details, therefore, are provided in appendices. We end the paper by considering a number of questions raised by some difficulties and shortcomings of the entropy model: which aspects of granular flow are amenable to simple, sweeping statistical approaches, and which aspects require closer attention to the details of grain interactions? What additional constraining information is most appropriate for these entropy methods? How can a model that is conceptually simple be both predictive and opaque?

2 Entropy Model

The model is developed in a two-dimensional setting of assemblies of disks that are poly-disperse but of a narrow size range. Because of the assumed small poly-dispersity, we can overlook the tendency of mono-disperse assemblies to develop ordered, crystallized (low entropy) particle arrangements. A narrow size range also permits characterizing the disk sizes by their mean diameter . The model is intended for durable and nearly rigid disks, in which contact indentations are small relative to particle size: the contacts are assumed rigid-frictional and without contact moments.

2.1 Contact quantities

The foundation of any entropy model is the set of chosen quantities that characterize the micro-states that comprise the full phase space of a physical system. The quantities that one chooses will largely determine the nature of the model and the phenomena that are addressed. For example, a phase space composed exclusively of particle arrangements will lead to a random isotropic state. A phase space composed exclusively of contact forces will lead to estimates of the static force distribution (e.g., Goddard, 2004). Our intention is to model conditions of movement, force, and arrangement during granular flow, and we choose a rather large, comprehensive phase space of certain contact quantities. These quantitites can be used to extract or to constrain the bulk deformation rate tensor, stress tensor, and fabric tensor, and these micro-quantities are all accessible from simulations of discrete particle systems.

A micro-state of a large two-dimensional assembly of disks is enumerated with six quantities at each of the assembly’s loading-bearing contacts, such that a micro-state of the entire assembly is a single point within a phase space of dimension . The six quantities, listed in Table 1, apply to a contact that is shared between two particles, and , which comprise an ordered pair (Fig. 1a).

Quantity Description
Compressive normal contact force,
Tangential contact force,
Orientation of contact normal vector
Orientation of branch vector
times contact slip, “
times rigid contact movement, “
Directional conventions are shown in Fig. 2.
and are also limited by Eq. (4).
Table 1: Contact quantities of the entropy model (Fig. 1)
Figure 1: Contact quantities of contact between two contacting disks . Direction conventions for biaxial symmetry are shown in Fig. 2.

Contacts are assumed to be enduring, non-collisional, and persistent over the brief interval for which bulk averages are extracted. Note that, henceforth, we will largely ignore the association of contacts with particles, thus overlooking the underlying topology of the contact network and relinquishing an ability to assure local equilibrium or kinematic compatibility. The compressive normal contact force acts at the contact “” upon particle “” (Fig. 1b). Because the average normal force in a two-dimensional assembly is roughly proportional to the mean diameter and to the mean stress , we normalize the individual normal forces as . The scalar tangential force on is , normalized as , and this force component is limited by friction, as characterized by coefficient : . The sign convention associated with and will be described further below.

Only tangential contact movements are modeled, since normal motions (those that alter contact indentations or cause the separation or creation of contacts) are known to contribute negligibly to bulk deformation during flow (characteristics A.1 and A.2 in the Introduction). During deformation, particles will both rotate and shift relative to each other, and we adopt two systems for describing these movements. In the first system, scalar represents an angular rate produced by translation of the center of relative to the center of — a rotational rate of the unit normal vector directed outward from (Fig. 1c). Note that the rate vector is orthogonal to . The relative motions of the particles at their contact are produced both by this tangential shifting of the particles’ centers (i.e., the velocity ) and by the two particles’ rotations, and (Fig. 1c). Motions in the first system are described by the quantities , , and . The second system more directly addresses the contact interactions, as expressed through three mechanisms: 1) slip (sliding) between the particles, 2) rolling movements that produce no slip, and 3) rigid rotations of the particle pair, also producing no slip. For equal-size disks of diameter , the mechanisms are as follows (Kuhn and Bagi, 2005):


all having the units time. Because all three mechanism usually occur simultaneously at any contact, we decompose the motions , , and as


forming an orthogonal separation of particle movements into three “” contact rates. Two quantities, and “”, are most relevant in the model, as they determine the bulk deformation and dissipation rates. Because both quantities can be expressed as linear combinations of and , the model does not require , and the reduction of three quantities to two greatly reduces computational demands.

As the final contact quantities, we introduce two orientation angles, and , with angle giving the direction of the contact normal , and angle representing the direction of the branch vector , of length , that joins the center of to that of (Fig. 1a). Although the two angles are identical for disks, they are distinguished because of their different roles in expressing stress and strain within the assembly. Both angles are measured from particle .

Although the model is quite general, it will be applied to the particular condition of biaxial compression, in which the width of the assembly is reduced while the height expands (Fig. 2a). We adopt a sign convention that takes advantage of the symmetry of these conditions and assists bookkeeping when the model is later compared with DEM results. The convention, shown in Fig. 2b, assigns a direction to the tangential unit vector corresponding to the sliding that would be expected if the particle motions roughly conformed to affine deformation during horizontal biaxial compression.

Figure 2: Conventions for unit vectors and and rotations when analyzing horizontal biaxial compression. “Forward” directions are shown.

This direction — clockwise or counterclockwise — alternates across quadrants of . Unit vector applies to the (scalar) tangential force on (i.e., and ) and to the relative contact movements (, , etc.). Positive values of , , etc., which are aligned with , are called “forward” forces and movements; whereas, negative values are in the “reverse” direction . Normal force is always compressive, and unit vector is directed outward from particle .

The domains of the six quantities are shown in Table 1. We also place unilateral and rigid-frictional constitutive limits upon the model,


further restricting the domains of , , and . These restrictions mean that the particles are unbonded (), and the tangential force is limited by friction to the range . Positive (forward) slip, , is only admissible with a corresponding positive tangential force, ; whereas, negative (reverse) slip is only possible for the opposite condition, . When the friction limit is not reached, , the slip is zero: . In adopting Eqs. (3) and (4), we forego a more elaborate elastic-frictional model, since critical state flow is insensitive to contact elasticity (characteristic A.2 of the Introduction). Note that the model is rate-independent, since forces depend only on the directions of and not on their magnitudes.

2.2 Probabilities and constraints

In the model, we do not individualize the six “” quantities among the contacts (a phase space of dimension ) but instead treat the quantities as continuous random variables with probability density :


dropping the subscripts. Although we write the complete phase space of possible micro-states as , all micro-states within this space that share the same probability density comprise a common macro-state, an element in a subspace of and written without parentheses: or . As will be seen, density contains comprehensive information about bulk quantities, such as stress and fabric, and about the correlations among contact orientation, force, and movement. However, by turning from individual “” contact quantities to their gross probability distribution, we forfeit the possibility of tracking individual contacts. In return, we gain a certain economy of expression for predicting overall distributions of the six quantities and the statistical relationships among them.

A wide range of the six quantities are observed in experiments and simulations, and correlations among the quantities are found during critical state flow, as enumerated in the Introduction. These correlations are expressions of order, which we consider a consequence of certain constraints that link the micro and bulk behaviors. We will restrict the admissible macro-states with five constraints: four fundamental constraints and one auxiliary constraint. These constraints bring a priori information that will bias the distribution . Each “”th constraint is expressed with a constraining function of the six random variables,


The expected value of the function is found through 6-fold integration across the full phase space, with each random variable integrated over the range given in Table 1:


and measure is implied in these integrations. The integrations are described in greater detail in A.

The expected values of functions will be constrained to certain average values,


with the prescribed averages, , developed below. That is, each restriction (8) is a moment constraint on density .

An essential “zeroth” constraint follows from the certainty that the six quantities lie within their full ranges:


such that integration of over its full domain equals 1.

The four fundamental constraints express general principles of isochoric dissipative flow at constant stress and are derived in Sections 2.32.6. The auxiliary constraint applies additional information gained from DEM simulations (Section 2.7). The constraints are summarized in Table 2.

Constraint Description Variables in Eq.
Fundamental 1 Mean stress (20)
2 Dissipation consistency (27)
3 Isochoric flow (38)
4 Biaxial rate, (43)
Auxiliary 5 Fraction of sliding contacts, (45)
Table 2: Summary of constraints in entropy model.

2.3 Fundamental constraint 1: Constant mean stress

We model steady state flow under constant mean stress , a condition commonly applied in geotechnical testing. The Cauchy stress in a large granular assembly is the volume average of the contact dyads ,


in which branch vector joins particle to , contact force acts upon , and is the number of contacts within the interior of the two-dimensional region of area (Rothenburg, 1980). The branch vector length is approximated as , so that


and the contact force is the sum of its tangential and (compressive) normal components:


In these and future expressions, and are unit orientation vectors associated with the contact normal and branch vectors:


For the particular sign convention shown in Fig. 2b, unit vector is


where the signum function alternates in sign between coordinate quadrants.

We replace the sum in Eq. (10) with an integration by substituting Eqs. (11) and (13), multiplying by probability and the number of contacts , and integrating across the full domain of the six quantities :


In this expression, we include a kernel function , which is simply the Dirac operator,


The Dirac kernel formally dictates the correspondence of the unit branch vector and the contact force vector for circular particles (by extension, unit vectors and ), as in the single sum of Eq. (10).

Equations (16)–(18) will later be used to extract the predicted stress tensor during critical state flow, but we now use them to enforce the first fundamental constraint on the probability density : the negative mean stress must equal the assigned pressure . The negative mean value of tensor function is


and equating its integration in Eq. (16) with leads to the first fundamental constraint


in which we have divided by the leading constant in Eq. (19). Equation (20) simply requires the dimensionless mean normal force (i.e., ) to equal the dimensionless quantity , thus enforcing pressure . The contact density  is about 1.5, a value that can be found with DEM simulations or estimated from characteristic A.4 of the Introduction.

2.4 Fundamental constraint 2: Dissipation consistency

During critical state flow at constant stress and volume, deformation is entirely plastic: all work done by the applied stress is dissipated through internal irreversible processes. We assume that the only available dissipation mechanism is frictional sliding between particles. This assumption is the basis of the second fundamental constraint on the probability density and provides the essential link among contact motions, contact forces, bulk deformation, and bulk stress. This constraint drives the anisotropies in force and movement.

The frictional dissipation rate at the contact of two disks is the product of its tangential contact force and the relative slip of the two particles’ surfaces at their contact: the “” rate in Eq. (1). Applying Eq. (2), this slip rate corresponds to as


Noting that and must conform to the constitutive rigid-frictional restrictions of Eqs. (3) and (4), the frictional dissipation rate per unit of area is


or, in terms of probability density ,

Dissipation (23)

where the Dirac kernel is applied again, enforcing the coincidence of the contact forces (associated with contact orientations ) and the contact motions (associated with branch vector orientations ) in the single sum of Eq. (22).

The internal work rate of the Cauchy stress is the inner product


where is the symmetric part of velocity gradient . Although Eq. (16) is an expression of , and an expression for is derived in the next section (as Eq. 35), directly substituting these expressions into the full inner product of Eq. (25) will lead to a non-linear constraint on probability , making its evaluation intractable. The situation is greatly simplified when deformation can be expressed with a single non-zero parameter. We consider the case of biaxial compression, for which is reduced to


and the work rate is .

Equating the rates of frictional dissipation and internal work in Eqs. (23) and (25) and substituting Eqs. (16), (17), and (26) lead to the second fundamental constraint on density for critical state flow in biaxial compression:


in which we have eliminated the leading term that appears in Eqs. (17) and (24). By substituting Eqs. (14), (15), (17), and (24), the function can be expanded as


which enforces consistency of work and dissipation.

2.5 Fundamental constraint 3: Isochoric flow

The third fundamental constraint (and, by far, the most difficult to derive) enforces the constant density (isochoric) condition of critical state flow. If tensor is the average, bulk velocity gradient within a granular assembly, its trace vanishes during steady state flow. To enforce this condition in an assembly with contacts, we require a means of estimating from the contact motions of particle pairs: from an -list of contact information . Several methods have been proposed for estimating (see Bagi, 2006 for a review), but most require additional data that is not supplied in the -list of the six quantities. Briefly, some methods involve a best-fit of particle motions to an affine velocity field, but these methods require the particle locations and velocities. Other methods use movements of boundary particles, but these methods also require knowledge of the particles’ locations as well as their velocities. Yet other methods use contact data alone, but require knowledge of the topological adjacencies within the contact network. Liao et al. (1997) have proposed an approximation of that requires only local contact data. For the current work, we use an alternative approach, which also requires only contact data (rather than particle data) but is suitable for integration in the form of Eq. (7). This estimate is derived from an analysis of the relative motions of contacting particle pairs that lie along the perimeter of a two-dimensional granular assembly.

The author has shown that the average velocity gradient within a two-dimensional region is exactly given by a double integral around its perimeter, in which perimeter locations are parameterized by the arc distances and (Kuhn, 2004a). Both distances are measured counterclockwise from a common, fixed point on the boundary (Fig. 3a).

Figure 3: Boundary chain of particles for finding the velocity gradient: (a) bounded assembly of particles, (b) chain of boundary branch vectors, and (c) two boundary particles.

The area-average velocity gradient is


where is the perimeter length; is the unit outward normal of the perimeter at position ; is the enclosed area; is the derivative of the velocity component along the perimeter, as point traverses this boundary; and kernel is a discontinuous function with domain and range given by


using the modulo function. Integral (30) is an objective rate and requires only two types of information around the boundary: the boundary normal and the relative motion along the perimeter, which is related to the contact motions. The enclosed area can also be expressed in a similar manner (Section 2.6). Equation (30) is a general and exact expression for the average spatial gradient within a closed two-dimensional region (bounded by a Jordan curve), provided that the boundary derivative is integrable and consistent with a field of boundary displacements, such that .

For an assembly of disks, we focus on the closed polygonal chain of segments (branch vectors) formed by the contacting disks along the assembly’s perimeter (Fig. 3b). The perimeter segment between particles and has outward unit normal which is perpendicular to the branch vector, , with arc length and rotation tensor that effects a counterclockwise rotation of :


In applying Eq. (30) to critical state flow, we assume that the particles are rigid and that bulk deformation is entirely due to the relative tangential motions between contacting particles. That is, we neglect any small changes in contact indentations during deformation and neglect the opening (extinction) and closing (creation) of contacts during small increments of bulk deformation. Such normal motions have a negligible contribution to the bulk deformation (characteristic A.1 in the Introduction). Considering only tangential motions, the velocity of particle relative to is equal to the angular rate multiplied by branch length and by the unit tangential direction , producing the rate vector (Fig. 3c). This relative velocity occurs along a branch vector of length , so that derivative . Rate can be expressed in terms of “” rates (see Eq. 2), as


The double integration in (30) distinguishes between “” and “” quantities: is associated with contact movement, whereas is associated with the outer normal . In a similar manner, we distinguish between the “c” and “” directions, and , and write (30) as the double sum:


where is the number of perimeter contacts, and and are measured counterclockwise from a single, common perimeter point. This double sum is an exact expression for average strain in a mono-disperse assembly of disks, provided that contact movements are limited to tangential displacement.

To develop a corresponding function that can be integrated as in Eq. (7), we idealize the material region in the limit of a very large assembly of disks in which the perimeter chain of contacts has an orientation that increases monotonically from one contact to the next, forming a convex, nearly smooth boundary. In making this idealization, we neglect the non-convex “inside corners” that would typically occur around a polygonal loop of branch vectors. Arc distances and can now be approximately parameterized with angles and . We also assume that the probability density of these angles is the same as that within the interior of the region: the density . The shape of the boundary chain will depend upon this density, which furnishes the relative numbers of contacts at different orientations. With these assumptions, Eqs. (30), (31), and (34) are written as


with function and kernel ,


and with and defined in Eq. (32) and defined in Eq. (33). Unlike the Dirac kernel ( in Eq. 18), the kernel of Eq. (37) effects the double sum of Eq. (34), in which the two arguments ( and , or and ) correspond to different points on the assembly boundary.

Although Eqs. (30) and (34) yield exact average gradients, Eqs. (35)–(37) are an estimate. The author’s DEM simulations show errors of 5%–25% for disk assemblies in critical state flow. Notwithstanding the small error, the approximation (35) is only intended as a modest constraint on the contact motions, holding them to a nearly isochoric condition. Errors result from a number of sources: (1) the non-tangential (normal) contact movements that can occur between non-rigid grains; (2) the approximation of in Eq. (30) with in Eq. (35), which is a truncation of the full expansion of ; (3) the non-monotonic variation of around a non-convex perimeter chain of branch vectors; and (4) subtle correlations among the branch vector lengths , contact movements , and orientations and .

The isochoric constraint on the probability density is the requirement of a zero expected value of trace , expressed as


in which the leading constant term in (36) is canceled. The dyad in Eq. (38) is


with trace


where is given by Eq. (33), and function of Eq. (15) enforces the biaxial symmetry shown in Fig. 2.

2.6 Fundamental constraint 4: Loading rate

We must specify the deformation rate that was included in the second (dissipation) constraint of Eqs. (27) and (28). Applying Eqs. (35) and (36) to the single deformation component of Eq. (26) yields the following constraint on probability density :


We can eliminate the leading coefficient by again considering the deformation gradient of Eqs. (30), (35), and (36), but for the special case of an ideal dilation, with and , which produces a known deformation rate with trace :


Dividing Eq. (41) by Eq. (42) and rearranging terms gives the fourth constraint on density :


which formally asserts the compression rate , assuring consistency of second and third constraints (Eqs. 27 and 38). The function is computed as


with given by Eq. (33).

2.7 Auxiliary constraint 5: Fraction of sliding contacts

The four constraints given above describe a model of critical state flow under biaxial loading. Although a simple four-constraint model will predict nearly all trends observed in experiments and numerical simulations, the qualitative agreement is, in some respects, in poorer accord: in particular, it over-predicts the fraction of sliding contacts and the activity of contact movements (a prediction of over 80% sliding contacts). These aspects can be improved with the supply of new information to the model — information of an entirely empirical origin. We include a single modest parcel of information. A four-constraint model predicts an excess in the fraction of sliding contacts (more than 80%), compared with the much smaller values noted in characteristic C.1 of the Introduction. The final constraint limits the fraction of sliding contacts :


where a value of is empirically derived.

2.8 Entropy

In the view of an experimentalist or simulator, a granular medium is a deterministic system — although one of bewildering complexity — in which the evolving conditions of each contact (its motions, forces, etc.) are uniquely determined by its own condition, by the conditions of all other contacts in the system, and by the boundary motions and forces. The probabilities are accessible, in this view, by frequent empirical observation of the contacts: from data in the form of the ” values. Jaynes (1957) describes this as an “objective” approach to probability, when probability is an expressed expectation based on observation. To Jaynes, statistical mechanics is based upon an alternative “subjective” school of thought in which probabilities are simply expressions of expectation based upon general information that is usually very limited. Any shortcomings of such predictions, when compared with experimental observation, are attributed by the “subjectivist” not to insufficient acuity but to insufficient information. Indeed, to the subjectivist, the notion of probability is irrelevant when all information of a deterministic system is available beforehand: such an exercise is not one of prediction but of certitude.

From a subjective viewpoint, this lack of information is synonymous with uncertainty, disorder, or “entropy” and is quantified, in our case, with the differential Shannon entropy


where must satisfy the constraints of the previous sections. Although restrictive, this small set of constraints does not, by itself, determine a unique macro-state, as many macro-states will satisfy the same constraints. The most likely macro-state — the most likely probability density — is the one that encompasses the greatest breadth of admissible micro-states among the space of possibilities (see Jaynes, 1957). For a large sample size (i.e. large ), this probability density maximizes while satisfying any available information that restricts the admissible micro-states: information in the form of our five moment constraints


If, instead, a granular system was to settle upon a macro-state of lower entropy, this occurence would indicate a bias toward greater order (a macro-state with a smaller breadth of micro-states), suggesting an influence of other unaccounted information (i.e., other constraints). This possibility is addressed in the conclusions.

Following the Jaynes formalism of maximizing (i.e., the “MaxEnt” principle), the condition


and the zeroth constraint of Eq. (9) lead to the most likely density


with normalizing (partition) function ,


and with five Lagrange multipliers that are computed so that the five moment constraints are satisfied (again, Jaynes, 1957). These multipliers are computed as the solutions of five non-linear equations, requiring a rather taxing evaluation of multiple multi-dimensional integrations, as described in B.

3 DEM Simulations and Model Verification

This section describes the small set of parameters that are required to implement and solve the model (Section 3.1), the DEM simulations that were used to evaluate the model (Section 3.2), and the results of this evaluation (Section 3.3).

3.1 Model implementation

The model requires three input parameters: the inter-particle friction coefficient , the fraction of sliding contacts , and the contact density coefficient . Mean stress and strain rate also appear, but only as scaling factors, and the forces and sliding rates are normalized with respect to these factors. Coefficient is largely a scaling factor that does not significantly affect results and was was assigned a value of 1.5 in all calculations (see the end of Section 2.3). Results depend primarily on the friction coefficient, and a range of values were used (, 0.3, 0.5, 0.7, and 0.9), with most results reported for . The fraction ranged from 5% to 37% in the simulations, depending on friction . An average value of 15% was used with the model.

Equations (47)–(50) define the density that corresponds to maximum disorder, a density expressed in terms of five multipliers . Once the are computed (B), we can extract meaningful information from the model by evaluating various expected values and marginal distributions of . These evaluations and an associated notational system are described in C.

3.2 DEM simulations

Statistics of contact force, movement, and orientation were gathered from simulations of 168 assemblies, each with 676 bi-disperse disks, which were sheared in horizontal biaxial compression. Details of the simulations are found in D. The bulk behavior is illustrated in Fig. 4, which shows the evolution of stress, fabric, and porosity. The Satake fabric tensor is the average of the contact orientation vectors , and the ratio is a fundamental measure of bulk anisotropy (Satake, 1982).

Figure 4: Stress, fabric, and porosity from DEM simulations (friction coefficient ).

A large initial stiffness causes the stress ratio to rise quickly from 1.0 to a peak condition at strain , and the critical state condition is reached at compressive strains beyond 16–18%, with a mean value of ratio of 1.80. While in the critical state, the bulk response is seen to fluctuate (with a standard deviation of 0.03), indicating the agitated, tumultuous nature of the underlying particle interactions.

By running 168 simulations with random initial configurations and by taking several snapshots during the subsequent critical state flow, we collected over 800,000 contact samples of orientations, forces, and movements. Once the critical state is reached, we found that any bulk average of the micro-scale data (stress, fabric, density, etc.) or any marginal distribution thereof attains a nearly constant, steady condition.

3.3 Comparison of model and DEM results

The Introduction recounts numerous micro-scale characteristics of granular flow. We will find strong qualitative agreement with the DEM results: with most characteristics described in the Introduction, the model and simulations exhibit the same trends. In presenting the model results, the following paragraphs are preceded with raised items (e.g., “B.7”) that refer to particular characteristics in the Introduction, and the corresponding references to the literature can be found there. Primary predictions of the entropy model are shown in Table 3, with references to the Introduction shown as raised items in the first column. The third column in Table 3 defines the various quantities, using the notation of C. Only averages are reported, although dispersions from the mean can also be extracted with the model. Unless otherwise noted, the results of the model and the simulations are for the single case .

In the following, the normal contact forces are divided by (i.e., the factor in Eq. 20), so that the mean normed value is 1 (row 6 of Table 3). In terms of an actual force , the normalized equals .

Description Definition (see C) DEM Model
1 Fabric tensor ratio, 1.30 1.42
2 Deviatoric stress ratio, 0.573 0.605
3 Tangential force contribution to the deviatoric stress ratio, 0.068 0.045
4 Normal force contribution to the deviatoric stress ratio, 0.505 0.560
5 Percent contribution of tangential forces to deviatoric stress 11.9% 7.5%
6 Mean normal contact force, normalized 1 1
7 Fraction of deviator stress from weak contacts 7.2% 11.0%
8 Fraction of mean stress from weak contacts 28.7% 26.7%
9 Fabric tensor ratio among weak contacts 1.05 1.17
10 Fabric tensor ratio among strong contacts 1.76 1.80
11 Fraction of sliding contacts, 11.2% 15.0%
12 Fraction of forward sliding contacts, 6.6% 7.7%
13 Fraction of reverse sliding contacts, 4.6% 7.3%
14 Ratio of forward and reverse sliding contacts 1.43 1.05
15 Mean normal contact force among non-sliding contacts 1.073 1.15
16 Mean normal contact force among forward sliding contacts 0.335 0.174
17 Mean normal contact force among reverse sliding contacts 0.306 0.159