Maximum disorder model for dense steady-state flow of granular materials
MMaximum disorder model for dense steady-state flowof granular materials
Matthew R. Kuhn a a Donald P. Shiley School of Engineering, University of Portland, 5000 N. Willamette Blvd.,Portland, Oregon, USA 97203, [email protected], (503)943-7361
Abstract
A flow model is developed for dense shear-driven granular flow. As described in thegeomechanics literature, a critical state condition is reached after sufficient shearingbeyond an initial static packing. During further shearing at the critical state, thestress, fabric, and density remain nearly constant, even as particles are being con-tinually rearranged. The paper proposes a predictive framework for critical stateflow, viewing it as a condition of maximum disorder at the micro-scale. The flowmodel is constructed in a two-dimensional setting from the probability density of themotions, forces, and orientations of inter-particle contacts. Constraints are appliedto this probability density: constant mean stress, constant volume, consistency of thecontact dissipation rate with the stress work, and the fraction of sliding contacts. Thedifferential 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 theunderlying phase space. The resulting distributions of contact force, movement, andorientation are compared with two-dimensional DEM simulations of biaxial compres-sion. The model favorably predicts anisotropies of the contact orientations, contactforces, contact movements, and the orientations of those contacts undergoing slip.The model also predicts the relationships between contact force magnitude andcontact motion. The model is an alternative to affine-field descriptions of granularflow.
Keywords: granular matter, entropy, critical state, fabric, anisotropy Introduction
The critical state concept in geomechanics holds that dense granular materi-als, when loaded beyond an initial static packing, eventually attain a steady statecondition of constant density, fabric, and stress (Schofield and Wroth, 1968). Thiscondition 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,
Preprint submitted to Elsevier December 20, 2018 a r X i v : . [ c ond - m a t . s o f t ] D ec nd physicists, who have devoted great effort in understanding the state’s underlyingmechanics. Density, fabric, and deviatoric stress at the critical state are known to de-pend upon the particles’ shapes and contact properties as well as on the mean stressand intermediate principal stress (Zhao and Guo, 2013). Even so, the eventual bulkcharacteristics for a given assembly are insensitive to the initial particle arrangementand to the stress path that ends in the critical state: for example, materials that areinitially either loose or dense eventually arrive at the same density condition aftersufficient shearing. This convergent characteristic resembles that of thermal systemsthat approach an equilibrium condition with sufficient passage of time.Another pervasive feature of critical state flow is the continual and intense activityof grains at the micro-scale, yet this local tumult produces a monotony in the bulkfabric, stress, and density. Micro-scale activity occurs in three ways: (a) statically,as alterations of the inter-particle contact forces, (b) geometrically, as changes inthe particles’ configuration and local density, and (c) topologically, as changes inthe load-bearing contact network among the particles. In these three respects, weview the critical state as the condition of maximum disorder that emerges duringsustained shearing. Early work by Brown et al. (2000) investigated disorder in thelocal density, and a recent paper by the author explored topological disorder at thecritical state (Kuhn, 2014). The current paper addresses statical disorder as expressedin a two-dimensional (2D) setting, focusing on anisotropies and distributions of thecontacts’ orientations, forces, and movements. The analysis applies to the criticalstate flow of dry unbonded frictional materials of sufficient density to develop aload-bearing (persistently jammed) network of contacts during slow (quasi-static,non-collisional) shearing. Particles are assumed durable (non-breaking) and nearlyrigid, such that deformations of the particles are small, even in the vicinity of theircontacts.During flow, granular materials have an internal organization of movement andforce, an organization that is in some respects pronounced but in others subtle. Webriefly review these characteristics of the critical state, as observed in laboratoryexperiments, numerical simulations, or both.A.1. The motions of individual particles do not conform to an affine, mean deforma-tion field, and fluctuations from the mean field are large and seemingly erratic(Kuhn, 2003). The rates at which contacting particles approach or withdrawfrom each other (i.e. contact movements in their normal directions) are gener-ally much smaller than those corresponding to an affine field. In contrast, thetransverse, tangential movements between contacting pairs are much largerthan those of affine deformation (Kuhn, 2003). As a result, bulk deformation isalmost entirely attributed to the tangential movements of particles (Kuhn andBagi, 2004).A.2. Strength, expressed as a ratio of the principal stresses, is insensitive to the con-tacts’ elastic stiffness and to the mean stress, such that simulations of either2oft 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 (frac-tion of sliding contacts, contact anisotropy, etc.) are also insensitive to contactstiffness (da Cruz et al., 2005).A.3. Particle rotations are large when compared with the bulk deformation rate (Kuhnand Bagi, 2004). In particular, the rolling motions between particles are muchlarger than the sliding movements (Kuhn, 2004b).A.4. The contact network is relatively sparse, in that the number of contacts withinthe load-bearing contact network is sufficient to produce a static indeterminacy(hyperstaticity) but with only a modest excess of contacts (Thornton, 2000). Theexcess very nearly corresponds to the number of contacts that are sliding (Kruytand Antony, 2007). This modest indeterminacy is consistent with observationsof intermittent, sudden reductions of stress, which result from periodic col-lapse events that are occasioned by fresh slip events or loss of contacts (Peñaet al., 2008). This condition of marginal hyperstaticity is referred to as “jammed”within the granular physics community.B.1. During critical state flow, the contact fabric is anisotropic, with the normalsof the contacts oriented predominantly in the direction of the major principalcompressive stress (Rothenburg and Bathurst, 1989).B.2. The normal contact forces are larger among those contacts oriented in the direc-tion of the major principal compressive stress; whereas, the averaged tangentialforces are larger for contact surfaces oblique to the principal stress directions(Rothenburg and Bathurst, 1989; Majmudar and Bhehringer, 2005).B.3. Deviatoric stress is primarily borne by the normal contact forces between parti-cles; whereas, the tangential contact forces make a much smaller contributionto the deviatoric stress (Thornton, 2000).B.4. When considering only the normal forces, deviatoric stress is primarily carried bythose 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).B.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).B.6. Many contacts slide in the “wrong direction” with respect to the direction thatcorresponds to an affine deformation (Kuhn, 2003).B.7. Compared with other orientations, contacts that are oriented in the directionof extension have a greater average magnitude of slip, but the mean slip veloc-ity is largest among contacts that are oriented obliquely to the directions ofcompression and extension (Kuhn and Bagi, 2004).B.8. Frictional sliding is more common among those contacts with a smaller-than-mean normal force (Radjai et al., 1998).3.9. The more mobile contacts — those with large sliding movements — tend to bethose that bear a smaller-than-average normal force (i.e., weak contacts) (Kruytand Antony, 2007).B.10. Deformation, when measured at the meso-scale of particle clusters, is relatedto contact orientation: contacts with branch vectors that are more aligned withthat of bulk compression tend to produce local dilation; whereas, contactsthat are more aligned with the direction of bulk extension tend to produce localcompression. These trends have been determined by studying the elongations ofvoids that are surrounded by rings of particles and their branch vectors (Nguyenet al., 2009).B.11. The probability density of the normal contact forces usually decreases expo-nentially for forces that are greater than the mean (Majmudar and Bhehringer,2005). With forces less than the mean, however, the density is more uniformthan exponential.B.12. Strength at the critical state increases with an increasing inter-granular frictioncoefficient, but the relation is non-linear, and little strength gain occurs whenthe 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).C.1. Relatively few contacts attain the frictional limit at any particular momentduring flow, and sliding typically occurs among only 8%–20% of the contactsin assemblies of disks and spheres (Radjai et al., 1998; Thornton, 2000). Thefraction of sliding contacts is reduced when the friction coefficient is increased(Thornton, 2000).C.2. For contacts that are not sliding, the probability density of the mobilized frictionis greatest for the nearly neutral condition of zero tangential force, and thedensity decreases with increasing mobilized friction (Majmudar and Bhehringer,2005).D.1. Internal force, movement, and deformation are highly heterogeneous and arespatially organized across scales of ten or more particles (Kuhn, 2003). Contactforces are patterned in “force chains” of highly loaded particles that are roughlyaligned with the major principal compressive stress (Majmudar and Bhehringer,2005). Deformation is localized into obliquely oriented microbands of thickness1–3 particle diameters (e.g., Kuhn, 1999) and into shear bands with a thicknessof 8–20 diameters (Desrues and Viggiani, 2004). Local stiffness (Tordesillaset al., 2011) and local dilation and contraction are often clustered (Kuhn, 1999),and chains of rapidly rotating particles are organized obliquely to the majorprincipal stress axes (Kuhn, 1999). The non-affine particle movements in 2Dassemblies are coordinated in large circulating vorticity cells that encompassseveral dozens of particles (Williams and Nabha, 1997). In 2D materials, the4opology of the contact network is also spatially organized: voids surroundedby 3-5 particles form elongated chains, and larger voids are often organized asclusters (Kuhn, 2014).These and other characteristics, as determined from new DEM results, will be dis-cussed and compared in relation to a proposed model of critical state flow. Whereasthe results listed above were gained from descriptive, observational studies, the cur-rent work is largely predictive. The above characteristics are arranged as follows: “A”characteristics serve as the basis in developing the model, “B” characteristics arefavorably predicted by the model, “C” characteristics are not well predicted or mustbe forced into agreement with an ad hoc intervention, and “D” characteristics areinaccessible, due to the nature of the model.The proposed model focuses on the movements, forces, and orientations ofinter-particle contacts during shear-driven flow. The model is premised on the con-jecture that the critical state is a condition of maximum micro-scale disorder amongthese contact characteristics. This maximally disorder state is extracted, and its re-sulting conditions are shown to exhibit many of the behaviors listed above. Withsome reluctance, we use the term “entropy” interchangeably with disorder, withentropy connoting a statistical measure of micro-scale unpredictability (i.e., as inShannon, 1948) rather than an extensive thermodynamic quantity having an inten-sive 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 forseveral decades. Brown et al. (2000) conducted experiments on two-dimensional as-semblies of spheres, and by applying a back-and-forth shearing, found that disorderin the local density and coordination number increased with each shearing cycle.The Jaynes maximum entropy (MaxEnt) formalism has typically been used to predictthe condition of maximum disorder (Jaynes, 1957). This approach has been appliedto the local fabric of disk assemblies by categorizing voids into several canonical types(Brown, 2000). Similar maximum entropy approaches have also been applied to localpacking density (Yoon and Giménez, 2012), to the contact forces (Chakraborty, 2010),to the contact orientations (Troadec et al., 2002), and to the contact displacementsand bulk elastic moduli (Rothenburg and Kruyt, 2009). Among these aspects, the dis-tribution of contact force has recently received the greatest attention. Early theoriesderived probability densities of force distribution by addressing the bulk mean stressbut without respecting the local force equilibrium (Goddard, 2004). More recentapproaches to the distribution of contact force have enforced the local equilibriumof particles and examined the statically admissible space of contact forces within agranular system (Chakraborty, 2010). These disorder theories address frictionless,non-flowing assemblies and do not address the anisotropies of force, fabric, andmovement that accompany flow. The paper not only admits flow but relies upon theshearing deformation to drive the contact movements and the contact forces that5enerate the bulk anisotropies of fabric and force.Our approach also differs from mean field theory, which is currently the dominantparadigm 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 estimatingby assuming that the particle motions or contact forces around a central particleconform to an affine, homogeneous kinematic (motion) or static (stress) field. Wewill not impose such affine restrictions, as doing so would imply a strong order inthe critical state condition — an order that simply does not exist (Kuhn, 2003). Mostmean-field approaches are also inherently deterministic insofar as they result in theunique response of a particle for a given description of its neighborhood. Rothenburgand Kruyt (2009) extended such methods by adopting probability densities for boththe surroundings and the response. The current work also adopts a probabilisticframework, but the approach is intended for granular flow in which the micro-scaleinteractions are dominated by tangential contact motions.The plan of the paper is as follows. We build a micro-scale flow model by firstidentifying an essential set of contact quantities, forming a rather comprehensivephase space of motion, force, and arrangement (Sec. 2.1). We regard these quantitiesas random variables with a bulk probability density distribution (Sec. 2.2) and applyconstraints, both rational and empirical, to this distribution (Secs. 2.3–2.7). By them-selves, these constraints do not describe a unique distribution. We assume that allsuch micro-states satisfying the constraints are equiprobable during critical stateflow and that the most likely macro-state (i.e., density distribution) is the one thatencompasses 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 pre-dicted distribution (Sec. 2.8). Section 3 describes the model’s solution and comparesit with past observations and with the results of new DEM simulations.Although the principles presented in these sections are fundamental, the set ofconstraints sparse, and the results promising, their evaluation is computationallytaxing. Complete computational details, therefore, are provided in appendices. Weend the paper by considering a number of questions raised by some difficulties andshortcomings of the entropy model: which aspects of granular flow are amenable tosimple, sweeping statistical approaches, and which aspects require closer attentionto the details of grain interactions? What additional constraining information is mostappropriate for these entropy methods? How can a model that is conceptually simplebe both predictive and opaque?
2. Entropy Model
The model is developed in a two-dimensional setting of assemblies of disks thatare poly-disperse but of a narrow size range. Because of the assumed small poly-6 able 1:
Contact quantities of the entropy model (Fig. 1)Quantity Description g n ∈ R + Compressive normal contact force, f n / ( p o (cid:96) ) g t a,b Tangential contact force, f t / ( p o (cid:96) ) θ c ∈ [0 , π ) a Orientation of contact normal vector n c θ (cid:96) ∈ [0 , π ) a Orientation of branch vector l (cid:96) ˙ φ slip ∈ R a,b (cid:112) / times contact slip, “slip k ” ˙ φ rigid ∈ R a √ times rigid contact movement, “rigid k ” a Directional conventions are shown in Fig. 2. b g t and ˙ φ slip are also limited by Eq. (4). dispersity, we can overlook the tendency of mono-disperse assemblies to developordered, crystallized (low entropy) particle arrangements. A narrow size range alsopermits characterizing the disk sizes by their mean diameter (cid:96) . The model is intendedfor durable and nearly rigid disks, in which contact indentations are small relative toparticle size: the contacts are assumed rigid-frictional and without contact moments. The foundation of any entropy model is the set of chosen quantities that charac-terize the micro-states that comprise the full phase space of a physical system. Thequantities that one chooses will largely determine the nature of the model and thephenomena that are addressed. For example, a phase space composed exclusively ofparticle arrangements will lead to a random isotropic state. A phase space composedexclusively 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, andarrangement during granular flow, and we choose a rather large, comprehensivephase space of certain contact quantities . These quantitites can be used to extractor to constrain the bulk deformation rate tensor, stress tensor, and fabric tensor,and these micro-quantities are all accessible from simulations of discrete particlesystems.A micro-state of a large two-dimensional assembly of disks is enumerated withsix quantities at each of the assembly’s M loading-bearing contacts, such that amicro-state of the entire assembly is a single point within a phase space of dimension M . The six quantities, listed in Table 1, apply to a contact k that is shared betweentwo particles, i and j , which comprise an ordered pair ( i, j ) k (Fig. 1a). Contactsare assumed to be enduring, non-collisional, and persistent over the brief intervalfor which bulk averages are extracted. Note that, henceforth, we will largely ignorethe association of contacts with particles, thus overlooking the underlying topologyof the contact network and relinquishing an ability to assure local equilibrium orkinematic compatibility. The compressive normal contact force f n k acts at the contact7 igure 1: Contact quantities of contact k between two contacting disks ( i, j ) k . Direction conventionsfor biaxial symmetry are shown in Fig. 2. “ k ” upon particle “ i ” (Fig. 1b). Because the average normal force in a two-dimensionalassembly is roughly proportional to the mean diameter (cid:96) and to the mean stress p o ,we normalize the individual normal forces as g n k = f n k / ( p o (cid:96) ) . The scalar tangentialforce on i is f t k , normalized as g t k = f t k / ( p o (cid:96) ) , and this force component is limited byfriction, as characterized by coefficient µ : | g t k | ≤ µg n k . The sign convention associatedwith f t and g t will be described further below.Only tangential contact movements are modeled, since normal motions (thosethat alter contact indentations or cause the separation or creation of contacts) areknown to contribute negligibly to bulk deformation during flow (characteristics A.1and A.2 in the Introduction). During deformation, particles will both rotate and shiftrelative to each other, and we adopt two systems for describing these movements. Inthe first system, scalar ˙ n k represents an angular rate produced by translation of thecenter of j relative to the center of i — a rotational rate of the unit normal vector n k directed outward from i (Fig. 1c). Note that the rate vector ˙ n k is orthogonal to n k . Therelative motions of the particles at their contact are produced both by this tangentialshifting of the particles’ centers (i.e., the velocity (cid:96) ˙ n k ) and by the two particles’ rota-tions, ˙ θ i and ˙ θ j (Fig. 1c). Motions in the first system are described by the quantities ˙ n k , ˙ θ i , and ˙ θ j . 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 (cid:96) , the mechanisms are asfollows (Kuhn and Bagi, 2005):slip k = ˙ n k − ( ˙ θ i + ˙ θ j ) roll k = ˙ θ i − ˙ θ j (1)rigid k = ( ˙ n + ˙ θ i + ˙ θ j ) − . Because all three mechanism usually occur simultane-ously at any contact, we decompose the motions ˙ n k , ˙ θ i , and ˙ θ j as ˙ φ slip ,k ˙ φ roll ,k ˙ φ rigid ,k = (cid:112) / − (cid:112) / − (cid:112) / (cid:112) / − (cid:112) / (cid:112) / (cid:112) / (cid:112) / ˙ n k ˙ θ i ˙ θ j (2)forming an orthogonal separation of particle movements into three “ ˙ φ ” contact rates.Two quantities, ˙ n k and “slip k ”, are most relevant in the model, as they determine thebulk deformation and dissipation rates. Because both quantities can be expressed aslinear combinations of ˙ φ slip ,k and ˙ φ rigid ,k , the model does not require ˙ φ roll ,k , and thereduction of three quantities to two greatly reduces computational demands.As the final contact quantities, we introduce two orientation angles, θ c and θ (cid:96) ,with angle θ c giving the direction of the contact normal n k , and angle θ (cid:96) representingthe direction of the branch vector l k , of length (cid:96) , that joins the center of i to that of j (Fig. 1a). Although the two angles are identical for disks, they are distinguishedbecause of their different roles in expressing stress and strain within the assembly.Both angles are measured from particle i .Although the model is quite general, it will be applied to the particular conditionof biaxial compression, in which the x width of the assembly is reduced while the x height expands (Fig. 2a). We adopt a sign convention that takes advantage ofthe symmetry of these conditions and assists bookkeeping when the model is latercompared with DEM results. The convention, shown in Fig. 2b, assigns a directionto the tangential unit vector t k corresponding to the sliding that would be expectedif the particle motions roughly conformed to affine deformation during horizontalbiaxial compression. This direction — clockwise or counterclockwise — alternatesacross quadrants of θ . Unit vector t k applies to the (scalar) tangential force on i (i.e., f t k and g t k ) and to the relative contact movements ( ˙ θ , ˙ φ slip ,k , etc.). Positive values of f t k , ˙ φ slip ,k , etc., which are aligned with t k , are called “forward” forces and movements;whereas, negative values are in the “reverse” direction − t k . Normal force f n k is alwayscompressive, and unit vector n k is directed outward from particle i .The domains of the six quantities are shown in Table 1. We also place unilateraland rigid-frictional constitutive limits upon the model, g n k ∈ R + (3) g t k ∈ − µg n k ⇔ ˙ φ slip ,k < − µg n k , µg n k ) ⇔ ˙ φ slip ,k = 0 µg n k ⇔ ˙ φ slip ,k > (4)further restricting the domains of g n k , g t k , and ˙ φ slip ,k . These restrictions mean that theparticles are unbonded ( f n k , g n k > ), and the tangential force f t k is limited by friction9 igure 2: Conventions for unit vectors n and t and rotations ˙ θ when analyzing horizontal biaxialcompression. “Forward” directions are shown. to the range [ − µf n k , µf n k ] . Positive (forward) slip, ˙ φ slip ,k > , is only admissible witha corresponding positive tangential force, f t k = µf n k ; whereas, negative (reverse) slipis only possible for the opposite condition, f t k = − µf n k . When the friction limit isnot reached, | f t k | < µf t k , the slip is zero: ˙ φ slip ,k = 0 . In adopting Eqs. (3) and (4), weforego a more elaborate elastic-frictional model, since critical state flow is insensitiveto contact elasticity (characteristic A.2 of the Introduction). Note that the model israte-independent, since forces f t k depend only on the directions of ˙ φ slip ,k and not ontheir magnitudes. In the model, we do not individualize the six “ k ” quantities among the M contacts(a phase space of dimension M ) but instead treat the quantities as continuousrandom variables with probability density p ( · · · ) : p ( · · · ) = p ( g n , g t , θ c , θ (cid:96) , ˙ φ slip , ˙ φ rigid ) (5)dropping the k subscripts. Although we write the complete phase space of possible micro-states as { g n k , g t , θ c k , θ (cid:96)k , ˙ φ slip ,k , ˙ φ rigid ,k } , all micro-states within this M spacethat share the same probability density p ( · · · ) comprise a common macro-state , an el-ement in a subspace of R and written without parentheses: ( g n , g t , θ c , θ (cid:96) , ˙ φ slip , ˙ φ rigid ) or ( · · · ) . As will be seen, density p contains comprehensive information about bulkquantities, such as stress and fabric, and about the correlations among contact ori-entation, force, and movement. However, by turning from individual “ k ” contactquantities to their gross probability distribution, we forfeit the possibility of tracking10ndividual contacts. In return, we gain a certain economy of expression for predict-ing overall distributions of the six quantities and the statistical relationships amongthem.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 enu-merated in the Introduction. These correlations are expressions of order, which weconsider a consequence of certain constraints that link the micro and bulk behaviors.We will restrict the admissible macro-states with five constraints: four fundamentalconstraints and one auxiliary constraint. These constraints bring a priori informa-tion that will bias the distribution p ( · · · ) . Each “ i ”th constraint is expressed with aconstraining function Γ i ( · · · ) of the six random variables, Γ i ( g n , g t , θ c , θ (cid:96) , ˙ φ slip , ˙ φ rigid ) or simply Γ i ( · · · ) (6)The expected value (cid:104) Γ i (cid:105) of the function is found through 6-fold integration acrossthe full phase space, with each random variable integrated over the range given inTable 1: (cid:104) Γ i ( · · · ) (cid:105) = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) Γ i ( · · · ) p ( · · · ) (7)and measure dg n dg t dθ c dθ (cid:96) d ˙ φ slip d ˙ φ rigid is implied in these integrations. The integra-tions are described in greater detail in Appendix A.The expected values of functions Γ i ( · · · ) will be constrained to certain averagevalues, Constraint i ⇒ (cid:104) Γ i ( · · · ) (cid:105) = Γ i (8)with the prescribed averages, Γ i , developed below. That is, each restriction (8) is amoment constraint on density p ( · · · ) .An essential “zeroth” constraint follows from the certainty that the six quantitieslie within their full ranges: (cid:104) Γ ( · · · ) (cid:105) = Γ , where Γ ( · · · ) = 1 and Γ = 1 (9)such that integration of p ( · · · ) over its full domain equals 1.The four fundamental constraints express general principles of isochoric dis-sipative flow at constant stress and are derived in Sections 2.3–2.6. The auxiliaryconstraint applies additional information gained from DEM simulations (Section 2.7).The constraints are summarized in Table 2. We model steady state flow under constant mean stress p o , a condition commonlyapplied in geotechnical testing. The Cauchy stress in a large granular assembly is the11 able 2: Summary of constraints in entropy model.Constraint i Description Variables in Γ i ( · · · ) Eq.Fundamental 1 Mean stress = p o ( g n , · , · , · , · , · ) (20)2 Dissipation consistency ( g n , g t , θ c , · , ˙ φ slip , · ) (27)3 Isochoric flow ( · , · , θ c , θ (cid:96) , ˙ φ slip , ˙ φ rigid ) (38)4 Biaxial rate, ˙ ε ( · , · , θ c , θ (cid:96) , ˙ φ slip , ˙ φ rigid ) (43)Auxiliary 5 Fraction of sliding contacts, η ( g n , g t , · , · , · , · ) (45) volume average of the contact dyads l (cid:96)k ⊗ f c k , σ = 1 A M int (cid:88) k =1 l (cid:96)k ⊗ f c k (10)in which branch vector l (cid:96)k joins particle i to j , contact force f c k acts upon i , and M int is the number of contacts within the interior of the two-dimensional region of area A (Rothenburg, 1980). The branch vector length is approximated as (cid:96) , so that l (cid:96)k ≈ (cid:96) n (cid:96)k (11)and the contact force is the sum of its tangential and (compressive) normal compo-nents: f c k = − f n k n c k + f t k t c k (12) = (cid:96) p o (cid:0) − g n k n c k + g t k t c k (cid:1) (13)In these and future expressions, n c and n (cid:96) are unit orientation vectors associatedwith the contact normal and branch vectors: n c = (cid:20) cos θ c sin θ c (cid:21) , n (cid:96) = (cid:20) cos θ (cid:96) sin θ (cid:96) (cid:21) (14)For the particular sign convention shown in Fig. 2b, unit vector t c is t c = K ( θ c ) (cid:20) − sin θ c cos θ c (cid:21) , K ( θ c ) = sgn(cos θ c sin θ c ) (15)where the signum function sgn () ∈ {− , } alternates in sign between coordinatequadrants.We replace the sum in Eq. (10) with an integration by substituting Eqs. (11)and (13), multiplying by probability p ( · · · ) and the number of contacts M int , andintegrating across the full domain of the six quantities ( · · · ) :12 = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) Γ σ ( · · · ) p ( · · · ) (16) Γ σ = (cid:16) M int (cid:96) p o /A (cid:17) K ( θ c , θ (cid:96) ) n (cid:96) ⊗ ( − g n n c + g t t c ) (17)In this expression, we include a kernel function K , which is simply the Dirac operator, K ( θ c , θ (cid:96) ) = δ ( θ c − θ (cid:96) ) (18)The Dirac kernel formally dictates the correspondence of the unit branch vector andthe contact force vector for circular particles (by extension, unit vectors n (cid:96) and n c ),as in the single sum of Eq. (10).Equations (16)–(18) will later be used to extract the predicted stress tensor duringcritical state flow, but we now use them to enforce the first fundamental constrainton the probability density p ( · · · ) : the negative mean stress − tr ( σ ) must equal theassigned pressure p o . The negative mean value of tensor function Γ σ is − tr ( Γ σ ) = 12 M int (cid:96) p o A K ( θ c , θ (cid:96) ) g n (19)and equating its integration in Eq. (16) with p o leads to the first fundamental con-straint (cid:104) Γ ( · · · ) (cid:105) = Γ , where Γ ( · · · ) = K ( θ c , θ (cid:96) ) g n and Γ = 2 (cid:16) M int (cid:96) /A (cid:17) − (20)in which we have divided by the leading constant in Eq. (19). Equation (20) simplyrequires the dimensionless mean normal force (cid:104) g n (cid:105) (i.e., (cid:104) f n (cid:105) /(cid:96)p o ) to equal the di-mensionless quantity Γ , thus enforcing pressure p o . The contact density M int (cid:96) /A is about 1.5, a value that can be found with DEM simulations or estimated fromcharacteristic A.4 of the Introduction. During critical state flow at constant stress and volume, deformation is entirelyplastic: all work done by the applied stress is dissipated through internal irreversibleprocesses. We assume that the only available dissipation mechanism is frictionalsliding between particles. This assumption is the basis of the second fundamentalconstraint on the probability density p ( · · · ) and provides the essential link amongcontact motions, contact forces, bulk deformation, and bulk stress. This constraintdrives the anisotropies in force and movement.The frictional dissipation rate at the contact of two disks is the product of itstangential contact force f t k = (cid:96)p o g t k and the relative slip of the two particles’ surfaces13t their contact: the “slip k ” rate in Eq. (1 ). Applying Eq. (2), this slip rate correspondsto ˙ φ slip as slip k = (cid:112) / φ slip ,k (21)Noting that g t and ˙ φ slip ,k must conform to the constitutive rigid-frictional restrictionsof Eqs. (3) and (4), the frictional dissipation rate per unit of area isDissipation = (cid:96) p o A M int (cid:88) k =1 (cid:112) / g t ˙ φ slip ,k (22)or, in terms of probability density p ( · · · ) ,Dissipation = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) Γ f ( · · · ) p ( · · · ) (23) Γ f = M int (cid:96) p o A (cid:112) / K ( θ c , θ (cid:96) ) g t ˙ φ slip (24)where the Dirac kernel K is applied again, enforcing the coincidence of the contactforces (associated with contact orientations θ c ) and the contact motions (associatedwith branch vector orientations θ (cid:96) ) in the single sum of Eq. (22).The internal work rate of the Cauchy stress σ is the inner productWork = σ : D (25)where D is the symmetric part of velocity gradient L . Although Eq. (16) is an expres-sion of σ , and an expression for L is derived in the next section (as Eq. 35), directlysubstituting these expressions into the full inner product of Eq. (25) will lead to anon-linear constraint on probability p ( · · · ) , making its evaluation intractable. Thesituation is greatly simplified when deformation can be expressed with a single non-zero parameter. We consider the case of biaxial compression, for which D is reducedto D = (cid:20) − ˙ ε
00 ˙ ε (cid:21) (26)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 constrainton density p ( · · · ) for critical state flow in biaxial compression: (cid:104) Γ ( · · · ) (cid:105) = Γ , where Γ ( · · · ) = Γ f + ˙ ε (Γ σ, − Γ σ, ) and Γ = 0 (27)14n which we have eliminated the leading term M int (cid:96) p o /A that appears in Eqs. (17)and (24). By substituting Eqs. (14), (15), (17), and (24), the function Γ ( · · · ) can beexpanded as Γ ( · · · ) = K ( θ c , θ (cid:96) ) (cid:104)(cid:112) / φ slip g t (28) + ˙ εg n ( − cos θ (cid:96) cos θ c + sin θ (cid:96) sin θ c ) − K ( θ c ) ˙ ε g t (cos θ (cid:96) sin θ c + sin θ (cid:96) cos θ c ) (cid:105) which enforces consistency of work and dissipation. The third fundamental constraint (and, by far, the most difficult to derive) en-forces the constant density (isochoric) condition of critical state flow. If tensor L ij = ∂v i /∂x j is the average, bulk velocity gradient within a granular assembly,its trace vanishes during steady state flow. To enforce this condition in an assemblywith M contacts, we require a means of estimating L from the contact motions ofparticle pairs: from an M -list of contact information { g n k , g t k , θ c k , θ (cid:96)k , ˙ φ slip ,k , ˙ φ rigid ,k } .Several methods have been proposed for estimating L (see Bagi, 2006 for a review),but most require additional data that is not supplied in the M -list of the six quanti-ties. Briefly, some methods involve a best-fit of particle motions to an affine velocityfield, but these methods require the particle locations and velocities. Other methodsuse movements of boundary particles, but these methods also require knowledgeof the particles’ locations as well as their velocities. Yet other methods use contactdata alone, but require knowledge of the topological adjacencies within the contactnetwork. Liao et al. (1997) have proposed an approximation of L that requires onlylocal contact data. For the current work, we use an alternative approach, which alsorequires only contact data (rather than particle data) but is suitable for integration inthe form of Eq. (7). This estimate is derived from an analysis of the relative motions ofcontacting particle pairs that lie along the perimeter of a two-dimensional granularassembly.The author has shown that the average velocity gradient L within a two-dimen-sional region is exactly given by a double integral around its perimeter, in whichperimeter locations are parameterized by the arc distances s p and s q (Kuhn, 2004a).Both distances are measured counterclockwise from a common, fixed point on theboundary (Fig. 3a). The area-average velocity gradient is15 igure 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. L = 1 A (cid:90) (cid:90) A ∂v i ∂x j dA (29) = 1 A (cid:90) S (cid:90) S dv i ds p Q ( s p , s q ) m qj ds p ds q (30)where S is the perimeter length; m q is the unit outward normal of the perimeter atposition s q ; A is the enclosed area; dv i /ds p is the derivative of the velocity compo-nent v i along the perimeter, as point s p traverses this boundary; and kernel Q is adiscontinuous function with domain [0 , S ) and range ( − / , / given by Q ( s p , s q ) = − S mod ( s q − s p , S ) , s q (cid:54) = s p , s q = s p (31)using the modulo mod ( ) function. Integral (30) is an objective rate and requires onlytwo types of information around the boundary: the boundary normal m q and the relative motion dv i /ds p along the perimeter, which is related to the contact motions.The enclosed area A can also be expressed in a similar manner (Section 2.6). Equa-tion (30) is a general and exact expression for the average spatial gradient withina closed two-dimensional region (bounded by a Jordan curve), provided that theboundary derivative dv i /ds is integrable and consistent with a field of boundarydisplacements, such that (cid:82) S ( dv i /ds ) ds = 0 .For an assembly of disks, we focus on the closed polygonal chain of M per segments(branch vectors) formed by the contacting disks along the assembly’s perimeter(Fig. 3b). The perimeter segment k between particles i and j has outward unit normal m qk which is perpendicular to the branch vector, m qk = − R · n (cid:96)k , with arc length16 ji ds q = (cid:96) and rotation tensor R that effects a counterclockwise rotation of π/ : R = (cid:20) −
11 0 (cid:21) (32)In applying Eq. (30) to critical state flow, we assume that the particles are rigidand that bulk deformation is entirely due to the relative tangential motions betweencontacting particles. That is, we neglect any small changes in contact indentationsduring deformation and neglect the opening (extinction) and closing (creation) ofcontacts during small increments of bulk deformation. Such normal motions have anegligible contribution to the bulk deformation (characteristic A.1 in the Introduc-tion). Considering only tangential motions, the velocity of particle j relative to i isequal to the angular rate ˙ n k multiplied by branch length (cid:96) and by the unit tangentialdirection t c k , producing the rate vector ˙ n k (cid:96) t c k (Fig. 3c). This relative velocity occursalong a branch vector of length (cid:96) , so that derivative dv i /ds p = ˙ n k t c k . Rate ˙ n k can beexpressed in terms of “ φ ” rates (see Eq. 2), as ˙ n k = (cid:16) φ slip ,k + √ φ rigid ,k (cid:17)(cid:46) √ (33)The double integration in (30) distinguishes between “ p ” and “ q ” quantities: s p isassociated with contact movement, whereas s q is associated with the outer normal m q . In a similar manner, we distinguish between the “c” and “ (cid:96) ” directions, t c k and n (cid:96)k , and write (30) as the double sum: L = − (cid:96) A M per (cid:88) p =1 M per (cid:88) q =1 ˙ n p Q ( s p , s q ) t c p ⊗ (cid:16) R · n (cid:96)q (cid:17) (34)where M per is the number of perimeter contacts, and s p and s q are measured coun-terclockwise from a single, common perimeter point. This double sum is an exactexpression for average strain in a mono-disperse assembly of disks, provided thatcontact 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 whichthe perimeter chain of contacts has an orientation θ that increases monotonicallyfrom one contact to the next, forming a convex, nearly smooth boundary. In makingthis idealization, we neglect the non-convex “inside corners” that would typicallyoccur around a polygonal loop of branch vectors. Arc distances s p and s q can nowbe approximately parameterized with angles θ c and θ (cid:96) . We also assume that theprobability density of these angles is the same as that within the interior of the region:the density p ( · , · , θ c , θ (cid:96) , · , · ) . The shape of the boundary chain will depend upon thisdensity, which furnishes the relative numbers of contacts at different orientations.17ith these assumptions, Eqs. (30), (31), and (34) are written as L ≈ (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) Γ L ( · · · ) p ( · · · ) (35)with function Γ L and kernel K ( θ c , θ (cid:96) ) , Γ L ( · · · ) = − (cid:0) M per (cid:96) (cid:1) A K ( θ c , θ (cid:96) ) ˙ n t c ⊗ ( R · n (cid:96) ) (36) K ( θ c , θ (cid:96) ) = − π mod( θ (cid:96) − θ c , π ) , θ (cid:96) (cid:54) = θ c , θ (cid:96) = θ c (37)and with n (cid:96) and R defined in Eq. (32) and ˙ n defined in Eq. (33). Unlike the Dirackernel ( K in Eq. 18), the kernel K of Eq. (37) effects the double sum of Eq. (34), inwhich the two arguments ( θ c and θ (cid:96) , or s p and s q ) correspond to different points onthe assembly boundary.Although Eqs. (30) and (34) yield exact average gradients, Eqs. (35)–(37) are anestimate. The author’s DEM simulations show errors of 5%–25% for disk assembliesin critical state flow. Notwithstanding the small error, the approximation (35) isonly intended as a modest constraint on the contact motions, holding them to anearly 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 ds in Eq. (30) with M per (cid:96)p ( · · · ) dθ in Eq. (35), which is atruncation of the full expansion of s ( θ ) ; (3) the non-monotonic variation of θ arounda non-convex perimeter chain of branch vectors; and (4) subtle correlations amongthe branch vector lengths (cid:96) , contact movements ˙ n , and orientations θ c and θ (cid:96) .The isochoric constraint on the probability density p ( · · · ) is the requirement of azero expected value of trace tr ( L ) , expressed as (cid:104) Γ ( · · · ) (cid:105) = Γ , where Γ ( · · · ) = K ( θ c , θ (cid:96) ) tr (cid:104) ( ˙ n t c ) ⊗ ( R · n (cid:96) ) (cid:105) and Γ = 0 (38)in which the leading constant term in (36) is canceled. The dyad in Eq. (38) is ˙ n t c ⊗ ( R · n (cid:96) ) = K ( θ c ) ˙ n (cid:20) sin θ c sin θ (cid:96) − sin θ c cos θ (cid:96) − cos θ c sin θ (cid:96) cos θ c cos θ (cid:96) (cid:21) (39)with trace tr (cid:104) ˙ n t c ⊗ ( R · n (cid:96) ) (cid:105) = K ( θ c ) ˙ n (cid:16) sin θ c sin θ (cid:96) + cos θ c cos θ (cid:96) (cid:17) (40)where ˙ n is given by Eq. (33), and function K ( θ c ) of Eq. (15) enforces the biaxialsymmetry shown in Fig. 2. 18 .6. Fundamental constraint 4: Loading rate We must specify the deformation rate ˙ ε that was included in the second (dis-sipation) constraint of Eqs. (27) and (28). Applying Eqs. (35) and (36) to the singledeformation component L = − ˙ ε of Eq. (26) yields the following constraint onprobability density p ( · · · ) : − (cid:0) M per (cid:96) (cid:1) A (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) K ( θ c , θ (cid:96) ) ˙ n (cid:104) t c ⊗ ( R · n (cid:96) ) (cid:105) , = − ˙ ε (41)We can eliminate the leading coefficient ( M per (cid:96) ) /A by again considering the defor-mation gradient of Eqs. (30), (35), and (36), but for the special case of an ideal dilation,with dv /ds p = n p and dv /ds p = n p , which produces a known deformation ratewith trace tr ( L ) = 2 : − (cid:0) M per (cid:96) (cid:1) A (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) K ( θ c , θ (cid:96) ) tr (cid:16) n c ⊗ ( R · n (cid:96) ) (cid:17) = 2 (42)Dividing Eq. (41) by Eq. (42) and rearranging terms gives the fourth constraint ondensity p ( · · · ) : (cid:104) Γ ( · · · ) (cid:105) = Γ , where Γ ( · · · ) = K ( θ c , θ (cid:96) ) (cid:26) ˙ ε tr (cid:104) n c ⊗ ( R · n (cid:96) ) (cid:105) + ˙ n (cid:104) t c ⊗ ( R · n (cid:96) ) (cid:105) , (cid:27) and Γ = 0 (43)which formally asserts the compression rate L = − ˙ ε , assuring consistency ofsecond and third constraints (Eqs. 27 and 38). The function Γ ( · · · ) is computed as Γ ( · · · ) = K ( θ c , θ (cid:96) ) (cid:18) ˙ ε K ( θ c )(sin θ c cos θ (cid:96) − cos θ c sin θ (cid:96) ) + ˙ n sin θ c sin θ (cid:96) (cid:19) (44)with ˙ n given by Eq. (33). The four constraints given above describe a model of critical state flow underbiaxial loading. Although a simple four-constraint model will predict nearly all trends observed in experiments and numerical simulations, the qualitative agreement is, insome respects, in poorer accord: in particular, it over-predicts the fraction of slidingcontacts and the activity of contact movements (a prediction of over 80% slidingcontacts). These aspects can be improved with the supply of new information to the19odel — information of an entirely empirical origin. We include a single modestparcel of information. A four-constraint model predicts an excess in the fraction ofsliding contacts (more than 80%), compared with the much smaller values notedin characteristic C.1 of the Introduction. The final constraint limits the fraction ofsliding contacts η : (cid:104) Γ ( · · · ) (cid:105) = Γ , where Γ = η and Γ ( · · · ) = (cid:40) , g t ∈ {− µg n , µg n } , g t ∈ ( − µg n , µg n ) (45)where a value of η is empirically derived. In the view of an experimentalist or simulator, a granular medium is a deter-ministic system — although one of bewildering complexity — in which the evolvingconditions of each contact (its motions, forces, etc.) are uniquely determined byits own condition, by the conditions of all other contacts in the system, and by theboundary motions and forces. The probabilities p ( · · · ) are accessible, in this view, byfrequent empirical observation of the contacts: from data in the form of the M “ k ”values. Jaynes (1957) describes this as an “objective” approach to probability, whenprobability is an expressed expectation based on observation. To Jaynes, statisticalmechanics is based upon an alternative “subjective” school of thought in which prob-abilities are simply expressions of expectation based upon general information that isusually very limited. Any shortcomings of such predictions, when compared with ex-perimental observation, are attributed by the “subjectivist” not to insufficient acuitybut to insufficient information. Indeed, to the subjectivist, the notion of probabilityis 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 un-certainty, disorder, or “entropy” and is quantified, in our case, with the differentialShannon entropy H ( p ( · · · )) = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) p ( · · · ) ln ( p ( · · · )) (46)where p ( · · · ) must satisfy the constraints of the previous sections. Although restric-tive, 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 p ( · · · ) — is the one that encompasses the greatestbreadth of admissible micro-states among the M space of possibilities (see Jaynes,1957). For a large sample size (i.e. large M ), this probability density maximizes H (cid:104) Γ i ( · · · ) (cid:105) = Γ i , i = 1 , , . . . , (47)If, instead, a granular system was to settle upon a macro-state of lower entropy, thisoccurence would indicate a bias toward greater order (a macro-state with a smallerbreadth 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 H (i.e., the “MaxEnt” principle),the condition ∂H ( p ( · · · )) ∂p ( · · · ) = 0 (48)and the zeroth constraint of Eq. (9) lead to the most likely density p ( · · · ) = 1 Z ( · · · ) exp (cid:32) − (cid:88) i =1 λ i Γ i ( · · · ) (cid:33) (49)with normalizing (partition) function Z , Z ( · · · ) = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) exp (cid:32) − (cid:88) i =1 λ i Γ i ( · · · ) (cid:33) (50)and with five Lagrange multipliers λ i that are computed so that the five momentconstraints are satisfied (again, Jaynes, 1957). These multipliers are computed asthe solutions of five non-linear equations, requiring a rather taxing evaluation ofmultiple multi-dimensional integrations, as described in Appendix B.
3. DEM Simulations and Model Verification
This section describes the small set of parameters that are required to implementand solve the model (Section 3.1), the DEM simulations that were used to evaluatethe model (Section 3.2), and the results of this evaluation (Section 3.3).
The model requires three input parameters: the inter-particle friction coefficient µ , the fraction of sliding contacts η , and the contact density coefficient ( M int (cid:96) /A ) − .Mean stress p o and strain rate ˙ ε also appear, but only as scaling factors, and theforces and sliding rates are normalized with respect to these factors. Coefficient ( M int (cid:96) /A ) − is largely a scaling factor that does not significantly affect results andwas was assigned a value of 1.5 in all calculations (see the end of Section 2.3). Results21 abricStress PorosityStrain, − ε P o r o s i t y S t r e ss r a t i o , σ / σ F a b r i c r a t i o , F / F Figure 4: Stress, fabric, and porosity from DEM simulations (friction coefficient µ = 0 . ). depend primarily on the friction coefficient, and a range of values were used ( µ = 0 . ,0.3, 0.5, 0.7, and 0.9), with most results reported for µ = 0 . . The fraction η rangedfrom 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 p ( · · · ) that corresponds to maximumdisorder, a density expressed in terms of five multipliers λ i . Once the λ i are computed(Appendix B), we can extract meaningful information from the model by evaluatingvarious expected values and marginal distributions of p ( · · · ) . These evaluations andan associated notational system are described in Appendix C. Statistics of contact force, movement, and orientation were gathered from sim-ulations of 168 assemblies, each with 676 bi-disperse disks, which were sheared inhorizontal biaxial compression. Details of the simulations are found in AppendixD. The bulk behavior is illustrated in Fig. 4, which shows the evolution of stress,fabric, and porosity. The Satake fabric tensor F ij is the average (cid:104) n c i n c j (cid:105) of the con-tact orientation vectors n c , and the ratio F /F is a fundamental measure of bulkanisotropy (Satake, 1982). A large initial stiffness causes the stress ratio σ /σ torise quickly from 1.0 to a peak condition at strain − ε = 2% , and the critical statecondition is reached at compressive strains beyond 16–18%, with a mean value ofratio σ /σ 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 theunderlying particle interactions.By running 168 simulations with random initial configurations and by takingseveral snapshots during the subsequent critical state flow, we collected over 800,000contact samples of orientations, forces, and movements. Once the critical state isreached, 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.22 odelDEM data x x Contact orientation, θ P r o b a b ili t y d e n s i t y o f θ ◦ ◦ ◦ ◦ Figure 5: Density distribution of contact orientation for horizontal biaxial compression in the x direc-tion, D θ c (cid:104) (cid:105) . Orientation θ = 0 corresponds to contact normals that are aligned with the compressiondirection ( x in Fig. 2). The Introduction recounts numerous micro-scale characteristics of granular flow.We will find strong qualitative agreement with the DEM results: with most characteris-tics described in the Introduction, the model and simulations exhibit the same trends.In presenting the model results, the following paragraphs are preceded with raiseditems (e.g., “B.7”) that refer to particular characteristics in the Introduction, and thecorresponding references to the literature can be found there. Primary predictions ofthe entropy model are shown in Table 3, with references to the Introduction shownas raised items in the first column. The third column in Table 3 defines the variousquantities, using the notation of Appendix C. Only averages are reported, althoughdispersions from the mean can also be extracted with the model. Unless otherwisenoted, the results of the model and the simulations are for the single case µ = 0 . .In the following, the normal contact forces g n are divided by M int (cid:96) /A ) − (i.e.,the Γ factor in Eq. 20), so that the mean normed value (cid:104) g n (cid:105) normed is 1 (row 6 of Table 3).In terms of an actual force f n , the normalized g nnormed equals f n ( M int (cid:96)/ p o A ) . ( B.1 ) As the most telling result, the model predicts an anisotropy of contact orien-tation that is consistent with the simulations: contacts are predominantly orientedin the direction of the major principal compressive stress, with a Satake fabric ratio F /F = (cid:104) n n (cid:105) / (cid:104) n n (cid:105) that is greater than one (row 1 of Table 3). The distributionof contact orientations is shown in Fig. 5, in which the [0 , π ) range has been foldedto [0 , π/ in accordance with the biaxial symmetry of loading. In this and otherfigures, the x direction (i.e., θ = 0 ) refers to contact normals that are oriented in the(horizontal) direction of compressive loading, and x is in the direction of extension(Fig. 2). Although the model and simulations display some differences, the modelcorrectly predicts the general trend of the orientation distribution. ( B.2 ) The model predicts that stress in the compression direction, − σ , is largerthan that in the extension direction, − σ , and the model predicts a deviator stress23 able 3: Results of DEM simulations and model (friction coefficient µ = 0 . ). Raised items in the firstcolumn refer to characteristics in the Introduction. ValuesDescription Definition(see Appendix C) DEM Model1 (B.1)
Fabric tensor ratio, F /F (cid:104) n n (cid:105) / (cid:104) n n (cid:105) (B.2) Deviatoric stress ratio, q/p o (cid:104) Γ σ, − Γ σ, (cid:105) /p o (B.3) Tangential force contribution tothe deviatoric stress ratio, q t /p o ∗ α (cid:104) g t ( n t − n t ) (cid:105) /p o (B.3) Normal force contribution to thedeviatoric stress ratio, q n /p o ∗ α (cid:104)− g n ( n n − n n ) (cid:105) /p o (B.3) Percent contribution of tangentialforces to deviatoric stress q t /q (cid:104) g n (cid:105) normed (B.4) Fraction of deviator stress fromweak contacts (cid:104) Γ σ, − Γ σ, (cid:105)| g n < (B.4) Fraction of mean stress from weakcontacts (cid:104) Γ σ, + Γ σ, (cid:105)| g n < (B.5) Fabric tensor ratio among weakcontacts (cid:104) n n (cid:105)| g n < / (cid:104) n n (cid:105)| g n < (B.5) Fabric tensor ratio among strongcontacts (cid:104) n n (cid:105)| g n > / (cid:104) n n (cid:105)| g n > (C.1) Fraction of sliding contacts, η (cid:104) (cid:105)| g t ∈{− µg n ,µg n } (B.6) Fraction of forward sliding con-tacts, η + (cid:104) (cid:105)| g t = µg n (B.6) Fraction of reverse sliding contacts, η − (cid:104) (cid:105)| g t = − µg n (B.6) Ratio of forward and reverse slidingcontacts η + /η − (B.8) Mean normal contact force amongnon-sliding contacts (cid:104) g n | g t ∈ ( − µg n , µg n ) (cid:105) normed (B.8) Mean normal contact force amongforward sliding contacts (cid:104) g n | g t = µg n (cid:105) normed (B.8) Mean normal contact force amongreverse sliding contacts (cid:104) g n | g t = − µg n (cid:105) normed ∗ α = M int (cid:96) p o /A odelDEM data x x Contact orientation, θ A v e r a g e n o r m a l f o r ce , g n , n o r m e d ◦ ◦ ◦ ◦ ModelDEM data x x Contact orientation, θ A v e r a g e t a n g e n t f o r ce , g t , n o r m e d ◦ ◦ ◦ ◦ Figure 6: Average contact force as a function of contact orientation: (a) average normal contact force, (cid:104) g n (cid:105) ( θ c ) , and (b) average tangential contact force, (cid:104) g t (cid:105) ( θ c ) . ratio, q/p o = − ( σ − σ ) /p o that is close to the value from the DEM simulations(Table 3, row 2). As in Eqs. (10), (16), and (17), deviatoric stress can result from threesources: from contact orientations biased in the compression direction, from largernormal contact forces in this same direction, and from the tangential contact forces.The predicted anisotropy of contact orientation was verified above (Fig. 5 and row 1).The model also predicts an anisotropy of the average normal forces, as is illustratedin Fig. 6a. The results of the simulations and model are similar: normal forces are, onaverage, larger in the direction of bulk compression. Anisotropy of tangential contactforce is shown in Fig. 6b. For both simulations and model, the averaged tangentialforce is largest at an angle of about 40 ◦ oblique to the compression direction. Althoughthe numerical values of simulations and theory differ, the model captures the primarytrend of anisotropy in the tangential force. ( B.3 ) The separate contributions of the normal and tangential forces to the fulldeviatoric stress can be extracted with the model, as these contributions are propor-tional to the expected values (cid:104)− g n ( n n − n n ) (cid:105) and (cid:104) g t ( n t − n t ) (cid:105) . For bothsimulations and model, the relative contribution of the tangential forces to the fulldeviator stress q is minor: a 11.9% contribution in the simulations and 8.9% from themodel (rows 3–5, Table 3). ( B.11 ) The model gives a fairly good prediction of the probability density of normalcontact forces, as seen in Fig. 7a. The model predicts an exponential tail for largenormal forces g n (although somewhat less steep than the DEM data), and it predictsa curved, flattened “shoulder”for normal forces g n less than 1.5. The model, however,does give a steep rise in the distribution for the smallest normal forces (those with g n < . ). ( C.2 ) The model poorly predicts the density of the tangential contact forces (seethe plot of g t /µg n in Fig. 7b). When considering only the non-sliding contacts (with | g t | < µg n ), the density of this ratio in the DEM simulations decreases with an25 odelDEM data Normal contact force, g n , normed P r o b a b ili t y d e n s i t y o f g n DEM data
Forward g t →← Reverse g t Mobilized contact friction, g t / ( µ g n ) P r o b a b ili t y d e n s i t y o f ψ Figure 7: Density distributions of contact force: (a) normal force distribution, D g n (cid:104) (cid:105) , and (b) distributionof the fraction of mobilized friction of tangential forces, D g t /µg n (cid:104) (cid:105) . increasing mobilized friction, consistent with experimental results (Majmudar et al.,2007). The model, however, gives a nearly uniform density for this same condition.This shortcoming is discussed in the conclusions. ( B.4,B.5 ) Contacts with a normal force less than the mean force (weak contacts)are known to operate differently than those with a greater-than-mean normal force(strong contacts). The orientation distribution of weak contacts is known to be nearlyisotropic, with only a small bias in the direction of loading. Consistent with thesimulations, the model predicts a Satake fabric ratio for weak contacts of only slightlygreater than one (Table 3, row 9); whereas the corresponding value for the strongcontacts is 1.8 (row 10): that is, the contact anisotropy of all contacts ( = q is derived from the weak contacts (row 7). The same weakcontacts have a more significant role in bearing the mean stress, as more than 26% of p o is attributed to these contacts (row 8). The model is consistent with the simulationsin all of these trends. ( B.6 ) The DEM simulations show that the directions of the contact movementsdo not always conform to those expected of an affine deformation field. As shownin Fig. 2, two particles oriented in the first quadrant would move in a conforming,“forward” direction if particle “ j ” moves upward and to the left over particle “ i ”.Although most contact slip does occur in this forward direction, the simulationsshow that about 40% of sliding contacts slip in the reverse direction (rows 12–13). Themodel predicts this same trend of a large fraction of contacts slipping in the reversedirection. ( B.7 ) The likelihood that a contact is sliding depends upon its orientation and also26 odelDEM dataAvg. η , model Avg. η , DEM x x Contact orientation, θ F r a c t i o n s li d i n g c o n t a c t s ◦ ◦ ◦ ◦ ModelDEM data Avg. η , DEM Normal contact force, g n , normed F r a c t i o n s li d i n g c o n t a c t s Figure 8: Fraction of sliding contacts as a function of (a) contact orientation, (cid:104) | g t ∈ {− µg n , µg n }(cid:105) ( θ c ) ,and (b) normal contact force, (cid:104) | g t ∈ {− µg n , µg n }(cid:105) ( g n ) . upon the magnitude of its normal force. Figure 8a shows that contact slip is mostlikely among contacts with normals n k that are oriented in the direction of extension(direction x for the biaxial loading conditions of this study, Fig. 2). Although themodel’s results in 8a are shifted because of a different η value, the same trend ispredicted. ( B.8 ) Figure 8b shows the relationship between the normal contact force g n andthe fraction of sliding contacts. With both model and simulations, the prevalence ofcontact sliding decreases with an increasing contact force: contact slip is far morelikely for lightly loaded contacts than for those bearing a large normal force. Thistrend is also captured by measuring the average normal forces among the slidingand non-sliding contacts (Table 3, rows 15–17): the average normal force amongsliding contacts is less than a third of that among non-sliding contacts. Moreover,the normal force among forward sliding contacts is slightly greater than the forceamong contacts sliding in the reverse direction. All of these trends are predicted withthe model. ( B.7,B.9 ) The sliding (slip) rate ˙ φ slip also depends upon contact orientation (Fig. 9).As with the DEM data, the model predicts a more vigorous magnitude of the slip rateamong contacts that are oriented in the direction of extension (Fig. 9a). Althoughthe model is consistent with the DEM data in this trend, the magnitudes of the sliprate in the model are considerably greater than those of the simulations: the modelgreatly over-predicts the vigor of contact activity. Figure 9b shows the mean sliprate as a function of contact orientation. Consistent with the DEM data, the modelpredicts that the mean slip rate is, on average, in the “forward” (positive) directionand that the mean is largest at orientations oblique to the directions of compressionand extension (at an angle of about 50 ◦ ). ( B.8 ) Figure 10 shows the relationship between the contact sliding rate and thecontact normal force. Although quantitative agreement is poor, the simulations and27 odel
DEM data x x M e a n m a g n i t ud e s li p r a t e , h | ˙ φ s li p | i / ˙ ε Contact orientation, θ ◦ ◦ ◦ ◦
210 (a)
ModelDEM data x x Contact orientation, θ M e a n s li p r a t e , h ˙ φ s li p i / ˙ ε ◦ ◦ ◦ ◦ Figure 9: Sliding rate ˙ φ slip and contact orientation: (a) average magnitude of sliding rate, (cid:104)| ˙ φ slip | / ˙ ε (cid:105) ( θ c ) ,and (b) mean sliding rate, (cid:104) ˙ φ slip / ˙ ε (cid:105) ( θ c ) . ModelDEM data
Normal contact force, g n , normed M e a n m a g n i t ud e s li p r a t e , h | ˙ φ s li p | i / ˙ ε Figure 10: Sliding rate ˙ φ slip and contact normal force, (cid:104)| ˙ φ slip | / ˙ ε (cid:105) ( g n ) . the model show a reduction in the magnitude of slip movements with increasingnormal force. ( B.10 ) Past DEM analyses of strains within particle clusters have shown that dila-tion is associated with clusters that are elongated in the direction of compression;whereas, clusters that are elongated in the direction of extension tend to contract(Nguyen et al., 2009). With the model, this aspect of granular deformation was stud-ied by extracting the dilation rate predicted when integration is restricted to differentcontact orientations θ c . The deformation function Γ L of Eq. (36) was analyzed withthe methods of Appendix C (see Eq. C.4) to compute the average dilation rate at-tributed to contacts of any given orientation θ c : the average (cid:104) tr ( Γ L ) / ˙ ε (cid:105) ( θ c ) . Althoughthe relationship between contact orientation and dilation was not computed forthe current DEM simulations, the model results, shown in Fig. 11, do confirm pastobservations: contacts oriented in the direction of compression (direction x ) are28 odel x x Contact orientation, θ V o l u m e t r i c r a t e C o n t r a c t i o n D il a t i o n ◦ ◦ ◦ ◦ Figure 11: Average volumetric deformation rate attributed to contacts with orientation θ c : the expecta-tion (cid:104) tr ( Γ L ) / ˙ ε (cid:105) ( θ c ) . ModelDEM data
Friction coefficient, µ S t r e n g t h , q / p o Figure 12: Effect of inter-particle friction coefficient µ on the deviatoric stress during critical state flow. associated with the dilation of an assembly; whereas those oriented in the directionof extension tend to produce contraction. ( B.12 ) As the final aspect of behavior, we consider the effect of the inter-particlefriction coefficient µ on bulk strength at the critical state. DEM simulations wereconducted with five coefficients and the model was solved with these same values.Figure 12 compares the model’s predictions of strength with the results of DEMsimulations, expressed as the deviator stress ratio q/p o . Although the model over-predicts strength at the critical state, its results follow a similar trend of the DEM data:strength increases with an increasing friction coefficient, but the increase is rathersmall for coefficients greater than 0.30.
4. Discussion
In the Introduction, we listed several micro-scale characteristics of critical stategranular flow that have been observed in past experiments and simulations. Theseobservations as well as new DEM simulations were used to evaluate a proposedmodel for granular flow. Nearly all of the model’s predictions are in qualitative (if not29uantitative) agreement with observed characteristics. The model favorably predictsanisotropies of the contact orientations, contact forces, and contact movementsand of the orientations of those contacts undergoing slip. The model also favorablypredicts relationships between the contact force magnitude, contact motion, andcontact slip direction. Although several pages were required for its derivation, and itssolution requires quite demanding computations, the model’s concept is fairly simple:the contact landscape of motion, force, and orientation is maximally disordered atthe critical state, after one accounts for the biases that arise from certain fundamentaland empirical aspects of flow.The model adopts individual contacts as its generic units, but without identifyingtheir locations or their affiliations with other contacts. As such, it is certainly unableto resolve the spatial localization and meso-scale patterning that is usually foundwithin granular flows (item D.1 in the Introduction). Past studies have shown thatmotions and forces are spatially coordinated (correlated) over distances of severalparticle diameters (Kuhn, 2003). Such patterning and correlation affirm a greaterorder (and lower entropy) than is accessible with the model.Biased by the loading direction in the second constraint, the model predicts amaximally disordered condition that is anisotropic. Although anisotropy, by itself,connotes greater order than isotropy, the small increase in the order of contact orien-tation is more than offset by an increase in the disorder associated with the contactforces and movements.In some respects, the model either poorly predicts behavior or requires an ad-ditional constraint to force agreement. These shortcomings include the density ofthe tangential contact forces among non-sliding contacts (Fig. 7b) and the fractionof sliding contacts (constraint 5). In other respects, the model is consistent withtrends in the simulation results but is in poor quantitative agreement (for example,in the vigor of the contact movements, Fig. 9a). Although its successes are promising,the model’s shortcoming are more revealing, raising this question: why are somebehaviors not amenable to a simple, contact-centric statistical treatment?An underlying assumption of the Jaynes maximum entropy (MaxEnt) approach isthat each micro-state is equally probable, provided that the micro-state is consistentwith available information (i.e., the imposed constraints). In the current setting,each of the “ k ” micro-states {· · · } in the dimension- M phase space belongs to aparticular macro-state ( · · · ) and is assumed equiprobable with the other micro-statesthat populate the same macro-state. In the model, the micro-states only incorporatecontact information and are devoid of content on the particles’ locations or theirtopologic arrangement. The equiprobable assumption ignores the fundamentalmeso-scale restrictions that are shared by the particles and contacts: the contactforces on each particle must be in equilibrium, and the contact movements mustbe compatible with a corresponding set of particle movements and rotations. Theseadditional restrictions operate at a meso-scale much larger than individual contacts,30nd will bias some micro-states and disallow others altogether, creating an unevenlandscape of the dimension- M phase space of micro-states. For example, withoutthe rather coersive restriction on the fraction of sliding contacts (constraint 5), themodel predicts a fraction greater than 80%, rather than the 11.2% measured in thesimulations. This result suggests that particle motions are coordinated (ordered)in a manner that greatly reduces the population of sliding contacts. Reasoning of afundamental nature, applied to meso-scale conditions, is palpably needed to replacebrute force, empirical constraints, such as our constraint on the fraction of slidingcontacts and other constraints that are approximate in nature (constraints 3 and 4).We consider the lack of meso-scale information the primary deficiency of the model,and a means of introducing such information is suggested in the following paragraphs.On the other hand, those micro-scale characteristics and trends that are favorablypredicted are apparently insensitive to meso-scale equilibrium and compatibilitylimitations, but are instead determined (or biased) by the rather broad constraints ofpressure, of dissipation consistency, and of the isochoric condition.The paper makes exclusive use of information that is formulated as momentconstraints on the density p ( · · · ) (Eqs. 7 and 8). An alternative form of informationwas introduced by the author in a paper on topologic entropy (Kuhn, 2014). In thiswork, a relative entropy principle (cross-entropy or Kullback-Leibler entropy) wasapplied, admitting information in the form of certain a priori inclinations (priors) of p ( · · · ) . Because these inclinations can have either an empirical or theoretical origin,this approach allows information — even imperfect information — gained fromone theory (for example, an affine field approach) to be included in a shear-drivenentropy model.We should also consider the entropy model in relation to mean-field upscalingmodels, which have been used with reasonable success in estimating the small-strain stiffness of dense granular materials. These models can be more descriptivethan predictive, as they require extensive additional information, in the form ofcontact stiffnesses, packing data, etc. Mean-field models, however, do provide aclearer, deterministic connection between the micro-scale and bulk behaviors: themicro and macro behaviors are unambiguously linked through contact rules. Thecorresponding connection is rather opaque in the paper’s approach. For example,why do the constraints lead to a greater number of forward sliding contacts thanreverse sliding contacts? What features of the model lead to smaller normal forcesamong sliding contacts than for non-sliding contacts? Although the model yieldssome convincing predictions, the model is certainly enigmatic in regard to its results,as the model’s predictions are extracted only through a painstaking integration ofdistribution averages.Finally, we suggest future extensions of this work. With some difficulty, the modelcould certainly be extended to three dimensions, which would require that the tan-gential forces and motions be expressed with Euler angles rather than the single θ θ c and θ (cid:96) (i.e., by relaxing the Dirac restriction).In the Introduction, we had noted that granular materials consistently convergetoward the critical state when loaded from differing initial conditions. Although thecurrent model addresses the terminal condition at the critical state, the disorder inthe paper’s contact quantitites (and other quantities, as well) could be tracked duringthe evolving condition of granular loading. An ideal entropy measure would increasein a consistent, monotonic manner during the loading process. This approach mightlead to an improved entropy definition based upon a more appropriate phase space. Acknowledgement
The author greatfully acknowledges productive discussions of this work with Dr.Niels Kruyt.
Dedication
The paper is dedicated to the memory of Dr. Masao Satake (1927–2013), whomade significant fundamental contributions to granular mechanics.
Appendix A. Integrations
Each constraint involves integrating a product Γ i ( · · · ) p ( · · · ) , as sketched inEq. (7). Any evaluation of these integrals must, in practice, reconcile the non-smoothnature of these integrands: both the rigid-frictional constitutive constraint of Eq. (4)and the kernel K () are discontinuous. In the dimension-6 phase space, large non-zero values of ˙ φ slip can exist alongside an enforced zero value at the juncture of slidingand non-sliding behaviors: at hyper-planes g t = ± µg n . This difficulty is resolved bysplitting the original integral into three parts, so that each part addresses a singlebranch of the three constitutive cases in Eq. (4):32 · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) = (cid:90) R + g n (cid:90) [0 , π ) θ c (cid:90) [0 , π ) θ (cid:96) (cid:90) R − ˙ φ slip (cid:90) R ˙ φ rigid with g t = − µg n and ˙ φ slip < (A.1) + (cid:90) R + g n (cid:90) ( − , g t (cid:90) [0 , π ) θ c (cid:90) [0 , π ) θ (cid:96) (cid:90) R ˙ φ rigid with g t ∈ ( − µg n , µg n ) and ˙ φ slip = 0+ (cid:90) R + g n (cid:90) [0 , π ) θ c (cid:90) [0 , π ) θ (cid:96) (cid:90) R + ˙ φ slip (cid:90) R ˙ φ rigid with g t = µg n and ˙ φ slip > The second integral on the right addresses the non-sliding branch, and the firstand third integrals address sliding in the “reverse” and “forward” directions (slidingdirections that are contrary to and consistent with those of affine deformation, asdescribed in Section 2.1). Note that the three integrals are 5-fold, compared withthe original 6-fold integrals (on the left of Eq. A.1), greatly reducing the complexityof numerical evaluations. Separating (A.1) as three integrals also permits a directapplication of constraint 5 (Eq. 45) and aids in computing separate statistics for thenon-sliding and sliding (both forward and reverse) contacts (Section 3.3).
Appendix B. Entropy numerics
The entropy H in Eqs. (46)–(50) is maximized by finding the five multipliers λ i that satisfy the constraints (cid:104) Γ i ( · · · ) (cid:105) = Γ i of Eq. 8. That is, the proper λ i are the rootsof five equations f i ( λ ) = A i ( λ ) − Γ i Z ( λ ) = 0 (B.1)where argument λ represents the list { λ , λ , . . . , λ } , and A i ( λ ) = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) Γ i ( · · · ) exp − (cid:88) j =1 λ j Γ j ( · · · ) (B.2)We used the minpack library to solve these equations by minimizing the sum of theirsquared residuals (specifically, function hybrj1 , which applies the Powell hybridmethod) (Moré et al., 1984). The solver requires evaluation of the partition function Z ( λ ) , the five functions f i ( λ ) , and the × Jacobian33 f i ∂λ j = − B ij + Γ i A j (B.3) B ij = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) Γ i ( · · · )Γ j ( · · · ) exp (cid:32) − (cid:88) k =1 λ k Γ k ( · · · ) (cid:33) (B.4)The total of integrands were evaluated within each iteration ofthe solver hybrj1 . Each symbolic integral is itself the sum of three integrals, each5-fold (see Eq. A.1). All integrals were evaluated with the CUBA library (specifically thefunction llcuhre , which applies adaptive polynomial cubature rules and permits avector of integrands to be evaluated with each call). Except for the variables g t , θ c ,and θ (cid:96) , which have ranges of [ − µg n , µg n ] or [0 , π ) , all integrations are improper withinfinite range. A change of variable, for example u = tan − g n , was applied in suchcases, and integrations were evaluated with the following support: g n ∈ [0 , , ˙ φ slip ∈ [ − , , and ˙ φ rigid ∈ [ − , . Two hundred million points were queried inevaluating each integrand within each of the three integrals within each iteration.The mean stress p o and strain rate ˙ ε were both set to 1.0, and the factor M int (cid:96) /A was1.50.After solving the multipliers λ i , meaningful results, such as those in Section 3.3must be extracted by evaluating integrals with appropriate integrands to computethe relevant (cid:104)·(cid:105) quantities. These calculations are described in the next appendix. Appendix C. Statistical calculations and notation
After determining the multipliers λ i , meaningful information can be extractedfrom the probability density p ( · · · ) by post-processing of various integrals in theform of Eq. (A.1). Statistical notations used in the paper are intended to suggest theirnumerical evaluation, although some differ from traditional notations. The expectedvalue of a function h ( · · · ) is written as (cid:104) h ( · · · ) (cid:105) or (cid:104) h (cid:105) and is computed, as in Eqs. (7)and (A.1), by integrating across the full domain ( g n , g t , θ c , θ (cid:96) , ˙ φ slip , ˙ φ rigid ) . When thedomain is restricted by a condition Ψ , we use the notation (cid:104) h ( · · · ) (cid:105)| Ψ , (cid:104) h ( · · · ) (cid:105) | Ψ = (cid:90) · · · (cid:90) ( g n ,g t ,θ c ,θ (cid:96) , ˙ φ slip , ˙ φ rigid ) h ( · · · ) p ( · · · ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Ψ (C.1)in which the integrand is zero outside of this restricted domain (i.e., where condition Ψ is not met). For example, the probability of condition Ψ is simply (cid:104) (cid:105)| Ψ , with h = 1 .34he conditional expectation of function h ( · · · ) subject to Ψ is written as (cid:104) h ( · · · ) | Ψ (cid:105) and is computed as (cid:104) h ( · · · ) | Ψ (cid:105) = (cid:104) h ( · · · ) (cid:105)| Ψ (cid:104) (cid:105)| Ψ (C.2)We also compare the statistics of contacts that are not sliding or are sliding in the“reverse” or “forward” directions. The three cases correspond to different Ψ conditions: ˙ φ slip = 0 , ˙ φ slip = < , and ˙ φ slip = > (or alternatively, g t ∈ ( − µg n , µg n ) , g t = − µg n ,and g t = µg n ), as explicitly isolated with the three integrals in Eq. (A.1).The marginal probability density of a single contact quantity will be written asa derivative, in keeping with the manner in which it is computed. For example, themarginal probability density of the contact orientation θ c is written as the derivative D θ c (cid:104) (cid:105) : D θ c (cid:104) (cid:105) = lim ∆ θ → (cid:104) (cid:105)| θ c <θ<θ c +∆ θ ∆ θ (C.3)As another example, the average expected normal force g n as a function of contactorientation θ c is written and computed as (cid:104) g n (cid:105) ( θ c ) = D θ c (cid:104) g n (cid:105)D θ c (cid:104) (cid:105) = lim ∆ θ → (cid:104) g n (cid:105)| θ c <θ<θ c +∆ θ (cid:104) (cid:105)| θ c <θ<θ c +∆ θ (C.4) Appendix D. DEM simulations, continued
Multiple discrete element (DEM) simulations of biaxial compression were con-ducted on bi-disperse assemblies of 676 particles. The two disk varieties have ratiosof 1.5:1 in size, 1:2.25 in number, and 1:1 in cumulative area (that is, 468 particlesof size . and 208 of size . ). The assemblies were small enough to prevent grossnon-homogeneity in the form shear bands (no bands were observed, by using themethods of Kuhn, 1999), yet large enough to capture the average, bulk material be-havior. To develop more robust statistics, 168 different assemblies were created bycompacting random sparse frictionless mixtures of the two disk sizes into denseisotropic packings within periodic boundaries until the average contact indentationwas . times the average radius. In the subsequent biaxial loading sequences,linear contact stiffnesses were applied between particles with equal tangential andnormal coefficients ( k t = k n ), and the base friction coefficient µ = 0 . was en-forced during the pair-wise particle interactions (da Cruz et al., 2005). Coefficients µ = 0 . , 0.3, 0.7, and 0.9 were used in separate simulations, all starting from the same168 assemblies. Using the standard DEM algorithm, the initially square assemblieswere horizontally compressed in increments ∆ (cid:15) = 1 × − while maintaining aconstant mean stress of × − k n (Cundall and Strack, 1979). The strain rate wassufficiently slow to maintain the quasi-static condition, with the inertial number I equal to . × − (i.e., the ratio of the shear time / ˙ ε to the inertial time of a35article (cid:112) m/p o , see da Cruz et al., 2005). These loading conditions coincide withthose described in Sections 2.5–2.6 and result in the stress, fabric, and volumetricbehavior that are shown in Fig. 4. A large initial stiffness causes the deviatoric stressto rise quickly from zero to a peak stress at strain − ε = 2% , and the critical statecondition is attained at compressive strains − ε of 16–18%. During subsequentsteady-state deformation, the contact conditions were interrogated at five strainsbetween − ε = 16% and (Fig. 4), and rates ˙ n and ˙ φ were determined withpairs of interrogations separated by strain ∆ ε = 1 × − . Applying the ergodicityprinciple, micro-state statistics were averaged across the five strains and the 168assemblies, involving 840 pairs of interrogations containing about 820,000 contacts.Another series of simulations were run with poly-disperse assemblies having particlesizes spanning a range of 3.0. The results were similar to those of the bi-disperseassemblies, and only the latter are reported herein. ReferencesReferences
Bagi, K., 2006. Analysis of microstructural strain tensors for granular assemblies. Int. J. Solids Struct.43 (10), 3166–3184.Blumenfeld, R., Edwards, S. F., 2009. On granular stress statistics: Compactivity, angoricity, and someopen issues. The Journal of Physical Chemistry B 113 (12), 3981–3987.Brown, C. B., 2000. Entropy and granular materials: model. J. Eng. Mech. 126 (6), 599–604.Brown, C. B., Elms, D. G., Hanson, M. T., Nikzad, K., Worden, R. E., 2000. Entropy and granular materials:experiments. J. Eng. Mech. 126 (6), 605–610.Chakraborty, B., 2010. Statistical ensemble approach to stress transmission in granular packings. SoftMatt. 6 (13), 2884–2893.Cundall, P. A., Strack, O. D. L., 1979. A discrete numerical model for granular assemblies. Géotechnique29 (1), 47–65.da Cruz, F., Emam, S., Prochnow, M., Roux, J.-N., Chevoir, F., 2005. Rheophysics of dense granularmaterials: Discrete simulation of plane shear flows. Phys. Rev. E 72 (2), 021309.Desrues, J., Viggiani, G., 2004. Strain localization in sand: an overview of the experimental resultsobtained in Grenoble using stereophotogrammetry. Int. J. Numer. and Anal. Methods in Geomech.28, 279–321.Goddard, J. D., 2004. On entropy estimates of contact forces in static granular assemblies. Int. J. SolidsStruct. 41 (21), 5851–5861.Härtl, J., Ooi, J. Y., 2008. Experiments and simulations of direct shear tests: porosity, contact frictionand bulk friction. Granul. Matter 10 (4), 263–271.Jaynes, E. T., 1957. Information theory and statistical mechanics. Phys. Rev. 106 (4), 620–630.Jenkins, J. T., Strack, O. D. L., 1993. Mean-field inelastic behavior of random arrays of identical spheres.Mech. of Mater. 16 (1), 25–33.Kruyt, N. P., Antony, S. J., 2007. Force, relative-displacement, and work networks in granular materialssubjected to quasistatic deformation. Phys. Rev. E 75 (5), 051308.Kruyt, N. P., Rothenburg, L., 2014. On micromechanical characteristics of the critical state of two-dimensional granular materials. Acta Mechanica 225 (8), 2301–2318.Kuhn, M. R., 1999. Structured deformation in granular materials. Mech. of Mater. 31 (6), 407–429.Kuhn, M. R., 2003. Heterogeneity and patterning in the quasi-static behavior of granular materials.Granul. Matter 4 (4), 155–166. uhn, M. R., 2004a. Boundary integral for gradient averaging in two dimensions: application to polygonalregions in granular materials. Int. J. Numer. Methods Eng. 59 (4), 559–576.Kuhn, M. R., 2004b. Rates of stress in dense unbonded frictional materials during slow loading. In:Antony, S. J., Hoyle, W., Ding, Y. (Eds.), Advances in Granular Materials: Fundamentals and Applica-tions. Royal Society of Chemistry, London, U.K., pp. 1–28.Kuhn, M. R., 2014. Dense granular flow at the critical state: maximum entropy and topological disorder.Granul. Matter 16 (4), 499–508.Kuhn, M. R., Bagi, K., 2004. Contact rolling and deformation in granular media. Int. J. Solids Struct.41 (21), 5793–5820.Kuhn, M. R., Bagi, K., 2005. On the relative motions of two rigid bodies at a compliant contact: applicationto granular media. Mech. Res. Comm. 32 (4), 463–480.Liao, C.-L., Chang, T.-P., Young, D.-H., 1997. Stress-strain relationship for granular materials based onthe hypothesis of best fit. Int. J. Solids Struct. 34 (31–32), 4087–4100.Majmudar, T. S., Bhehringer, R. P., 2005. Contact force measurements and stress-induced anisotropy ingranular materials. Nature 435 (1079), 1079–1082.Majmudar, T. S., Sperl, M., Luding, S., Behringer, R. P., 2007. Jamming transition in granular systems.Phys. Rev. Lett. 98 (5), 058001.Moré, J. J., Sorensen, D. C., Hillstrom, K. E., Garbow, B. S., 1984. The MINPACK project. In: Cowell, W. R.(Ed.), Sources and Development of Mathematical Software. Prentice-Hall, Englewood Cliffs, N.J., pp.88–111.Nguyen, N.-S., Magoariec, H., Cambou, B., Danescu, A., 2009. Analysis of structure and strain at themeso-scale in 2d granular materials. Int. J. Solids Struct. 46 (17), 3257–3271.Peña, A. A., Lizcano, A., Alonso-Marroquin, F., Herrmann, H. J., 2008. Biaxial test simulations using apacking of polygonal particles. Int. J. Numer. and Anal. Methods in Geomech. 32 (2), 143–160.Radjai, F., Wolf, D. E., Jean, M., Moreau, J.-J., 1998. Bimodal character of stress transmission in granularpackings. Phys. Rev. Lett. 80 (1), 61–64.Rothenburg, L., 1980. Micromechanics of idealized granular systems. Ph.D. thesis, Carleton University,Ottawa, Ontario, Canada.Rothenburg, L., Bathurst, R., 1989. Analytical study of induced anisotropy in idealized granular materials.Géotechnique 39 (4), 601–614.Rothenburg, L., Kruyt, N. P., 2009. Micromechanical definition of an entropy for quasi-static deformationof granular materials. J. Mech. Phys. Solids 57 (3), 634–655.Satake, M., 1982. Fabric tensor in granular materials. In: Vermeer, P. A., Luger, H. J. (Eds.), Proc. IUTAMSymp. on Deformation and Failure of Granular Materials. A.A. Balkema, Rotterdam, pp. 63–68.Schofield, A. N., Wroth, P., 1968. Critical state soil mechanics. McGraw-Hill, New York.Shannon, C. E., 1948. A mathematical theory of communication. Bell System Technical Journal 27 (3),379–423.Thornton, C., 2000. Numerical simulations of deviatoric shear deformation of granular media. Géotech-nique 50 (1), 43–53.Tordesillas, A., Lin, Q., Zhang, J., Behringer, R. P., Shi, J., 2011. Structural stability and jamming of self-organized cluster conformations in dense granular materials. J. Mech. Phys. Solids 59 (2), 265–296.Troadec, H., Radjai, F., Roux, S., Charmet, J. C., 2002. Model for granular texture with steric exclusion.Phys. Rev. E 66, 041305.Williams, J. R., Nabha, R., 1997. Coherent vortex structures in deforming granular materials. Mech.Cohesive-Frict. Mater. 2 (3), 223–236.Yoon, S. W., Giménez, D., 2012. Entropy characterization of soil pore systems derived from soil-waterretention curves. Soil Sci. 177 (6), 361–368.Zhao, J., Guo, N., 2013. A new definition on critical state of granular media accounting for fabric aniso-tropy. In: AIP Conference Proceedings. Vol. 1542. pp. 229–232.uhn, M. R., 2004a. Boundary integral for gradient averaging in two dimensions: application to polygonalregions in granular materials. Int. J. Numer. Methods Eng. 59 (4), 559–576.Kuhn, M. R., 2004b. Rates of stress in dense unbonded frictional materials during slow loading. In:Antony, S. J., Hoyle, W., Ding, Y. (Eds.), Advances in Granular Materials: Fundamentals and Applica-tions. Royal Society of Chemistry, London, U.K., pp. 1–28.Kuhn, M. R., 2014. Dense granular flow at the critical state: maximum entropy and topological disorder.Granul. Matter 16 (4), 499–508.Kuhn, M. R., Bagi, K., 2004. Contact rolling and deformation in granular media. Int. J. Solids Struct.41 (21), 5793–5820.Kuhn, M. R., Bagi, K., 2005. On the relative motions of two rigid bodies at a compliant contact: applicationto granular media. Mech. Res. Comm. 32 (4), 463–480.Liao, C.-L., Chang, T.-P., Young, D.-H., 1997. Stress-strain relationship for granular materials based onthe hypothesis of best fit. Int. J. Solids Struct. 34 (31–32), 4087–4100.Majmudar, T. S., Bhehringer, R. P., 2005. Contact force measurements and stress-induced anisotropy ingranular materials. Nature 435 (1079), 1079–1082.Majmudar, T. S., Sperl, M., Luding, S., Behringer, R. P., 2007. Jamming transition in granular systems.Phys. Rev. Lett. 98 (5), 058001.Moré, J. J., Sorensen, D. C., Hillstrom, K. E., Garbow, B. S., 1984. The MINPACK project. In: Cowell, W. R.(Ed.), Sources and Development of Mathematical Software. Prentice-Hall, Englewood Cliffs, N.J., pp.88–111.Nguyen, N.-S., Magoariec, H., Cambou, B., Danescu, A., 2009. Analysis of structure and strain at themeso-scale in 2d granular materials. Int. J. Solids Struct. 46 (17), 3257–3271.Peña, A. A., Lizcano, A., Alonso-Marroquin, F., Herrmann, H. J., 2008. Biaxial test simulations using apacking of polygonal particles. Int. J. Numer. and Anal. Methods in Geomech. 32 (2), 143–160.Radjai, F., Wolf, D. E., Jean, M., Moreau, J.-J., 1998. Bimodal character of stress transmission in granularpackings. Phys. Rev. Lett. 80 (1), 61–64.Rothenburg, L., 1980. Micromechanics of idealized granular systems. Ph.D. thesis, Carleton University,Ottawa, Ontario, Canada.Rothenburg, L., Bathurst, R., 1989. Analytical study of induced anisotropy in idealized granular materials.Géotechnique 39 (4), 601–614.Rothenburg, L., Kruyt, N. P., 2009. Micromechanical definition of an entropy for quasi-static deformationof granular materials. J. Mech. Phys. Solids 57 (3), 634–655.Satake, M., 1982. Fabric tensor in granular materials. In: Vermeer, P. A., Luger, H. J. (Eds.), Proc. IUTAMSymp. on Deformation and Failure of Granular Materials. A.A. Balkema, Rotterdam, pp. 63–68.Schofield, A. N., Wroth, P., 1968. Critical state soil mechanics. McGraw-Hill, New York.Shannon, C. E., 1948. A mathematical theory of communication. Bell System Technical Journal 27 (3),379–423.Thornton, C., 2000. Numerical simulations of deviatoric shear deformation of granular media. Géotech-nique 50 (1), 43–53.Tordesillas, A., Lin, Q., Zhang, J., Behringer, R. P., Shi, J., 2011. Structural stability and jamming of self-organized cluster conformations in dense granular materials. J. Mech. Phys. Solids 59 (2), 265–296.Troadec, H., Radjai, F., Roux, S., Charmet, J. C., 2002. Model for granular texture with steric exclusion.Phys. Rev. E 66, 041305.Williams, J. R., Nabha, R., 1997. Coherent vortex structures in deforming granular materials. Mech.Cohesive-Frict. Mater. 2 (3), 223–236.Yoon, S. W., Giménez, D., 2012. Entropy characterization of soil pore systems derived from soil-waterretention curves. Soil Sci. 177 (6), 361–368.Zhao, J., Guo, N., 2013. A new definition on critical state of granular media accounting for fabric aniso-tropy. In: AIP Conference Proceedings. Vol. 1542. pp. 229–232.