UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Continuum limits of granular systems DeGiuli, Eric 2012

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2013_fall_degiuli_eric.pdf [ 3.1MB ]
JSON: 24-1.0071936.json
JSON-LD: 24-1.0071936-ld.json
RDF/XML (Pretty): 24-1.0071936-rdf.xml
RDF/JSON: 24-1.0071936-rdf.json
Turtle: 24-1.0071936-turtle.txt
N-Triples: 24-1.0071936-rdf-ntriples.txt
Original Record: 24-1.0071936-source.json
Full Text

Full Text

Continuum Limits of Granular Systems by Eric DeGiuli  B.Sc. (Hons.), University of Toronto, 2006 M.Sc., University of British Columbia, 2009  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Mathematics)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2013 c Eric DeGiuli 2012  Abstract Despite a century of study, the macroscopic behaviour of quasistatic granular materials remains poorly understood. In particular, we lack a fundamental system of continuum equations, comparable to the Navier-Stokes equations for a Newtonian fluid. In this thesis, we derive continuum models for two-dimensional granular materials directly from the grain scale, using tools of discrete calculus, which we develop. To make this objective precise, we pose the canonical isostatic problem: a marginally stable granular material in the plane has 4 components of the stress tensor σ ˆ , but only 3 continuum equations in Newton’s laws ∇ · σ ˆ = 0 and σ ˆ=σ ˆT . At isostaticity, there is a missing stress-geometry equation, arising from Newton’s laws at the grain scale, which is not present in their conventional continuum form. We first show that a discrete potential ψ can be defined such that the stress tensor is written as σ ˆ = ∇ × ∇ × ψ, where the derivatives are given an exact meaning at the grain scale, and converge to their continuum counterpart in an appropriate limit. The introduction of ψ allows us to understand how force and torque balance couple neighbouring grains, and thus to understand where the stress-geometry equation is hidden. Using this formulation, we derive the missing stress-geometry equation ∆(Fˆ : ∇∇ψ) = 0, introducing a fabric tensor Fˆ which characterizes the geometry. We show that the equation imposes granularity in a literal sense, and that on a homogeneous fabric, the equation reduces to a particular form of anisotropic elasticity. We then discuss the deformation of rigid granular materials, and derive the mean-field phase diagram for quasistatic flow. We find that isostatic states are fluid states, existing between solid and gaseous phases. The appearance of isostaticity is linked to the saturation of steric exclusion and Coulomb inequalities. Finally, we present a model for the fluctuations of contact forces using tools of statistical mechanics. We find that force chains, the filamentary networks of contact forces ubiquitously observed in experiments, arise from an entropic instability which favours localization of contact forces.  ii  Preface A version of Chapter 2 was published in Physical Review E, Oct. 2011 [DM11]. All research was conducted by the author, except the following: Jim McElwaine helped find the gauge group of ψ in three dimensions, and helped organize the paper. Throughout this thesis, use has been made of Jim McElwaine’s Discrete Element Method (DEM) software for simulating granular materials. Figures 2.1, A.1, 3.3, 3.4, 3.5, 3.6, 3.7, and 6.2 have been produced using output from this program.  iii  Table of Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  List of Tables  vii  Preface  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Symbols - Chapters 1-4  . . . . . . . . . . . . . . . . . . . . . .  List of Symbols - Chapters 5-6  . . . . . . . . . . . . . . . . . . . . . . xiii  Acknowledgments  x  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xv  1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Outline of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1 8 9  2 Discrete Calculus . . . . . . . . 2.1 Intrinsic formulation . . . . . 2.2 Extrinsic formulation . . . . 2.3 Discussion and extensions . . 2.3.1 Physical interpretation 2.3.2 Polydisperse packings 2.3.3 Force law . . . . . . . 2.4 Conclusion . . . . . . . . . .  . . . .  . . . . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  11 14 17 20 20 23 23 24  3 Isostaticity . . . . . . . . . . . 3.1 Geometric formulation . . 3.2 Discrete calculus equations 3.3 Final equation . . . . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  . . . .  26 28 36 40  . . . .  . . . .  iv  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  42 43 43 45 47 49  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  51 51 53 56 57 60 62 65 69 70 71 72 72  5 Hyperstaticity . . . . . . . . . . . . . . . . . . . . . . 5.1 Phase space . . . . . . . . . . . . . . . . . . . . . 5.2 Measure . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Edwards’ approach . . . . . . . . . . . . . 5.2.2 Force Network Ensemble . . . . . . . . . . 5.3 Properties of the enlarged Force Network Ensemble 5.3.1 Contact potential . . . . . . . . . . . . . . 5.3.2 Variational principles . . . . . . . . . . . . 5.3.3 Microcanonical v.s. Canonical Ensembles . 5.3.4 Scaling . . . . . . . . . . . . . . . . . . . . 5.4 A note on percolation . . . . . . . . . . . . . . . .  . . . . . . . . . .  . . . . . . . . . . .  . . . . . . . . . . .  . . . . . . . . . . .  . . . . . . . . . . .  . . . . . . . . . . .  . . . . . . . . . . .  . . . . . . . . . . .  . . . . . . . . . . .  73 74 76 76 78 79 79 80 81 82 83  6 Mean-field Theory of Hyperstatic Packings 6.1 Gauge fields . . . . . . . . . . . . . . . . . 6.2 Isostaticity . . . . . . . . . . . . . . . . . . 6.3 Isostatic cluster partition function . . . . . 6.3.1 Mean-field approximations . . . . . 6.3.2 Computing the partition function . 6.4 Global partition function . . . . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  85 87 89 92 93 95 97  3.4  3.3.1 Scaling . . . . . . . . . . . . . . . . . 3.3.2 Solutions . . . . . . . . . . . . . . . . 3.3.3 Homogeneous and isotropic fabric . . 3.3.4 Homogeneous and anisotropic fabric . 3.3.5 Inhomogeneous and anisotropic fabric Discussion . . . . . . . . . . . . . . . . . . .  4 Deformation . . . . . . . . . . . 4.1 Kinematic variables . . . . . 4.2 Constitutive laws . . . . . . 4.3 Geometric compatibility . . . 4.3.1 Floppy modes . . . . 4.3.2 General deformations 4.4 Fabric . . . . . . . . . . . . . 4.5 Summary of Chapters 2-4 . . 4.6 Homogeneous solutions . . . 4.6.1 Simple shear flow . . 4.6.2 Pure shear flow . . . 4.6.3 Isotropic expansion . 4.6.4 Conclusion . . . . . .  . . . . . . . . . . . . .  v  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  6.5  6.6  6.4.1 Mean-field contact forces 6.4.2 Free-energy minimization 6.4.3 Saddle-point equations . Free energy . . . . . . . . . . . . 6.5.1 Equation of state . . . . 6.5.2 Fluctuations . . . . . . . Conclusion . . . . . . . . . . . .  7 Conclusion  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  98 99 100 102 103 105 109  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114  Appendices A Torque Balance in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . 127 B Discrete Calculus in 3D . . . . . . . . . . . . . . . . . . . . . . . . . 128 C Convexity of psi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 D Numerical Simulations  . . . . . . . . . . . . . . . . . . . . . . . . . 133  E Differential Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . 134 F Renormalization of rho  . . . . . . . . . . . . . . . . . . . . . . . . . 136  G Green’s Function for a Half-plane  . . . . . . . . . . . . . . . . . . 139  H Fabric Change due to Contact Opening . . . . . . . . . . . . . . . 141 I  Jacobians I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143  J Jacobians II  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146  K Mean-field Relation Between rho’s . . . . . . . . . . . . . . . . . . 148 L Contact-angle Anisotropy . . . . . . . . . . . . . . . . . . . . . . . . 151  vi  List of Tables 4.1  Deformation-stress duality . . . . . . . . . . . . . . . . . . . . . . .  vii  58  List of Figures 1.1  Sketch of 2D packing of rigid grains, showing the position of a contact, r c , and the position of a grain, r g . The hatched grains are rattlers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4  Delaunay triangulation and Voronoi tesselation . . . . . . . . . . . Maxwell-Cremona Diagram . . . . . . . . . . . . . . . . . . . . . .  16 21  Contact-network tessellation and its dual . . . . . . . . . . . . . . Notation convention for c , tc , and sg vectors. We also have sc = r − r = tc − tc . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Geometric interpretation of ψ. . . . . . . . . . . . . . . . . . . . . 3.4 Averaging cells Ωm . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Simple-shear simulation snapshot. . . . . . . . . . . . . . . . . . . ˆ . . . . . . . . . . . . . . . . . 3.6 Convergence of fabric tensor g ˆ to δ. 3.7 Convergence of fabric tensor Fˆ . . . . . . . . . . . . . . . . . . . . . 3.8 Green’s function to a normal forcing on an isotropic geometry . . . 3.9 Green’s function to a normal forcing on an anisotropic geometry . 3.10 Green’s function to a shear forcing on an isotropic geometry . . . . 3.11 Green’s function to a shear forcing on an anisotropic geometry . .  27  2.1 2.2 3.1 3.2  28 32 33 34 35 41 45 48 48 49  4.1  Mean-field phase diagram for rigid, quasistatic, 2D granular materials. 67  6.1 6.2 6.3 6.4 6.5 6.6  Simple-shear simulation snapshot. . Isostatic cluster decomposition . . . √ Simulation results: P vs. z and ∆z Proportion of rattlers . . . . . . . . . Prediction of 1/α . . . . . . . . . . . Ratio of α predictions . . . . . . . .  . . . . . . . . vs φ. . . . . . . . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  . . . . . .  85 91 107 108 108 109  A.1 Local geometry around a grain . . . . . . . . . . . . . . . . . . . . 127 B.1 Subset of Delaunay triangulation and Voronoi tesselation . . . . . 128  viii  K.1 Cluster homotopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 L.1 Probability distribution of γ c . . . . . . . . . . . . . . . . . . . . . . 152 c . . . . . . . . . . . . . . . . . . . . . 153 L.2 Probability distribution of fN  ix  List of Symbols - Chapters 1-4 Symbol Ac A Cg d D ˆ D Fg Fˆ Fˆ fgc c fN fTc ˆc G g g ˆ H H I Lc c c g  ˆ± M NL NRC NRG NS ns P qgc  Definition sc  c  × Maxwell-Cremona area Set of contacts belonging to grain g Dimension of space Mean linear dimension of grains 1 T 2 (∇w + (∇w) ) Body force exerted on grain g Fabric tensor on loop Fabric tensor Contact force exerted on grain g at contact c Normal component of contact force at c Tangential component of contact force at c Contact normal fabric tensor at c Gravity vector Fabric tensor Generalized hyperstaticity Hyperstaticity Inertial number Set of loops adjacent to contact c Contact vector Contact vector δˆ ± µ1 εˆ Number of loops Number of real (force-bearing) contacts Number of real (force-bearing) grains Number of sliding contacts Number of sliding contacts per grain Pressure Weight  x  Equation Reference (3.18) (2.19) (1.5) (1.7) (1.2) (1.5) (3.32) (3.39) (1.5) (4.8) (4.8) (4.12) (4.55) (3.11) (2.5) (2.5) (1.2) (2.9) Figure 3 Figure 2.1 (4.17) (2.6) (2.5) (2.5) (1.7) (1.7) (2.28) (4.22) Continued on next page  Symbol [z] rc rg r sc stg tcg tc tr u V Vc vc v vg w wgc z¯ z¯L ziso z ˆ Γ ˆc Γ ∆¯ z εˆ κ µ ρ ρg ρt σ ˆ φ ϕ ϕg ϕ ψ  Table 1 – continued from previous page Definition Equation Reference Real part of z (3.48) Position of contact c (1.5) Centre-of-mass of grain g (1.5) Position of loop (2.11),(3.9) Cross-contact vector Figure 3 Vector pointing across triangle t (2.13) rc − rg (4.1) c r −r (3.3) Tensor trace (1.12) Continuum velocity field (4.45) Volume of the packing (1.8) Relative velocity of g (c) with respect to g(c) at c (4.2) Weight (3.5) Continuum grain velocity (4.24) Grain linear velocity (4.1) Continuum contact velocity (4.24) Contact velocity of g at c (4.1) Contact number (1.7) Loop contact number (2.7) Isostatic contact number (1.7) Loop contact number (3.8) Continuum Cosserat strain rate tensor (4.43) Cosserat strain rate tensor at c (4.3) z¯ − ziso Permutation tensor (Levi-Civita symbol) (1.11) Stiffness number (1.1) Coulomb friction coefficient (1.6) Continuum loop force (3.27) Loop force on loop , oriented to act on grain g (2.9), (3.1) Loop force on triangle t (2.14), (2.15) Stress tensor (1.8) Area fraction (4.53) Continuum torque potential (3.27) Torque potential on loop , oriented to act on grain g (2.11), (3.2) Torque potential on loop (3.8) Airy stress function (2.4),(3.27) Continued on next page  xi  Symbol ψc ψc ξ ω ωg · (∇ϕ)c (∇v)c (∆ϕ) (∆v)g  Table 1 – continued from previous page Definition Equation Reference Airy stress function on contact c (2.16) Airy stress function on contact c in loop (3.7) trFˆ (3.40) Continuum grain angular velocity (4.24) Grain angular velocity (4.1) Coarse-grained average (3.13) Discrete gradient of ϕ at c (3.18) Discrete gradient of v g at c (4.31) Discrete Laplacian of ϕ at (3.17) Discrete Laplacian of v g at g (4.30)  xii  List of Symbols - Chapters 5-6 Symbol ˆ A C D F [P] f¯ f¯c N  H H Ji c 0 (c) NF C NGC NRC NRG NS NST NV C ns P q S[P] T T X Z Zi ZM F z¯  Definition Angoricity (P − P )(T − T ) Multiple integration Free-energy functional Mean-field contact force − c · f¯c , normal component of f¯c Generalized hyperstaticity Hyperstaticity Cluster Jacobian Contact vector Loop adjacent to c, exterior to cluster Number of contacts on force chains Number of potential contacts Number of real (force-bearing) contacts Number of real (force-bearing) grains Number of sliding contacts Number of spanning trees Number of virtual contacts Number of sliding contacts per grain Pressure log(NST )/NRC Entropy functional ¯ Shear part of σ Topology of contact network Compactivity Partition function Isostatic cluster partition function Mean-field partition function Contact number  xiii  Equation Reference (5.4) (5.7) (5.11)  (5.1) (5.2) (6.16) Figure 3 (6.20) (6.12) (5.2) (5.2) (5.1) (6.8) (6.11) (5.1) (6.44) (6.8) (5.10) (6.44) (5.7) (5.4) (5.7) (6.14) (6.13) (5.2) Continued on next page  Symbol ziso α αc α ˆ β γc ∆¯ z δ[f |V C ] δf δρ δψ η λ ν ρ¯ ρi ρ ρt σ ˆ ψ¯ ζ ω ·  Table 2 – continued from previous page Definition Equation Reference Isostatic contact number, d + 1 (1.7) Isotropic part of α ˆ (6.46) c c α ˆ: (6.26) Inverse of angoricity (5.4) c ) αc (f¯N (6.36) 0 Weight (6.27) z¯ − ziso Virtual contact constraint δ-function (6.8) ¯ f −f ρ − ρ¯ ψ − ψ¯ Shear part of α ˆ (6.46) Contact potential (5.5) Principal axis angle of α ˆ (6.46) Mean-field loop force Cluster loop force (6.22) Loop force on loop (6.7) Loop force on triangle t (6.5), (6.6) Stress tensor (1.8) Mean-field stress function exp( log αc ) (6.35) ¯ Principal axis angle of σ (6.44) Ensemble average (5.14)  xiv  Acknowledgments This thesis would not have been possible without frequent discussions with Neil Balmforth, Christian Schoof, and Ian Hewitt. For hosting me in Cambridge in 2009, and immeasurable help with numerical simulations, I’m indebted to Jim McElwaine. At conferences, I’ve benefited from discussions with Dapeng Bi, Rafi Blumenfeld, Bulbul Chakraborty, Corentin Coulais, Olivier Dauchot, Silke Henkes, Mitch Mailman, Lisa Manning, Farhang Radjai, Brian Tighe, and Alexei Tkachenko.  xv  Chapter 1  Introduction Pour sand on a table and it forms a pile; tilt the table and the sand flows. On the one hand, we have a solid; on the other, a fluid. What are the stresses exerted on the table? To what angle do we need to tilt the table to cause flow? How long is the runout from the flow? Though timeless and trivial to state, these questions still lack reliable quantitative answers. Sand is one of many granular materials, that is, collections of grains sufficiently large that their thermal motion is negligible. Ubiquitous at home and in industry, granular materials pose outstanding challenges in physics, fluid mechanics, and engineering [JNB96, dG99, vH10, FP08, Yu02]. Common to these disciplines is the overarching problem of predicting macroscopic behaviour from knowledge of the grain-scale physics. To restate our elementary questions in more refined terms, we would like to know what equation governs stress transmission within a granular solid, and how large a shear stress need be applied to induce flow. Together, we call these questions the fundamental problem of granular solids. We will only partially address the difficult problem of what happens once flow has commenced. For simplicity, we consider only dry, convex, cohesionless, macroscopic granular materials. Dry means that the grains are immersed in a fluid of sufficiently low density that it can be neglected. Convex means that two grains can touch at, at most, one contact. Cohesionless means that there are no attractive forces between grains, for example due to small liquid films [REYR06]. Macroscopic means that the grains are sufficiently large that electrostatic forces are negligible. Typical examples of such granular materials are sand grains, salt grains, glass beads, and mustard seed. It is a crucial point that under typical experimental conditions the Young’s modulus of the grains, E, is much larger than the external pressure imposed on the system, P , so that the stiffness number [AR07] κ≡  2/3  E P  1.  (1.1)  For example, for sand grains with E ∼ 30GPa and a density of ρg ∼ 3×103 kg/m3 , supporting 20cm of grains under gravity, we have κ ∼ 5 × 104 . For elastic grains of mean linear dimension D, the deflection of a typical contact is D/κ. The largeness 1  of κ suggests that grains can be treated as rigid [MH06]; we will not always do so, but we are always interested in the limit of large κ. A static packing is physically equivalent to one deforming quasi-statically; this can be quantified by considering the inertial number [MiD01] I=  γD ˙ , P/ρg  (1.2)  where γ˙ is the shear rate, P is the pressure, and ρg is the grain density. I measures the importance of inertia; it is a ratio of a microscopic time scale D/ P/ρg to a macroscopic one 1/γ. ˙ For example, if a 20cm thick layer of sand grains, with D = 2mm, is subject to simple shear under gravity with the upper wall moving at 1 cm/s, we have I ∼ 10−4 along the centre-line. We will exclusively consider the quasistatic limit I 1. In the fluid mechanics literature, it is customary to measure the importance of inertia with the Reynolds number Re = U L/ν, where U is a velocity scale, L is a length scale, and ν is the kinematic velocity of the fluid. We can similarly define a granular Reynolds number Re, where for simplicity we will consider a simple shear geometry. We take L = D and ν = τ /(ρg γ), ˙ where τ is the shear stress. For rigid cohesionless grains, the only force scale comes from the pressure, hence τ = µe P for some µe = µe (I). In simple shear flow, the function µe is known empirically and tends to a constant at small I [MiD01]. Finally, there are two velocity scales, deriving from the aforementioned microscopic and macroscopic time scales. It √ is natural to define a velocity scale from the velocity fluctuations U = δV 2 = V 2 − V 2 . The latter exhibit a nontrivial scaling δV 2 = (C0 γD) ˙ 2 I α with α ≈ −1 and C0 a dimensionless constant [MiD01]. We find Re =  UD I 2+α/2 ≈ C0 ν µe (I)  (1.3)  Because of the weak dependence of µe on I, this is monotonic with I, i.e., Re → ∞ as I → ∞ and Re → 0 as I → 0. However, because it involves the empirical functions δV (I) and µe (I), it is more convenient to use I to measure the importance of inertia. Owing to their practical importance, scientific study of granular media goes back at least two centuries1 . Charles Coulomb proposed that along a failure surface, for example the surface of a sandpile, the shear stress τ is related to the 1  We note, however, that the basic law of solid friction on which the Mohr-Coulomb relation is based was discovered by Leonardo da Vinci in the 15th century, and rediscovered by Guillaume Amontons in the 18th century [dV80, Amo06].  2  normal stress σn by τ = σn tan φf + c,  (1.4)  where φf and c are material constants, due to friction and cohesion, respectively 2 [Cou21]. In the engineering literature, (1.4) is known as the Mohr-Coulomb relation, after Christian Otto Mohr, and Charles Coulomb [Ned92]. Before failure occurs, granular materials are typically modelled as isotropic elastic media. The usual equations of elasticity, together with (1.4), then provide a first solution to our fundamental problem. In the last century, a huge body of theoretical, experimental, and computational work, mostly performed in engineering departments, has led to sophisticated refinements of this picture. The latter emphasize the history dependence of stresses, not accounted for in classical elasticity [Woo90]. This program has been largely successful in certain disciplines, in particular soil mechanics [Woo90]. However, more elaborate industrial manipulations of granular materials are often inefficient and dangerous, owing presumably to incomplete understanding of our fundamental problem [dG99]. To advance fundamental understanding, some twenty years ago a new generation of physicists began performing controlled experiments on granular materials [MOB96, VHBC99, GHL+ 01, GLBH01, GRCB03]. In particular, several experiments measured the normal stress under a sand pile, with the surprising conclusion that the stress profile depends on how the sand pile was created [VHBC99, SRC+ 01, GHL+ 01, GLBH01, GRCB03, ABG+ 01]. When created by pouring from an extended source, one finds that the stress is maximal directly below the centre of the pile, as one would naively expect. However, when the pile is created by pouring from a point source, the stress is a local minimum below the centre of the pile, and is instead maximal on a ring some distance from the centre. To explain this phenomenon, a lively debate has ensued on the form of equations governing stress transmission within a granular solid, so far without a satisfactory resolution. One issue brought to the forefront by these experiments is the importance of the microscopic arrangement of the grains to stress transmission within the packing [Oda72a, Oda72b, RB89, CWBC98, dG99, EG99, EG01, BB02, ABG+ 01]. With regard to the stress under a sandpile, it was determined that pouring from a point source causes a dominant anisotropy in the arrangement of grain-grain contacts, known as the packing fabric. This occurs because of continual avalanching of grains from the pile centre during pile formation [ABG+ 01]. In contrast, when poured from an extended source, avalanching occurs in random directions and anisotropy 2  We have tan φf = µe (I = 0) in terms of µe introduced above.  3  is greatly reduced. Hence a fully microscopic theory of stress transmission in a granular solid must incorporate, at least statistically, the packing fabric. From the perspective of physics, the way forward is to reconcile empirical laws with theoretical understanding under the simplest possible circumstances. In this thesis, we attempt to derive equations governing stress transmission in granular materials directly from Newton’s laws. To make this objective precise, we pose the isostatic problem [Ale98, Mou98, EG99, EG01, BB02].  rc rg  Figure 1.1: Sketch of 2D packing of rigid grains, showing the position of a contact, r c , and the position of a grain, r g . The hatched grains are rattlers. Consider a quasistatic frictional packing of N rigid grains in d dimensions (Figure 1.1). Of the N grains, a fraction x0 may be trapped within the packing but not contribute to mechanical stability; these are known as rattlers and will be discussed below. Mechanical equilibrium requires that each of the remaining NRG ≡ (1 − x0 )N grains g is subject to Newton’s laws of force and torque balance: Fg = −  fgc ,  (r c − r g ) × fgc .  0=  c∈C g  c∈C g  4  (1.5)  Here fgc is the contact force exerted on grain g at contact c, r c is the position of c, r g is position of the centre-of-mass of g, and F g is the body force exerted on the grain, typically due to gravity. The real (non-rattling) grains touch at NRC real (force-bearing) contacts. Each contact force is subject to the Coulomb friction inequality |fT | ≤ µfN ,  (1.6)  where fN and fT are the normal and tangential components of the force, and µ > 0 is the microscopic friction coefficient. Some finite number of contacts, NS , will saturate this inequality; these correspond to sliding contacts. For rigid grains, we can consider the grain positions as fixed, and the contact forces as degrees-of-freedom which are subject to the constraints (1.5) and (1.6). Force and torque balance give dNRG and d(d − 1)NRG /2 scalar equations, respectively. For frictional grains we have dNRC − NS degrees-of-freedom in the contact forces. The number of linearly-independent constraints must not exceed the number of degrees-of-freedom, hence assuming linear independence of Newton’s laws and the sliding constraints we must have dNRC − NS ≥ d(d + 1)NRG /2, or z¯ ≥ d1 ns + ziso ,  (1.7)  where z¯ = 2NRC /NRG is the average number of real contacts per real grain, known as the mean contact number3 , ns = 2NS /NRG is the number of sliding contacts per real grain, and ziso = d + 1. A packing with z¯ = ns /d + ziso is said to be generalized isostatic [SvHvS07] 4 . Such a packing is marginally stable in the sense that removing a single contact will lead to grain rearrangement. The isostatic state is thus the boundary between solidlike and fluidlike behaviour in a granular material, and has been likened to a critical point [LN98, OSLN03, LN10]. It occupies a special role in the study of stress transmission, for at isostaticity Newton’s laws can be solved for the contact forces. Given the external forces and the grain positions, there is a unique solution for the contact forces which does not depend on the deformation history, which is in general needed to know the frictional forces. Of course, the constitutive laws at contacts still hold; but at isostaticity all the relevant information is encoded in the geometry: changing the stiffness of a contact will not change the contact force there. This implies that these stiffnesses cannot appear in macroscopic equations. 3  z¯ is also called the mean coordination number. Friction can also be modelled by geometrical asperities on otherwise spherical particles, interacting with normal, repulsive forces, without Coulomb friction. Such packings are found to be exactly isostatic, obviating the need to consider sliding contacts [POS12]. 4  5  To see why this presents a problem, let us take a macroscopic point of view. The macroscopic object of interest is the stress tensor [RS81, KR96] σ ˆ=−  1 V  (r c − r g )fgc ,  (1.8)  g∈G c∈C g  where V is the volume of the packing. If a continuum description exists, then mechanical equilibrium requires that the stress tensor satisfies [LL86, CL00] F =∇·σ ˆ, σ ˆ=σ ˆT ,  (1.9)  where F is the continuum body force. These equations are the macroscopic analog of (1.5), and give d(d + 1)/2 equations for the d2 components of the stress tensor. To solve for the stress tensor, we need a further d(d − 1)/2 equations. For an elastic solid, these are provided by Saint Venant’s compatibility condition, which ensures that the strain can be derived from a displacement field [LL86]. Hooke’s law is then used to write this as an equation for the stress tensor. However, in the isostatic state, we appear to have no Hooke’s law, hence this program fails. It is believed that the isostatic state is characterized by d(d − 1)/2 stress-geometry equations which relate the external forces, the stress tensor, and the packing fabric [EG99, EG01, BB02, Blu04]. Determination of these missing equations is the isostatic problem. The isostatic problem is particularly attractive because, at least in principle, everything is known at the microscopic level. The difficulty only arises in coarsegraining this microscopic information to produce a continuum description. This is in contrast to a general hyperstatic packing (¯ z > ns /d + ziso ) for which the entire deformation history may be relevant. Isostaticity has thus become a central, but controversial, paradigm in the study of granular materials [EG99, Mou98, Rou00, MTDT02, BB02]. Moukarzel and Roux have shown, and it has been confirmed by numerical simulations, that frictionless5 packings of spheres become exactly isostatic in the rigid limit κ → ∞ [Mou98, Rou00, MTDT02, OSLN03]. Numerical studies have shown that frictional packings of spheres generally remain hyperstatic in the rigid limit [AR07, SvHvS07]. However, Kruyt recently performed numerical simulations of extremely large packings of 50 000 polydisperse disks under uniaxial compression [Kru10]. He found that after strains of ∼ 5%, the bound in (1.7) is approached to within 2%, which is on the order of the softness 1/κ. Hence, for frictional spheres in the 5  For frictionless packings, we must modify the above counting, since the contact forces then have only NRC degrees of freedom. For spherical grains, the torque balance conditions become redundant, and hence ziso = 2d. For non-spherical grains, we find ziso = d(d + 1).  6  double limit N > κ → ∞, isostaticity may be attained under typical strains as an attractor of the dynamics. Realistic grains are never perfectly spherical. It has been reported that numerical simulations of frictionless ellipsoids violate the rigidity inequality analogous to (1.7) [vH10]. In the frictionless limit µ → 0, the Coulomb inequality can only be met at real contacts (i.e. those with a non-zero force) by letting |fT | → 0. This forces all d − 1 components of fT to vanish, which amounts to letting NS → (d − 1)NRC . Rearranging (1.7) then gives z¯ ≥ d(d + 1). This can be violated if the constraints are no longer linearly independent. Indeed, vanishing of tangential forces is equivalent to f c = nc nc · f c at each contact, with nc the contact normal. In the spherical limit, contact normals are collinear with contact vectors, and the torque balance equations become redundant; equivalently, the frictionless constraints together with the torque balance equations have NRG linear dependences. This reflects the fact that frictionless spheres can roll without affecting other grains. For aspherical grains, one can no longer see this linear dependence explicitly, but simulations indicate that M linear dependences remain, with M → NRG in the spherical limit [MSOC09, SMCO12]. These dependences are associated with extended rotational modes of deformation. If constraint counting is modified to account for these dependencies, then these packings can be considered generalized isostatic. Unfortunately, so far no theoretical argument can predict M for a given grain shape [vH10, SMCO12]. In this work we exclusively consider frictional granular materials, for which rotation of spheres is no longer a trivial mode and (1.7) is always satisfied. The physical argument for isostaticity is extremely simple. Imagine a static packing formed with grains of some finite stiffness, and consider what happens as the stiffness is increased; this is equivalent to reducing the applied pressure. Because all contact forces are repulsive and there is no cohesion, this process will break unnecessary contacts, i.e., those not needed to support the load. If an essential contact is broken, the packing will rearrange to find a new configuration which supports the load. As decompression ensues, the only way an unnecessary contact can be maintained is for the adjacent grains to move precisely such that the relative velocity at the contact vanishes; mathematically, grain positions must stay on a surface in configuration space of zero measure. For real grains with some amount of geometric randomness, no matter how small, it is believed that this never happens [Mou98, Rou00]. The delicate point of the argument is, of course, the final one. A simple illustrative example is that of a local triangular lattice, which can only be maintained under decompression if grain radii satisfy a certain linear relation. However, in the presence of friction, it is not obvious that such linear relationships will not be maintained, because friction tends to make grains move collectively in clusters. 7  Hence the degree to which isostaticity is applicable to realistic granular materials is still controversial. However, it remains a compelling starting point from which to achieve fundamental progress. In this thesis, we therefore consider granular materials which are isostatic or weakly hyperstatic. Several issues mentioned above deserve elaboration. First, we have neglected from our analysis the rattlers, grains which do not contribute to mechanical stability. When only statics are considered, this definition is tautological; it is only well-defined when dynamics are considered. Operationally, one usually defines rattlers as those grains having fewer than 2 contacts6 , or all of whose contact forces are below some threshold reflecting numerical precision [SvHvS07, AR07, vH10]. In the presence of gravity, this needs to be adjusted [AR07]. A systematic study of such definitions, as well as prediction of the rattler concentration x0 , is still lacking. Second, in deriving the stability inequality (1.7), we have assumed that grains are rigid. This allows us to explicitly separate the geometry from the forces. As discussed above, realistic granular materials are nearly rigid, but we do expect corrections due to finite stiffness. These should be negligible compared to other uncertainties in the problem. Finally, an important but seldom mentioned point is that the rigid limit is singular with respect to ambient thermal fluctuations. Consider an experiment which undergoes fluctuations in the ambient temperature of order ∆T , for example due to diurnal cycles. The importance of these fluctuations is measured by the ratio of linear expansion of a grain to the typical deflection of a contact, ∆δ αL D∆T = , δ D/κ  (1.10)  where αL is the coefficient of linear thermal expansion; αL ≈ 8 × 10−6 K −1 for glass. For κ = 5 × 104 and ∆T = 1K, we have ∆δ/δ ∼ 0.4, which is not at all negligible. This effect is presumably responsible for long-time creep of granular materials [NDBC11]. We will not consider such time-dependent effects.  1.1  Outline of this thesis  In the next chapter, we develop analytic tools, under the rubric of discrete calculus, which allow us to rewrite Newton’s laws (1.5) in a framework in which a continuum limit can be taken. In particular, we show that, subject to some constraints, 6  For frictionless spheres, one can make sharper statements. Hilbert showed that in d dimensions, d + 1 contacts are needed for mechanical stability, hence grains with fewer than d + 1 contacts are rattlers [HCV52].  8  Newton’s laws can be solved by a discrete potential analogous to the Airy stress function of elasticity. This allows us to reformulate the isostatic problem. In Chapter 3, we use discrete calculus to derive a general stress-geometry equation for granular materials in two dimensions. In Chapter 4, we consider the deformation of granular materials and derive equations relating the fabric, deformation, and stress. At the end of Chapter 4, we summarize the results of Chapters 2-4 and present the mean-field phase diagram for quasistatic granular materials. In Chapter 5, we move away from continuum models and instead investigate statistical-mechanical models which take into account local fluctuations. We first give an overview of statistical-mechanical approaches to understanding static granular materials. We present arguments in favour of a statistical mechanics based on the Force Network Ensemble (FNE), the ensemble of mechanically equilibrated force configurations on a given geometrical arrangement of grains. In Chapter 6, a detailed mean-field theory is presented for the FNE, making use of the tools of discrete calculus. Finally, in Chapter 7 we conclude and discuss our results in a broad context.  1.2  Notation  In this thesis we deal with scalars, vectors, tensors, and their pseudo counterparts, defined on grains g, contacts c, loops , triangles t, and tetrahedra v. When unambiguous, these are indicated by a single superscript, but it is often necessary to add subscripts to indicate orientation. These are explained as they appear. The sets of all real grains, virtual grains, real contacts, virtual contacts, triangles, and tetrahedra in the packing are denoted by RG, V G, RC, V C, T , and V , respectively. Together, the real and virtual grains form the set of grains G = RG ∪ V G, and similarly, the real and virtual contacts together form the set of contacts C = RC ∪ V C. We also consider local sets, for example C g denotes the set of contacts surrounding a grain g, and likewise for other neighbouring quantities. The sets of boundary contacts, triangles, and tetrahedra are denoted BC, BT , and BV , respectively. ˆ Outer Tensors are denoted with a hat. The identity tensor is written δ. (dyadic) products are simply denoted with a space, hence for vectors a and b, ˆ = ab is a tensor. For vectors and 2-tensors, operators always act on adjacent A indices, in the same order as they are normally written, and all contractions are indicated with a dot, the number of dots indicating the number of contracted ˆ : ∇(B ˆ (∇ · C)) ˆ is indices; e.g., the kmth component of A i,j,l Aij ∂i (Bjk (∂l Clm )). The 3D cross-product is defined with the Levi-Civita symbol. We take the free ˆ B) ˆ ijm = index to be the first, e.g., (A× k,l Aik εjkl Blm . To conform to convention 9  we make an exception to the above rules for the tensor curl in 3D, letting (∇ × ˆ ij = B) k εikl ∂k Bjl . ˆ × B) ˆ il = We make repeated use of the 2D cross product (A j,k Aij εjk Bkl , with εˆ the 2D Levi-Civita symbol, ε12 = −ε21 = 1, ε11 = ε22 = 0. We frequently ˆ The 2D curl is (∇ × u)ij = use the properties εˆT = −ˆ ε and εˆ· εˆT = δ. k εik ∂k uj , i.e., the 2D curl only differs from the gradient by multiplication by εˆ: ∇× = εˆ· ∇. ˆ be symmetric can be written In both 2D and 3D, the condition that a 2-tensor A ˆ=A ˆT ⇐⇒ εˆ : A ˆ = 0. A  (1.11)  Note that the trace of a 2-tensor can be written ˆ = δˆ : A. ˆ tr(A)  (1.12)  Contours are always oriented anticlockwise, and we use +/− to indicate geometric elements anticlockwise/clockwise from given elements, around grains, loops, etc. For example, in 2D, c+ = c+ (g, ) is the contact anticlockwise from loop when looking from grain g.  10  Chapter 2  Discrete Calculus In a granular solid, mechanical equilibrium requires a delicate balance of forces and geometry. To understand how macroscopic rigidity can emerge in this amorphous solid, it is crucial that we understand how Newton’s laws pass from the disordered grain scale to the laboratory scale. In this chapter, we introduce an exact discrete calculus, in which Newton’s laws appear as differential relations at the scale of a single grain. Using this calculus, we introduce gauge variables which describe identically force- and torque-balanced configurations. In a first, intrinsic formulation, we use the topology of the contact network, but not its geometry. In a second, extrinsic formulation, we introduce geometry with the Delaunay triangulation. These formulations show, with exact methods, how topology and geometry in a disordered medium are related by constraints. In particular, we derive Airy’s expression for a divergence-free, symmetric stress tensor in two and three dimensions. In this chapter we consider a static, frictional packing Ω in dimensions d = 2 and d = 3, in the absence of body forces. The packing is kept rigid by the application of external forces at NEC external contacts. For simplicity, we consider identical disks (d = 2) or spheres (d = 3). Each of the NRG non-rattling grains g is subject to force and torque balance: fgc ,  0= c∈C g  (r c − r g ) × fgc .  0=  (2.1)  c∈C g  As above, fgc is the contact force exerted on grain g at contact c, r c is the position of c, and r g is position of the centre of g. Each contact force must also satisfy the Coulomb inequality (1.6). To relate contact forces to macroscopic objects, we need to consider the stress tensor 1 σ ˆ=− (r c − r g )fgc , (2.2) V g g∈G c∈C  where V is the volume of the packing. In what follows we will decompose this c g c g into contributions from each grain g, σ ˆ g ≡ 1/V g c∈C g (r − r )fg , where V is a volume associated to a grain, discussed below. The notation is described in 1.2. 11  In this Chapter we will be interested in the continuum limit, in which the entire discrete packing Ω shrinks to a point r in the continuum, and (2.2) becomes the value of the field σ ˆ at r. This limit can be taken in many ways, for example by including in (2.2) a smoothing function which is centred at the point r and decays to zero at infinity. It is theoretically convenient to omit such a smoothing function, and instead consider a neighbourhood of r defined topologically. We discuss this more in Chapter 3. As in the Introduction, in the continuum the stress tensor satisfies 0=∇·σ ˆ, σ ˆ=σ ˆT ,  (2.3)  the continuum analog of (2.1). In two dimensions (2D), in the absence of body forces, it implies that σ ˆ can be written as [Mus63, HC09] σ ˆ = ∇ × ∇ × ψ,  (2.4)  where ψ is a scalar, known as the Airy stress function. Here the 2D curl operator is (∇ × ψ)i = ij ∂j ψ, in analogy with the 3D curl operator (∇ × ψ) = ijk ∂j ψk , with ij and ijk the 2D and 3D Levi-Civita tensors, respectively. Hence σij = ik ∂k jl ∂l ψ. In component form, σ ˆ=  ∂yy ψ −∂xy ψ . −∂xy ψ ∂xx ψ  ˆ known In three dimensions (3D), (2.4) also holds, with ψ replaced by a tensor ψ, as the Beltrami stress tensor. A σ ˆ of this form is identically divergence-free and symmetric, so ψ automatically yields stress configurations in force and torque balance. The extension of (2.4) to granular materials was initiated by the seminal works of Satake [Sat86, Sat92, Sat93, Sat97, Sat04] and Ball and Blumenfeld [BB02]. These authors found “loop forces” ρ such that σ ˆ = ∇ × ρ, as well as partial analogs of ψ. However, they were not able to derive ρ from a discrete version of ψ which gives identically torque-balanced configurations, so that the analogy with (2.4) is incomplete1 . In this chapter we extend their results by deriving an exact, discrete analog of (2.4) for frictional sphere packings in two and three dimensions. We first present an intrinsic formulation based solely on the topology of the contact network. We then present a complementary extrinsic formulation, by building the contact network 1  We note that for strictly frictionless packings in 2D, a discrete Airy stress function was introduced by Henkes [Hen09, HC09].  12  into a triangulation of space. The added structure of the triangulation depends heavily on the geometry of the contact network. The new variables ρ and ψ define changes of coordinates in phase space. In this chapter, we consider the rigid limit κ → ∞ in which the deformation of grains can be ignored. This allows us to fix the grain positions and consider the Force Network Ensemble of force configurations on a fixed packing geometry [SVvHvS04, TSVvH10]. The geometry consists of NRG real (force-bearing) grains as well as NV G virtual grains, or rattlers, which are trapped in the packing, but do not contribute to mechanical stability. The real grains touch at NRC force-bearing contacts2 , so phase space is spanned by the NRC contact forces and has dimension dNRC . Force and torque balance, as well as the sliding constraints, restrict force configurations to a subset of phase space, of dimension H = dNRC − NS −  d(d + 1) NRG , 2  (2.5)  which is a measure of hyperstaticity. The sliding constraints will play no role in what follows, so it is convenient to write H = H + NS . Since H depends only on the geometry, it is a topological invariant under any coordinate change in force space. The intrinsic formulation uses the contact network, which assigns a vertex to each grain centre and a link to each contact. Mechanical stability requires that real grains form closed loops . This property implies a combinatorial identity between the number of independent loops NL , real grains NRG , and real contacts NRC , which we can use to rewrite H . To see this identity, we consider an arbitrary grain g and inductively build the entire packing, one loop at a time. Initially, we have one grain, no contacts, and no loops. We choose an arbitrary contact belonging to g and trace out a shortest3 loop back to g, containing, say, k new grains. In doing so we have added k grains, k + 1 contacts, and 1 loop, so NRG − NRC + NL is preserved at its original value, 1. We can continue in this way, tracing out paths from existing grains, through new grains, and back to the existing grains, each time adding a loop when the path rejoins the existing contact network. In doing so, always beginning loops on existing grains, we eventually trace out the entire interior contact network, so that NL = NRC − NRG + 1. 2  (2.6)  These include NEC external contacts at the boundary at which external forces are applied; we assume that each of the NBG boundary grains is subject to one external force. 3 ’Shortest’ means fewest number of contacts.  13  This is a version of Euler’s formula4 for the contact network. We call the set of independent loops chosen by this process a net. In 2D, the net is unique, but in higher dimensions, many nets are possible5 . Letting z¯ = 2NRC /NRG be the average number of contacts per grain, and z¯L = 2NRC /NL be the average number of contacts around a loop, Euler’s formula implies 2/¯ zL = 1 − 2/¯ z + 1/NRC , or 1 1 1 1 = + , + z¯ z¯L 2 2NRC  (2.7)  which explicitly exhibits the duality between grains and loops.  2.1  Intrinsic formulation  The relation (2.6) allows us to write H in two distinguished ways, by eliminating either NRC , or NRG . The first, H = dNL −  d(d − 1) NRG − d, 2  (2.8)  indicates that phase space is spanned by a vector field defined on the loops, providing dNL degrees of freedom, subject to torque balance, providing 12 d(d−1)NRG constraints, and a single extra vector constraint. To exhibit this formulation explicitly, we assign to each loop a fixed orientation, and a pseudovector loop force ρ . Setting ρg = +ρ if the oriented loop enters the grain g, and ρg = −ρ if it exits g, we now write each contact force as fgc =  ρg ,  (2.9)  ∈Lc  where Lc is the set of loops adjacent to the contact c. This definition ensures that Newton’s 3rd law is satisfied. Moreover, since each loop going through g enters at precisely one contact, and exits at precisely one contact, when we sum the contact forces incident on a grain, each loop force appears in equal and opposite pairs, so that force balance is satisfied identically. 4  Euler’s formula states that V − E + F = χ for a graph with V vertices, E edges, and F faces on a surface with Euler characteristic χ. For a finite planar graph, we can embed the graph on the surface of a sphere, identifying the outside of the graph with a face that wraps around the sphere. This amounts to taking χ = 1 in this formula. 5 In graph theory, what we call loops are known as cycles, and NL is called the cyclomatic number, or nullity. Nets are closely related to spanning trees, which can be explicitly counted by Kirchoff’s matrix-tree theorem. See [Bol79].  14  In 2D, each contact is adjacent to two loops, and the loops can all be taken to have anticlockwise orientation. This makes (2.9) a simple difference of loop forces, so that adding a constant to each ρ leaves invariant the contact forces. In this case, there is a vector gauge freedom and so a vector constraint is needed to fix a gauge. However, in higher dimensions, the net is not unique, and there is no natural orientation of the loops. It is not clear in this case whether a gauge freedom exists, or whether the extra constraint is a net-dependent consistency requirement. We may also write H =  d(d + 1) d(d − 1) d(d + 1) NL − NRC − , 2 2 2  (2.10)  in which the number of grains no longer appears. This formulation is achieved by supplementing ρ with a field ϕ , a pseudoscalar in 2D and a pseudovector in 3D, defined so that r c × fgc = ϕg + r × ρg , (2.11) ∈Lc  where r is the position of , defined as the average of the position of the contacts belonging to . As with ρg , ϕg = +ϕ if the oriented loop enters the grain g, and ϕg = −ϕ if it exits g. Summing the torques applied to a grain, it is easily seen that torque balance is satisfied identically. To ensure that torques result from tangential contact forces, the torques inferred from ρ must equal those determined from ϕ , which requires ϕg + r − r c × ρg .  0=  (2.12)  ∈Lc  These are the contact constraints appearing in (2.10). As with the ρ formulation, the residual constraint may depend on the net. In two dimensions (2.9) and (2.11) reduce to Satake’s formulation [Sat93]. The intrinsic formulation can be put into a differential form by defining topological divergence and curl operators. The former, defined for vector fields on c the contacts, is (d∗f )g ≡ c∈C g fg , while the latter, defined for pseudovector and pseudoscalar fields on the loops, is (d∗ρ)cg ≡ ∈Lc ρg . These operators satisfy the identity d∗d∗ = 0. Dropping subscripts and superscripts, we can write Newton’s laws as d∗f = 0 and d∗(r × f ) = 0, (2.9) as f = d∗ρ, and (2.11) as r × f = d∗(ϕ + r × ρ). Continuing in this fashion, one can define a complete topological calculus, with analogs of Stokes theorem and Poincar´e’s lemma. There is a difficulty, however; vector fields, such as f , have no direct continuum meaning, since they change sign 15  with the orientation of the contact. They are, properly, 1-forms. This difficulty may be overcome by pairing f with a vector field; this is precisely what is accomplished by the stress tensor, (1.8). This introduces geometry into the formulation. In 2D the loops tile the plane, so their geometry is simple, but in 3D, they may connect in a complex way. Moreover, in 3D nets are not uniquely defined, so given a contact c, there is no way of knowing how many loops are adjacent to c without some knowledge of the whole net. Rather than proceeding in this way, we therefore present a complementary extrinsic formulation based on a triangulation of space. By using a triangulation, the topology is locally regular and simply related to the geometry. Around any grain, we can determine the topology and geometry of the triangulation knowing only the positions of the grain’s neighbours. In particular, we use the Delaunay triangulation and its dual, the Voronoi tessellation6 , shown in Figure 2.1. 26.5  26.5  26  26  c g  25.5  25.5  25  c  t  24.5  stg  24  scg  25 24.5  g  g  24  23.5  23.5  23  23  22.5  22.5  22  22  21.5 18  19  20  21  22  23  24  21.5 18 25  19  20  21  22  23  24  25  Figure 2.1: Delaunay triangulation (left) and Voronoi tesselation (right) in 2D. In the Delaunay triangulation, the lines are contact vectors, for both real (solid) and virtual (dotted) contacts. The contact vectors cg circulate anticlockwise around the Delaunay triangle t. The vector stg connects contacts belonging to the grain g. In the Voronoi tessellation, the vectors scg circulate anticlockwise around grains g. Although in reality soft spheres deform on contact, for simplicity real contacts are shown with overlaps. The Delaunay triangulation is formed by adding to the contact network links between adjacent grains which do not touch, so that space is partitioned into simplices: triangles in 2-space and tetrahedra in 3-space7 . Physically, the new 6  Also known as the Dirichlet tessellation This is done in 2-space (3-space) in such a way that the circumcircle (circumsphere) formed from each elementary triangle (tetrahedron) contains no other grains. Hence, the vertices of the 7  16  “virtual” contacts correspond to contacts that could be created under a small deformation8 . The Voronoi tesselation is the dual of the Delaunay triangulation; each Voronoi cell is centred on a grain and contains the points nearer to that grain than any other. We define the volume associated to a grain, V g , as the volume of its Voronoi cell. Our general approach to derive an analog of (2.4) is to define an exact discrete calculus on these graphs in which (1.5) appear as differential relations9 . By Poincar´e’s Lemma, ∇ · F = 0 ⇒ F = ∇ × G, these suggest natural forms for ρ and ψ. Routine calculations then show that (1.5) holds identically in the new variables.  2.2  Extrinsic formulation  In the plane, contacts correspond to edges of Voronoi cells; we assign vectors scg along these edges, circulating anticlockwise around grains (Figure 2.1). We also assign vectors cg pointing from the center of g to the center of its neighbour at c. These allow a natural definition for the (signed) area of a contact, Ac = cg × scg . For a tensor field defined on the contacts, we define its divergence by analogy with a tensor version of Gauss’ theorem, defining the resulting line integral with the s vectors: dA ∇ · σ ˆ ≡ Ag (∇ · σ ˆ )g ≡ − g  c∈C g  scg × σ ˆc ≡ −  dr × σ ˆ. ∂g  We emphasize that the underlying fields and definitions are discrete. Because the discrete fields satisfy a similar algebra as in the continuum, the continuum notation is useful for intuition and calculations. We first rewrite (1.8) as a sum over contacts, such that the contribution from each interior contact is Ac σ ˆ c . One finds σ ˆ c = −(1/Ac ) c f c . From this definition we see that scg × σ ˆ c = f c , so (∇ · σ ˆ )g = 0 is equivalent to force balance. By triangles (tetrahedra) correspond to the nearest grain centres. 8 In fact, the triangulation provides a means of defining “small deformation” rigorously, as a deformation which preserves the topology of the triangulation. This partitions the phase space of grain positions into equivalence classes, the boundaries of which correspond to the swapping of virtual contacts across a loop. This will be important in Chapter 7. 9 Elements of a discrete calculus have been reinvented many times in the literature. Our elementary approach emphasizes the connection with classical derivatives, whereas other works use the language of differential forms. Whereas we allow vector- and tensor-valued quantities on each graph element, others always have scalar-valued quantities, interpreting these as vectors or pseudovectors according to the graph element on which they are defined. See [ID91, SMGS99]. For a more mathematical approach, written in the modern language of differential forms, see [DHLM05].  17  inspection of (1.8), and using (1.11), we see that σ ˆ g = (ˆ σ g )T is equivalent to torque balance, and hence we reproduce the expected continuum equations (1.9) in the discrete calculus, at the grain scale. We will define the pseudovector ρ on the triangles, which correspond to vertices of Voronoi cells; we assign vectors stg circulating anticlockwise around grains, connecting the Voronoi vertices. Using a form of Green’s theorem, for a vector field defined on the triangles, we may define its curl as dA ∇ × ρ ≡ Ag (∇ × ρ)g ≡ − g  Noting that 2stg =  stg ρt ≡ −  c+ g  −  c− , g  dr ρ.  (2.13)  ∂g  t∈T g  we see that σ ˆ g = (∇ × ρ)g if −  +  fgc = ρt − ρt ,  (2.14)  where t+ = t+ (g, c) is the triangle anticlockwise from contact c. Summing the contact forces incident on a grain, all loop forces appear in equal and opposite pairs, so force balance is identically satisfied. By working with the Delaunay triangulation, we treat real and virtual contacts on the same footing. Therefore, equation (2.14) applies to both real and virtual contacts, and hence generically leads to virtual contact forces. To obtain physical force configurations, we must explicitly impose that these vanish: if c ∈ V C g , 0 = − + fgc = ρt − ρt . Since virtual contacts only occur in the interior of loops, this shows that ρt must be constant on loops; this also shows that we recover the formulation of Satake [Sat93, Sat97, Sat04] and Ball and Blumenfeld [BB02] if we sum over the virtual contacts. An immediate consequence of σ ˆ g = (∇ × ρ)g is that the stress tensor for the packing can be written as a boundary sum [KR96, HOC07] st ρ t ≡ −  dA σ ˆ ≡ Aˆ σ=− Ω  t∈BT  dr ρ, ∂Ω  where the s’s connect boundary contacts in an anticlockwise manner. We may check that H is conserved under this coordinate change by counting triangles. Applying Euler’s formula to a single loop , NV G − NV C + NT = 1. Adding this relation over the whole packing, we find NV G − NV C + NT = NL , so that, using (2.5), H = 2NT − (2NV C + NRG − 2NV G + 2). This indicates that the loop forces are constrained by the virtual contact constraints and by torque balance, and that 2NV G constraints become redundant. This redundancy is a consequence of the definition of loop forces. Indeed, since loop forces automatically describe grains in force balance, for a virtual grain with z virtual contacts, it is 18  sufficient to impose z − 1 virtual contact contraints to guarantee that all z virtual contact forces vanish. Finally, the loop forces have a gauge freedom ρ → ρ + ∆ρ, with ∆ρ any constant, as is clear from their definition. The torque balance constraint can be rewritten in terms of ρ by first defining dA ∇ · ρ ≡ Ag (∇ · ρ)g ≡ −  stg × ρt ≡ −  g  dr × ρ. ∂g  t∈T g  Comparing this equation with (2.13), and using σ ˆ g = (∇ × ρ)g , we see by inspecg g T g tion that σ ˆ = (ˆ σ ) is equivalent to (∇ · ρ) = 0. This motivates a search for ψ c such that ρt = ∇ × ψ c . We define dA ∇ × ψ ≡ At (∇ × ψ)t ≡ − t  c t  ψc ≡ −  dr ψ. ∂t  c∈C t  It is not obvious from this definition, but as shown in Appendix A, this yields force configurations that identically satisfy torque balance. Since each triangle is surrounded by three contacts, and contacts interior to a loop are shared by two triangles, 3NT = 2NV C + NRC . This yields the global counting relation 3NT = 2NC + NBC , in which BC is the set of triangle edges that circulate around the boundary. We can therefore write H = (NC +NBC )−(2NV C −3NV G +3). The stress function ψ is defined on the real and virtual contacts, and constrained on the virtual contacts by the requirement of no virtual contact force. Because rattlers automatically satisfy force and torque balance, 3NV G constraints are redundant. Finally, ψ has a three dimensional gauge freedom ψ → ψ + ∆ψ, with (∆ψ)c = A · r c + B any linear function of position, which follows from the computation ∇ × (∆ψ)c = εˆ · A. We have therefore succeeded in writing σ ˆ g = ∇ × ρt = ∇ × ∇ × ψ c . It is possible to invert this transformation, up to gauge freedom. For a contour of adjacent triangles (t0 , t1 , . . . , tn ), separated by contacts (c1 , c2 , . . . , cn ), we find n  ρtn − ρt0 =  n  f ci = i=1  tn  sc i × σ ˆ ci ≡  dr × σ ˆ,  (2.15)  t0  i=1  where the forces are exerted from the left side of the contour to the right side. Similarly, n−1  ψ cn − ψ c0 = −  sti × ρti ≡ −  cn  dr × ρ.  (2.16)  c0  i=0  Force and torque balance ensure that the line integrals in (2.15) and (2.16), respectively, are path-independent. 19  Equation (2.16) allows us to relate ψ c to ϕ in the intrinsic formulation. We first allow ψ to take two values on each real contact, one per adjacent loop ± . We define a mean ψ within a loop by z ψ¯ ≡ c∈RC ψ c , where z is the number of contacts around the loop. Then (2.16) implies ψ c = ψ¯ + (r − r c ) × ρ , where z r ≡ c∈RC r c . Comparing with (3.3) we see that ϕ = ψ¯ and the consistency c = ψ c . In the next constraints in the intrinsic formulation are equivalent to ψ+ − chapter, we will exploit this observation. With the Voronoi and Delaunay constructions, it is possible to build an extrinsic formulation in 3D as well. This is presented in Appendix B.  2.3 2.3.1  Discussion and extensions Physical interpretation  As gauge variables, ρ and ψ do not admit a unique physical interpretation. However, in 2D, the inversion formulae suggest a macroscopic interpretation of potential differences. In particular, equation (2.15) indicates that a difference of loop forces gives the flux of stress across any contour connecting two triangles. Similarly, equation (2.16) indicates that a difference of stress function gives the flux of loop force across any contour connecting two contacts. A linear change in stress function gives a constant flux of loop force, and hence no stress. Therefore, stresses correspond to curvature of the stress function. We prove in Appendix C that in 2D, in regions where macroscopic normal stresses are repulsive, the stress function is convex. By considering how a single contact force is written in terms of ψ, it is also possible to obtain a microscopic interpretation. In 2D, an arbitrary contact force fgc can be decomposed as fgc =  1 1 + + A A−  c c gψ  1  + c ∈C c  A  (c )  c  ψc ,  (2.17)  where C c is the set of contacts surrounding the loops adjacent to c, and the c are oriented in the same direction as cg . Considering rigid grains so that the geometry is fixed, we see from this expression that an increase in ψ c directly increases the normal force at c. It also propagates to neighbouring contacts, adding to their contact forces a component in the direction of cg . If the neighbouring loops contain virtual contacts, the disturbance further propagates to the surrounding contacts in such a way that force and torque balance are preserved. We can therefore associate a change in ψ c with a change in the normal force at c, as well as the necessary  20  changes in normal and tangential forces at neighbours of c so as to preserve forceand torque-balance. ρ¯gext ρ  +  c ρ  −  ρ¯g  Figure 2.2: Subset of a Maxwell-Cremona reciprocal diagram. Each polygon corresponds to a single grain. Also shown is the subdivision into triangles for one particular grain, and a triangle corresponding to an external contact. In 2D, one can also obtain some geometrical intuition for these quantities by plotting the loop forces as points in force space, and drawing lines connecting the loop forces of adjacent loops, which are the contact forces. The Maxwell-Cremona reciprocal diagram thus obtained [Max64, Sat92, TvEV08] is a tiling of polygons, each tile corresponding to a single grain (Figure 2.2). The total area of the tiling depends only on the boundary loop forces, and hence only on the boundary contact forces. It can be estimated for large packings, as follows. We define nominal tile centres for each grain, ρ¯g , and divide each tile into − triangles, one for each contact. The signed area of each triangle is acg = 21 (ρ − + ρ¯g )×(ρ − ρ¯g ). For interior contacts, adding the contributions from each adjacent grain, we have +  −  ac ≡ acg+ + acg− = 12 fgc− × (ρ¯g − ρ¯g ).  (2.18)  The total area of the tiling is ac −  A= c∈RC  acgext , c∈EC  21  (2.19)  where the second term sums over the NEC external contacts at the boundary, to correct for overcounting in the first sum. It is clear from Figure 2.2 that A is independent of the choice of ρ¯g , but the relative size of the terms in (2.19) is sensitive to this choice. For example, the trivial choice ρ¯g = 0 makes all ac = 0, and the entire contribution comes from the boundary term. This work has shown that σ ˆ = ∇ × ρ, which implies that the mean-field ¯ variation of ρ is given by ρ(r) = σ ˆ × r in the continuum, up to an irrelevant constant of integration. If we choose ρ¯g = σ ˆ × r g , then, on average, each ρ¯g is in the middle of the tile corresponding to g, as shown in Figure 2.2. The sums √ in (2.19) therefore scale with the number of terms, i.e, as NRC ∼ N , and NEC ∼ N , respectively, and hence for large enough packings, the boundary term is negligible. In fact, with this choice of ρ¯g , the interior term is exactly ac = 21 (ˆ ε·σ ˆ · εˆ) : c∈RC  fgc  c g  (2.20)  c∈RC  = − 12 (ˆ ε·σ ˆ · εˆ) : Aˆ σT  (2.21)  = A det(ˆ σ ),  (2.22)  √ so that A = A det(ˆ σ ) + O( N ). For periodic packings, there are no external contacts EC and hence this result is exact [TV10]. The above choice of ρ¯g does not guarantee that each ac is positive, which might be desirable in applications [TvEV08]. If we choose, instead, ρ¯g = det(ˆ σ )/P εˆ·r g , c where P = 1/2 tr(ˆ σ ) is the pressure, then we find a = − det(ˆ σ )/(2P ) fgc− · cg− . In this representation, repulsion of the contact forces is equivalent to positivity of the areas. Moreover, c∈RC ac = A det(ˆ σ ) as before, so the boundary term must again be negligible for large enough packings. Since the tiling area depends only on the boundary forces, it is invariant under force-balance preserving changes of interior forces. In our formulation, localized force rearrangements can be generated by changing the stress function ψ c at a particular contact, which changes the loop forces at adjacent loops. This construction can therefore be used to explore the Force Network Ensemble, and serves the same purpose as the “wheel moves” of Tighe and collaborators [TSS+ 05, TV10, TV11]. Attempts to extend the reciprocal diagram construction to 3D have not yet been successful [TV10]. In fact, since the loop forces are proper pseudovectors in 3D, they cannot simply be plotted as points in space and therefore cannot define the vertices of any polyhedra. Hence, if an analog of the reciprocal diagram exists, it is not obviously related to the loop forces defined in this work.  22  2.3.2  Polydisperse packings  For clarity we have considered only packings of identical d-spheres, but many of our results extend to polydisperse packings. Indeed, because it is topological in nature, the intrinsic formulation (2.9), (2.11), and (3.3) holds for packings of arbitrary convex grains, without modification. To extend the extrinsic formulation to polydisperse spheres, it is most convenient to use the radical Voronoi tesselation [GF82], which ensures that the edges of Voronoi cells pass through contacts. Discrete derivatives can be defined exactly as in this paper, and the definition of ρ is easily made, so that σ ˆ = ∇ × ρ. However, ψ as defined in this chapter does not describe identically torque-balanced configurations.  2.3.3  Force law  For simplicity, the discussion above was limited to hard disks and spheres, interacting without cohesion. In fact, these assumptions are inessential. We considered hard disks and spheres, and framed this chapter as a discussion of force configurations on a fixed geometry, to show the power of constraint counting. However, the definitions of ρ and ϕ hold for disks and spheres of arbitrary softness. Provided contacts appear at the midpoint of contact vectors, the definition of ψ also holds. We also considered purely repulsive packings. In a repulsive packing, all real grains form closed loops. When grains admit cohesive forces, in addition to the loop-forming grains which bear the external load, there may be tree-like subsets of grains which extend into loops. Starting from the leaves of the tree, it is easy to see that each contact force on each of these grains must vanish. Therefore, provided the grains are treated as virtual, all of the results of this paper hold without modification. The physical novelty of cohesion is that the vanishing contact forces may be composed of an adhesive component, for example due to hydrodynamic forces, as well as an elastic component. This allows the grains to remain in place under an infinitesmal disturbance, unlike the rattlers. More generally, this discussion shows that our results have little to do with the form of the force law and the repulsive constraints at the contacts, although the latter are an essential feature of granular physics [Mou98, Rou00]. Rather, the key physical requirement is that grains interact locally, so that a well-defined contact network exists. Our results then follow from the general form of Newton’s laws, (1.5), acting on this network. Our results may therefore be useful in other problems involving local balance constraints, e.g. rigidity percolation [JT96] and metabolic networks [SVC02, AKV+ 04].  23  2.4  Conclusion  We have shown that Newton’s laws (1.5) can be interpreted as differential relations in a discrete calculus. In an intrinsic topological formulation, these take the form d∗f = 0,  d∗(r × f ) = 0.  (2.23)  These fields can be written in terms of gauge fields ρ and ϕ as f = d∗ρ,  r × f = d∗(ϕ + r × ρ),  (2.24)  which are required to satisfy, for consistency, d∗(ϕ + r × ρ) = r × d∗ρ.  (2.25)  As discussed above, the stress tensor does not appear in this formulation, so the relation (2.4) is absent. Equations (2.23), (2.24) and (2.25) are exact, and hold for polydisperse, cohesive or non-cohesive, soft d-sphere packings. In this paper we considered d = 2 (disks) and d = 3 (spheres), but their extension to higher dimensions is straightforward. To make contact with continuum equations, we also constructed an extrinsic formulation using the Delaunay triangulation and Voronoi tesselation. Using the triangulation, we constructed a discrete calculus in which Newton’s laws (1.5) reproduce their continuum equivalent (1.9) at the scale of a single grain, i.e., (∇ · σ ˆ )g = 0,  σ ˆ g = (ˆ σ g )T .  (2.26)  In 2D, we introduced gauge fields ρ and ψ, defined so that σ ˆ g = (∇ × ρ)g ,  ρt = (∇ × ψ)t .  (2.27)  ˆ To obtain physical force The same relations hold in 3D with ρ → ρˆ and ψ → ψ. configurations, these must be subject to the geometry-dependent virtual contact constraint σ ˆ c |c∈V C = 0. Equations (2.26), and (2.27) are exact, and hold for monodisperse, cohesive or non-cohesive, packings of soft disks or spheres. Together, (2.27) imply that σ ˆ g = (∇ × ∇ × ψ)g , which is an exact, discrete representation of (2.4). On a homogeneous packing geometry, in the continuum ˆ in 3D. From this expression we limit, σ ˆ will satisfy σ ˆ = ∇ × ∇ × ψ, with ψ → ψ see that the pressure P ≡  1 tr(ˆ σ) = d  1 2 2∇ ψ 1 2 ˆ 3 ∇ tr(ψ)  24  −  1 3∇  ˆ · (∇ · ψ)  in 2D in 3D.  (2.28)  We also proved that, in 2D, if macroscopic stresses are repulsive, then ψ is a convex function. These relations are insufficient to completely specify the continuum limit of this problem, as discussed in the Introduction. The discrete representation of ψ developed in this work allows insight into the missing stress-geometry equations [TW99, EG99, BB02]. Indeed, at the microscopic level, ψ is only required to satisfy the virtual contact constraints, and the repulsive constraints. The former are present in all packing geometries but a very special one: the triangular lattice in 2D. On this lattice, the grains are locally and globally as densely packed as possible. In a sense, the geometry is trivial. If we assume that the repulsive constraint requires only that macroscopic stresses are repulsive, then ψ can be any convex function. In general, however, virtual contact constraints exist and couple the ψ field throughout the packing. Their distribution is intimately related to the size and shape of loops. For hyperstatic packings, Newton’s laws are insufficient to fully specify the stress state at the microscopic level, so we expect a family of solutions at the macroscopic level as well. These could, in principle, depend on the microscopic force law and the packing history. In the simplest case, they would depend only on the distribution of virtual contacts. Only at isostaticity is the microscopic stress state fully specified by the geometry. In the next chapter, we use discrete calculus to present a solution to the isostatic problem.  25  Chapter 3  Isostaticity In this chapter we consider the problem of finding the missing stress-geometry equations in an isostatic granular medium. We fix the load, a collection of boundary forces, and the geometry, the positions and orientations of a set of rigid grains. Generalizing from the previous chapter, the grains can be of arbitrary convex shape. As discussed in the Introduction, at isostaticity we expect a unique solution for all the interior contact forces which depends only on the geometry and the load. In terms of the macroscopic stress tensor, in d dimensions d(d − 1)/2 equations are missing, which arise microscopically from Newton’s laws. Our starting point is the intrinsic formulation of Chapter 2. We consider a static packing in the absence of gravity, initially in either dimension d = 2 or d = 3. We define fields ρ and ϕ such that fgc =  ρg ,  (3.1)  ∈Lc  r c × fgc =  ϕg + r × ρg .  (3.2)  ∈Lc  Writing the contact forces and torques in this way, force and torque balance are identically satisfied. However, we must require that the torques computed from ρ are consistent with those computed from ϕ. Writing tc = r c − r , this imposes d(d − 1)NRC /2 constraints ϕg − tc × ρg ,  0=  (3.3)  ∈Lc  d(d − 1)/2 for each contact c. In the limit of large friction coefficient, these are precisely the missing stress-geometry equations, since all other constraints have been satisfied. At finite friction, there are also constraints due to sliding contacts, which hold in addition to (3.3). Hence as a starting point we consider first (3.3). Our goal is to rewrite these equations in such a way that a continuum limit may be taken. In Chapter 2, we saw that continuum expressions are naturally related to sums of discrete expressions around closed contours. Since NRC = NL + NRG − 1, it is natural to sum these equations around grains and loops to form an equivalent set of constraints that are more easily interpreted as continuum equations. 26  We will define discrete integrals over the loops, exactly analogous to those over Delaunay triangles considered in Chapter 2. In fact, each loop integral can be considered as a sum over Delaunay triangles within a loop. The loops have the advantage of duality with respect to the grains in d = 2, and they lead to generalizations of the results of Chapter 2 to grains of arbitrary convex shape (Figure 3.1).  sc  c  g  Figure 3.1: Contact-network tesselation (left) and its dual (right) in 2D. In the contact-network tessellation, the lines are contact vectors, c , circulating anticlockwise around loops . In its dual, the vectors scg circulate anticlockwise around grains g. Summing around a grain g, we find rc ×  0= c∈C g  ρg ,  (3.4)  ∈Lc  since each loop containing g shares exactly 2 contacts with g, which contribute ρ and ϕ in equal and opposite pairs. Recalling (3.1), this is exactly the original equation of torque balance. However, by performing a weighted sum around a loop with weights v c , we find vc  0= c∈C  ϕg − tc × ρg  .  (3.5)  ∈Lc  vc  We will specify below. Euler’s identity NRC = NL + NRG − 1 indicates that these equations must have a degeneracy. Although we have no explicit expression 27  for this degeneracy, it is unimportant since any solution of (3.4) and (3.5) certainly gives a solution to (3.3). These equations hold in both d = 2 and d = 3, and can easily be extended to arbitrary dimension. However, as mentioned in Chapter 2, only in d = 2 can all the loops be oriented consistently1 . To give (3.3) a continuum meaning in d ≥ 3, it appears necessary to sum around loops with antisymmetric tensorial coefficients, which have the interpretation of loop orientations2 . We leave this complication for future work and proceed in d = 2. In this chapter, we will frequently write vectors c without subscripts which indicate how the contact is oriented. All such vectors will be oriented as shown in Figure 3. We will often use the identities A × B = A · εˆ · B and εˆT = −ˆ ε = εˆ−1 .  g tc  c c  sg  Figure 3.2: Notation convention for r − r = tc − tc .  g  c , tc ,  and sg vectors. We also have sc =  In the plane, we can orient all loops anticlockwise; (3.3) becomes 0 = ϕ − ϕ − tc × ρ + tc × ρ .  3.1  (3.6)  Geometric formulation  Let us first establish the relationship between ρ and ϕ. We introduce auxiliary variables ψ c ≡ ϕ − tc × ρ , 1  (3.7)  This is closely related to MacLane’s planarity criterion, a result in graph theory [Bol79, Mac37, O’N73] 2 In the plane, all antisymmetric tensors are multiples of εˆ, hence this is superfluous.  28  which, we will see, are the values of the Airy stress function. Our original equations (3.6) are equivalent to continuity of ψ c at a contact, ψ c+ (c) = ψ c− (c) . The definition (3.7) has a simple geometric interpretation (Figure 3.3). In (r, ψ) space the plane with normal (ˆ ε · ρ , +1) which passes through (r , ϕ ) is described by the equation 0 = (r −r , ψ −ϕ )·(ˆ ε ·ρ , +1) = (r −r )×ρ +ψ −ϕ . Recalling that tc = r c −r , the definition (3.7) says that for each loop, we create a facet of a surface which passes through (r c , ψ c ) for each contact. Continuity of ψ across a contact is equivalent to continuity of the surface at that contact. Hence, in the continuum, ψ is a continuous function. Further properties of the geometric interpretation are discussed in Appendix E. In this discussion, we began with the fields ρ and ϕ and defined ψ from them. Suppose, instead, that we are given ψ. The necessary and sufficient condition that ψ forms a valid discrete stress function is that it varies linearly across loops. Since stress corresponds to curvature of the stress function, this is equivalent to requiring that all of the stress is concentrated on the grains. This is the physical meaning of the stress-geometry equation. We return to this point below. If we sum (3.7) around a loop, we see that ϕ =  1 z  ψc  (3.8)  rc ,  (3.9)  c∈C  provided r =  1 z  c∈C  where z is the number of contacts around a loop. Hence ϕ is nothing but a loop average of ψ. Again summing (3.7) around a loop, but with coefficients c , we find  −1 c c  ρ = −  t · εˆ  c∈C  = (ˆ g )−1 ·  c  ·  ψc  c∈C  −1 A  dr ψ ∂  = (ˆ g )−1 · (∇ × ψ)  (3.10)  where g ˆ = =  1 A 1 A  c c  t · εˆ  c∈C c c  r · εˆ  c∈C  29  (3.11)  is a fabric tensor, the first of several that we will encounter. To understand the behaviour of g ˆ , we use the remarkable geometric identity A δˆ =  dA ∇r = εˆ ·  c1 g 2 (r  dr r = εˆ · ∂  + r g ),  (3.12)  c∈C  which is valid for any closed loop , for which the area enclosed by the contour is A . Here the integrals are in the continuum and we exactly discretize them by considering the contour formed by the contact vectors. We see that if r c = 1 g g ˆ Defining grain centres ˆ· g ˆ · εˆT = δˆ and hence g ˆ = δ. 2 (r + r ) for each c, then ε by centre-of-mass, the condition that contacts lie midway between adjacent grain centres is always satisfied by monodisperse disks. Hence g ˆ − δˆ measures local deviations from a monodisperse disk packing. The averaging properties of fabric tensors are closely related to the symmetries of the problem. At an individual contact, there is a reflection symmetry between the two participating grains, which we call contact symmetry. This implies that the average of any quantity which is antisymmetric at contacts must vanish up to boundary terms. Deviations of the fabric from a monodisperse disk packing satisfy contact symmetry and, moreover, are intensive. Hence averages of g ˆ over ˆ a sufficiently large area always yield δ, even if the fabric is anisotropic. For any field Ψ defined on loops, we define its average Ψ  ≡  1 AΩ  AΨ  (3.13)  ∈Ω  over some domain Ω with area AΩ . To compare with numerical simulations, for each loop we will let Ω = Ωm be the mth coordination shell of , i.e., the set of loops that can be reached from by crossing m or fewer contacts. This is illustrated in Figure 3.4. To illustrate the behaviour of fabric tensors, we performed numerical simulations of highly polydisperse disks undergoing simple shear (Figure 3.5)3 . We construct fabric tensors from the frame shown in Figure 3.5, for which z¯ = 3.1. It is obvious from the picture that the forces are anisotropic; as will be shown below, the fabric is anisotropic as well: there are more contacts in the compressive direction than the dilatant one. However, Figure 3.6 shows that g ˆ converges to ˆ δ as averaging size increases, as anticipated above. 3  Details of the numerics are discussed in Appendix D.  30  Note that with discrete calculus, g ˆ can be written as dr r · εˆ  g ˆ = ∂  = (∇ × r) · εˆT = εˆ · (∇r) · εˆT .  (3.14)  Since ∇r = δˆ in the continuum, convergence of g ˆ to δˆ is consistent with this expression; indeed, it could have been anticipated without reference to the powerful identity (3.12).  31  0  30  Figure 3.3: Geometric interpretation of ψ. The horizontal coordinates are (x, y) while the vertical coordinate is ψ. The vertical scaling is arbitrary. Each facet corresponds to a loop, its colour indicating the loop size, from hot (large loop) to cold (small loop). Dark blue corresponds to triangular loops (z = 3), through lime green for z = 6, to dark red for z = 11. The vertices of the facets are plotted at (r c , ψ c ). Each hole in the surface corresponds to a grain. 32  50  48  46  44  42  40  38  36  34  32  30  −5  0  5  10  15  20  25  30  Figure 3.4: Averaging cells Ωm for m = 0 (red), m = 1 (red+yellow), m = 2 (red+yellow+cyan), and m = 3 (red+yellow+cyan+blue). The mth averaging cell includes those loops that can be reached from by crossing m or fewer contacts.  33  35  50  40  30  20  10  0  Figure 3.5:0 Frame of10a highly polydisperse packing undergoing simple shear. 60Line 20 30 40 50 thickness is proportional to the normal force.  34  1.2 1 0.8  gij  0.6 0.4 0.2 0 −0.2 −0.5  0  0.5  1  1.5 m  2  2.5  3  3.5  ˆ From simulations of highly Figure 3.6: Convergence of fabric tensor g ˆ to δ. polydisperse disks, we compute averages of the fabric tensor g ˆ over neighbourhoods of of m = 0, 1, 2, 3 neighbours; i.e., for m = 2 we compute an average over all loops which can be reached from by crossing 2 or fewer contacts. The plot shows the mean and standard deviation of the individual components of g ˆ (gxx (red,x), gxy (green,o), gyx (blue,+), gyy (black, )). We see that, on average, g ˆ is the identity, and that its fluctuations across the packing decrease with increasing m. For clarity, gyx and gyy have been shifted.  35  3.2  Discrete calculus equations  Now let us return to the consistency equations (3.4) and (3.5). Equation (3.4) becomes dr × ρ ≡ −Ag (∇ · ρ)g ,  sg × ρ ≡  0= +  (3.15)  ∂g  ∈Lg −  where sg = r c − r c circulates around the grain (Figure 3). This already has a natural continuum meaning as (∇ · ρ)g = 0. The loop equations become vc ϕ − ϕ  0=  v c tc × ρ − tc × ρ  −  c∈C  ,  (3.16)  c∈C  which will take some work to be put into a continuum form. We begin with the ϕ term. There is a special choice of v c that makes vc ϕ − ϕ  A (∆ϕ) ≡ c∈C  into a discrete Laplacian4 . To see this, consider the continuum Green’s theorem in the form dA ∆ϕ = −  dr × ∇ϕ.  Ω  (3.17)  ∂Ω  The left-hand side is naturally defined by dA ∆ϕ ≡ Ω  A (∆ϕ) ∈L  vc ϕ − ϕ  =  ,  c∈∂Ω  where = (c) is the loop to the right of c, orienting the sum anticlockwise. The second equality follows since all the contributions from interior contacts have cancelled. We aim to define ∇ϕ and v c so that (3.17) holds identically. A natural definition of ∇ϕ on the contacts is (∇ϕ)c = 4  1 εˆ · Ac  c  ϕ −ϕ  ,  (3.18)  For a detailed discussion of possible weights for discrete Laplacians, see [WMKG07, BS01]  36  with sc = r − r and Ac = sc × c , since this obeys a discrete gradient theorem5 . Indeed, for a contour of adjacent loops ( 0 , 1 , . . . , n ), separated by contacts (c1 , c2 , . . . , cn ), we find n  ϕ  n  −ϕ  0  n  sci · (∇ϕ)ci ≡  =  dr · (∇ϕ).  (3.19)  0  i=1  With this definition, c  dr × ∇ϕ ≡ ∂Ω  × (∇ϕ)c = −  c∈∂C  c∈∂C  | c |2 ϕ −ϕ Ac  .  To enforce (3.17), we should take v c = | c |2 /Ac . Similarly, we can write v c tc × ρ − tc × ρ  vc rc × ρ − ρ  =  − r ×ρ −r ×ρ  c∈C  c∈C  c  =  × (∇ρ)c × r c + (∇(r × ρ))c  c∈C  dr × (∇ρ)c × r c + (∇(r × ρ))c  = ∂  = −A ∇ · (∇ρ)c × r c + (∇(r × ρ))c .  (3.20)  Hence we obtain discrete calculus equations 0 = (∇ · ρ)g ,  (3.21)  0 = (∆ϕ) + (∇ · ((∇ρ) × r + (∇(r × ρ)))) .  (3.22)  As shown in Chapter 2, we also have σ ˆ g = (∇ × ρ)g .  (3.23)  The Coulomb inequality provides additional constraints on solutions of (3.21) c can be written as the pair of and (3.22). Indeed, at each contact, |fTc | ≤ µfN inequalities (δˆ ± µ1 εˆ) : nc f c ≤ 0,  (3.24)  5 Note that (∇ϕ)c can be defined with εˆ · c /Ac replaced by ac /(sc · ac ), for any vector ac , and (3.19) will still hold. The choice (3.18) is fixed by requiring that εˆ · (∇r)c · εˆT gives g ˆ upon appropriate averaging.  37  where nc is normal to the contact. We can write this as (δˆ ± µ1 εˆ) : (nc  c  ·σ ˆ c ) ≥ 0.  (3.25)  By contact symmetry, contact normals and contact vectors are collinear, up to fluctuations. Hence on taking a continuum average (or by considering disks), this can more simply be written σ ˆ c : (δˆ ± µ1 εˆ) ≥ 0. On taking an average, σ ˆ is symmetric, and hence this reduces to positivity of pressure P ≥ 0, or ∆ψ ≥ 0.  (3.26)  We will refine this condition in the next chapter. In the geometric interpretation, the surface ψ has positive mean curvature (see Appendix E). Note that the Coulomb constraint has continuum consequences beyond (3.26) in setting the initial fabric, and boundary conditions; for example, in packings created by deposition it is usually assumed that the Coulomb constraint is saturated at the surface [WCCB96]. Equations (3.21) and (4.33) are exact and hold at each grain g, and loop , respectively. To obtain continuum equations, we need to take averages · , as described above. The averaging cell must be large enough to suppress grainscale fluctuations, but smaller than macroscopic length scales; if there is a wide separation between these scales, we can assume that the packing is homogeneous on the averaging scale. We let ρ= ρ ,  ϕ= ϕ ,  ψ = ψc  (3.27)  and have, as discussed above, g ˆ  ˆ = δ.  (3.28)  The average ρ can be identified with a continuous function evaluated at the point r. With data from numerical simulations, this averaging can be explicitly implemented. However, it is analytically intractable. Hence to proceed theoretically, we will need to make approximations on some averages. Recalling that ρ = (ˆ g )−1 · (∇ × ψ) , this is already in a discrete calculus form ρ = (ˆ g )−1 · (∇ × ψ) .  (3.29)  Since fluctuations of g ˆ are small, it is natural to consider a mean-field approximation ρ = g ˆ −1 · (∇ × ψ) = ∇ × ψ. Likewise, since ϕ = z1 c∈C ψ c , the simplest mean-field approximation is ϕ = ψ. Under these assumptions, the torque balance equation ∇ · ρ = 0 is 38  identically satisfied. Because the packing is homogeneous on the averaging scale, ∇ · (∇ρ) × r = ∇ · (∇ρ) × r and ∇ · (∇(r × ρ)) = ∆(r × ρ); then the loop equation becomes 0 = ∆ϕ+∇· (∇ρ)×r +∇(r ×ρ) = ∆ϕ+∇·(ˆ ε ·ρ) = ∆(ϕ−ψ). Both equations are identically satisfied, hence these mean-field assumptions fail to give a continuum equation [BB02]. We can understand the failure of these mean-field assumptions through the geometric interpretation of ψ. Indeed, consider a smooth surface obtained by interpolating ψ across the grains, as easily imagined in Figure (3.3). The surface is alternately flat and curved on loops and grains, respectively. In passing to a continuum, we forget which parts of space were loops, and which were grains; hence we implicitly replace the original patchwork surface with a coarse-grained smooth surface which is not flat on loops. Each loop shrinks to a point, and the discrete loop equation (3.22) will become an equation valid at the point r . ¯ Denoting this coarse-grained surface by ψ(r), we have, by construction, ψ = ψ¯ in the continuum. The mean-field assumption ϕ(r) = ψ(r) is equivalent to ϕ = ¯ ) in terms of this smooth surface. But positivity of pressure implies that the ψ(r surface has positive mean curvature; in many cases the surface will be convex, as apparent in Figure 3.3. This means that ϕ will systematically deviate from ¯ ), and hence ϕ = ψ in the continuum. ψ(r Discrete calculus allows us to obtain self-consistent corrections to the mean¯ Since the coarse-grained surface ψ(r) ¯ field assumption ϕ = ψ. is assumed smooth, c ¯ ) around ψ(r ¯ ), and compute ϕ = 1 ¯ c we can Taylor expand ψ(r c∈C ψ(r ) z exactly, introducing fabric tensors which characterize the local geometry. Here ¯ we will fit ψ(r) to a quadratic polynomial; this is physically equivalent to homogeneously distributing stress across a loop. Later we will discuss the effect of higher-order terms6 . We have ¯ ), ¯ ¯ ) + (r − r ) · ∇ψ(r ¯ ) + 1 (r − r )2 : ∇∇ψ(r ψ(r) = ψ(r 2 6  (3.30)  For a typical loop, a quadratic polynomial is the lowest order polynomial that fits through all n ¯ ψ . To see this, consider an mth degree polynomial ψ(r) = m ˆ n :: · · · :: (r − r0 )n , in which n=0 a n (r − r0 )n is the n-fold outer product of r − r0 and :: · · · :: is shorthand for n contractions. Each th coefficient a ˆ n is a symmetric n -order tensor, with n+1 degrees-of-freedom, so the polynomial has j(m) = 21 (m + 1)(m + 2) degrees-of-freedom. Within a loop, the requirement that ψ¯ goes through all (r c , ψ c ) gives z constraints. The identity 12 = 1/¯ z +1/ z implies that at isostaticity, z = 6. Since j(1) = 3, j(2) = 6 and j(3) = 15, the lowest-order polynomial which passes through all ψ c has m = 2; i.e., it is quadratic. Higher-order polynomials include effects of neighbouring loops. c  39  and hence ¯ )+ 1 ϕ = ψ(r 2z  ¯ ) tc tc : ∇∇ψ(r c∈C  ¯ ) + 1 Fˆ : ∇∇ψ(r ¯ ), = ψ(r 2  (3.31)  defining a fabric tensor 1 Fˆ = z  tc tc .  (3.32)  c∈C  Note that Fˆ can be rewritten 1 Fˆ = z  tc r c c∈C   1  = z   rc rc − z r r  .  (3.33)  c∈C  Hence, using a mean-field approximation on the term Fˆ : ∇∇ψ , in the continuum we have ϕ − ψ = 12 Fˆ : ∇∇ψ.  (3.34)  Positivity of mean curvature means that the right-hand side is typically positive. ¯ ) with ρ . In Appendix F we show that up to We may likewise compare ρ(r correlations on sub-averaging scales, ρ = ρ¯ to all orders in a Taylor expansion.  3.3  Final equation  Up to correlations on sub-averaging scales, we find ρ=∇×ψ ϕ = ψ + 21 Fˆ : ∇∇ψ.  (3.35)  Hence 0 = ∇ · ρ is identically satisfied. However, the loop equation gives the leading-order stress-geometry equation 0 = ∆ϕ − ∆ψ = 21 ∆ Fˆ : ∇∇ψ 40  (3.36)  If desired, one can use ∇∇ψ = tr(ˆ σ )δˆ − σ ˆ to write this directly in terms of the stress tensor. In view of the identity P = 12 ∆ψ and equation (3.34), the stress-geometry equation has a simple physical interpretation: in a homogeneous continuum which is macroscopically equivalent to the original discrete packing, the pressure in the loops, 21 ∆ϕ, must equal the pressure in the grains, 12 ∆ψ; hence the stress-geometry equation imposes granularity in a literal sense. It is natural to define the average of Fˆ over Ω as Fˆ =  1 2NRC  tc tc .  (3.37)  ∈Ω c∈C  This choice of weighting ensures that, when NRC is fixed, Fˆ = tc tc .  (3.38)  The behaviour of Fˆ under averaging is shown in Figure 3.7, as computed ˆ but from the packing displayed in Figure 3.5. We notice that Fˆ is close to δ, has a systematic deviation which expresses fabric anisotropy. We also notice that fluctuations are much larger than in g ˆ, since there is no exact cancellation on averaging. 0.4  F  ij  0.3 0.2 0.1 0 −0.1 0  1  2 m  3  4  Figure 3.7: Convergence of fabric tensor Fˆ under averaging. As in Figure 3.6, we compute Fˆ over neighbourhoods of m = 0, 1, 2, 3 neighbours. The plot shows the mean and standard deviation of the individual components of Fˆ (Fxx (red,x), Fxy (green,o), Fyx (blue,+), Fyy (black, )). For clarity, Fyx and Fyy have been shifted. A dashed line is plotted at 21 |tc |2 . Our fabric tensor Fˆ differs from fabric tensors commonly appearing in the literature, which are constructed with contact vectors [Oda72a, Oda72b, MNN83, 41  RR05, RDAR12]. A common choice is Fˆ = Since tc /|tc | = εˆ ·  c /| c |  1 NRC  c c  .  c∈RC  up to fluctuations, over a sufficiently large area we have Fˆ ∝ εˆ · Fˆ · εˆT ,  (3.39)  with the constant of proportionality reflecting different normalizations. This implies that over a sufficiently large area Fˆ and Fˆ share principal axes.  3.3.1  Scaling  ¯ In deriving (3.36) we neglected higher-order terms in the Taylor expansion of ψ. By a scaling analysis, here we show that these terms should be irrelevant. For rigid grains there is a single stress scale given by the pressure P . However, there are two length scales. The grain geometry gives a microscopic one, ξ, which we can take as |tc |2 . We also have a macroscopic length scale L coming from boundary conditions. This could be the size of the domain, or the size of the region to which a force is applied. The continuum equation (3.36) can only apply when δ ≡ ξ/L 1. Moreover, in a continuum model all continuum variables must vary only over macroscopic length scales L. Variables and derivatives naturally scale as ψ ∼ L2 P,  Fˆ ∼ ξ 2 ,  ∇∼  1 . L  (3.40)  Consider the Taylor expansion ¯ c ) = ψ(r ¯ )+ ψ(r n≥1  1 c n n ¯ ), (t ) :: · · · :: ∇n ψ(r n!  (3.41)  n  with (tc )n the n-fold outer product of (tc ), and :: · · · :: shorthand for n contractions. The nth term will give rise to an nth -order fabric tensor Fˆn , which scales as ξ n . Hence the nth order term will scale as δ n L2 P . In the macroscopic limit δ → 0, we need only keep leading order terms in δ. Noting that ∆ψ ∼ P,  ∆(Fˆ : ∇∇ψ) ∼ δ 2 P,  (3.42)  we see that the simplest mean-field approximation corresponds to keeping terms up to n = 1 in (3.41). This fails to yield a non-vanishing continuum equation, 42  but we have seen that it is sufficient to work to quadratic order n = 2. Hence the failure of the simplest mean-field approximation is intimately related to the microscopic character of stress localization, as announced above. In the language of statistical field theory, stress localization renormalizes the relationships between ψ, ϕ, and ρ. This is necessary to obtain a nontrivial continuum equation, just as renormalization is used in quantum field theory to remove infinities [Par88, ID91].  3.3.2  Solutions  The stress geometry equation (3.36) is a closed equation for the stress function ψ, given the specification of the geometry via Fˆ and appropriate boundary conditions. Together with the Coulomb constraint, we have ∆ Fˆ : ∇∇ψ = 0,  ∆ψ ≥ 0.  (3.43)  Using ∇∇ψ = 2P δˆ − σ ˆ we can also write ∆ Fˆ : (2P δˆ − σ ˆ ) = 0,  P ≥ 0.  (3.44)  Note that we have derived (3.43) under the assumption that body forces are absent. The latter are easily added at the continuum level, because of linearity. For example, a particular solution to a gravitational body force ∇ · σ ˆ 1 = g is ˆ The Coulomb condition P ≥ 0 applies to the total stress, hence if σ ˆ 1 = g · r δ. σ ˆ=σ ˆ0 + σ ˆ 1 , we have P0 + g · r ≥ 0. Here we consider solutions to (3.43) under simple boundary conditions.  3.3.3  Homogeneous and isotropic fabric  ˆ In this The simplest geometry is an isotropic and homogeneous fabric, Fˆ ∝ δ. case we find ∆P = 0 or ∆∆ψ = 0,  (3.45)  which is exactly the equation satisfied by ψ for an isotropic elastic medium [Mus63]. It is noteworthy that we derive isotropic elasticity without the introduction of any strain variable. Some authors have claimed that isotropic elasticity yields solutions in agreement with experiments, but for values of the Poisson ratio violating the thermodynamic bound ν < 12 [SRC+ 01]. Our result may resolve this paradox, since we obtain ∆∆ψ = 0 without reference to strain, Hooke’s law, or energy7 . 7  See also the force chain splitting model of Bouchaud et al, who derive isotropic elasticity by different means [BCLO19].  43  By introducing complex coordinates, we can write down the general solution to ∆∆ψ = 0 [Mus63]. We let z = x + iy, ∂z =  1 2 (∂x  z¯ = x − iy, ∂z¯ = 21 (∂x + i∂y ).  − i∂y ),  (3.46)  It is easily checked that ∆ = 4∂z ∂z¯,  (3.47)  which allows immediate integration of (3.45). The latter introduces 4 unknown functions, which are reduced to 2 by the condition that ψ(x, y) is real; the result is that ψ(z, z¯) =  z¯f (z) + g(z)  (3.48)  for arbitrary functions f and g, with [f (z)] = 21 (f (z) + f (z)). The form of the potentials f and g is known for many boundary conditions. An important example is the Green’s function for a semi-infinite half-plane y < 0. The boundary conditions are σ ˆ ·y  y=0  = 2πA δ(x),  (3.49)  where the vector A indicates the direction of forcing, and the factor of 2π is there for later convenience. This amounts to 4 scalar equations: on y = 0, x = 0, we have σ ˆ · y = 0, and the traction on an infinitesimal curve around the origin is 2πA. Since ∂yy ψ −∂xy ψ σ ˆ= −∂xy ψ ∂xx ψ these are equivalent to ∂x2 ψ  y=0  = +2πAy δ(x),  ∂x ∂y ψ|y=0 = −2πAx δ(x).  (3.50)  For these boundary conditions, the potentials take the form dg(z) = A¯ log z, dz  f (z) = −A log z,  (3.51)  where A = Ax + iAy , A¯ = Ax − iAy [Mus63]. The stresses described by these potentials are given explicitly in Appendix G, and are plotted in Figure 3.8 for a normal forcing A = (0, 1). 44  14 14  0  0 12  12 −0.2  −0.2  10  10 −0.4 −0.6  8 −0.4  8  6  6 −0.6 4  4 −0.8  −0.8  2  2 −1 −1  −0.5  0  0.5  −1 −1  1  −0.5  0  0.5  1 15  10 0  0  −0.2  5 −0.2  −0.4  −0.4  10 0  −0.6  −0.6 −5  −0.8 −1 −1  5  −0.8  −0.5  0  0.5  1  −10 −1 −1  0 −0.5  0  0.5  1  Figure 3.8: Stresses resulting from a normal point forcing at the origin of a semi-infinite half plane, with isotropic fabric. (σyy (top left), σxx (top right), σxy (bottom left), P (bottom right))  3.3.4  Homogeneous and anisotropic fabric  In the more general case of a homogeneous but anisotropic fabric, we find Fˆ : ∇∇∆ψ = 0.  (3.52)  In the principal axes of Fˆ , this is Fx ∂x2 ∆ψ + Fy ∂y2 ∆ψ = 0, with Fx and Fy the eigenvalues of Fˆ . Empirically, deviations of the fabric Fˆ from δˆ are typically on the order of 10-20%, so we expect that both eigenvalues Fx and Fy are positive [RDAR12]. So long as this is true, this is an elliptic equation for the pressure, and solutions will be qualitatively similar to those from isotropic elasticity. In fact, (3.52) is equivalent to the equation satisfied by ψ in anisotropic elasticity [Sad09]. Because (3.52) only involves a 2nd-order tensor, it is not the most general form possible in anisotropic elasticity: (3.52), which can be written in any dimension, describes materials with orthotropic symmetry, i.e. d perpendicular hyperplanes of symmetry in d dimensions. In d = 2 this amounts to two perpendicular lines of symmetry. Physically, these symmetries represent those of 45  simple shear: reflection symmetry between the positive and negative shear axes, and reflection symmetry about the positive shear axis. Remarkably, the general solution to (3.52) can also be written down [Sad09]. As above, we want to factor the operator Fˆ : ∇∇∆. Let us write Fˆ = F0 (δˆ + η ˆ),  (3.53)  with tr(ˆ η ) = 0. A general η ˆ can be written η ˆ=  a b , b −a  (3.54)  so that − det(ˆ η ) = a2 + b2 > 0. The latter measures the strength of fabric anisotropy and is expected to be much smaller than 1. To factor the operator Fˆ : ∇∇, we look for solutions of the form ψ(x, y) = f (x + µ± y). This implies 0 = 1 + a + 2bµ± + (1 − a)µ2± 1 + µ2± f (x + µ± y) = 0.  (3.55)  Solutions to f (x + µ± y) = 0 give stresses at ∞ and will not interest us; however, we find nontrivial solutions when µ± is a root of the characteristic equation 0 = 1 + a + 2bµ± + (1 − a)µ2± 1 + µ2± . The roots µ± = ± i give solutions depending on z and z¯, as earlier; but we also find µ± =  1 −b ± 1−a  a2 + b2 − 1 .  (3.56)  For a2 + b2 < 1 the roots are complex conjugates, which, along with the condition that ψ is real, leads to the general solution ψ(z, w) =  f (z) + g(w) ,  (3.57)  with w = x + µy, µ = µ+ , and f and g arbitrary functions. We notice two differences from the general isotropic solution (3.48): first, the new complex coordinate w is a more complicated mixture of x and y than z or z¯, involving a rotation whose magnitude depends on the anisotropy; and second, that the general form of the solution is actually simpler. We now have a sum of two functions of different variables, unlike (3.48); the mathematical degeneracy of the latter, obviously related to the root degeneracy of the characteristic equation, represents the physical degeneracy that isotropy implies infinitely many planes of symmetry.  46  The Green’s function for (3.49) is known [Sad09]. We find df (z) = A1 log z, dz  dg(w) = A2 log w, dw  (3.58)  with A1 =  Ax + µAy , π(1 + iµ)  A2 = −  Ax + iAy π(1 + iµ)  (3.59)  In Figure 3.9, we show the stresses described by (3.58) for a normal forcing Ax = 0, Ay = 1, when a = 0.1 and b = 0.3. We notice that stress contours are stretched and rotated by fabric anisotropy, as observed in experiments [SRC+ 01, GHL+ 01, GRCB03]. We also notice that a region is produced in which P < 0, shown in light colours. The Coulomb condition requires that Ptot ≥ 0, where Ptot is the total pressure. If the material were not pre-stressed or subject to gravity, then an instability is produced in these regions, since we require P ≥ 0 everywhere. In general, regions where (3.58) has P < 0 are potentially unstable. In Figures 3.10 and 3.11 we show the stresses described by (3.58) for a shear forcing Ax = 1, Ay = 0, on the same isotropic and anisotropic fabrics. We note a very large region of potential instability.  3.3.5  Inhomogeneous and anisotropic fabric  For a general inhomogeneous fabric, we can expand (3.43) into (∆Fˆ ) : ∇∇ψ + 2(∇Fˆ ) :·∇∇∇ψ + Fˆ : ∇∇∆ψ = 0.  (3.60)  If the fabric varies on a mesoscopic scale L which is smaller than macroscopic scales L, then this equation can admit solutions which differ strongly from those of elasticity. Let us consider an extreme case (∆Fˆ ) : ∇∇ψ = 0.  47  (3.61)  15 0  0  −0.2  15  −0.2 10  −0.4  −0.4  −0.6  5 −0.6  −0.8  −0.8 0  −1 −1  −0.5  0  0.5  10  5  0  −1 −1  1  −0.5  0  0.5  1 20  0  10 0  −0.2  10  −0.2 5  −0.4  −0.4 0  0 −0.6  −0.6  −0.8  −5 −0.8  −1 −1  −0.5  0  0.5  −10  −10−1 −1  1  −20 −0.5  0  0.5  1  Figure 3.9: Stresses resulting from a normal point forcing at the origin of a semiinfinite half plane, with anisotropic fabric. (σyy (top left), σxx (top right), σxy (bottom left), P (bottom right)) 0  0  −0.2  −0.2  −0.4  −0.4  −0.6  −0.6  −0.8  −0.8  −1 −1 0  −0.8  −0.6  −0.4  −0.2  0  0.2  0.4  0.6  0.8  1  −1 −1  −0.2  −0.2  −0.4  −0.4  −0.6  −0.6  −0.8  −0.8  −1 −1  −0.8  −0.6  −0.4  −0.2  0  0.2  0.4  0.6  0.8  1  0  −0.8  −0.6  −0.4  −0.2  0  0.2  0.4  0.6  0.8  −1 1 −1  −0.8  −0.6  −0.4  −0.2  0  0.2  0.4  0.6  0.8  1  Figure 3.10: Stresses resulting from a shear point forcing in the +x direction at the origin of a semi-infinite half plane, with isotropic fabric. (σyy (top left), σxx (top right), σxy (bottom left), P (bottom right))  48  0  20  0  10  15 −0.2  −0.2 10  5 −0.4  5  −0.4 0  0  −0.6  −0.6  −0.8  −5 −0.8  −5  −1 −1  −0.5  0  0.5  −1 −1  1  −10 −15 −0.5  0  0.5  1 20  0  15 0  −0.2  −0.2  −0.4  10 −0.4  −0.6  −0.6  10  0  5 −0.8 −1 −1  −10 −0.8  −0.5  0  0.5  1  0 −1 −1  −20 −0.5  0  0.5  1  Figure 3.11: Stresses resulting from a shear point forcing in the +x direction at the origin of a semi-infinite half plane, with anisotropic fabric. (σyy (top left), σxx (top right), σxy (bottom left), P (bottom right)) This is an algebraic equation in the stress tensor; it says that σ ˆ has principal axes which are rotated π/4 from those of ∆Fˆ . If ∆Fˆ has eigenvalues of opposite sign, it is a wave equation for ψ and describes propagating stresses. It has been argued in the literature that an equation of this form applies to isostatic granular materials [WCCB96, TW99, BB02, Blu04].  3.4  Discussion  At isostaticity, the grain-scale equations exactly determine all contact forces, and hence it is natural to expect that a continuum equation exists which replaces these grain-scale equations. This is to be compared with a hyperstatic granular material, for which static Newton’s laws are insufficient to determine all contact forces, even at the microscopic level: the history of deformation is needed to know the tangential forces. In deriving our central result, the stress-geometry equation (3.36), we never explicitly used isostaticity. Moreover, we made no reference to the locations of contacts which are sliding. However, we have seen that stresses described by the stress-geometry equation are equivalent to those in anisotropic elasticity. On a homogeneous fabric and for general stress boundary conditions, we expect a unique stress solution. Hence in passing from the grain scale to the 49  continuum, the role of hyperstaticity and sliding is relegated to three areas: first, in the validity of the mean-field assumptions made in taking the continuum limit; second, in the parameters of the theory, i.e., the fabric tensor Fˆ ; and third, in boundary conditions. Regarding mean-field assumptions, we will see in Chapter 7 that force correlation lengths decrease away from isostaticity, and hence mean-field assumptions should become better away from isostaticity. Of course, it is possible that correlations alter the picture presented here; we leave this for Chapter 7. The fabric certainly contains information about hyperstaticity. For example, its trace is ξ 2 ≡ tr(Fˆ ) = |tc |2 ≈  2Ag A = , πNL π(¯ z − 2)φ(1 − x0 )  (3.62)  where have used NL πξ 2 ≈ A and A = Ag NRG /((1−x0 )φ) with φ the area fraction, Ag the average area per grain, and x0 the fraction of rattlers. At fixed volume and rattler fraction, this manifestly depends on the contact number of the packing: hyperstaticity decreases the size of loops and hence decreases ξ. Sliding may affect the rattler fraction and the contact number. Indirectly, hyperstaticity and sliding may also affect the orientational degrees of freedom in Fˆ . This leaves boundary conditions. Since we have not derived Hooke’s law or any variant thereof, we can consider only problems with boundary conditions on the stress. It is commonly assumed that when a packing is formed by slow deposition, as a sandpile, that the grains at the surface came to rest from sliding [WCCB96]; if this is the case, then the Coulomb condition will be saturated at the surface and sliding will directly influence the stress boundary conditions. These points call attention to the evolution of the fabric, and the relation of the latter to deformation. We discuss this in the next Chapter.  50  Chapter 4  Deformation In the previous chapter, we derived the stress-geometry equation, which allows the determination of stresses in a quasistatic granular material, if the fabric is known. In the simplest cases, the latter can be assumed homogeneous and its orientation obtained from boundary conditions. However, in general, we need to know how the fabric evolves under flow. Therefore, in this chapter we consider the deformation of quasistatic granular materials. We first introduce the basic kinematic variables involved in deformation, and show how these are related to stress, in the case of rigid grains. The kinematic variables are subject to geometric compatibility conditions, expressing the grain-scale relationships between grain movement, grain rotation, and contact movement. After discussing allowable deformations, we then derive a differential equation for the evolution of the fabric.  4.1  Kinematic variables  Consider a deformation defined by grain velocities v g and angular velocities ω g in the plane. In general, there are two velocities at each contact, corresponding to movement of what was, in the previous instant, the contact point r c , on grains g (c) and g(c). Under deformation, the contact point c on grain g moves at a rate wgc = v g − ω g εˆ · tcg ,  (4.1)  where tcg = r c − r g . The relative velocity at a contact is V c = wgc − wgc = v g − ω g εˆ · tcg − v g + ω g εˆ · tcg .  (4.2)  At the discrete level, v, ω, w, and V have distinct physical meanings. For any grain g, we can move into a Lagrangian frame in which v g = 0 and ω g = 0. Hence the local deformation around g is specified by V c at contacts belonging to g. When V c = 0, the contact c stays fixed under deformation; the local deformation is a combination of rolling and rotation of the contact normal nc . 51  When VNc = nc · V c < 0, grains g and g are moving closer together; when VNc > 0, the grains are separating. Finally, sliding deformation is described by VTc = −nc × V c = 0. To pass from the grain scale to the continuum, these fields need to be given continuum meanings. v, ω, and w have well-defined continuum averages, but V c is antisymmetric at contacts and does not have a continuum meaning. Just as the stress tensor encodes the information from contact forces, which are antisymmetric at contacts, we can create a tensor ˆ c = − 1 εˆ · sc V c , Γ Ac  (4.3)  ˆ the Cosserat strain which has units of strain rate1 . Following Kruyt, we will call Γ rate [Kru03]. As we will see, it differs from the usual strain rate because of grain rotation. ˆ is related to V in a similar way as the stress tensor is related to contact Γ forces. For a contour of grains (g0 , g1 , . . . , gn ) separated by contacts (c1 , c2 , . . . , cn ), we have gn  n  ˆ=− dr · Γ  g0  ci  ×  i=1  1 ci ci s V Ac  n  V ci ,  =  (4.4)  i=1  which should be compared with n  n  sc i ×  dr × σ ˆ=− 0  i=1  1 Aci  ci  f ci  n  f ci ,  =  (4.5)  i=1  for a contour of loops ( 0 , 1 , . . . , n ) separated by contacts (c1 , c2 , . . . , cn ). This ˆ is to V as σ shows that εˆT · Γ ˆ is to f . When an open contour is taken which shrinks to a point in the continuum, ˆ = nc V , where the unit vector n is directed from g0 to gn , (4.4) becomes n · Γ nc is the contact density in the direction n, and V is the mean relative velocity across contacts oriented with contact normal n; this is the analog of Cauchy’s stress theorem. 1ˆ  Γ is essentially the same as Satake’s void strain γ and Kruyt’s Cosserat strain [Kru03, Sat04].  52  ˆ ∼ +ˆ When sliding is predominantly anticlockwise, we have V c ∼ −sc and Γ ε. ˆ Likewise, if sliding is predominantly clockwise, we have Γ ∼ −ˆ ε. Under affine ˆ ˆ ∼ −δ. compression, V c ∼ εˆ · sc and Γ ˆ is further illustrated by energy considerations. Under The significance of Γ the action of external forces F g and torques τ g , the rate at which work is done in a quasistatic deformation is dW = dt  (F g · v g + τ g ω g ) g∈Ω  fgc · v g + (tcg × fgc )ω g  =− g∈Ω c∈C g  fc · V c  =  (4.6)  c∈Ω  1 Ac  =− c∈Ω  c c  f : εˆ · sc V c  ˆc Ac σ ˆc : Γ  =− c∈Ω  ˆ dA σ ˆ : Γ,  ≡−  (4.7)  Ω  the continuum version of which is attributed to Hill [Hil63]. Equation (4.6) is the classical theorem of virtual work, which expresses energy conservation [Sat93, ˆ describes strain which costs energy. Rou00, Kru03, Sat04]. It indicates that Γ  4.2  Constitutive laws  ˆ is directly related to deformation, it is the kinematic variable that is diBecause Γ ˆ and σ rectly related to stress. The relation between Γ ˆ depends on the constitutive laws at contacts; here we consider the simplest case, that of rigid grains. When grains are rigid, contact forces are Lagrange multipliers which are determined by the constraints of grain incompressibility (steric exclusion), repulsion of contact forces, and Coulomb friction. We decompose the relative velocity and contact force at c into their normal and tangential parts V c = +nc VNc + εˆ · nc VTc c  c  f = −n  c fN  c  − εˆ · n  fTc ,  (4.8) (4.9)  where nc is a unit contact normal at c, with the same sign convention as c . Note that V c is the relative velocity of g with respect to g, while f c is the force exerted by g on g. 53  At each contact c, stability with respect to normal movements requires that c > 0 and the normal relative velocity either the contact bears a repulsive force fN c vanishes, VN = 0, or the normal relative velocity is repulsive, VNc > 0, and the c = 0. This condition is called the Signorini condition normal force vanishes, fN and can be written [RBR96] c fN ≥ 0,  VNc ≥ 0,  c fN VNc = 0.  (4.10)  c , then no sliding can occur, Similarly, Coulomb friction requires that if |fTc | < µfN c VT = 0, while if sliding is occurring, then the contact force must be at the boundary of the Coulomb cone, and friction must oppose motion. This can be expressed as c |fTc | ≤ µfN ,  fTc VTc ≥ 0,  c |fTc | − µfN VTc = 0.  (4.11)  It is easy to see that these laws imply dW/dt ≥ 0, i.e., quasistatic flow is dissipative. The equalities in the discrete Signorini and Coulomb laws imply that the f and V fields are strongly correlated at the discrete level, since at each contact, either c = 0 or V c = 0, and similarly for the tangential components. Nevertheless, fN N they can be put into a continuum form with mean-field assumptions, as shown below. We comment on the validity of these assumptions in the conclusion of this chapter. We can write the Signorini and Coulomb conditions in a continuum form by introducing a fabric tensor related to contact normals, ˆ c = nc nc , G  (4.12)  so that Ac ˆ c ˆ c G :Γ sc × nc Ac ˆ c · εˆT ) : Γ ˆc. VTc = c (G s × nc  Ac ˆ c G :σ ˆc nc · c Ac ˆ c · εˆT ) : σ fTc = c c (G ˆc n ·  VNc =  c fN =  For convex grains, we have  c  (4.13) (4.14)  · nc > 0, and we expect that sc × nc > 0 so that  c fN ≥0  ⇐⇒  VNc  ⇐⇒  ≥0  ˆc : σ G ˆc ≥ 0 ˆc : Γ ˆ c ≥ 0. G  (4.16)  ˆc · M ˆ± : σ G ˆ c ≥ 0.  (4.17)  (4.15)  ˆ ± = δˆ ± 1 εˆ. Then we also have Let M µ c fN ± µ1 fTc ≥ 0  ⇐⇒  54  c ≥ 0. Since (ˆ ˆ c = −f c V c /Ac , the Note that these inequalities imply fN σ c )T · Γ equalities of the Signorini and Coulomb conditions can be written  ˆc : G ˆ c, 0 = (ˆ σ c )T · Γ ˆc : M ˆ+ · G ˆ c · εˆ (ˆ ˆc : M ˆ− · G ˆ c · εˆ . 0 = (ˆ σ c )T · Γ σ c )T · Γ The second is nonlinear because of the inherent nonlinearity of |fTc |. Under the simplest mean-field assumptions the above conditions become ˆ:σ G ˆ≥0 ˆ:Γ ˆ≥0 G ˆ·M ˆ± : σ G ˆ ≥ 0,  (4.18)  ˆ : G, ˆ 0= σ ˆ ·Γ ˆ : M ˆ+ · G ˆ · εˆ σ ˆ : M ˆ− · G ˆ · εˆ . 0= σ ˆ ·Γ ˆ ·Γ  (4.19)  and  c V c = 0 could These forms are not unique: for example, the discrete equation fN N ˆ ˆ ˆ also give a continuum equation (G : σ ˆ )(G : Γ) = 0, using (4.13). This equation differs from (4.19) in the mean-field assumptions implicit in writing it. We comment further on this below. ˆ and G ˆ vary smoothly in the continuum, then the domain will be If σ ˆ , Γ, partitioned into regions in which only one direction of sliding is permitted. In this case the nonlinearity in the Coulomb condition can be avoided, because in each domain  ˆ : M ˆ± · G ˆ · εˆ 0= σ ˆ ·Γ  (4.20)  will hold for one of M+ and M− . At the discrete level, it is clear that, in principle, nearby grains may slide in an opposing sense. Hence it is implicitly assumed in the mean-field approximation that in an averaging domain all grains slide with the same sense. Because nc = sc · εˆ/|sc | up to fluctuations, over a sufficiently large packing we have 1 ˆ≡ 1 G nc nc = εˆT · Fˆ · εˆ. (4.21) NRC ξ c∈RC  This can be used to rewrite the constitutive equations in terms of Fˆ , if desired. Equations (4.19) and (4.20) are linear and algebraic and can be solved with long, but straightforward computations. Before delving into this solution, we derive continuum equations which represent the relationship between grain kineˆ matics v and ω, and contact kinematics w and Γ. 55  4.3  Geometric compatibility  ˆ are related through discrete equations (4.1) and (4.2); The fields v, ω, w, and Γ admissible deformations must satisfy these equations of geometric compatibility. To understand admissible deformations in the continuum, these equations must be given a continuum meaning. Let us first show that v g and ω g can be recovered from wgc . Consider qgc wgc = v g  qgc  c∈C g  − ω g εˆ ·  c∈C g  qgc tcg . c∈C g  For convex particles, we can always choose positive coefficients qgc so that the vector c∈C g qgc tcg = 0. Then we have vg = with q g =  c c∈C g qg ,  1 qg  qgc wgc ,  (4.22)  c∈C g  so that v is a weighted grain average of w. Similarly, scg εˆ · tcg ,  scg wgc = −ω g c∈C g  c∈C g  so that ωg =  1 2Ag  1 = 2Ag  scg · wgc c∈C g  dr · w ∂g  = − 12 tr (∇ × w)g ,  (4.23)  where Ag = 12 c∈C g tcg × scg is an area associated to a grain. Under the simplest mean-field closures on ω = ω g , w = wc , and v = v g , these become ω = − 12 tr(∇ × w) v = w.  (4.24)  Note that the relation between ω and w can be rewritten as ε ω εˆ = − 12 tr(∇ × w)ˆ = − 12 εˆεˆ : (∇w)T = 12 (∇w − (∇w)T ), 56  (4.25)  hence grain rotation is equivalent to macroscopic vorticity in the mean-field approximation2 . This equivalence has been observed in two-dimensional experiments [CCL97]. Having related v and ω to w, we need only relate w to V . The form of these geometric compatibility equations will depend on the type of deformation; to see what is possible we reexamine the theorem of virtual work.  4.3.1  Floppy modes  Equation (4.6) indicates that any modes satisfying V c = 0 at all contacts cost no energy, the so-called floppy modes [vH10]. If an imposed boundary deformation is compatible with a floppy mode, then deformation can proceed without any external work being done on the system. We therefore ask, by counting equations, when floppy modes exist. Let us first generalize floppy modes in the presence of sliding contacts. With a finite coefficient of friction, a certain number of contacts, NS , will be sliding, in a direction ucS . McNamara and Herrmann have shown that we can assume that the list of contacts which are sliding is fixed during infinitesimal deformations and only changes when the Coulomb inequality is saturated at a contact [MH06]. The natural generalization of a floppy mode is a mode for which the non-sliding contacts satisfy V c = 0, while those which are sliding satisfy the weaker condition V c ∝ ucS . Assuming linear independence of the constraints, the number of such generalized floppy modes is d(d + 1) NRG − d(NRC − NS ) + (d − 1)NS 2 d = NRG ziso + d1 ns − z¯ , 2  K=  (4.26)  where z¯ = 2NRC /NRG and ns = 2NS /NRG as before. We find that generalized floppy modes exist when z¯ < ziso + ns /d,  (4.27)  hence generalized isostaticity marks the arrival of these modes. Consider an isostatic granular material subject to a slowly varying load. If a single contact which was not previously sliding saturates the Coulomb inequality, and hence begins sliding, then there is a unique floppy mode which governs the ensuing infinitesimal deformation. This property of isostaticity has been discussed extensively by Combe and Roux [CR00, Rou00, RC02]. 2  Here we have used the identity  ij kl  = δik δjl − δil δjk .  57  Now let us find continuum equations for floppy modes. Initially we will consider a strict floppy mode for which V c = 0 at each contact. We will exploit a remarkable similarity between relative velocities and Newton’s laws. Because V c = wgc − wgc , the floppy mode condition V c = 0 is equivalent to continuity of w across a contact. There is a clear analogy between wgc and ψ c in the preceding chapter, laid out in detail in Table 4.1. Although this duality between deformation and stress has been noted by several authors [Sat92, Rou00, Kru03, Sat04], we take it further by introduction of the field ψ. In the present context, duality allows us to derive a deformation-geometry equation analogous to the stress-geometry equation of the preceding chapter. Our strategy, as before, is to rewrite contact constraints as grain and loop equations, and then to put them into a continuum form. grains  loops  vg  ϕ  ωg  ρ  wgc  ψc  wgc = wgc  ψc = ψc  (∇ × ω) = 0  (∇ · ρ)g = 0  (∆v)g = · · ·  (∆ϕ) = · · ·  Table 4.1: Deformation-stress duality Summing V c = 0 around a loop, we find +  −  tcg − tcg  0 = εˆ ·  ωg  g∈G  = εˆ ·  dr ω ∂  = εˆT · A (∇ × ω) ,  (4.28)  so 0 = εˆT · (∇ × ω) = (∇ω) . Summing around grains, with coefficients z c , we find z c v g − v g − εˆ ·  0= c∈C g  z c tcg ω g − tcg ω g . c∈C g  58  (4.29)  As before, we write Ag (∆v)g ≡ −  sc × (∇v)c ,  dr × (∇v) = ∂g  (4.30)  c∈C g  with sc oriented clockwise around the grain, and (∇v)c ≡ −  1 εˆ · sc v g − v g . Ac  (4.31)  Green’s theorem indicates that we should take z c = |sc |2 /Ac . We have  c∈C g g  c∈C g  sc × (∇(rω))c  r c sc × (∇ω)c −  z c tcg ω g − tcg ω g = =A  ∇ · (∇ω)r  g  c∈C g g  −A  ∇ · ∇ ωr  g  (4.32)  which bears an obvious similarity to (3.20). Hence we find exact discrete calculus equations 0 = (∇ω)  (4.33)  0 = (∆v)g + ∇ · (∇ω)r  g  · εˆ − ∇ · ∇ ωr  g  · εˆ,  (4.34)  which are analogous to (3.21) and (3.22). Assuming that the packing is homogeneous on the averaging scale, we can ignore correlations between r and ω, and we find 0 = ∇ω = − 12 ∇tr(∇ × w), ˆ · εˆ 0 = ∆w − ∇ · (ω δ)  (4.35)  The loop equation can be written 0 = ∇ · (∇w · εˆT )T , which is solved by writing (∇w · εˆT )T = ∇ × u for some velocity field u. This can be rearranged to ∇w = (∇u)T .  (4.36)  Separating this into symmetric and antisymmetric parts, we see that ˆ ≡ 1 (∇w + (∇w)T ) = + 1 (∇u + (∇u)T ) D 2 2  (4.37)  ω εˆ = 12 (∇w − (∇w)T ) = − 21 (∇u − (∇u)T ),  (4.38)  ˆ described by u is the same as that described by w, but hence the strain rate D their vorticities are inverted. The grain equation gives 0 = ∆w = ∇ · (∇u)T = ∇(∇ · u). 59  We recognize the Navier-Cauchy equation of linear elasticity for the velocity field u, in the special case of zero shear modulus. Since boundary conditions on w can be converted into conditions on u by inversion of shear directions, we infer that a granular material supporting floppy modes has no resistance to shear. Above, we found that a floppy mode exists when z¯ < ziso , hence we predict that the shear modulus vanishes as z¯ → ziso . This is precisely what is found in simulations of frictional disks [vH10]. We can check that these floppy modes involve grain rotation. Indeed, let us count how many floppy modes there are for which grain rotation is a constant ω g = ω0 . The earlier counting is modified by two factors: there are N − 1 fewer degreesof-freedom in (r, ω), but also 2NL fewer constraints, since the loop equations (4.33) are identically satisfied. The number of such translational modes is then KT = 2(NRG − NRC + NL ) + 1 = 3,  (4.39)  using the identity NRC = NL + NRG − 1. The only modes are the rigid-body motions, and they require NS = 0. Hence the floppy modes described by u involve variable grain rotation, i.e., rolling. At first this appears to contradict the fact that ∇ω = 0, i.e., there is no macroscopic vorticity. The resolution to this apparent paradox is that, topologically, rolling exactly cancels around a loop: the sense of rotation of grains must alternate across each contact, and every loop must have an even number of grains3 . This implies that contact networks supporting floppy modes do not contain triangles of three grains in contact: the latter frustrate rolling. The importance of these microscopic properties of rolling has been stressed by Rivier [Riv06, RARC+ 09]. They appear here in mean-field as the condition ∇ω = 0. The special state in which all loops are even and no contacts are sliding has been called by Rivier the dry fluid, since a granular material in such a state flows without dissipation.  4.3.2  General deformations  For the general situation in which some contacts are sliding, we can impose a relative velocity V c at each contact, which may be zero at many contacts. We then sum the equations of geometric compatibility relating V c to ω g and wc around loops and grains. The right-hand sides of these equations are the same as (4.33) and (4.34) considered above, so we need only sum V c . Summing around a 3 Mathematically, the condition that all loops have an even number of grains means that the contact network is bipartite; grains can alternately be coloured red and blue, such that red grains are neighbours of blue grains and vice-versa.  60  loop, we find c  Vc=  × −  c∈C  c∈C  1 c c sV Ac  ˆ dr × εˆ · Γ  =− ∂  ˆ ), = A (∇ · εˆ · Γ  (4.40)  while summing around a grain with coefficients z c = |sc |2 /Ac , we find zcV c = c∈C g  sc × − c∈C g  1 εˆ · sc V c Ac  ˆ dr × Γ  =− ∂g  ˆ g. = Ag (∇ · Γ)  (4.41)  ˆ modify (4.34) and (4.33) to The terms involving Γ ˆ = (∇ω) (∇ · (ˆ ε · Γ)) ˆ g = (∆v)g + ∇ · (∇ω)r (∇ · Γ)  g  · εˆ − ∇ · ∇ ωr  g  · εˆ,  (4.42)  which are the exact discrete calculus equations for arbitrary deformations of granular materials in the plane. Under the mean-field closures v = w and ω = − 12 tr(∇ × w), we find ˆ + ∇ω 0 = −∇ · εˆ · Γ ˆ + ∇w − ω εˆ . 0=∇· −Γ  (4.43) (4.44)  ˆ hence Using (4.25) we find ∇w − ω εˆ = 12 (∇w + (∇w)T ) = D, ˆ+D ˆ . 0=∇· −Γ We see that the grain equation constrains the strain while the loop equation conˆ + 1 (∇w · εˆT )T ), strains the vorticity. The latter can be written 0 = ∇ · (−ˆ ε·Γ 2 which implies ˆ = 1 (∇w)T − 1 ∇u Γ 2 2  61  (4.45)  for a velocity field u. The grain equation becomes ˆ + 1 ∇w + 1 (∇w)T 0=∇· −Γ 2 2 T 1 ˆ ˆ = ∇ · − Γ + Γ + 2 (∇u)T + 21 (∇w)T ˆ−Γ ˆ T + 1 ∇(∇ · u) + 1 ∇(∇ · w) = −∇ · Γ 2  2  ˆ + 1 ∇(∇ · u), ˆ−Γ ˆ T + 1 ∇(∇ · u) + ∇trΓ = −∇ · Γ 2 2 so that u must solve ˆ−Γ ˆ T − ∇tr(Γ), ˆ ∇(∇ · u) = ∇ · Γ  (4.46)  which is the analog of the equation which, for true floppy modes, resembled linear ˆ By similar manipulations, it elasticity. Here we see no resemblance for general Γ. can be shown that this equation is equivalent to ∆u = ∆w,  (4.47)  ˆ and which may be more convenient if w is known. The 6 degrees-of-freedom in Γ w are reduced to 4 in w and u.  4.4  Fabric  So far we have related stress to fabric, stress to contact deformation, and contact deformation to grain motion. It remains to establish the relationship between the fabric and the deformation. During quasistatic flow, there will be periods of time during which the contact network is maintained, punctuated by instantaneous events of contact opening and contact closing. During the continuous part of deformation, we can simply differentiate d ˆ 1 d F = dt z dt 1 d = z dt  tc tc c∈C  tc r c c∈C    =  =  1 d  z dt 1 z  c∈C   rc rc − z r r  c∈C  dr c c dr c r + rc dt dt 62  −  dr dr r −r , dt dt  where we have used the fact that z is constant during the continuous part of deformation. We need an expression for dr c /dt. For strict floppy modes, the contact velocity w is well-defined at contacts. However, when contacts slide, at each contact we have two velocities, wgc and wgc , corresponding to the movement of each grain ¯ c , depends on participating in the contact. The true velocity of the contact c, w the local curvature of the grains. We can proceed as in the previous Chapter by ¯ c ) about the loop centre. It will be Taylor expanding the continuum field w(r sufficient to expand to linear order. Hence we let dr c ¯ c ) = w(r ) + tc · ∇w(r ) = w(r dt and dr /dt = w(r ), which implies 1 d ˆ F = (∇w(r ))T · dt z  tc r c + c∈C  1 z  r c tc · ∇w(r ) c∈C  = (∇w(r )) · Fˆ + Fˆ · ∇w(r ). T  It is clear that retaining higher order terms in a Taylor expansion will lead to terms involving higher order fabric tensors; as in the previous Chapter, these are expected to be irrelevant for large systems. Under the simplest mean-field assumptions, this gives the continuum equation dFˆ = (∇w)T · Fˆ + Fˆ · (∇w). dt The right hand side describes advection of fabric by grain movement. However, this is only valid so long as the contact network is maintained. When contacts open and break, the fabric tensor changes discontinuously. In Appendix H, we show that if loops and merge into by opening of their common contact c, then the change in NRC Fˆ is δ Fˆ− =  zz c 2z 2z c c (t − tc )(tc − tc ) − tc tc − t t , z z z  (4.48)  where z = z , z = z , and z = z . We can obtain a mean-field form by making replacements z, z → z¯L , c  t − tc → 2tc tc tc → Fˆ− tc tc → Fˆ− , 63  where Fˆ− is the mean fabric which is destroyed by contact opening, and z¯L is the mean loop contact number. Using z = z + z − 2, this results in δ Fˆ− ≈ 2¯ zL Fˆ−  (4.49)  Under the same assumptions, contact creation will lead to a change in NRC Fˆ of −2¯ zL Fˆ+ , where Fˆ+ is the mean fabric which is created by contact closing. These processes modify the fabric evolution equation to d zL J− Fˆ− − J+ Fˆ+ , (4.50) (NRC Fˆ ) = NRC (∇w)T · Fˆ + NRC Fˆ · (∇w) + 2¯ dt where J± are the rates of contact opening and closing, which must satisfy J+ − J− =  dNRC . dt  (4.51)  ˆ by consideration of The quantities J± and Fˆ± should be related to the strain rate D grain-scale deformation. The latter exhibits rich behaviour, including localization of deformation into slip bands [Kuh99]. Without entering into these details, a simple mean-field model which reproduces expected behaviour under compression and shear is   1 ˆ, Fˆ± = 21 ξ 2 ˆ D (4.52) δ∓ ˆ | det D| where ξ 2 = |tc |2 = tr(Fˆ ). The rates of contact opening and closing are determined by the proximity of neighbouring grains and the strength of velocity fluctuations; these are already complicated in isotropic packings, and barely studied in anisotropic situations [MN07, CCPZ12]. However, one simplification is afforded by quasistatic flow: the only time scale comes from the strain, hence we expect   ˆ tr(D)  ˆ J± = NRG I± z¯, φ, | det D|, (4.53) ˆ | det D| for some functions I± , where φ is area fraction. We will not propose a form for these functions; this is left for future work. We can likewise differentiate the relation A = Ag N/φ to see that dφ φ dA =− dt A dt = −φ∇ · v = −φ∇ · w,  (4.54)  using the fact that A∇ · v = dA/dt, which follows from the definition of ∇v. 64  4.5  Summary of Chapters 2-4  Using tools of discrete calculus, in Chapters 2-4 we have developed a nearly complete mean-field theory of rigid, frictional, quasistatic granular materials in the plane. The basic variables in the theory are the stress tensor σ ˆ , the fabric tensors ˆ the contact number z¯, the area fraction φ, the contact velocity w, and Fˆ and G, ˆ the Cosserat strain rate Γ. By linearity of Newton’s laws with stress, we can reintroduce gravity at the continuum level. We write g for the gravity vector, normalized to include density. The above variables are then related by Newton’s laws g =∇·σ ˆ  (4.55)  0 = εˆ : σ ˆ,  (4.56)  0 = ∆ Fˆ : (2P δˆ − σ ˆ) ,  (4.57)  the stress-geometry equation  the frictional rigid-contact constitutive laws ˆ :G ˆ 0= σ ˆ ·Γ ˆ : M ˆ± · G ˆ · εˆ , 0= σ ˆ ·Γ  (4.58) (4.59)  the equations of geometric compatibility ˆ − ∇ω 0 = ∇ · εˆ · Γ ˆ + ∇w − ω εˆ , 0=∇· −Γ  (4.60) (4.61)  geometry evolution equations dFˆ dFˆ ˆ = [F , w, z¯], dt dt d¯ z d¯ z = [w, z¯, Fˆ , φ] dt dt dφ = −φ∇ · w dt and miscellaneous constitutive laws ˆ = 1 εˆT · Fˆ · εˆ G ξ ω = − 12 tr(∇ × w)  (4.62) (4.63) (4.64)  (4.65) (4.66)  tr(Fˆ )  (4.67)  ˆ ± = δˆ ± 1 εˆ. M µ  (4.68)  ξ=  65  Admissible solutions are also constrained by inequalities ˆ:Γ ˆ≥0 G ˆ·M ˆ± : σ G ˆ ≥ 0.  (4.69)  ˆ in terms of potentials ψ and These equations can be simplified by writing σ ˆ and Γ u as σ ˆ = (g · r)δˆ + ∇ × ∇ × ψ ˆ = 1 (∇w)T − 1 ∇u. Γ 2  2  (4.70) (4.71)  ˆ identically solve Written in this way, σ ˆ and Γ g =∇·σ ˆ  (4.72)  0 = εˆ : σ ˆ  (4.73)  ˆ − 1 (∇w)T − 1 ∇w 0=∇· Γ 2 2  (4.74)  and the remaining equation of geometric compatibility becomes ∆(w − u) = 0.  (4.75)  The primary physical ingredient which is presently missing is understanding of the contact opening and contact closing processes, which appear in dFˆ /dt and d¯ z /dt. As a result, the system is not closed. However, even if we were to specify phenomenological closures, it is not obvious how one would go about solving these equations, or whether it is possible to solve all equations simultaneously. To examine these issues, we synthesize the results of degree-of-freedom counting, presented in this Chapter and in the Introduction. We will assume that no further relations exist between the variables besides those implied by the above. In the Introduction, we found that a granular material can support a load in mechanical equilibrium if ∆¯ z ≡ z¯ −ziso −ns /2 ≥ 0. If this is satisfied, we can write σ ˆ in terms of ψ and the stress-geometry equation (4.57) must hold. If ∆¯ z < 0, then mechanical equilibrium cannot hold on a generic geometry: if it holds at all, relations must exist between the grain positions and orientations, for example requiring that these lie on a lattice. Since we here assume that no such relations exist, we infer that the packing cannot support a nontrivial load. In the presence of inertia, we could have a kinetic stress σ ˆ = 0, but in the quasistatic limit, we have σ ˆ = 0. Thus, when ∆¯ z < 0 the rigid-contact constitutive laws are trivially ˆ Note also that the Coulomb inequality is saturated. satisfied, independent of Γ. 66  ∆z¯  =0  )  In this Chapter, we counted generalized floppy modes, which are deformations that preserve the steric exclusions. For rigid grains, these are the only deformations possible. We found that the equations of geometric compatibility can be ˆ can be written in terms of w and u, solved if ∆¯ z ≤ 0. If this is satisfied, Γ and (4.75) must hold. When ∆¯ z > 0, deformations are not possible on a generic geometry: any deformation would involve elastic deformation of grains, which is ˆ = 0 and w = 0, so that the contact forbidden with rigid grains. Thus we expect Γ constitutive laws are trivially satisfied, independent of σ ˆ . In this case the steric exclusion inequality is saturated. We also counted strict floppy modes, which are deformations which cost no energy. We found that these exist when z¯ < ziso . This indicates that when z¯ < ziso , a deformation which costs no energy could release any stress, and hence we expect σ ˆ = 0, even if the packing is inertial. The (generalized) isostatic line ∆¯ z = 0 is marginal with respect both to mechanical equilibrium and geometric compatibility; indeed, these are closely related by duality, as discussed extensively by Moukarzel, Roux, and collaborators [Mou98, Rou00, CR00, Mou25, MTDT02]. Along the isostatic line, the packing can deform while supporting a load: it is a quasistatic fluid. For z¯ > ziso , deformation requires sliding, which dissipates energy: the fluid is dissipative. At the special isostatic point z¯ = ziso , ns = 0, deformation is along strict floppy modes: this is Rivier’s dry fluid [Riv06, RARC+ 09].  ˆ = 0) solid (Γ  flui  d(  ns  ˆ = 0) gas (ˆ σ=Γ  I ziso  z¯  Figure 4.1: Phase diagram for rigid, quasistatic granular materials in the plane. The isostatic point is indicated by I. In the white region, we expect only inertial phases. These results are summarized in the phase diagram of Figure 4.1. We find solid and gas phases, separated by the fluid/isostatic line and the isostatic point I. Note that when σ ˆ = 0, sliding contacts do not give constraints, hence the gas 67  phase is along the line ns = 0. ˆ and σ In the fluid phase, both Γ ˆ are non-zero, in general, and all of the above equations need to be solved together. The 5 degrees-of-freedom in ψ, w, and u are constrained by the 3 scalar equations (4.57) and (4.58), and the vector equation (4.75). We will present some explicit solutions below, modulo the aforementioned uncertainty in the geometric evolution. ˆ = 0, In the solid phase, ∆¯ z > 0, the packing is static and the fabric is fixed: Γ ˆ w = 0, and dF /dt = 0. In this regime, we only need to solve the stress-geometry equation, subject to the Coulomb inequality. If the solution to a given boundaryvalue problem violates the Coulomb inequality, then the packing fails along the ˆ ·M ˆ± : σ line G ˆ = 0, for one of M+ and M− . In the failed region, we expect that the packing is isostatic; here the full fluid equations need to be solved, subject to stress boundary conditions specified by the failure condition. In the gas phase, ∆¯ z < 0, the packing flows without stress. Since sliding ˆ = 0. In this phase contacts pose no constraints in this case, ns = 0 and Γ we need only solve the geometric compatibility condition (4.75), subject to the steric exclusion inequality. As mentioned above, the former is equivalent to the Navier-Cauchy equation of linear elasticity, for a material with zero shear modulus. If the solution to a given boundary-value problem violates the steric exclusion ˆ : Γ ˆ = 0. In inequality, then a jamming transition occurs along the curve G the condensed region, we expect that the packing is isostatic, and the full fluid equations need to be solved, subject to deformation boundary conditions specified by steric-exclusion. Note that the unjamming and jamming transitions need not be the reverse of each other; in fact, hysteresis is expected because the transitions correspond to the saturation of different constraints. Isostaticity is notable for several reasons: first, only at isostaticity are the contact constitutive laws nontrivial, and second, both the steric exclusion inequality and the Coulomb inequality are marginal when approaching isostaticity from the solid and gas phases, respectively. This suggests that the present model, like the critical state model of soil mechanics, can be put into the language of plasticity theory, where similar behaviour is found [Hil98, Woo90]. Finally, the absence of deformations in the solid phase suggests that this phase may be inaccessible to experiment. However, because contacts can be created by infinitesimal deformation, it is possible that real packings slowly creep into this phase without any visible macroscopic deformation; indeed, the finite stability of real frictional packings suggest that this is the case. Also, the finite stiffness of real grains allows deformations which are not floppy modes.  68  4.6  Homogeneous solutions  We now examine some solutions in the fluid phase, where the full system of equations needs to be solved. It is convenient to use Pauli matrices       0 1 0 −i 1 0 , , , ςˆ1 ≡  ςˆ2 ≡  ςˆ3 ≡  1 0 i 0 0 −1 which satisfy ςˆa · ςˆb = δab δˆ + i  abc  ςˆc ,  where abc is the 3D Levi-Civita symbol. We also have det(ˆ ς a ) = −1, tr(ˆ ς a ) = 0. We can decompose any 2-tensor in the plane in terms of Pauli matrices and the identity. We let ˆ = ΓV δˆ + γa ςˆa Γ σ ˆ = P (δˆ + τa ςˆa ) ˆ= G  1ˆ 2 δ + ga ςˆa ,  (4.76) (4.77) (4.78)  ˆ implies that τ2 = g2 = 0. with implied sums over a. Symmetry of σ ˆ and G Substituting these relations into the contact constitutive laws (4.19) and (4.20) gives, after some algebra, 0 = ΓV 1 + 2τa ga + 2iγ2 τ1 g3 − τ3 g1 + γ1 2g1 + τ1 + γ3 2g3 + τ3  (4.79)  and 2ΓV τ1 g3 − τ3 g1 + γ1 2g3 + τ3 + iγ2 1 − 2τ1 g1 − 2τ3 g3 − γ3 2g1 + τ1 2 = ± [ΓV + γ1 τ1 + γ3 τ3 ] . (4.80) µ ˆ · εˆ = A ˆT − δˆ tr(A) ˆ and (4.65) we can check that Using the identity εˆ · A Fˆ : (2P δˆ − σ ˆ ) = Fˆ : (ˆ ε·σ ˆ · εˆT ) = (ˆ εT · Fˆ · εˆ) : σ ˆ ˆ:σ = ξG ˆ, hence the stress-geometry equation is equivalent to ˆ:σ 0 = ∆ ξG ˆ = 2∆ P ξ( 12 + g1 τ1 + g3 τ3 ) . 69  (4.81)  We can also write the Signorini and Coulomb inequalities in these variables. They are  ˆ·M ˆ± : σ G ˆ = 2P  1 2  ˆ:σ G ˆ = 2P 21 + g1 τ1 + g3 τ3 ≥ 0 ˆ:Γ ˆ = 2 1 ΓV + g1 γ1 + g3 γ3 ≥ 0 G 2  (4.82)  + τ1 (g1 ± µ1 g3 ) + τ3 (g3 ∓ µ1 g1 ) ≥ 0.  (4.84)  (4.83)  We now want to write the equalities as partial differential equations in u and w. To do so, from (4.45) we simply compute ˆ = 1 ∇ · (w − u) 2ΓV = tr(Γ) 2 T ˆ 2γ1 = ςˆ1 : Γ = 21 ∂y (wx − ux ) + 12 ∂x (wy − uy ) ˆ = i ∂y (wx + ux ) − i ∂x (wy + uy ) 2γ2 = ςˆ2T : Γ  (4.85)  ˆ = 1 ∂x (wx − ux ) − 1 ∂y (wy − uy ), :Γ 2 2  (4.88)  2  2γ3 =  ςˆ3T  2  (4.86) (4.87)  omitting some intermediate steps. Substituting these into (4.79), (4.80), and (4.81) then gives 3 PDE’s in u, w, and ψ, which are to be solved together with ∆(w − u) = 0 and the inequalities (4.82). Rather than pursuing an analytical general solution, we will look at the simplest case: the homogeneous solutions with σ ˆ, ˆ ˆ G, ∇w and ∇u constant. We will assume that G is fixed and look for solutions ˆ is varied. Note that for homogeneous stress and fabric, the stress-geometry as G equation is identically satisfied.  4.6.1  Simple shear flow  A natural flow problem is homogeneous shear flow w = (w0 y, 0). It is easy to see ˆ we can find a function u = (ax + by, cx + dy) that satisfies that for any constant Γ, ˆ ∆(w − u) = 0 and Γ = 21 (∇w)T − 21 u. Hence w is irrelevant in determining the stresses: it only features in the equation dFˆ /dt = 0, which we are ignoring. This ˆ implies that we need boundary conditions on Γ. c ˆ by considering Boundary conditions on slip velocity V can be applied to Γ ˆ n · Γ = nc V , which relates the mean relative velocity V across contacts oriented ˆ with nc the contact density in the direction n. Under with contact normal n to Γ, homogeneous shear flow we expect that when n = (0, 1), we have V = V0 (1, 0), and when n = (1, 0), we have V = V1 (0, −1); i.e. sliding is clockwise. We will assume that nc has been absorbed into V0 and V1 . The slip velocities V0 and V1 will be similar in magnitude, but it is not obvious how they are precisely related; here we let V1 = βV0 , where 0 ≤ β ≤ 1. Then ΓV = γ3 = 0, γ1 = 21 V0 (1 − β), and iγ2 = 21 V0 (1 + β). In steady simple shear flow, 70  it is observed that the fabric has principal axes at π/4 to the shear direction, so we let g3 = 0 [OK74]. The solution to the constitutive laws is then 4(β 2 + 1)µg1 , −µ + 2µβ − µβ 2 + 4g12 µ + 8µg12 β + 4µβ 2 g12 + 4g1 − 4β 2 g1 µβ 2 − 4β 2 g1 − µ − 4g12 µ + 8βg1 + 4µβ 2 g12 − 4g1 . τ3 = − −µ + 2µβ − µβ 2 + 4g12 µ + 8µg12 β + 4µβ 2 g12 + 4g1 − 4β 2 g1  τ1 =  (4.89) (4.90)  Symmetry requires that the principal stress direction is at π/4 to the shear direction, and hence τ3 = 0, which is achieved for generic µ only when β = 1. The expression for τ1 then reduces to τ1 =  1 . 2g1  (4.91)  Surprisingly, this is independent of µ. We notice that the solution saturates the steric exclusion inequality; this was previously noticed in simulations, and has been discussed in detail by Radjai and collaborators [RDAR12]. In fact, given the earlier discussion of the dynamics in the gaseous phase, it should not be surprising: a fluid phase appears precisely when this inequality is saturated. We also notice that the shear stress appear to diverge as the fabric becomes isotropic. Although we have no explanation for this phenomenon, we note that homogeneous shear flows are always observed with an anisotropic fabric, so these states may not be accessible [RDAR12].  4.6.2  Pure shear flow  Explicit numerical results on stress and fabric are available for pure shear flow, with w = (vx x, −vy y) [AR10, Kru10]. As above, we can always satisfy the geometˆ so that the latter must be constrained ric compatibility equations relating w to Γ, ˆ must take the form separately. By reflection symmetry of pure shear flow, Γ   1 0 ˆ = V0  , Γ (4.92) 0 −β where V0 and β are constants. Again, by symmetry, we must have g1 = τ1 = 0. The constitutive laws then become 0 = (1 − β)(1 + 2τ3 g3 ) + (1 + β)(2g3 + τ3 )  (4.93)  0 = (1 − β) + (1 + β)τ3 .  (4.94)  71  These imply τ3 = (β − 1)/(β + 1) and 0 = g3 (1 − τ32 ). The latter equation indicates that either the fabric is isotropic, g3 = 0, or one of the normal stresses vanishes, τ32 = 1. Neither of these conclusions are consistent with simulations: in pure shear, fabric anisotropy develops between compressive and expansive axes, and normal stresses are nonzero in both directions [AR10, Kru10]. This implies that either homogeneous solutions do not exist for pure shear, or that the mean-field assumptions leading to the contact constitutive laws are invalid, at least in this case. Future work should probe these questions.  4.6.3  Isotropic expansion  Another simple flow is isotropic expansion or compression, γ1 = γ2 = γ3 = 0. We find that − 12 = τ1 g1 + τ2 g3 ,  τ1 g3 = τ3 g1 ±  1 . µ  The steric exclusion inequality requires that ΓV ≥ 0, i.e., only expansion is allowed. This is an expression of Reynolds’ dilatancy and a direct consequence of grain rigidity. We find that the Coulomb inequality becomes ±  P ≥ 0, µ2  which must hold for both ± = + and ± = −. This is only possible if P = 0, which means the solution is trivial. Physically, if we try to homogeneously expand a packing of rigid grains, it loses its ability to support a stress immediately; because grains lack cohesive forces, this is not surprising.  4.6.4  Conclusion  The simple solutions presented here satisfy the stress-geometry equation trivially, ˆ they offer and because w can always be found which is compatible with the given Γ, little insight into the deformation equations. However, they show the importance of the steric exclusion and Coulomb inequalities in constraining solutions. The results qualitatively match with expectations from experiment, for simple shear and isotropic expansion, but not for pure shear. ˆ w, Fˆ , As discussed above, the relation between the geometric variables Γ, ˆ z¯, and ns is incomplete. Further investigation into the contact opening and G, contact closing processes is clearly an important direction for future work.  72  Chapter 5  Hyperstaticity In the rigid limit, frictional granular materials are generically hyperstatic: for any given arrangement of grains, there is a family of solutions of contact forces which satisfy the constraints of mechanical equilibrium. As before, we consider a packing of NRG force-bearing grains, touching at NRC = NRG z¯/2 force-bearing contacts. Force and torque balance on each grain constrain contact forces to lie in a space of dimension d(d + 1) H = dNRC − NRG 2 d z − ziso ), (5.1) = NRG (¯ 2 with ziso = d + 1. Sliding constraints further reduce the dimension of available contact force solutions to d(d + 1) H = dNRC − NS − NRG 2 d (5.2) = NRG (¯ z − d1 ns − ziso ), 2 where ns = 2NS /NRG . In the Introduction we argued, on physical grounds, that granular materials should be close to generalized isostatic; equivalently, the degree of generalized hyperstaticity H should be small. But so long as z¯ is bounded away from d1 ns +ziso , the dimension of available force configurations is extensive, and hence diverges in the macroscopic limit. It is not obvious that fluctuations in this large space do not play a role. In the mean-field models of the previous Chapter, we did not take into account these fluctuations. Here and in Chapter 6, we will investigate the role of fluctuations using statistical mechanics. These Chapters are independent of Chapters 3 and 4, though we will use the tools developed in Chapter 2. If fluctuations are important, macroscopic properties need to be determined by an ensemble average over all microscopically admissible states [Cha10, CRST12]. This requires knowing the appropriate phase space, and the frequency with which each state is visited. We address these points in turn. For simplicity, in Chapters 5 and 6 we will consider only packings of spheres or disks. 73  5.1  Phase space  Subject to a transient external forcing, a frictional granular medium rapidly dissipates energy through friction and inelastic collisions and forms a static packing in metastable mechanical equilibrium. Hence the dynamics can be viewed as a process of energy dissipation. Due to grain stiffness, there are two widely separated energy scales associated with quasistatic dynamics. Recall that for elastic grains of Young’s modulus E, subject to a pressure P , we define the stiffness number κ = (E/P )2/3 1, such that the typical deflection of a contact is δ = D/κ, where D is grain diameter. The energy needed to change a normal force by an amount on the order of its mean is ED3 (δ/D)5/2 = ED3 κ−5/2 . At each contact, tangential forces are bound to normal forces by the Coulomb constraint |fT | ≤ µfN . Normal force rearrangements will often lead to contact mobilization, and subsequent grain rotation will in turn adjust tangential forces. For spheres, grain rotation is sharply distinguished from centre-of-mass movement1 . As discussed in the previous chapter, grain movement near jamming is primarily perpendicular to contacts, which therefore need to be mobilized. Hence to move a grain by a distance on the order of its diameter requires an energy |fT |D ∼ µP D3 in a jammed, frictional packing2 . The ratio of these energy scales is 1/(µκ) 1. Since both grain rearrangement and force rearrangement are effected by grain movement, force configurations will be explored with a much greater frequency than geometrical configurations. Hence there are two time scales: a fast time scale, on which nearest-neighbour relationships are preserved and grains are fixed on grain-size length scales, and a slow time scale, on which the geometry evolves. On the fast time scale, the evolution is in an ensemble in which contacts can open and break, and contact forces can change [Bou02, SVvHvS04]. This separation of time scales has been explicitly observed in experiments [LDBB08, CBD12]. Hence the dynamical variables are the contact forces f = {f c } and the status (open/closed) of each contact. For each choice of closed contacts, T , the phase space is O(T ), the space of contact force configurations which satisfies Newton’s laws on each force-bearing grain and the Coulomb inequality at each force-bearing contact. O has dimension H . We will consider separately the dynamics of opening and closing contacts, and that of contact force evolution on a fixed contact network. Initially, we consider 1 Near jamming, non-spherical grains have rotational modes primarily involving grain rolling, whose energies are intermediate between translational modes and normal force rearrangement. The present argument needs refinement in this case. 2 In a frictionless packing, grain movement is hindered by random inelastic collisions, but the energy needed to overcome these depends strongly on the local geometry.  74  the latter. For each c, the Coulomb inequality for c defines a convex set V c in d dimensions. O is the intersection of this family with the subspace satisfying Newton’s laws. The geometry of O provides basic constraints on the dynamics. As shown by Unger et al [UKW05], O is convex: Consider any two solutions to Newton’s laws and the Coulomb inequality, f1 = {f1c } and f2 = {f2c }. O is convex if every point on the straight line between them is also a solution. Let f = f1 + λ(f2 − f1 ), with 0 ≤ λ ≤ 1. By linearity, f solves Newton’s laws. For each c, V c is convex, hence each f c lies in these convex sets. The field f thus lies in the intersection of all the V c , hence f lies in O. Convexity is important because it is impossible to get lost in a convex space. Roughly speaking, if states within a convex space have a similar a priori probability of being visited, then one expects that the states are actually visited with this probability, even in the limit of an infinite system; i.e., no glass transition takes place in which ergodicity is broken and the dynamics becomes stuck in a subset of phase space. However, we have so far neglected the presence of sliding contacts. The latter correspond to saturation of the Coulomb inequality. Since the Coulomb inequality delineates the boundary of O, if NS contacts are sliding, then configurations lie on a subset of the boundary of O, with dimension H = H − NS . It is easy to see that the boundary of a single Coulomb cone is not convex, hence a glass transition is possible if a certain number of contacts are forced to be sliding. However, under the rigid-grain dynamics discussed in Chapter 4, (4.11), the direction of sliding is fixed once sliding commences. The subset of the Coulomb cone with a given sliding direction is convex, hence if we specify the direction of sliding at each contact, then phase space is convex. Consider a horizontal layer of compacted granular material, initially with no contacts sliding. We slowly tilt the layer under gravity, inducing a shear stress on the packing. We assume the aforementioned separation of time scales so that the forces can be considered evolving on a fixed geometry. Then as the shear stress is increased, the force network is explored, with tangential forces increasing to support the load. This will lead to saturation of the Coulomb constraint, and hence sliding, at some contacts. Eventually the packing fails and flow is induced. This process can be viewed as a hierarchical shrinking of phase space: each time a Coulomb constraint is saturated, force networks are restricted to a subset of one fewer dimension. Since the direction of sliding at each contact is not expected to change during this process, phase space remains convex throughout. One then expects failure when the phase space becomes zero-dimensional, i.e., when generalized isostaticity is reached. The proposal that failure corresponds to generalized isostaticity was previously advanced by Wyart, who used the hypothesis to deduce the dependence of a horizontal layer’s avalanche angle on its 75  thickness, finding good agreement with experiments [Wya09]. We conclude that the relevant phase spaces for quasistatic dynamics are T and O(T ), with dynamical variables f and the status (open/closed) of each contact. Moreover, so long as sliding directions are maintained during the dynamics, then contact forces evolve within a convex space, and we expect force configurations to be ergodically explored.  5.2  Measure  For each T , quasistatic force configurations necessarily lie in the space O(T ), and in fact they will generically be on the boundary of O(T ). However, we do not know the frequency with which each state is visited. This frequency, or equivalently, the measure on phase space with which we take ensemble averages, must either be postulated phenomenologically, or derived from the governing equations of motion [Kur00]. It is helpful to compare the situation with usual Hamiltonian mechanics. There, the dynamical variables q and p evolve under a dynamics that conserves energy E. The Hamiltonian structure of the equations of motion can be used to show that the dynamics preserves the canonical measure exp(−βH(q, p)), where H is the Hamiltonian function, and β a constant3 (Liouville’s theorem). For a general dissipative dynamics, as applies to granular materials, there is no known measure which is preserved by the dynamics [Kur00]. Recent work has set constraints on possible measures of nonequilibrium steady states in dissipative systems, but these are yet to be applied to granular materials [HS01, Kur05, PEKK12]. So we must proceed phenomenologically.  5.2.1  Edwards’ approach  Let us first consider approaches which appear in the literature. More than 20 years ago, Edwards proposed that static granular materials can be described by an ensemble over all mechanically stable configurations of rigid grains that fill a certain volume V [EO89]. Later, it was realized that the macroscopic stress σ ˆ can also be imposed independently of V , hence the full microcanonical Edwards ensemble considers all configurations of grain positions and orientations that fill V , as well as contact forces that support σ ˆ [Edw05]. Edwards proposed that each such configuration in mechanical equilibrium be considered with equal probability. 3  Of course, β = 1/(kB T ) in terms of Boltzmann’s constant kB and temperature T .  76  Hence the microcanonical Edwards probability density is Pmc (r, ω, f ) =  1 δ V − V [r, ω] δ σ ˆ −σ ˆ [r, f ] Θ[r, ω, f ], Zmc  (5.3)  where δ(·) is a Dirac delta-function, Θ is a function which evaluates to unity if the configuration is in mechanical equilibrium, and zero otherwise, and Zmc is a normalization constant. In Edwards’ approach, the configuration space is much larger than O, since it considers all possible geometric configurations in addition to force configurations. The expected value of any macroscopic observable A(V, σ ˆ ) can be computed within this ensemble as Dr  A =  Dω  Df A[r, ω, f ] Pmc (r, ω, f ),  where Dr is shorthand for integration over all r g . Inserting 1 = and 1 = dV¯ δ(V¯ − V ) into this expression, we find A =  1 Zmc  Dr  Dω  ¯ dσ  ¯ σ ¯ −σ dσδ( ˆ)  ¯ ¯ ¯ eS(V ,σ) dV¯ A(V¯ , σ) ,  with ¯ = log S(V¯ , σ)  ¯ −σ Df δ(σ ˆ [r, f ])δ(V¯ − V [r, ω])Θ[r, ω, f ]  the microcanonical entropy. As in usual Hamiltonian mechanics, the microcanonical entropy is extremely difficult to compute, because of the δ-functions imposing the macroscopic constraints. Moreover, this is exacerbated by the complexity of Θ, to which we will return later. It is therefore natural to pass to a canonical ensemble in which fixed macroscopic observables V and σ ˆ are replaced by fixed temperature-like variables com−1 ˆ = V (∂S/∂ σ pactivity X = (∂S/∂V ) and angoricity A ˆ )−1 [Edw05]. Note that angoricity is a tensor. Given the validity of a microcanonical ensemble, the justification for the canonical ensemble is standard and discussed in, for example, [BDD06, HC09]. The canonical Edwards probability density is Pc (r, ω, f ) =  1 −V [r,ω]/X−V α:ˆ ˆ σ [r,f ] e Θ[r, ω, f ], Zc  (5.4)  ˆ−1 . where α ˆ=A The canonical form of the probability density implies relationships between the mean values and fluctuations of V and σ ˆ , which can be used to test the 77  validity of (5.4) [MK02, HOC07, MRDR+ 09]. Numerical simulations suggest that for frictionless packings, fluctuations of stress indeed have the canonical form indicated in (5.4), but tests of the volume fluctuations have been less conclusive [MK02, HOC07, MRDR+ 09]. Recently, both volume and stress dependence of (5.4) has been indirectly tested with two-dimensional experiments [PD12]. These authors found that the stress fluctuations follow (5.4), but the volume fluctuations do not. For small systems, it is possible to explicitly enumerate all the possible mechanically stable microstates and measure the frequency with which these states are accessed, for a given preparation protocol. This has been attempted for packings of frictionless disks, for N ≤ 100, in constant volume experiments and simulations [ABOS12, GBO06, GBOS09]. These authors found that microstates were accessed with a frequency that was far from uniform, suggesting that (5.3) is incorrect: more variables are needed to characterize the macroscopic properties of these granular systems. Because of these conflicting results, and because (5.4) lacks a microscopic justification, it has remained controversial. Despite obvious successes [MK02, HOC07], strict tests have also shown violations of (5.4), particularly in the behaviour of geometric observables [BKLS00, PD12, PF12]. Moreover, analytical calculations which could clarify these issues are intractable due to the complexity of (5.4)4 . Here we therefore consider a simpler ensemble, motivated by the aforementioned energy argument.  5.2.2  Force Network Ensemble  We have argued above that for nearly rigid grains, there is a separation of time scales between evolution of the forces and evolution of the geometry. Hence it is not entirely surprising that stress may equilibrate, in the sense of following (5.4), while volume does not. Here we therefore consider a restricted Edwards ensemble in which the geometry, and hence V , is fixed. The restriction of the Edwards ensemble to force configurations on a strictly fixed contact network is known as the Force Network Ensemble [SVvHvS04, TV10, TvEV08, TSVvH10]. However, because contacts can open and close with an infinitesimal amount of energy, it is more physically relevant to consider an enlarged Force Network Ensemble (eFNE) in which a contact network of potential contacts is fixed, and one considers all the possible configurations of closed contacts, T . The natural 4 However, see [HC05] where the partition function is calculated at isostaticity in the related stress ensemble, obtained from (5.4) in the limit X → ∞.  78  generalization of the canonical measure (5.4) to this case is Peq (T , f ) =  1 NRC −V α:ˆ λ e ˆ σ[f ] Θ[T , f ], Z  (5.5)  where λ is a contact potential which governs the strength of fluctuations in contact number, analogous to a chemical potential in thermal statistical mechanics. Here we suppress dependence on the geometry, which is hidden in Θ and σ ˆ . Explicitly, fgc  δ  Θ[T , f ] =  c∈C g  g  r c × fgc  δ c∈C g  c Θ fN − µ1 |fTc | c  ≡ δF B [f ] δT B [f ] ΘC [T , f ],  (5.6)  where in the final step we have collected contributions from force balance, torque balance, and the Coulomb inequality. The normalization constant in (5.5) is the partition function λNRC  Z=  ˆ σ [f ] Df e−V α:ˆ Θ[T , f ],  (5.7)  T  where Df ≡ c∈RC df c , and T is a sum over all topologies of the contact network, given the potential contact network. Equivalently, this sums over all choices of real contacts. In the next Chapter, we will present a model to compute Z, from which we can determine the thermodynamic properties of (5.5). We now establish some basic properties of (5.5).  5.3 5.3.1  Properties of the enlarged Force Network Ensemble Contact potential  It is clear from the expression for Z that λ is necessary for dimensional homogeneity. By defining an elementary phase space volume, it plays a role very similar to in classical statistical mechanics. Indeed, we can write √ √ √ ˆ −V √α :ˆ σ [ λf ] λ D( λf ) e Θ[ λf ],  √ 3NRG λ Z=  (5.8)  T  where√all forces in the right-hand √ side now appear in the dimensionless combination λf . One can consider λ as the dimensional weight assigned to each real contact, per unit force. It has units of inverse force. 79  Note that equation (5.8) implies the scaling property Z(α, ˆ λ) = √  5.3.2  1 λ  3NRG  √ Z(α/ ˆ λ, 1).  (5.9)  Variational principles  The canonical probability distribution (5.5) has the Gibbs form and therefore satisfies the variational principles associated with this distribution [Par88]. For any probability distribution P, we define its configurational entropy S[P] = − log(P) P , where P denotes an average taken with respect to P. Then, of all normalized P ¯RC , we have ¯ and NRC P = N restricted to O5 , for which σ ˆ P=σ S[Peq ] ≥ S[P],  (5.10)  where the constraints on σ ˆ and NRC fix the Lagrange multipliers α ˆ and λ. The Legendre transform of S with respect to α ˆ and − log λ is the free energy functional F [P] = − log(λ) NRC + α ˆ : Vσ ˆ  P  − S[P],  (5.11)  which satisfies the variational principle F [Peq ] ≤ F [P],  (5.12)  where the variation is taken over all normalized P restricted to O. One can easily compute the standard relations F ≡ F [Peq ] = − log(Z)  (5.13)  and 1 ∂F , V ∂α ˆ ∂F = −λ , ∂λ  1 ∂S V ∂ σ ˆ ∂S − log λ = ∂ NRC  σ ˆ = NRC where S = S[Peq ], and 5  α ˆ=  (5.14)  denotes an ensemble average taken with respect to Peq .  Meaning that, for each T , force configurations lie in O(T ).  80  5.3.3  Microcanonical v.s. Canonical Ensembles  Let us establish the relationship between the canonical and microcanonical en¯ σ − σ) ¯ and 1 = N¯RC δN¯RC ,NRC (T ) into (5.7), we sembles. Inserting 1 = dσδ(ˆ find ¯  λNRC  Z=  ¯ ˆσ ¯ −V α: dσe  ¯RC N ¯  λNRC  =  δN¯RC ,NRC (T )  ¯ Df δ(ˆ σ [f ] − σ)Θ[f ]  T ¯  ¯ Smc (σ, ¯ NRC ) ˆσ ¯ −V α: dσe e ,  (5.15)  ¯RC N  where ¯RC ) = log ¯ N Smc (σ,  δN¯RC ,NRC (T )  ¯ Df δ(ˆ σ [f ] − σ)Θ[f ]  (5.16)  T  ¯ and is the entropy of the micro-canonical eFNE corresponding to fixed stress σ ¯ fixed contact number NRC . In the thermodynamic limit, (5.15) is dominated by its saddle-point contribution (ˆ σ0 , N0 ), satisfying α ˆ=  1 ∂Smc ¯ V ∂σ  (5.17) ¯RC =N0 ¯ σ 0 ,N σ=ˆ  and − log λ =  ∂Smc ¯RC ∂N  (5.18) ¯RC =N0 ¯ σ0 , N σ=ˆ  In this limit, F = − log Z = −N0 log λ + α ˆ :Vσ ˆ 0 − Smc (ˆ σ0 , N0 ) + O(log Smc ),  (5.19)  hence σ ˆ0 =  1 ∂F = σ ˆ , V ∂α ˆ  N0 = −  ∂F = NRC ∂ log λ  (5.20)  and Smc = S[Peq ],  (5.21)  establishing equivalence of the canonical and micro-canonical ensembles in the thermodynamic limit. The notable point of the above manipulations is (5.16), which shows that the entropy is simply the logarithm of the volume of configuration ¯RC . Hence, computing Z is equivalent ¯ NRC (T ) = N space on the shell σ ˆ [f ] = σ, ¯RC are varied. Before presenting a ¯ and N to knowing the geometry of O(T ) as σ model for Z, we therefore examine some simple properties of O, and hence of Smc . 81  5.3.4  Scaling  In the special case of constant contact number, the linearity of Newton’s laws allows us to extract the dependence of S on the pressure P . Consider ¯ = log SF N E (σ)  ¯ Df δ(ˆ σ [f ] − σ)Θ[f ] ,  (5.22)  the entropy of the Force Network Ensemble on a strictly fixed contact network. ¯ = kσ Making a change of coordinates in (5.22) to σ ˆ , it is easy to show that SF N E (k σ ˆ ) = ( 12 NRG d∆z − d2 ) log(k) + SF N E (ˆ σ ),  (5.23)  where ∆z = z¯ − (d + 1). If we set k = P then σ ˆ has unit pressure. This can also be used to obtain the isotropic equation of state in the FNE. We write α ˆ = α(δˆ + η ˆ),  (5.24)  where η ˆ is symmetric and Tr η ˆ = 0. Using (5.23) with k = 1/α, we find ZF N E (α, η ˆ) = α−NRG d∆z/2  ˆ  ¯ SF N E (σ) ¯ ˆ σ ¯ −(δ+η):V dσe e ,  (5.25)  hence 1 Tr σ ˆ d 1 ∂FF N E 1 = + dV ∂α α  P =  i  ∂FF N E ∂ηii  .  (5.26)  The tensor ∂FF N E /∂ η ˆ is function of η ˆ. On an isotropic geometry, we can only ˆ form tensors from δ and η ˆ, hence its most general form is ∂FF N E /∂ η ˆ = f (I1 , I2 , . . . , Id ) δˆ + g(I1 , I2 , . . . , Id ) η ˆ,  (5.27)  where In is the nth invariant of the tensor η ˆ, and f and g are unknown functions. Under reflection symmetry r → −r, we have F → F , In → (−1)n In , and η ˆ→η ˆ, which implies that odd-n invariants must appear in ratios in f and g, e.g. I1 /I3 , etc. Since I1 = Tr η ˆ = 0, in d = 2 and d = 3 this implies that f and g are functions of I2 only. When subject to isotropic forcing (ˆ η = 0), there is no physical difference between shear in +v and −v directions, and hence under η ˆ → −ˆ η , we have F → F . This implies f (I2 )|η=0 = f (0) = 0. Hence ∂FF N E /∂ η ˆ|η=0 = 0. ˆ ˆ  82  Using F = − log Z we find αP =  NRG ∆z , 2V  (5.28)  which is the isotropic equation of state in the Force Network Ensemble [TV11]. In physical systems, where contacts can open and close, the equation of state is no longer a consequence of scale invariance. In the next Chapter, we will generalize (5.28).  5.4  A note on percolation  The statistical mechanical approaches to understanding rigidity in a granular material discussed in this Chapter derive from the work of Edwards, and are not closely connected to models commonly appearing in the statistical mechanics literature. A separate thread of work has attempted to understand the rigidity of amorphous solids as a percolation phenomenon, aiming to draw on the large body of knowledge on the latter [GRH+ 90, MD95, JT95, JT96, Mou98]. More precisely, links have been drawn between granular materials and rigidity percolation, in which one considers a regular lattice of springs, diluted at random with a probability p, i.e., each spring is present with probability 1 − p [JT95, MD95, JT96]. In the thermodynamic limit, at a critical p = pc , the network of springs gains rigidity. The critical network is a fractal object, with structure at all scales [MD95]. However, in a granular material at the onset of rigidity, the contact network is not a fractal object [Mou98]. Moukarzel has suggested that the critical difference between networks of springs and granular materials is that the former contain a large number of redundant bonds [Mou98]. As discussed in the Introduction, in granular materials, the constraint that contact forces are repulsive leads to the removal of redundant contacts as the pressure decreases towards zero, i.e., as the onset of rigidity is approached. Hence the contact network of a granular material at the onset of rigidity does not resemble a critical percolation cluster. More recently, a link to percolation has been drawn using not the contact network itself, but the network formed by considering contacts carrying normal forces larger than a threshold f [OSN06]. In general, the resulting network consists of a number of disconnected clusters, separated by ’weak’ contacts carrying forces smaller than f . There is a critical threshold fc such that when f > fc , no cluster spans the system, but when f < fc , almost all of the selected contacts belong to the same cluster. Around the critical f = fc , the statistics of force network clusters display system-size scaling, exactly as in conventional percolation. Moreover, universality is found: variations in friction, polydispersity, and crystallization of the contact network do not affect the associated scaling exponents. 83  Although these scaling properties of the force networks are related to our statistical mechanical approach, we will not pursue this connection further. We leave this for future work.  84  Chapter 6  Mean-field Theory of Hyperstatic Packings  50  40  30  20  10  0  Figure 6.1: Frame10of a polydisperse packing undergoing simple shear. 60Line 0 20 30 40 50 thickness is proportional to the normal force. As discussed in the previous Chapter, in the rigid limit there is a separation of time scales between evolution of the geometry and evolution of the contact forces. This motivates a statistical mechanics in which the geometry is fixed on grain-size length scale, and only the contact forces fluctuate. In this Chapter we 85  present a mean-field theory for such a scenario, in two dimensions. Our main goal is the computation of the enlarged Force Network Ensemble partition function Z, from which we can derive the macroscopic properties of the system. The complexity of this task means that numerous mean-field approximations will be necessary; although physically motivated these are often crude. Nevertheless, we will be able to explain one ubiquitous phenomenon of stress distribution in granular packings: at mesoscopic scales, it is highly heterogeneous. A glance at the contact forces in a numerical simulation, reproduced in Figure 6.1, reveals the existence of filamentary force-carrying structures, known as force chains. These are observed in numerical simulations and experiments, but so far are without an accepted theoretical explanation [RWJM98, MB05, MSLB07]. We will see that force chains can be explained as a result of an entropic instability: there are simply more force configurations with a heterogeneous distribution of stress than without. We consider as fixed a contact network of potential contacts, called the slow geometry. The NGC potential contacts are what are measured geometrically: they correspond to grains which are touching, or grains which could touch by an infinitesimal movement. We ignore the geometrical change associated with contact opening and closing, since this is infinitesimal. In general, the set of NRC contacts which bear a non-zero force is a subset of the potential contacts1 . All lengths and areas are fixed, and the dynamics explores different topologies of the potential contact network, as well as the space of contact force configurations on each contact network. In this Chapter, we do not deal with the evolution of the slow geometry, but instead suppose that its statistical properties are known. To obtain macroscopic stress-transmission properties, we integrate over all mechanically equilibrated microscopic force configurations consistent with macroscopic constraints. As discussed in the previous Chapter, we consider a canonical distribution Peq (f ) =  1 NRC −Aα:ˆ λ e ˆ σ[f ] Θ[f ], Z  (6.1)  where volume V has been replaced by area A in two dimensions, and fgc  δ  Θ[f ] = g  c∈C g  r c × fgc  δ c∈C g  ≡ δF B [f ] δT B [f ] ΘC [f ].  c Θ fN − µ1 |fTc | c  (6.2)  1  In the work of Makse and collaborators, the contact number measured from the potential contacts is called the geometrical contact number, while that measured from the real contacts is called the mechanical contact number [SWM08].  86  The normalization constant in (6.1) is the partition function λNRC  Z=  ˆ σ [f ] Df e−Aα:ˆ Θ[f ],  (6.3)  T  where Df ≡ c∈RC df c , and T is a sum over all topologies of the contact network, given the potential contact network. Because (6.1) has the canonical form, macroscopic properties are computed from derivatives of the free energy F = − log Z. Hence our main goal is computation of Z, in the mean-field approximation. The reader who is not interested in the details of the computation can skip to section 6.5, where we present the result for F .  6.1  Gauge fields  From the perspective of statistical mechanics, force and torque balance are interactions that are implicitly specified by Newton’s laws. In their original form, they make evaluation of (6.3) extremely difficult. In Chapter 2, we showed how these constraints can be satisfied identically by the introduction of ρ = {ρt } and ψ = {ψ c }. For each grain g, the loop forces ρ and stress function ψ satisfy, in an exact discrete calculus, the relations σ ˆ g = (∇ × ρ)g = (∇ × ∇ × ψ)g .  (6.4)  Here we will use the extrinsic formulation of Chapter 2, based on the Voronoi tesselation and the Delaunay triangulation (Figure 2.1). The latter divides the packing into triangles t, whose edges are the contact vectors c . In addition to real contacts, these include virtual contacts, which connect grains which are geometrical neighbours, but not touching. In what follows, the real contacts RC and virtual contacts V C will play complementary roles. Explicitly, the contact force at c exerted on grain g can be written −  +  fgc = ρt − ρt ,  (6.5)  in terms of the loop forces on the Delaunay triangles adjacent to c. Each loop force ρt is written in terms of the stress function at neighbouring contacts as ρt = (∇ × ψ)t = −  1 At  c t  ψc,  (6.6)  c∈C t  where At is the area of the triangle t, and the ct circulate anticlockwise around the triangle. In this formulation, we treat real and virtual contacts on the same 87  footing, so to obtain physical force configurations we must explicitly require that virtual contacts carry no force, which is equivalent to requiring that loop forces are constant within the loops of the packing. Therefore, for each loop , ρ = (∇ × ψ) = −  1 A  c  ψc,  (6.7)  c∈RC  where the sum is over the real contacts around the loop. The expressions for f [ρ] and f [ψ] indicate explicitly how force and torque balance cause ρ and ψ to interact. In the new variables Newton’s laws do not appear in Θ. In two dimensions, Θ is a product of three terms. The first term, δG , chooses a gauge. The second term is δ (f c ) ,  δ[f |V C ] ≡ c∈V C  which ensures that virtual contacts carry no force. In terms of ψ, δ[f |V C ] requires that ψ varies linearly across loops. Finally, ΘC is the Coulomb constraint at the contacts, discussed above. In terms of these variables, the partition function can be written λNRC J  Z=  ˆσ Dψ e−Aα:ˆ δG δ[f |V C ]ΘC .  (6.8)  T  The Jacobian of the change of variables, J, depends only on the geometry and hence factors out of the integrals. As shown in Appendix I, up to an irrelevant 2 , the square of factor which depends only on the slow geometry, it is exactly NST the number of spanning trees in the contact network. The latter is exponential in the system size, NST = eqNRC , with q dependent on the connectivity of the contact network. We will suppose that q is a constant in the thermodynamic limit. Recalling the geometric formulation of Chapter 3, the partition function (6.8) has a simple geometric interpretation: it sums over all continuous triangulated surfaces which are flat on the interior of loops and satisfy the Coulomb inequality ˆ σ weights these surfaces by their curat each vertex. The Boltzmann factor e−Aα:ˆ vature, while the contact potential λNRC weights them by the number of vertices where curvature is allowed. Recalling that the normal to the surface (ψ, r) on loop is (1, ρ ), we see that the difference in normals across a vertex is exactly c is a constraint on the curvature (0, f c ). Hence the Coulomb inequality |fTc | ≤ µfN of the surface. In the frictionless limit, it corresponds exactly to convexity of the surface. With a finite coefficient of friction, a finite amount of twisting curvature is allowed at each vertex. 88  Sums over random surfaces have been studied by theoretical physicists since the early 1980’s, as models for quantum gravity [Dav85, ADF85]. If (6.8) did not contain the Coulomb inequality, it would be equivalent to one of these theories, with the dependence of Z on λ related to the so-called string susceptibility exponent. It was found that in the thermodynamic limit, the sum over surfaces is dominated by surfaces with infinite fractal (Hausdorff) dimension, and hence the model was declared unphysical. When subject to fixed total curvature and a convexity constraint at each vertex, it is not possible to crinkle the surface in such a way as to obtain an infinite fractal dimension. Hence our random surface theory differs essentially from these theories by microscopic constraints. In fact, we will see that convergence of Z requires the Coulomb inequality. Incidentally, this suggests that one could obtain a finite theory of quantum gravity by including a local convexity constraint, which is related to assuming that mass is always positive. To the author’s knowledge, this has not been attempted2 . In (6.8) we also have a gauge fixing constraint. Since stress corresponds to curvature, we can apply any rigid-body motion to the (ψ, r) surface without affecting the stress. We can fix a gauge by choosing three contacts in the interior of the packing and setting their ψ c = 0, which amounts to pinning the surface there. Having done so, this constraint will be omitted in the sequel. In what follows, we measure all lengths by the grain diameter D, so that any contact vector c between touching grains has unit length.  6.2  Isostaticity  The gauge and virtual contact constraints in Z restrict the integral to O, which has dimension H = 2NRC − 3N = NRG (¯ z − ziso ) ≥ 0,  (6.9)  where z¯ = 2NRC /NRG is the mean contact number and ziso = 3 in the plane. A packing is isostatic when H vanishes, i.e., when z¯ = ziso . When friction is present, the fast dynamics will, in general, cause contacts to slide, in which case the components of their contact forces are related by µfN = |fT |. This creates additional constraints to go in (6.8) and (6.9), and introduces an additional complexity which will be reported elsewhere: here we restrict ourselves 2  After the crinkling problem was discovered, some authors added terms to the Einstein-Hilbert action which depended on higher powers of the curvature. These were found to not alter the conclusion that typical surfaces had infinite Hausdorff dimension. However, a strict inequality such as the Coulomb inequality is only equivalent to adding infinitely many such terms to the action. Apparently this has not been attempted [BV95].  89  to the case of perfect friction, µ → ∞. In this limit the Coulomb condition reduces to repulsion of normal forces c fN ≥ 0.  (6.10)  To maintain continuity, we will continue to call this the Coulomb inequality. In general, the contact number of a frictional packing depends on its preparation history [SEmcG+ 02, AR07]. However, in the realistic range of pressures for macroscopic grains, the force indeterminacy per grain is small, i.e., ∆z = H /NRG = z¯ − ziso  1.  These relations neglect boundary terms, which are irrelevant to the global force indeterminacy. However, we can sharpen the concept of force indeterminacy by considering a cluster of loops i, and measuring how determinate its interior degreesof-freedom when those on its boundary are fixed. In particular, we consider a finite connected cluster of loops, i, with some prescribed values of ψ on its boundary ∂i, and ask: when do the virtual contact constraints in the cluster determine the values of all ψ c inside the cluster? The answer depends on the value of the local hyperstaticity H i = NCi − 2NVi C i = 3 + 2NRC − 3N i − z i ,  (6.11)  i where the equality is a consequence of Euler’s relation. Here, NRC and NVi C are i the number of real and virtual contacts strictly inside the cluster, z is the number of contacts around the cluster boundary, N i is the number of non-rattling grains i inside the cluster, and NCi = NRC + NVi C . The crucial difference between (6.11) and (6.9) is in the boundary terms. For a single loop, H i = 3 − z i ≤ 0, because virtual contact constraints couple the ψ field at the real contacts surrounding the loops. However, for a large enough cluster, H i → H > 0. Hence, there is a critical size ξI ∼ 1/δz at which H i = O(1), and the cluster is isostatic. By construction, ξI is a point-to-set force correlation length [MC10, MC11]. The above argument is closely related to other bulk-surface arguments related to the jamming transition, but here we use ψ, which is less constrained than f [TW99, WNW05, MC10, MC11]. Any packing can be decomposed into isostatic clusters, separated by a percolating network of real contacts, which we call force chains. Justification for this name will follow from the analysis. An example of such a decomposition, where each H i = 0, is shown in Figure 6.2. We observe that the clusters have a broad range of sizes, and with our algorithm, which forces H i = 0, they occasionally  90  45  40  35  30  25  20  15  10  5  10  15  20  25  30  35  40  45  50  Figure 6.2: Isostatic cluster decomposition on a monodisperse packing undergoing simple shear flow with z¯ = 3.4. The colours have no meaning. Multiple colours on a single triangle indicate that the corresponding loop belongs to overlapping isostatic clusters. overlap: typically, 1% of loops belong to multiple clusters. We will assume that these overlaps can be neglected. The clusters, and hence the chains, are not unique; this is because the grains which participate in force chains are equivalent to the grains which form the isostatic clusters. However, for any cluster decomposition, we have NF C = N ∆z + N∂C ,  (6.12)  a consequence of assuming H i = 0√for each cluster. 1, there will be many, large isostatic For systems large enough that N ∆z clusters. In this work we seek a simple model for such systems, in particular when they are homogeneous but perhaps anisotropic. As a mean-field approximation, 91  we suppose that each of these isostatic clusters is independent, and to ensure that the results do not depend on the cluster decomposition, we include a flat average over all decompositions. The partition function then partially factorizes as λNRC e2qNRC  Z ≈ ZM F = T  1  Zi  ZF C  i  ,  (6.13)  C  where · C denotes an average over cluster decompositions, Z i is the partition function for cluster i, and ZF C will correct for overcounting, as discussed below. Here, and in the following, the subscript M F indicates that a mean-field approximation has been employed.  6.3  Isostatic cluster partition function  In the following we denote by Li the loops in i, G∂i the grains on the cluster boundary, Gi the grains in the interior of the cluster, C ∂i the contacts on the cluster boundary, RC i the real contacts in the interior of i, V C i the virtual contacts in the interior of i, and C i = RC i ∪ V C i . For each isostatic cluster, we have a partition function Zi =  i  ˆ σ) Dψ|C i e−α:(Aˆ ΘiC δVi C .  Dψ|C ∂i  (6.14)  In making the mean-field approximation, the Coulomb constraints on force chains will need to be replaced by appropriate mean-field expressions. For definiteness, we suppose that each H i = 0. Then in each cluster, each interior ψ is exactly determined by those on the boundary. These interior integrations disappear, giving a Jacobian Dψ|C i δVi C .  Ji =  (6.15)  As shown in Appendix J J i can be estimated for a class of geometries, with the result Ji ≈  A 8N  NVi C  .  (6.16)  Hence we have the cluster partition function Zi = J i  i  ˆ σ) Dψ|C ∂i e−α:(Aˆ ΘiC .  92  (6.17)  Note that Z i is holographic: the properties of the cluster are determined by the values of the gauge field ψ on its boundary. In the geometric interpretation, Z i is a sum over continuous surfaces, bounded by a simple closed curve, whose interior is fixed by the boundary values and the sliding constraints. By writing Z as a product over Z i , instead of integrating over surfaces, we need only to integrate over curves. The difficulty is that the shape of the surface in the interior of the cluster is determined by constraints which are only implicitly specified. To each choice of contact network and isostatic cluster decomposition, there is an ensemble mean variation of the fields ψ, ρ, and f , which must respect force ¯ and torque balance. Because of linearity, it is convenient to introduce δψ = ψ − ψ, ¯ and δf = f − f¯, the deviations from ensemble means. δρ = ρ − ρ,  6.3.1  Mean-field approximations  ˆ i in terms To compute Z i , we need to express all quantities appearing in ΘiC and σ ∂i of the values of ψ on the boundary C . The stress tensor for the cluster is (Aˆ σ )i = −  c c c∈C i  c∈C i  dA σ ˆ=  =  Ac σ ˆc  f =  i c  =−  ρ  dA ∇ × ρ = − i 0 (c)  dr ρ ∂i  (6.18)  c∈C ∂i  where 0 (c) is the loop adjacent to c exterior to the cluster, and we have used discrete calculus to write this as a boundary sum. Ignoring intraloop coupling outside of the cluster, the best approximation for ρ 0 (c) is 1 (c) ρM0 F = ρ¯ 0 (c) + ¯ δψ c c , AL  (6.19)  where, for simplicity, we assume that the external loop has the mean loop area A¯L = A/NL . Hence, as a consequence of ∇ · σ ˆ = 0, the stress tensor is easily written in terms of boundary ψ. We find ¯ i− (Aˆ σ )iM F = (Aσ) c∈C ∂i  ¯ i=− with (Aσ)  c∈C ∂i  cρ ¯ 0 (c)  =−  c∈C i  93  1 δψ c ¯ AL  c f¯c .  c c  ,  (6.20)  More difficult are the Coulomb constraints ΘiC , which involve each contact force in the cluster. In principle, the interior constraints can be solved to find f c = f c [ψC ∂i ],  (6.21)  however in practice this is intractable. We will be interested in the case where clusters are large, for which the stress-geometry equation should hold. Then, in the mean-field approximation, the solution (6.21) will resemble the solution to the corresponding PDE. As we have seen, the Green’s function decays slowly as 1/|r|. A tractable mean-field approximation is to ignore this decay altogether, and replace occurrences of ρ with ρi =  1 Ai  A ρ =− ∈Li  1 Ai  c  ψc,  (6.22)  c∈C ∂i  the average loop force across the cluster, which depends only on boundary ψ. In the geometrical interpretation, we want to know how the surface normal (1, ρ ) changes across the cluster as the boundary curve fluctuates. Since (1, ρi ) is the average surface normal across the cluster, the simplest approximation is that these normals fluctuate together, i.e., for each ∈ Li , we let δρM F = λ δρi ,  (6.23)  where λ is a parameter, unrelated to the chemical potential λ. In Appendix K we show that, in the simplest mean-field approximation, λ = 1. For c on the cluster boundary, f c = ρ 0 (c) −ρ (c) , where (c) is the loop adjacent to c interior to the cluster. Using f¯c = ρ¯ 0 (c) − ρ¯ (c) , the mean-field expression for each f on the cluster boundary is then 1 c i c c ¯c fM . F = f − δρ + ¯ δψ AL  (6.24)  Equivalently, we can write c δfM F =  Gcc  c  δψ c ,  c ∈C ∂i  with Gcc = 1/Ai +  1/A¯L 0  if c = c if c = c  c and δf = f − f¯. For interior contacts c ∈ RC i , δfM F = 0. Hence each interior contact in i has its contact force set entirely by its mean value f¯c , and the fluctuations from ψ affect only the contacts on C ∂i . This has the fortunate consequence  94  that the Coulomb constraints in the interior of i only constrain the mean values f¯. It is worth considering the physical implications of this prescription for the contact forces. Since all contact forces are written in terms of loop forces, all grains in the cluster, including those on the boundary, are identically in force balance. Provided that the mean values f¯ are chosen to satisfy torque balance, then all interior grains will be identically in torque balance. If we consider a grain g on the cluster boundary, then the net torque exerted on g is τ g = − 12  c−  ×  1 δψ c− A¯L  c−  − δρi +  1 c+ 2  1 δψ c+ A¯L  ×  c+  − δρi ,  where c+ = c+ (g) is the contact anticlockwise around the cluster boundary from g, and c− the contact clockwise around the boundary from g. We see that the direct contributions from δψ c± vanish identically, leaving only the indirect contribution τ g = 12 (  c−  −  c+  ) × δρi .  Hence the grains on the cluster boundary are not identically in torque balance, but receive small contributions from fluctuations of ψ. However, since g∈∂Gi τ g = 0, the entire cluster is in torque balance. c ) ≥ 0 is equivalent to It is easily checked that the Coulomb condition Θc (fM F c c δψ c ≤ ψM ≡ A¯L f¯N +  6.3.2  c  · δρi .  (6.25)  Computing the partition function  ¯ i ), we find Extracting the mean-stress contribution Z¯ i = exp(−α ˆ : (Aσ) 1  eSi ≡ Z i /Z¯ i = J i  dψ c e A¯L  αc δψ c  ΘiC ,  (6.26)  c∈C ∂i  with αc = α ˆ : c c . This integral is nontrivial because each Coulomb constraint couples all the boundary contacts through δρi . Nevertheless, it can be computed exactly using the Fourier transform. In particular, we impose the coupling with a δ-function, i.e., we use 1= =  dρi δ(ρi − ρi [ψ]) Ai A¯L  2  Ai i i i i dρi i A¯L (ik +ν )·(ρ −ρ [ψ]) dk e , 2 (2π)  95  where we have introduced a Fourier representation of the δ-function. The vector ν i will need to be chosen in a particular way to ensure convergence of the integrals. Substituting in Z i , changing the order of integration, and performing a trivial change of variables to δψ and δρ, we find Si  e  =J  i  dk  i d(δρ  i)  (2π)2  e  c ψM  Ai i i i ¯ (ik +ν )·δρ A L  1  d(δψ c )e A¯L  (αc +γ c )δψ c  ,  −∞  c∈C ∂i  with γ c = (iki + ν i ) ·  c  (6.27)  and J i = J i (Ai /A¯L )2 . The ψ integrals now decouple into a product over the boundary contacts, each one being straightforward. Provided αc + ν i · c > 0 for each contact c, we find eSi = J  i  =J  i  1 c d(δρi ) AA¯ i (iki +ν i )·δρi 1 (αc +γ c )ψM e L e A¯L 2 c c (2π) α +γ c∈C ∂i   i 1 c (αc +γ c )f¯N  d(δρ ) eai ·δρi , dki  e αc + γ c (2π)2 ∂i  dki  (6.28)  c∈C  with  Ai ai = ¯ (iki + ν i ) + AL  (αc + γ c )  c  c∈C ∂i  zi  and J i = J i (A¯L ) . Convergence of the δρi integral requires that Re[ai ] = 0, and hence Ai (αc + ν i · c ) c . (6.29) 0 = ¯ νi + AL c∈C ∂i  The δρi integral then becomes δ(Im[ai ]), so that Ai 0 = ¯ ki + ki · AL  c c  ˆ. = ki · M  (6.30)  c∈C ∂i  ˆ Hence ˆ is non-singular and approximately equal to (Ai /A¯L + 1 z i )δ. The matrix M 2 i i 3 the only solution for k is the zero wavenumber k = 0 . We find   1 c i c c ¯ eSi = J i  (6.31) e(α +ν · )fN  , c + νi · c α ∂i c∈C  3 This is a consequence of the neglect of inter-cluster coupling. Future work should explore the effect of the latter on k.  96  ˆ ) ≈ J i /(Ai /A¯L + 1 z i )2 . Note with ν i evaluated with (6.29) and J i = J i / det(M 2 that ν i is a linear functional of the αc ’s, which represents an additive renormalization of αc from coupling with the cluster. Solving (6.29), −1    A¯L ν i = − ˆ δ+ i A  c c c∈C ∂i  ·  A¯L Ai  αc c . c∈C ∂i  Since c∈C ∂i c = 0, under an isotropic forcing αc = α we have ν i = 0. One can also show that when the cluster is isotropic, ν i = 0 even for anisotropic αc . Hence ν i = 0 only from coupling between stress anisotropy and geometrical anisotropy. In this work we consider the physically relevant case of large clusters, z i 1, for which ν i is negligible independent of α. ˆ Indeed, Ai ∼ (z i )2 , hence ν i ∼ 1/z i 1. Neglecting ν i , we have our result eSi = J i c∈C ∂i  A¯L αc f¯c e N. αc  (6.32)  We saw earlier that convergence of the integrals requires that αc + ν i · c > 0 for each contact c. Hence αc > 0 for each contact, which requires that α ˆ be a positive definite tensor. When we assemble the cluster partition functions, their product will count each force chain contact twice, except contacts at the boundary of the domain4 . To correct for over counting, we set ZF C = (A¯L )NF C c∈F C  6.4  1 −α:A c c c ¯c e ˆ σ¯ eα fN . c α  Global partition function  Hence we obtain the global mean-field partition function λNRC e2qNRC  ZM F = T  A 8N  NV C  A NL  4  NF C  ¯ ˆ σ e−α:A c∈F C  1 αc f¯c e N αc  . C  Our mean-field approach, which breaks up the domain into overlapping subdomains, the isostatic clusters, is a type of Bethe-ansatz. See [ID91].  97  6.4.1  Mean-field contact forces  We have yet to specify the value of the mean field f¯. Implicitly, we have assumed that it satisfies force and torque balance, as well as the Coulomb constraints. The ¯ can be obtained by basic variation of f¯ needed to sustain an applied stress σ ¯ in the continuum. The solving σ ˆ = ∇ × ρ = ∇ × ∇ × ψ for a constant σ ˆ =σ result is ¯ ¯ × r, ψ(r) = − 12 r × σ ¯ ¯ ρ(r) = −r × σ where we have set to zero the integration constants, which are gauge freedoms. These expressions indicate the basic secular variation of ρ and ψ needed to support ¯ c ) and the applied stress. The discrete mean fields are then obtained by ψ¯c ≡ ψ(r ¯ ), for each real contact c and loop , where r is the centroid of the loop. ρ¯ ≡ ρ(r The basic mean-field variation of contact forces is then − + ¯ f¯0c = ρ¯ − ρ¯ = −sc × σ.  (6.33)  This expression, derived from continuum considerations, gives the local variation of contact forces needed to support a constant applied stress. It accounts for the local geometry of loops, but it ignores the Coulomb constraints as well as the structure of isostatic clusters. Away from the frictionless limit, the Coulomb constraint is rarely violated by (6.33), and we do not consider it further5 . The model treats contacts in the interior of isostatic clusters differently from those on the force chains, since all the fluctuations are concentrated on the latter. As a result, there is no a priori reason to expect that force chain contacts will experience the same mean stress as those in the interior of isostatic clusters. A powerful consequence of the canonical probability form is that we can consider a family of models, with different expressions for their mean stress, and determine the model which is closest to the true (non mean-field) form using the variational property of the free energy, (5.10). A mean-field theory using this property is often called variational mean-field theory [CL00]. We therefore consider a mean-field contact force ¯ f¯c = (1 + kφc )f¯0c = (1 + kφc )sc × σ,  (6.34)  where k is a variational constant, and φc is a function such that φc = +1 on FC, φc averages to −NF C /(NRC − NF C ) on RC - FC, and varies smoothly in between, To satisfy the Coulomb constraint, one could project f¯ inside the Coulomb cone. We leave this for future work. 5  98  to preserve force and torque balance. For k > 0, the effect of this modulation is to increase the mean contact forces on FC at the expense of those on RC-FC. When k = 0, the mean contact forces depend only on the local geometry and are otherwise the same in the interior and boundary of isostatic clusters. In order to have repulsive forces everywhere, we must have −1 ≤ k ≤ (NRC − NF C )/NF C . This modulation can be considered the first in a series of corrections to (6.33) which accounts for correlations in the discrete geometry, going from the continuum scale downward. We assume that these are the most important corrections, and that the geometry is otherwise self-averaging, i.e., any extensive subset of the packing tends to its expected value in the thermodynamic limit. Applying the cluster modulation, the resulting mean stress is ¯ =σ ¯0 + σ ¯0 =σ  k A  c ¯c f0 c∈F C  −  k NF C A NRC − NF C  c ¯c f0 c∈RC−F C  NRC − NF C NF C NF C 1−k +k NRC NRC − NF C NRC  ¯ 0, =σ ¯ 0 NF C /N . where we have used self-averaging to replace, e.g., − c∈F C c f¯0c by Aσ By construction, the total stress is not affected by the cluster modulation.  6.4.2  Free-energy minimization  As discussed above, the minimization property of the free energy, (5.10) can be used to determine the optimal value of the parameter k without an ad-hoc prescription: we simply minimize F over all admissible values of k; at the local minimum, the model is closest to the true canonical form, which minimizes F over all admissible probability distributions [CL00, ID91, Par88]. The crucial point in c term comes ¯ in Z comes from all contacts, the f¯N this minimization is that while σ c = (1 + k)(f¯c ) for contacts on force only from contacts on force chains. Hence f¯N N 0 chains, and λNRC e2qNRC  ZM F = T  A 8N  −NRC  A NL  NF C  ¯0 ˆ σ e−α:A c∈F C  1 (1+k)αc (f¯c )0 N e αc  , C  where we have ignored an irrelevant constant term. For each cluster decomposition, we assume that the number of isostatic clusters and the number of contacts on force chains will diverge in the thermodynamic limit. In this case, Boltzmann weights can be replaced by average values log ζ = log αc , 99  (6.35)  and c β ≡ αc (f¯N )0 = α ˆ:  c c  c (f¯N )0 ≥ 0.  (6.36)  Assuming that cluster decompositions are self-averaging, the average over decompositions gives −NRC  A 8N  λNRC e2qNRC  ZM F = T  ¯ 0 β(1+k)NF C ˆ σ e−α:A e  A ζNL  NF C  .  The k-dependence is via a term exp(kβ). Since β ≥ 0, exp(β) ≥ 1. For any topology T , this minimizes the free energy (maximizes Z) when k takes its maximal value (NRC − NF C )/NF C . This minimization indicates that a homogeneous stress distribution k = 0 is entropically unstable: the system can increase its entropy by concentrating stress on the force chains. The only limit to this process is that all real contacts must have repulsive forces. In particular, this means that contacts at the interior of the isostatic clusters are expected to have very small forces, while those on force chains may sustain very large stresses. Note that, in principle, this need not have happened: we could have found that the free energy is minimized when k = 0, which would favour a homogeneous distribution of stress, or k < 0, which would favour large forces in the interior of isostatic clusters.  6.4.3  Saddle-point equations  The sum over all topologies is now purely combinatorial: NGC  = T  NRC =0  NGC , NRC  where NGC is the number of potential contacts. In the thermodynamic limit, the only extensive terms in the free energy come from saddle point contributions: we ¯0 ˆ σ S0 and seek the maximum value of S over N write ZM F = e−α:A 0 RC . NRC e In the thermodynamic limit, the sum is equivalent to an integral dNRC . The exponent that we need to extremize is S0 = log  NGC NRC  + NRC log  8N λe2q eβ A  + (2NRC − 3N ) log  A ζNL  , (6.37)  where we have used the identity (6.12) to eliminate NF C , and ignored sub-extensive boundary terms, irrelevant in the thermodynamic limit. We seek ∂S0 /∂NRC = 0.  100  Applying Stirling’s approximation log N ! = N log(N/e) + N → ∞, it is easy to see that log  N K  1 2  log N + O(1/N ) as  = N log N − K log K − (N − K) log(N − K) + O(log N ).  The extremal condition is A NGC − NRC 8N λe2q eβ NRC A ζNL 2q β Z − z¯ 32Aλe e −2∆z/(¯z −2) = e z¯ N ζ 2 (¯ z − 2)2 32Aλe2q eβ = g(¯ z ; Z) , N ζ2  1=  2  e−(2NRC −3N )/NL  (6.38)  where Z = 2NGC /NRG is the geometrical contact number, and we have grouped the z-dependence into g(¯ z ; Z). This equation determines NRC , or equivalently, z¯, as a function of λ and α. ˆ To leading order in ∆z, its solution is z¯ − ziso =  1−  13 ˜ 13 + h  Z − ziso −  3 ˜ h  ,  where 8π e2q eβ λ 32A e2q eβ λ ˜ α) = , h(λ, ˆ = N ζ2 φ ζ2 where φ is area fraction. Here we have used A = Ag N/φ with Ag = π/4 the average area per grain. ¯ S0 ˆ σ The final expression for the partition function is ZM F = e−α:A e , with S0 evaluated using the saddle point value of NRC . Comparing with (5.11) and (5.13), we thus identify the entropy S as S = S0 − NRC log λ. Since we have neglected some constants above, our expression for S cannot be used to obtain numerical values for the entropy; however, it suffices to extract the parametric dependence of S, as needed to obtain macroscopic properties.  101  6.5  Free energy  The dimensionless free energy is F (λ, α) ˆ = − log(ZM F ) ¯ − S0 =α ˆ : Aσ ¯ − NGC log NGC + NRC log NRC + (NGC − NRC ) log(NGC − NRC ) =α ˆ : Aσ − NRC log(8N λeβ e2q /A) − NRG ∆z log(A/(ζNL )).  (6.39)  ¯ q, A, NRG , and NGC are constants; NRC = NRC (λ, α), In this expression, σ, ˆ c as specified implicitly by (6.38), ζ = ζ(α) ˆ = exp( log α ), and β = β(α) ˆ = c ) ; and z ¯ = 2NRC /N , NL = NRC − N + 1, αc = c · α ˆ · c . From F αc (f¯N 0 we can determine the equation of state, the contact potential equation, and the fluctuations and correlations in σ ˆ and z¯. To understand the physical meaning of (6.39), we recall the identification of entropy S as the logarithm of the number of states. The expectation value of a macroscopic variable Q = Q(ˆ σ , z¯) is Q =  d¯ z  ˆσ dˆ σ λNRC e−Aα:ˆ Ω(ˆ σ , z¯)Q(ˆ σ , z¯),  where Ω = eS is the density of states. We can write Ω as Ω = eS = eS1 eS2 eS3 , where S1 = NGC log NGC − NRC log NRC − (NGC − NRC ) log(NGC − NRC ), S2 = NRC log(8N eβ e2q /A), and S3 = N ∆z log(A/(ζNL )). The prefactors in these terms reflect their origin. The combinatorial terms in S1 correspond to the entropy of the contact network and count the number of ways to choose the real contacts from the potential contacts. The terms in S2 are contributed from each real contact: NRC log(8N e2q /A) arise from Jacobians and play a minor role. However, the term NRC β is the entropy of the Coulomb inequality: in appropriate units, β is the difference between the average normal force at a contact and the minimum c = 0). In a model with a finite coefficient of friction, allowable normal force (fN c − 1 |f¯c |) , with the same interpretation. this would be generalized to β = αc (f¯N µ T Hence β measures the size of phase space without regard to force or torque balance. Finally, the terms in S3 are contributed only from contacts on force chains. These correspond to the entropy of force configurations which satisfy Newton’s laws.  102  6.5.1  Equation of state  ¯ Using (5.13), The mean-field equation of state is σ ˆ = σ. ¯= σ  1 ∂F NRC ∂β N ∆z ∂ log ζ ¯− =σ + , A ∂α ˆ A ∂α ˆ A ∂α ˆ  since the terms from ∂NRC /∂ α ˆ vanish by virtue of ∂S0 /∂NRC = 0. This simplifies to 2∆z z¯  c c c  ·α ˆ·  c c ¯c fN  =  c  ,  (6.40)  where we have dropped the subscript on f¯N . We will resolve this into scalar equations below. Contracting along α ˆ this implies 2∆z =α ˆ: z¯  c c ¯c fN  = β.  (6.41)  Substituting this into (6.38) and solving in the limit ∆z → 0, the leading-order contact potential equation is z¯ − ziso =  1−  11 11 + h  Z − ziso −  3 h  ,  (6.42)  where h(λ, α) ˆ =  8πe2q λ , φ ζ2  where φ is area fraction. The contact potential equation implies that at isostaticity, the geometrical contact number is Ziso = ziso +  3φ ζ 2 . 8πe2q λ  (6.43)  This implies that even when z¯ = ziso , a range of geometrical contact numbers is possible, owing to fluctuations. It is convenient to put the equation of state in scalar form. We will take axes so that the stress tensor has the form   P + T 0  · Rω ¯ = R−ω ·  σ (6.44) 0 P −T   P + T cos(2ω) −T sin(2ω) , (6.45) = −T sin(2ω) P − T cos(2ω) 103  where Rω is an anticlockwise rotation matrix with angle ω. In this representation, P is the pressure, and T is the shear part of the stress tensor. By allowing T to take either sign, we can require 0 ≤ ω < π/2. The vectors sc connect the centres of neighbouring loops, and are, on average, c perpendicular to c . We let sc /|sc | = − c ×Rγ , so that γ c represents the deviation from orthogonality. When the Delaunay triangulation and Voronoi tesselation are divided into areas for each contact, we find Ac = 21 sc × c = 12 |sc | cos(γ c ). Writing c /|sc | = − c × Rγ c × σ ¯· c= contact vectors as c = (cos(θc ), sin(θc )), we have f¯N c c c cos(γ )P + cos(2θ + 2ω − γ )T . For simplicity, here we will consider only macroscopically isotropic geometries for which the distribution of contact angles is P(θc ) = 1/(2π); the anisotropic case is discussed in Appendix L. Furthermore, we will let γ c = 0; in Appendix L, we show that considering a distribution of γ c does not change the results in the isotropic case. Finally, as shown in Appendix L, Ac is uncorrelated with θc in simple shear flow, which simplifies calculations. To express the equation of state in terms of scalar equations, it is convenient c = M ˆ into M11 + M22 = f¯c , M11 − M22 = to decompose the matrix c c f¯N N c , and M + M c ¯c cos(2θc )f¯N 12 21 = sin(2θ )fN . Then with the above prescriptions 2A f¯N = P NRC AT cos(2θc )f¯N = cos(2ω) NRC AT sin(2θc )f¯N = − sin(2ω), NRC which gives the right-hand-side of the tensorial equation of state. To compute the left-hand-side, we likewise take coordinates so that   α+η 0  · Rν , α ˆ = R−ν ·  (6.46) 0 α−η where Rν is an anticlockwise rotation matrix with angle ν. We saw above that α ˆ must be positive definite; this requires α > |η| ≥ 0. By allowing η to take either sign, we can require 0 ≤ ν < π/2. In these coordinates, we have αc = c · α ˆ· c= c α + η cos(2θ + 2ν). Writing η/α = sin(2ξ)  104  (6.47)  with |ξ| ≤ π/4, we find the isotropic-geometry equations of state 1 2AP = α cos(2ξ) ∆zNRG tan(ξ) cos(2ν) AT − cos(2ω) = α cos(2ξ) ∆zNRG tan(ξ) sin(2ν) AT sin(2ω). − = α cos(2ξ) ∆zNRG These can be rewritten more simply as NRG ∆z 1 , 2A α cos(2ξ) T /P = −2 tan(ξ), P =  (6.48) (6.49)  ω = ν.  (6.50)  Not surprisingly, we find that the principal axes of the stress must coincide with those of α. ˆ The stress ratio T /P is found to be a function of the analogous ratio η/α in α. ˆ We also note that the isotropic equation of state generalizes the scaling result of the previous Chapter: when stress fluctuations are isotropic, ξ = 0 and we recover the earlier result. However, we notice that contact fluctuations do not affect these equations. Because we have no notion of a granular thermometer with which to measure α, ˆ the equations of state alone cannot be tested empirically. Just as in thermal statistical mechanics, we need additional relations which relate the temperature to the fluctuations. Recall that for an ideal gas, the equation of state relates the pressure to the temperature. This is only useful if we know how to measure T ; we do so through the additional relation mv 2 = kB T . It is by relating velocity fluctuations to temperature that we understand the meaning of temperature, and thus understand the meaning of the equation of state. Likewise, we can derive equations from F which relate α ˆ to the stress fluctuations, and later eliminate α ˆ to relate the mean stress to its fluctuations.  6.5.2  Fluctuations  The free energy also determines the fluctuations in stress and contact number, and their correlations. The stress fluctuations are determined by the 4-tensor 1 ∂2F A2 ∂ α∂ ˆ α ˆ c c c c NF C = A2 ( c·α ˆ · c )2  σ ˆσ ˆ − σ ˆ σ ˆ =−  105  .  (6.51)  By symmetry, this reduces to 5 scalar equations. Rather than present the full system of 5 stress covariances, we concentrate on the variance and covariance of ¯ = 2αP + 2ηT cos(2ω − 2ν), one can verify directly from P and T . Using α ˆ :σ (6.3) that −1 ∂ 2 F , 4A2 ∂α2 −1 ∂ 2 F (δT )2 ≡ (T − T )2 = , 4A2 ∂η 2 −1 ∂ 2 F C ≡ (P − P )(T − T ) = . 4A2 ∂α∂η (δP )2 ≡ (P − P )2 =  After some algebra we find 1 NRG 18g(¯ z) ∆z − 2 cos2 (ξ) α2 cos3 (2ξ) 4A2 z¯ g (¯ z) 1 NRG 2 sin2 (ξ) 18g(¯ z ) cos(2ξ) (δT )2 = 2 ∆z 1 − − 2 2 3 2 α cos (2ξ) 4A z¯ g (¯ z ) cos2 (ξ) tan (2ξ) 1 NRG 18g(¯ z) C=− 2 ∆z sin(2ξ) + 2 tan(ξ) cos(2ξ) 3 2 α cos (2ξ) 4A z¯ g (¯ z)  (δP )2 =  (6.52) (6.53) (6.54)  In each of these covariances, the term with ∆z comes directly from stress fluctuations through α, ˆ while the term with g(¯ z ) comes indirectly through the contact fluctuations. As ∆z → 0, the former vanish, because isostaticity implies that there is a single force configuration in mechanical equilibrium. However, stress fluctuations can be induced by fluctuations in the contact network. Similarly, the fluctuations in z¯ are determined by F . We have (δ¯ z )2 ≡ z¯2 − z¯  2  4 ∂2F 2 ∂(log λ)2 NRG 4 ∂NRC = 2 NRG ∂(log λ) 2 g(¯ z) =− . NRG g (¯ z) =−  (6.55)  A trivial consequence of this expression is that NRG δ¯ z depends only on z¯ and Z; this has been observed in numerical simulations [MLRJ+ 08]. These relations can be used to test the theory, as announced above. We have prepared 1000 homogeneous and isotropic packings at large friction µ = 10, over a 106  range of area fractions (0.798 < φ < 0.870) and system sizes (1000 < N < 6800). The packings are formed by randomly placing grains in a box at initial area fraction φ0 = 0.5, and then quasistatically, homogeneously, and isotropically compressing the gas to the desired final area fraction. Details of the numerics are discussed in Appendix D. Basic results of the simulations are shown in Figure 6.3. We reproduce the √ expected scalings P ∼ (∆z)2 and ∆z ∼ φ − φ0 obtained in earlier simulations [OSLN03, vH10]6 . ×106  3000 1.2  2500 1  2000 (∆ z)2  P1/2  0.8  1500 1000  0.4  500 0 3  0.6  0.2 0  3.5  4  4.50.8  z  0.82  0.84  0.86  φ  Figure 6.3: Results of numerical simulations showing the scaling of and (∆z)2 with φ.  √  P with z  To test the equation of state, we solve both the isotropic equation of state and the pressure fluctuation equation for α, and check that we get the same result. From numerical simulations, it is straightforward to compute everything needed to do so, except the geometrical coordination number Z. A simple physical hypothesis is that new contacts are formed by contact of a rattler with a forcebearing grain, and hence NGC = NRC +x0 N , where x0 is the proportion of rattlers. x0 is shown as a function of P in Figure 6.4. To compute fluctuations numerically, recall that the canonical ensemble we consider has a single value of angoricity and contact potential across the system, which leads to fluctuations in the stress and contact number in subsets of the system. Any such subsystem must be sufficiently large for the canonical distribution to hold. Numerical simulations of frictionless packings indicate that the angoricity 6  Note that, when comparing P to φ or ∆z, one must take into account the contact force law. Our code uses linear springs, which implies P ∼ φ − φ0 . Hertzian contacts would have P ∼ (φ − φ0 )3/2 . See [AR07, vH10].  107  is well-defined down to subsystems of three grains [HOC07]. In our simulations we compute fluctuations over subsystems of approximately 25 grains. 0.2  x0  0.15  0.1  0.05  0 0  2  4  6  8  P  10 6  x 10  Figure 6.4: Proportion of rattlers as a function of P , as obtained from numerical simulations.  7  7  x 10  6  6  5  5  4  4 1/α  1/α  x 10  3  3  2  2  1  1  0 0  0.2  0.4  0.6 ∆z  0.8  1  0 0  0.2  0.4  0.6 ∆z  0.8  1  Figure 6.5: Prediction of 1/α from equation of state (black, ×) and pressure fluctuations (red, 0). Each mark corresponds to a different simulation, over a range of N and φ. The left assumes Z = z¯ + 2x0 /(1 − x0 ), while the right assumes Z = 6. The resulting α comparison is shown in Figure 6.5. On the left we use Z = z¯ + 2x0 /(1 − x0 ), as suggested above, while on the right we consider an extreme 108  case Z = 6, meaning that all virtual contacts in the Delaunay triangulation could be closed by an infinitesimal movement. In both cases, the prediction of 1/α from the pressure exceeds the prediction from the fluctuations. As shown in Figure 6.6, the ratio of these predictions has a nontrivial scaling with ∆z. Disagreement between these predictions for α is a failure of the theory. The dynamics of the contact network appears to play a minor role. 14 12  αF/αE  10 8 6 4 2 0  0.2  0.4  0.6 ∆z  0.8  1  Figure 6.6: Ratio of α as predicted from fluctuations (αF ) to α as predicted from the equation of state (αE ), using Z = 6.  6.6  Conclusion  In this Chapter we have presented a statistical-mechanical theory of the enlarged Force Network Ensemble. Using tools developed in Chapter 2, we have been able to sum a large number of force configurations which are exactly in force balance, and macroscopically in torque balance. The weights given to these configurations depend on the total stress supported by the packing as well as the number of force-bearing contacts. Our first step in computation of the partition function is to divide the packing into isostatic clusters. These are clusters of loops such that when the stress function ψ is specified on the boundary of the cluster, its value everywhere in the interior in uniquely determined by constraints of mechanical equilibrium. The isostatic cluster construction is closely related to similar constructions in the literature, which count constraints in the interior and boundary of a domain [TW99, WNW05, MC10, MC11]. However, unlike the Tkachenko-Witten 109  and Wyart et al constructions, we consider hyperstatic, rather than hypostatic (¯ z < ziso ) packings. Mailman and Chakraborty consider regions in which the boundary contact forces are fixed, and the interior contact forces allowed to vary. The latter are more strongly correlated than the stress function ψ, since they are more strongly constrained. By working with the stress function ψ rather than contact forces, we are able to create clusters which are exactly isostatic, whose size can be precisely enumerated using Euler’s formula, i.e., by setting the local hyperstaticity H i = 0. In contrast, the Mailman-Chakraborty construction yields only scaling arguments for the size of isostatic regions. An example of an isostatic cluster decomposition is shown in Figure 6.2. We call the boundaries of isostatic clusters force chains. An isostatic cluster decomposition indicates the spatial extent to which forces are correlated by mechanical equilibrium. Indeed, suppose we have a configuration in force and torque balance, and now perturb ψ c at one contact which is on a force chain. Then, by construction, a new configuration in force and torque balance can be found by changing the field ψ in the two isostatic clusters adjacent to c, and nowhere else 7 . This explicitly indicates that the point-to-set correlation length for the field ψ, ξI , is the linear size of isostatic clusters [MC10, MC11]. Assuming that clusters are circular leads to a simple estimate of ξI . Indeed, the area and circumference of a typical cluster will then be πξI2 and 2πξI . If there are NI isostatic clusters, we find NI πξI2 = A and NI 2πξI = NF C , which implies 2A 2A 1 ∼ , NF C NRG ∆z √ where we have used NF C = NRG ∆z + O( NRG ) and assumed that the packing is large. Hence the size of typical isostatic clusters diverges as isostaticity is approached, though one must treat boundary terms carefully [WNW05, MC10, MC11]. Because ξI is the correlation length of forces induced by mechanical equilibrium, it is natural to suppose that the isostatic clusters are statistically independent. This is our first approximation in the mean-field theory; it reduces computation of the partition function to computation of the partition function of a single isostatic cluster. Within a cluster, we need only integrate the field ψ around its boundary: mechanical equilibrium then uniquely fixes the value of ψ in the cluster interior. At this point, the computation is still intractable; we therefore approximate the ξI =  7 This precise identification of the minimal changes between configurations in mechanical equilibrium on disordered packings is novel; however, see [TvEV08, TV11].  110  interior using knowledge of the stress-geometry equation, developed in Chapter 3. Beyond this point, the computation is straightforward. The mean-field approximation requires that we specify how ensemble-mean contact forces f¯c depend on the local geometry. Using a variational property of the canonical distribution (6.1), we have shown that it is entropically favourable for mean contact forces to be very large on the force chains and small in the interior of isostatic clusters. This provides a plausible explanation for the ubiquitous observation of filamentary networks of large forces, surrounded by regions of small forces, in experiment and simulation [RWJM98, MB05, MSLB07]. However, there is an important caveat: the isostatic cluster decomposition is not unique. This implies that in real systems, the observed forces will be a superposition of the forces obtained from different isostatic cluster decompositions. The present theory says nothing about how the weights of this superposition are chosen. The phenomenology of isostatic clusters is reminiscent of the decomposition into strong and weak force networks considered by Radjai and collaborators [RWJM98, SR05]. These authors have shown that strong contacts, that is, contacts bearing normal forces above the mean, are structurally different from the weak contacts. The strong contacts form a filamentary network which spans the system and carries nearly all of the shear stress. In contrast, the weak contacts contribute only a pressure, and thus play the role of an interstitial fluid. In the isostatic cluster framework, because of entropic instability, the weak contacts are precisely those in the interior of isostatic clusters. We thus attribute the existence of these distinct force networks as a consequence of entropic instability. Some of the predictions of the theory have been tested with numerical simulations. We find quantitative disagreement, indicating that the approximations employed need to be scrutinized. We leave this for future work.  111  Chapter 7  Conclusion In this thesis we have considered the quasistatic flow of two-dimensional granular materials. In Chapter 2, we developed tools of discrete calculus, which allow discrete equations to be written in a form which converges to their continuum counterpart when averaged over a sufficiently large volume. The main application of these tools is the motivation for and natural definition of stress and torque potentials, which can be used to satisfy the equations of mechanical equilibrium identically. In this way, we have recovered earlier results of Satake and Ball and Blumenfeld, in particular the definition of loop forces which satisfy force balance identically [BB02, Sat92, Sat93]. However, we have been able to extend these results by introducing a discrete Airy stress function, which plays a crucial role in later chapters. Moreover, we have extended these constructions to 3D, as presented in Appendix B. In Chapters 3 and 4, we have derived, directly from first-principles, the missing stress-geometry equation, as well as a mean-field phase diagram for quasistatic flow. These results are new, and are summarized in section 4.5. If the fabric of a packing is known, then the equation allows the complete determination of stresses, given stress boundary conditions. To illustrate this, we have derived the Green’s function for a point forcing at the edge of a semi-infinite plane. We have not attempted quantitative comparison with experiment, which is a definite goal of future research. One issue which arises in comparison with real data is that real geometries are finite, which makes solving the equation much more difficult. Numerical solutions to the stress-geometry PDE may be necessary. As discussed in Chapter 3, the stress-geometry equation is equivalent to the equation satisfied by σ ˆ in a particular flavour of anisotropic elasticity. Several researchers have claimed that stress distribution in granular materials should be described by anisotropic elasticity, but these arguments have all been heuristic [dG99, GG02, GG04]. By deriving the equation from first principles, we have been able to put these investigations on a sounder footing. Moreover, the derivation makes a precise and testable connection between the anisotropy tensor and fabric anisotropy [Oda72a, Oda72b, OKNN82]. In a comparison of isotropic elasticity with experimental data of stresses below a sandpile, Serero and collaborators found qualitative, but not quantitative, agreement; they suggested that anisotropic elas112  ticity would be a natural refinement, but could not test this hypothesis because the large number of fitting parameters in the anisotropy tensor would need to be constrained [SRC+ 01]. By showing the relation of the latter to fabric anisotropy, we remove this difficulty. In Chapter 4, we extended this discussion to include the deformation of granular materials, which shed some light on the development of fabric, and the jamming and unjamming transitions between solid and gaseous phases. Owing to incomplete understanding of contact opening and contact closing processes, the deformation equations are incomplete: we can solve them only by ignoring the relation between fabric and deformation. The results qualitatively agree with experiment, but a closer look into contact opening and closing is necessary to complete the model. The mean-field model thus goes a long way towards resolving the fundamental problem we posed in the Introduction: we have derived an equation for stress transmission in a granular material, and we have established conditions under which the material is predicted to fail. Future work is needed to close the gaps in the theory, and to establish the limits of the present approach with numerical and experimental tests. In Chapter 5, we shifted perspective from mean-field models which ignore fluctuations to statistical mechanics proper, where these fluctuations can be taken into account. Motivated by energy arguments, we considered a statistical mechanics of contact forces on a fixed geometry. We generalize earlier constructions by allowing contacts to open and close, defining the enlarged Force Network Ensemble. Before presenting a detailed mean-field model for the partition function of this ensemble, we derive some simple consequences of the assumed probability distribution. Finally, in Chapter 6, we use tools from Chapter 2 as well as mean-field approximations to construct the partition function for the enlarged Force Network Ensemble. In doing so, we found that the existence of force chains can be explained by an entropic instability. Although quantitative agreement with numerical simulations was not obtained, this theory provides a plausible scenario by which to construct a statistical mechanics for granular materials. Each approximation in the theory can be scrutinized with numerical simulations, a goal of future work. To overcome the difficulties of discrete models and naturally explore the continuum limit, it is customary in theoretical physics to pose statistical mechanics problems as computation of a functional integral, where one integrates over entire functions rather than a discrete collection of variables [CL00, ID91, Par88, HC05]. Because the results of Chapter 2 and 3 transform the discrete equations of force and torque balance into continuum equations, it is possible to formulate such a functional integral for granular materials. This is a promising direction for future work. 113  Bibliography [ABG+ 01] A. P. F. Atman, P. Brunet, J. Geng, G. Reydellet, P. Claudin, R. P. Behringer, and E. Cl´ement. From the stress response function (back) to the sand pile “dip”. The European Physical Journal E: Soft Matter and Biological Physics, 17(1):93–100, 2005-05-01. [ABOS12] SS Ashwin, Jerzy Blawzdziewicz, Corey S O’Hern, and Mark D Shattuck. Calculations of the structure of basin volumes for mechanically stable packings. Physical Review E, 85(6):061307, 2012. [ADF85] J. Ambjørn, B. Durhuus, and J. Fr¨ohlich. Diseases of triangulated random surface models, and possible cures. Nuclear Physics B, 257:433–449 [AKV+ 04] E. Almaas, B. Kov´acs, T. Vicsek, Z. N. Oltvai, and A.-L. Barab´asi. Global organization of metabolic fluxes in the bacterium Escherichia coli. Nature, 427:839–843, 2004. [Ale98] S. Alexander. Amorphous solids: their structure, lattice dynamics and elasticity. Physics Reports, 296(2-4):65–236, 1998. [Amo06] G. Amontons. Histoire de l’Academie Royale des Sciences avec les Memoires de Mathematique et de Physique, 1699-1708. Chez Gerald Kuyper, Amsterdam, 1706. [AR07] Ivana Agnolin and Jean-No¨el Roux. Internal states of model isotropic granular packings. i. assembling process, geometry, and contact networks. Phys. Rev. E, 76:061302, Dec 2007. [AR10] Emilien Az´ema and Farhang Radja¨ı. Stress-strain behavior and geometrical properties of packings of elongated particles. Physical Review E, 81(5):051304, 2010. [BB02] R. C. Ball and R. Blumenfeld. Stress field in granular systems: Loop forces and potential formulation. Phys. Rev. Lett., 88:115505, 2002.  114  [BCLO19] J. P. Bouchaud, P. Claudin, D. Levine, and M. Otto. Force chain splitting in granular materials: A mechanism for large-scale pseudoelastic behaviour. The European Physical Journal E: Soft Matter and Biological Physics, 4(4):451–457, 2001-04-19. [BDD06] E. Bertin, O. Dauchot, and M. Droz. Definition and relevance of nonequilibrium intensive thermodynamic parameters. Physical review letters, 96(12):120601 [BKLS00] A. Barrat, J. Kurchan, V. Loreto, and M. Sellitto. Edwards’ measures for powders and glasses. Physical review letters, 85(24):5034– 5037, 2000. [Blu04] Raphael Blumenfeld. Stresses in isostatic granular systems and emergence of force chains. Physical Review Letters, 93(10):108301–, 08 2004. [BM76] J.A. Bondy and U.S.R. Murty. Graph theory with applications, volume 290. Macmillan London, 1976. [Bol79] B. Bollobas. Graph Theory. Springer-Verlag, New York, U.S.A., 1979. [Bou02] J. Bouchaud. Granular media: some ideas from statistical physics. eprint arXiv:cond-mat/0211196, November 2002. [BS01] Alexander Bobenko and Boris Springborn. A discrete Laplace– Beltrami operator for simplicial surfaces. Discrete & Computational Geometry, 38(4):740–756, 2007-12-01. [BV95] W. Bock and J.C. Vink. Failure of the Regge approach in two dimensional quantum gravity. Nuclear Physics B, 438(1):320–344, 1995. [CBD12] C Coulais, RP Behringer, and O Dauchot. Dynamics of the contacts reveals widom lines for jamming. EPL (Europhysics Letters), 100(4):44005, 2012. [CCL97] F. Calvetti, G. Combe, and J. Lanier. Experimental micromechanical analysis of a 2d granular material: relation between structure evolution and loading path. Mechanics of Cohesive-frictional Materials, 2(2):121–163, 1997. [CCPZ12] Patrick Charbonneau, Eric I. Corwin, Giorgio Parisi, and Francesco Zamponi. Universal microstructure and mechanical stability of 115  jammed packings. 2012.  Physical Review Letters, 109(20):205501–, 11  [Cha10] B. Chakraborty. Statistical ensemble approach to stress transmission in granular packings. Soft Matter, 6:2884–2893, 2010. [CL00] P. M. Chaikin and T. C. Lubensky. Principles of Condensed Matter Physics. Cambridge University Press, Cambridge, U.K., October 2000. [Cou21] C. Coulomb. Theorie des machines simples, en ayant egard au frottement de leurs parties et A la roideur des cordages. Bachelier, Paris, 1821. [CR00] Ga¨el Combe and Jean-No¨el Roux. Strain versus stress in a model granular material: A devil’s staircase. Physical Review Letters, 85(17):3628–3631, 10 2000. [CRST12] M.P. Ciamarra, P. Richard, M. Schr¨oter, and B.P. Tighe. Statistical mechanics for static granular media: open questions. Soft Matter, 2012. [CWBC98] M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, and P. Claudin. Jamming, force chains, and fragile matter. Phys. Rev. Lett., 81(9):1841– 1844, Aug 1998. [Dav85] F. David. A model of random surfaces with non-trivial critical behaviour. Nuclear Physics B, 257:543–576 [dG99] P. G. de Gennes. Granular matter: a tentative view. Reviews of Modern Physics, 71(2):S374–S382, 03 1999. [DHLM05] M. Desbrun, A. N. Hirani, M. Leok, and J. E. Marsden. Discrete Exterior Calculus. 2005. [DM11] Eric DeGiuli and Jim McElwaine. Laws of granular solids: Geometry and topology. Phys. Rev. E, 84:041310, Oct 2011. [dV80] L. da Vinci. Static measurements of sliding and rolling friction. In Codex Arundel, volume 40. British Library, 1480. [Edw05] S.F. Edwards. The full canonical ensemble of a granular system. Physica A, 353:114–118, 2005.  116  [EG99] S. F. Edwards and D. V. Grinev. Statistical mechanics of stress transmission in disordered granular arrays. Phys. Rev. Lett., 82:5397, 1999. [EG01] S. F. Edwards and D. V. Grinev. Transmission of stress in granular materials as a problem of statistical mechanics. Physica A: Statistical Mechanics and its Applications, 302(1–4):162–186, 12 2001. [EO89] S.F. Edwards and R.B.S. Oakeshott. Theory of powders. Physica A, 157:1080–1090, 1989. [FP08] Yo¨el Forterre and Olivier Pouliquen. Flows of dense granular media. Annual Review of Fluid Mechanics, 40(1):1–24, 2008. [GBO06] Guo-Jie Gao, Jerzy Blawzdziewicz, and Corey S O’Hern. Frequency distribution of mechanically stable disk packings. Physical Review E, 74(6):061304, 2006. [GBOS09] Guo-Jie Gao, Jerzy Blawzdziewicz, Corey S O’Hern, and Mark Shattuck. Experimental demonstration of nonuniform frequency distributions of granular packings. Physical Review E, 80(6):061304, 2009. [GF82] B.J. Gellatly and J.L. Finney. Characterisation of models of multicomponent amorphous metals: the radical alternative to the Voronoi polyhedron. J. Non-Crys. Sol., 50:313–329, 1982. [GG02] C. Goldenberg and I. Goldhirsch. Force chains, microelasticity, and macroelasticity. Physical review letters, 89(8):84302, 2002. [GG04] C. Goldenberg and I. Goldhirsch. Small and large scale granular statics. Granular Matter, 6(2):87–96, 2004. [GHL+ 01] Junfei Geng, D. Howell, E. Longhi, R. P. Behringer, G. Reydellet, L. Vanel, E. Cl´ement, and S. Luding. Footprints in sand: The response of a granular material to local perturbations. Physical Review Letters, 87(3):035506–, 07 2001. [GLBH01] Junfei Geng, Emily Longhi, R. P. Behringer, and D. W. Howell. Memory in two-dimensional heap experiments. Physical Review E, 64(6):060301–, 11 2001. [GRCB03] Junfei Geng, G. Reydellet, E. Cl´ement, and R. P. Behringer. Green’s function measurements of force transmission in 2d granular materials. Physica D: Nonlinear Phenomena, 182(3–4):274–303, 8 2003. 117  [GRH+ 90] Etienne Guyon, St´ephane Roux, Alex Hansen, D Bideau, J-P Troadec, and H Crapo. Non-local and non-linear problems in the mechanics of disordered systems: application to granular media and rigidity problems. Reports on Progress in Physics, 53(4):373, 1990. [HC05] S. Henkes and B. Chakraborty. Jamming as a critical phenomenon: a field theory of zero-temperature grain packings. Physical review letters, 95(19):198002, 2005. [HC09] S. Henkes and B Chakraborty. Statistical mechanics framework for static granular matter. Phys. Rev. E, 79:061301, 2009. [HCV52] David Hilbert and Stephan Cohn-Vossen. Geometry and the imagination. Chelsea, New York, NY, 1952. [Hen09] S. Henkes. A statistical mechanics framework for static granular matter. PhD thesis, Brandeis U., 2009. [Hil63] R. Hill. Elastic properties of reinforced solids: Some theoretical principles. Journal of the Mechanics and Physics of Solids, 11(5):357– 372, 9 1963. [Hil98] Rodney Hill. The mathematical theory of plasticity, volume 11. Oxford university press, 1998. [HOC07] Silke Henkes, Corey S. O’Hern, and Bulbul Chakraborty. Entropy and temperature of a static granular assembly: An Ab Initio approach. Phys. Rev. Lett., 99:038002, Jul 2007. [HS01] T. Hatano and S. Sasa. Steady-state thermodynamics of langevin systems. Physical review letters, 86(16):3463, 2001. [ID91] C. Itzykson and J.-M. Drouffe. Statistical Field Theory. Cambridge University Press, Cambridge, U.K., 1991. [JNB96] H.M. Jaeger, S.R. Nagel, and R.P. Behringer. Granular solids, liquids, and gases. Rev. Mod. Phys., 68:1259–73, 1996. [JT95] Donald J Jacobs and Michael F Thorpe. Generic rigidity percolation: the pebble game. Physical review letters, 75(22):4051–4054, 1995. [JT96] D. J. Jacobs and M. F. Thorpe. Generic rigidity percolation in two dimensions. Phys. Rev. E, 53(4):3682–3693, Apr 1996.  118  [KR96] NP Kruyt and L. Rothenburg. Micromechanical definition of the strain tensor for granular materials. Journal of applied mechanics, 63(3):706–711, 1996. [Kru03] N.P. Kruyt. Statics and kinematics of discrete cosserat-type granular materials. International Journal of Solids and Structures, 40(3):511 – 534, 2003. [Kru10] Niels P. Kruyt. Micromechanical study of plasticity of granular materials. Comptes Rendus M´ecanique, 338(10–11):596–603, 2010/11// 2010. [Kuh99] M.R. Kuhn. Structured deformation in granular materials. Mechanics of materials, 31(6):407–429, 1999. [Kur00] J. Kurchan. Emergence of macroscopic temperatures in systems that are not thermodynamical microscopically: towards a thermodynamical description of slow granular rheology. Journal of Physics: Condensed Matter, 12(29):6611 0953–8984, 2000. [Kur05] J. Kurchan. In and out of equilibrium. Nature, 433(7023):222–225 [LDBB08] F. Lechenault, O. Dauchot, G. Biroli, and J. P. Bouchaud. Critical scaling and heterogeneous superdiffusion across the jamming/rigidity transition of a granular glass. EPL (Europhysics Letters), 83(4):46003, 2008. [LL86] L.D. Landau and E.M. Lifshitz. Theory of Elasticity, volume 7 of Course of Theoretical Physics. Pergamon Press, Oxford, U.K., 1986. [LN98] A.J. Liu and S.R. Nagel. Jamming is not just cool any more. Nature, 396, 1998. [LN10] Andrea J. Liu and Sidney R. Nagel. The jamming transition and the marginally jammed solid. Annual Review of Condensed Matter Physics, 1(1):347–369, 2010. [Mac37] S. MacLane. A combinatorial condition for planar graphs. Fundamenta Mathematicae, 28:22–32, 1937. [Max64] J.C. Maxwell. On reciprocal figures and diagrams of force. Philos. Mag., 27:250–261, 1864.  119  [MB05] T.S. Majmudar and R.P. Behringer. Contact force measurements and stress-induced anisotropy in granular materials. Nature, 435(7045):1079–1082, 2005. [MC10] M. Mailman and B. Chakraborty. Unjamming of Granular Packings as a Constraint Satisfaction Problem: Evidence for a Growing Static Length Scale in Frictionless Packings. ArXiv e-prints, July 2010. [MC11] M Mailman and B Chakraborty. A signature of a thermodynamic phase transition in jammed granular packings: growing correlations in force space. Journal of Statistical Mechanics: Theory and Experiment, 2011(07):L07002, 2011. [MD95] C. Moukarzel and P. M. Duxbury. Stressed backbone and elasticity of random central-force systems. Physical Review Letters, 75(22):4055–4058, 11 1995. [MH06] Sean C. McNamara and Hans J. Herrmann. Quasirigidity: Some uniqueness issues. Phys. Rev. E, 74:061303, Dec 2006. [MiD01] GDR MiDi. On dense granular flows. The European Physical Journal E: Soft Matter and Biological Physics, 14(4):341–365, 2004-08-01. [MK02] H.A. Makse and J. Kurchan. Testing the thermodynamic approach to granular matter with a numerical model of a decisive experiment. Nature, 415(6872):614–617, 2002. [MLRJ+ 08] V. Magnanimo, L. La Ragione, JT Jenkins, P. Wang, and HA Makse. Characterizing the shear and bulk moduli of an idealized granular material. EPL (Europhysics Letters), 81(3):34006, 2008. [MN07] Namiko Mitarai and Hiizu Nakanishi. Velocity correlations in dense granular shear flows: Effects on energy dissipation and normal stress. Physical Review E, 75(3):031305–, 03 2007. [MNN83] M. M. Mehrabadi and S. Nemat-Nasser. Stress, dilatancy and fabric in granular materials. Mechanics of Materials, 2(2):155–161, 8 1983. [MOB96] Brian Miller, Corey O’Hern, and R. P. Behringer. Stress fluctuations for continuously sheared granular materials. Physical Review Letters, 77(15):3110–3113, 10 1996. [Mou25] Cristian F. Moukarzel. Isostaticity in granular matter. Granular Matter, 3(1):41–52, 2001-01-25. 120  [Mou98] Cristian F. Moukarzel. Isostatic phase transition and instability in stiff granular materials. Physical Review Letters, 81(8):1634–1637, 08 1998. [MRDR+ 09] S. McNamara, P. Richard, S.K. De Richter, G. Le Ca¨er, and R. Delannay. Measurement of granular entropy. Physical Review E, 80(3):031301, 2009. [MSLB07] T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer. Jamming transition in granular systems. Phys. Rev. Lett., 98(5):058001, Jan 2007. [MSOC09] Mitch Mailman, Carl F Schreck, Corey S O’Hern, and Bulbul Chakraborty. Jamming in systems composed of frictionless ellipseshaped particles. Physical review letters, 102(25):255501, 2009. [MTDT02] Cristian Moukarzel, M. F. Thorpe, P. M. Duxbury, and M. F. Thorpe. Granular Matter Instability: A Structural Rigidity Point of View, pages 125–142. Springer US, 2002. [Mus63] N.I. Muskhelishvili. Some basic problems of the mathematical theory of elasticity. P. Noordhoff, 1963. [NDBC11] V.B. Nguyen, T. Darnige, A. Bruand, and E. Clement. Creep and fluidity of a real granular packing near jamming. Physical review letters, 107(13):138303 [Ned92] R. Nedderman. Statics and Kinematics of Granular Materials. Cambridge, Cambridge, U.K., 1992. [Oda72a] Masanobu Oda. Initial fabrics and their relations to mechanical properties of granular material. Soils and Foundations, 12(1):17–36, 1972. [Oda72b] Masanobu Oda. The mechanism of fabric changes during compressional deformation of sand. Soils and Foundations, 12(2):1–18, 1972. [OK74] Masanobu Oda and Junichi Konishi. Microscopic deformation mechanism of granular material in simple shear. Soils and Foundations, 14(4):25–38, 1974. [OKNN82] Masanobu Oda, Junichi Konishi, and Siavouche Nemat-Nasser. Experimental micromechanical evaluation of strength of granular materials: Effects of particle rolling. Mechanics of Materials, 1(4):269– 283, 12 1982. 121  [O’N73] P. V. O’Neil. A short proof of Mac Lane’s planarity theorem. Proc. Amer. Math. Soc., 37:617–618, 1973. [OSLN03] Corey S. O’Hern, Leonardo E. Silbert, Andrea J. Liu, and Sidney R. Nagel. Jamming at zero temperature and zero applied stress: The epitome of disorder. Phys. Rev. E, 68:011306, Jul 2003. [OSN06] Srdjan Ostojic, Ell´ak Somfai, and Bernard Nienhuis. Scale invariance and universality of force networks in static granular matter. Nature, 439(7078):828–830, 2006. [Par88] G. Parisi. Statistical field theory. Addison-Wesley, 1988. [PD12] J.G. Puckett and K.E. Daniels. Equilibrating temperaturelike variables in jammed granular subsystems. arXiv preprint arXiv:1207.7349, 2012. [PEKK12] C. P´erez-Espigares, A.B. Kolton, and J. Kurchan. Infinite family of second-law-like inequalities. Physical Review E, 85(3):031135, 2012. [PF12] F. Paillusson and D. Frenkel. Probing ergodicity in granular matter. Physical Review Letters, 109(20):208001, 2012. [POS12] S. Papanikolaou, C. S. O’Hern, and M. D. Shattuck. Isostaticity at frictional jamming. ArXiv e-prints, 1207.6010, 2012. [RARC+ 09] Nicolas Rivier, C´ecile Appert-Rolland, Fran¸cois Chevoir, Philippe Gondret, Sylvain Lassarre, Jean-Patrick Lebacque, and Michael Schreckenberg. Stability and Jamming Transition in Hard Granular Materials: Algebraic Graph Theory, pages 535–544. Springer Berlin Heidelberg, 2009. [RB89] L. Rothenburg and R.J. Bathurst. Analytical study of induced anisotropy in idealized granular materials. Geotechnique, 39(4):601– 614, 1989. [RBR96] Farhang Radjai, Lothar Brendel, and St´ephane Roux. Nonsmoothness, indeterminacy, and friction in two-dimensional arrays of rigid particles. Physical Review E, 54(1):861–873, 07 1996. [RC02] Jean-No¨el Roux and Ga¨el Combe. Quasistatic rheology and the origins of strain. Comptes Rendus Physique, 3(2):131–140, 2002.  122  [RDAR12] F. Radjai, J. Y. Delenne, E. Az´ema, and S. Roux. Fabric evolution and accessible geometrical states in granular materials. Granular Matter, 14(2):259–264, 2012. [REYR06] Vincent Richefeu, Moulay Sa¨ıd El Youssoufi, and Farhang Radja¨ı. Shear strength properties of wet granular materials. Phys. Rev. E, 73(5):051304, May 2006. [Riv06] N. Rivier. Extended constraints, arches and soft modes in granular materials. Journal of non-crystalline solids, 352(42):4505–4508 0022–3093, 2006. [Rou00] Jean-No¨el Roux. Geometric origin of mechanical properties of granular materials. Phys. Rev. E, 61:6802–6836, Jun 2000. [RR05] Farhang Radja¨ı and St´ephane Roux. Contact Dynamics Study of 2D Granular Media: Critical States and Relevant Internal Variables, pages 165–187. Wiley, 2005. [RS81] L. Rothenburg and A.P.S. Selvadurai. A micro-mechanical definition of the cauchy stress for particulate media. In AP.S. Selvadurai, editor, Proceedings of the International Symposium on Mechanical Behaviour of Structured Media, pages 469–486, 1981. [RWJM98] F. Radjai, D.E. Wolf, M. Jean, and J.J. Moreau. Bimodal character of stress transmission in granular packings. Physical review letters, 80(1):61–64, 1998. [Sad09] Martin H. Sadd. Elasticity - Theory, Applications, and Numerics (2nd Edition). Elsevier, 2009. [Sat86] M. Satake. Compatibility conditions for granular assembly. In G. Ishizaka, editor, Science on Form: Proceedings from the First International Symposium on Science on Form, pages 191–199, 1986. [Sat92] M. Satake. A discrete-mechanical approach to granular materials. International Journal of Engineering Science, 30(10):1525–1533, 10 1992. [Sat93] M. Satake. New formulation of graph-theoretical approach in the mechanics of granular materials. Mechanics of Materials, 16(1-2):65 – 72, 1993. Special Issue on Mechanics of Granular Materials.  123  [Sat97] M. Satake. Three-dimensional discrete mechanics of granular materials. In IUTAM Symp. on Mech. of Gran. Mat., pages 193 – 202, Dordrecht, The Netherlands, 1997. Kluwer Academic. [Sat04] M. Satake. Tensorial form definitions of discrete-mechanical quantities for granular assemblies. Int. J. of Sol. and Struc., 41(21):5775 – 5791, 2004. [SEmcG+ 02] Leonardo E. Silbert, Deniz Erta¸s, Gary S. Grest, Thomas C. Halsey, and Dov Levine. Geometry of frictionless and frictional sphere packings. Phys. Rev. E, 65:031304, Feb 2002. [SMCO12] Carl F Schreck, Mitch Mailman, Bulbul Chakraborty, and Corey S O’Hern. Constraints and vibrations in static packings of ellipsoidal particles. Physical Review E, 85(6):061305, 2012. [SMGS99] W. Schwalm, B. Moritz, M. Giona, and M. Schwalm. Vector difference calculus for physical lattice models. Phys. Rev. E, 59(1):1217– 1233, Jan 1999. [SR05] Lydie Staron and Farhang Radjai. Friction versus texture at the approach of a granular avalanche. Phys. Rev. E, 72:041308, Oct 2005. ´ Cl´ement, and D. Levine. [SRC+ 01] D. Serero, G. Reydellet, P. Claudin, E. Stress response function of a granular layer: Quantitative comparison between experiments and isotropic elasticity. The European Physical Journal E: Soft Matter and Biological Physics, 6(2):169– 179, 2001. [SVC02] D. Segr´e, D. Vitkup, and G.M. Church. Analysis of optimality in natural and perturbed metabolic networks. Proc. Natl. Acad. Sci., 99:15112–15117, 2002. [SvHvS07] Kostya Shundyak, Martin van Hecke, and Wim van Saarloos. Force mobilization and generalized isostaticity in jammed packings of frictional grains. Phys. Rev. E, 75(1):010301, Jan 2007. [SVvHvS04] Jacco H. Snoeijer, Thijs J. H. Vlugt, Martin van Hecke, and Wim van Saarloos. Force network ensemble: A new approach to static granular matter. Phys. Rev. Lett., 92(5):054302, Feb 2004.  124  [SW00] R. Shrock and F.Y. Wu. Spanning trees on graphs and lattices in d dimensions. Journal of Physics A: Mathematical and General, 33(21):3881 0305–4470, 2000. [SWM08] C. Song, P. Wang, and H.A. Makse. A phase diagram for jammed matter. Nature, 453(7195):629–632, 2008. [TSS+ 05] Brian P. Tighe, Joshua E. S. Socolar, David G. Schaeffer, W. Garrett Mitchener, and Mark L. Huber. Force distributions in a triangular lattice of rigid bars. Phys. Rev. E, 72:031306, Sep 2005. [TSVvH10] Brian P. Tighe, Jacco H. Snoeijer, Thijs J. H. Vlugt, and Martin van Hecke. The force network ensemble for granular packings. Soft Matter, 6:2908–2917, 2010. [TV10] Brian P Tighe and Thijs J H Vlugt. Force balance in canonical ensembles of static granular packings. J. Stat. Mech., 2010(01):P01015, 2010. [TV11] Brian P Tighe and Thijs J H Vlugt. Stress fluctuations in granular force networks. J. Stat. Mech., 2011(04):P04002, 2011. [TvEV08] Brian P. Tighe, Adrianne R. T. van Eerd, and Thijs J. H. Vlugt. Entropy maximization in the force network ensemble for granular solids. Phys. Rev. Lett., 100:238001, Jun 2008. [TW99] Alexei V. Tkachenko and Thomas A. Witten. Stress propagation through frictionless granular material. Phys. Rev. E, 60(1):687–696, Jul 1999. [UKW05] Tam´ as Unger, J´anos Kert´esz, and Dietrich E. Wolf. Force indeterminacy in the jammed state of hard disks. Phys. Rev. Lett., 94:178001, May 2005. [vH10] M van Hecke. Jamming of soft particles: geometry, mechanics, scaling and isostaticity. Journal of Physics: Condensed Matter, 22(3):033101, 2010. [VHBC99] L. Vanel, D. Howell, R. P. Behringer, and E. Clement. Memories in sand: Experimental tests of construction history on stress distributions under sandpiles. Phys. Rev. E, 60:R5040–R5043, 1999. [VRDEY09] C. Voivret, F. Radjai, J.Y. Delenne, and MS El Youssoufi. Multiscale force networks in highly polydisperse granular media. Physical review letters, 102(17):178001 125  [WCCB96] J. P. Wittmer, P. Claudin, M. E. Cates, and J. P. Bouchaud. An explanation for the central stress minimum in sand piles. Nature, 382(6589):336–338, 07 1996. [WMKG07] Max Wardetzky, Saurabh Mathur, Felix Kalberer, and Eitan Grinspun. Discrete Laplace operators: no free lunch. In Proceedings of the fifth Eurographics symposium on Geometry processing, pages 33–37, Barcelona, Spain, 2007. Eurographics Association. [WNW05] M. Wyart, S. R. Nagel, and T. A. Witten. Geometric origin of excess low-frequency vibrational modes in weakly connected amorphous solids. EPL (Europhysics Letters), 72(3):486, 2005. [Woo90] D. M. Wood. Soil Behaviour and Critical State Soil Mechanics. Cambridge, Cambridge, U.K., first edition, 1990. [Wu01] FY Wu. Number of spanning trees on a lattice. Journal of Physics A: Mathematical and General, 10(6):L113 0305–4470, 2001. [Wya09] Matthieu Wyart. On the dependence of the avalanche angle on the granular layer thickness. EPL (Europhysics Letters), 85(2):24003, 2009. [Yu02] M. Yu. Advances in strength theories for materials under complex stress state in the 20th century. Appl. Mech. Rev., 55(3):169–218, 2002.  126  Appendix A  Torque Balance in 2D 17 16.5  0  ℓc  16  +  ℓc  15.5  s tg  15  −  ℓc  14.5  g 14 13.5 13 12.5  Figure A.1: Geometry used in torque balance computation. 12 8  9  10  11  12  13  14  To verify that torque balance is satisfied identically when ρ = ∇ × ψ, we need to show that (∇ · ρ)g = 0. Labelling the contacts and triangles around a grain as 0 in Figure A.1, we note that stg is parallel to c (g,t) , so that Ag (∇ · ρ)g =  stg × ρt t∈T g  stg ×  = t∈T g  = t∈T g  = t∈T g  1 ( 2  −1 − At  c+  −  c−  c−  )×  −  ψc + −1 − At  1 − + At ψ c − At ψ c t A  =0 on summation around the grain.  127  c+  +  ψc +  c−  −  ψc +  c0  ψc  c+  0  +  ψc  Appendix B  Discrete Calculus in 3D  c g  c stgc g  Figure B.1: Subset of Delaunay triangulation and Voronoi tesselation in 3D. Real contacts are indicated by small dots, with their associated contact vectors solid. Virtual contacts are shown as unfilled circles, with their associated contact vectors dashed. In three dimensions, we can proceed analogously. The Voronoi tesselation assigns to each grain a convex polyhedron (Figure B.1). Each face corresponds to a contact, each edge corresponds to a Delaunay triangle, and each vertex corresponds to a Delaunay tetrahedron. As before, we can define a discrete divergence with  128  Gauss’ theorem: dV ∇ · σ ˆ ≡ V g (∇ · σ ˆ )g ≡  Ac ˜cg · σ ˆc ≡  g  dS · σ ˆ, ∂g  c∈C g  where ˜ = /| | and the natural area for a contact, Ac , is the area of the contact face. We rewrite (1.8) as a sum over contacts, such that each interior contact contributes V c σ ˆ c = − c f c . In analogy with 2D, the natural choice of V c is the product of the contact face area and the contact vector length, i.e., V c = Ac | c |. We then see that Ac ˜cg · σ ˆ c = − ˜cg · ˜c f c = −fgc , so again (∇ · σ ˆ )g = 0 is equivalent to force balance. Again σ ˆ g = (ˆ σ g )T is equivalent to torque balance, and hence we reproduce the expected continuum equations (1.9) in the discrete calculus, at the grain scale. We now seek to define ρˆ so that σ ˆ = ∇ × ρ. ˆ By Gauss’ theorem dV ∇ × ρˆ ≡ V g (∇ × ρ) ˆg≡ g  dSgt × (ρˆt )T ≡  dS × ρˆT . ∂g  t∈T g  The natural way to define dSgt is to decompose the contact faces on either side + − of t into triangles of area Atc± and set dSgt = Atc+ ˜cg + Atc− ˜cg ≡ c∈Cgt Atc ˜cg . It is easily seen from Figure B.1 that At ˜c = 1 (r t − r c ) × st , where r t is the c g  gc  2  intersection of the triangle t with the corresponding edge of the Voronoi cell of g, and stgc circulates anticlockwise around contacts. We can now rewrite (∇ × ρ) ˆ g as V g (∇ × ρ) ˆ g = − 21  (r c − r t ) × stgc × (ρˆt )T t∈T g c∈Cgt  (r c − r g ) × stgc × (ρˆt )T  = − 12 t∈T g c∈Cgt  c g  = − 14 c∈C g  where we have used the identity (1.8), we now see that σ ˆg = ∇ × equivalent to  × stgc × (ρˆt )T ,  t∈T c  t t c∈Cgt sgc = 0 to exchange r ρˆt if cg fgc = 12 t∈T c cg × stgc  fgc = − 21  ρˆt · stgc ≡ − 12 t∈T c  ρˆ · dr,  with r g . Using × (ρˆt )T . This is (B.1)  (∂c)g  t ˆt · c . The latter equation is satisfied identically if we let and 0 = g t∈T c sgc ρ ρˆt = − |s2t |2 ρtgc stgc . This allows (B.1) to be rewritten fgc = t∈T c ρtgc , which shows explicitly that the contact forces depend only on a vectorial quantity, albeit one  129  with an intrinsic orientation. It also shows reduction to the form of the intrinsic formulation. Summing the contact forces incident on a grain, all loop forces appear in equal and opposite pairs, so force balance is identically satisfied. As in 2D, we must explicitly require that there be no virtual contact forces. A new feature in 3D is that ρtgc has a nontrivial gauge freedom ρtgc → ρtgc + (∆ρ)tgc with (∆ρ)tgc = + − B v − B v , for any vector field B v defined on the Delaunay tetrahedra, which are dual to Voronoi vertices. Again, σ ˆ g = (∇ × ρ)g implies that the stress tensor for the packing can be written as a boundary sum [KR96]  Ω  dS × ρT ,  Ac lc × (ρt )T ≡  dV σ ˆ ≡Vσ ˆ=  ∂Ω  c∈BC  where the l’s are oriented outwards. ˆc such that ρˆt = ∇ × ψ ˆc . Using Stokes’ theorem, we set We now seek ψ ˆ · dS ≡ At ∇ × ψ ˆ (∇ × ψ)  t  ˆc · ψ  · s˜t ≡  t  c∈C t  c  ˆ · dr, ψ  ≡ ∂t  where s˜ = s/|s| and the contour is oriented anticlockwise around s˜t . A natural ˆc is to set ψ ˆc = 1 ˜c ˜c ψ c ; then we have At ρtgc = choice to guarantee ρˆt = ∇ × ψ 4 − 18 |st | c∈C t cgt ψ c , in close analogy with the 2D case. The computation that this leads to torque-balanced force configurations is almost identical to the 2D case, and for the same reasons, works only for identical spheres. We can count degrees-of-freedom to check these formulations. Applying Euler’s formula to the Voronoi polyhedron of grain g, we have NVg − NTg + NCg = 2. Each vertex is the confluence of three edges of the polyhedron, so that 3NVg = 2NTg . Summing over the entire packing, 4NV −NBV −3NT +NBT +2NC −NBC = 2N and 12NV −3NBV = 6NT −2NBT . Since the entire packing is itself a convex polyhedron with trivalent vertices, we have NBV − NBT + NBC = 2 and 3NBV = 2NBT , so that 4NV − 3NT + 2NC = 2N + 2 and 2NV = NT . With these relations we can write H = 3NT − 3NRG − 3(NV C − NV G ) − 3(NV − 1). This shows that the loop forces ρ are constrained by torque balance and the virtual contact constraints, NV G of which are redundant. The gauge freedom has a dimension of NV − 1 because adding a constant to B v does not change ∆ρ. For ψ, we write M3 = NC − 3NV C + 6NV G − 3(NT − NC ) − (NC − N ) − (6NC − 6 − 12NV + N ). All but the first three terms correspond to gauge freedoms. First, we have gauge transformations of the form ∆ψ c = r c · B, which lead to (∆ρ)tgc = −stgc ×B. Here B is a constant, but it may be derived from a fluctuating 130  field on the triangles, viz. B = t∈T c B t . This gives 3NT unknowns constrained + − by 3NC equations. Second, we can have ∆ψ c = | |−1 ˜c · (∆ψ g − ∆ψ g ), with ∆ψ g = r g B, where again B is a constant. In this case it can be derived from B = c∈C g B c , which gives NC unknowns constrained by N equations. Gauge transformations of this type leave invariant ρ. There remains a gauge subspace of dimension 6NC − 6 − 12NV + N , the significance of which is unclear.  131  Appendix C  Convexity of ψ The force acting on a plane with unit normal n and area A is An · σ ˆ , with normal component An · σ ˆ · n. This is always positive, and hence repulsive, if and only if σ ˆ is a positive-definite tensor. In 2D, σ ˆ = ∇ × ∇ × ψ can be inverted to write ∇∇ψ = −ˆ ε·σ ˆ · εˆ. Using a matrix identity, this can be written ∇∇ψ = tr(ˆ σ )δˆ − σ ˆ,  (C.1)  where δˆ is the identity tensor. ˆ ≡ ∇∇ψ, convexity of ψ is equivalent to positive definiteness of H. ˆ Writing H ˆ The latter is equivalent to the statement that H has positive eigenvalues. Writing λi and ui for the eigenvalues and eigenvectors of σ ˆ , we have tr(ˆ σ ) = i λi . From ˆ ˆ (C.1), H · uj = i=j λi uj , so uj is an eigenvector of H with eigenvalue i=j λi . This is positive if σ ˆ is positive definite.  132  Appendix D  Numerical Simulations We have performed numerical simulations using DEM, a Discrete Element Method code written by Jim McElwaine. The software models granular materials as a collection of soft disks or spheres, interacting via linear contact laws c fN = kN δ c  dfTc /dt = kT VTc , where kN and kT are spring constants, δ c is the overlap at c, and VT is the tangential component of the relative velocity at c. The contact forces are subject to the Coulomb inequality c |fTc | ≤ µfN ,  such that contacts begin sliding when this inequality is saturated. The code directly solves Newton’s laws mg  d2 r g = dt2  Ig  dω g dt  fgc c∈C g  r c × fgc ,  = c∈C g  where mg and I g are the mass and moment of inertia of grain g. For such interactions, the appropriate definition of the stiffness number is κ = kN /P. Simulations discussed in Chapter 3 were performed with κ > 103 , while those discussed in Chapter 6 have κ > 102 .  133  Appendix E  Differential Geometry In the exact geometrical formulation of ψ, it is natural to ask how physical quantities relate to mathematical quantities, particularly those for which differential geometry offers applicable theorems. Here we discuss the physical interpretation of curvature and the Gauss-Bonnet Theorem. A first, trivial point is that since ψ and r have different units, we can rescale ψ by any positive factor. By taking a factor which is arbitrarily small, the surface becomes asymptotically flat; we assume that such a rescaling has already been done. This allows us to keep only leading-order terms in ψ. We will work in the continuum, although these quantities have exact analogs for the original discrete surface. The most basic mathematical quantities are the metric tensor1 g ˆψ = δˆ + (∇ψ)(∇ψ),  (E.1)  which measures distances along the surface, and the shape tensor2 I=  1 1 + |∇ψ|2  ∇∇ψ,  which measures curvature. The Gaussian curvature is det I K= = det(∇∇ψ) + O(ψ 4 ) det g ˆψ  (E.2)  (E.3)  Using ∇∇ψ = εˆ · σ ˆ · εˆT , we have K = det(ˆ σ ) + O(ψ 4 ),  (E.4)  while the mean curvature is 1 tr g ˆψ · εˆ · I · εˆT det g ˆψ 1 tr σ ˆ + O(ψ 2 ) = 1 + O(ψ 2 ) = 2P + O(ψ 2 ).  H=  1 2  Also called the first fundamental form. Also called the second fundamental form.  134  (E.5)  The celebrated Gauss-Bonnet theorem states that kg ds + 2π,  dA K =  (E.6)  ∂Ω  Ω  where Ω is the surface and kg is the geodesic curvature of the boundary. The physical interpretation of this theorem is that dA det(ˆ σ)  A=  (E.7)  Ω  depends only on boundary quantities. In fact, A is exactly the area of the MaxwellCremona reciprocal tiling discussed in Chapter 2; compare equation (2.20).  135  Appendix F  Renormalization of ρ. ¯ ) with ρ . Throughout, we will ignore correlations between Here we compare ρ(r r and other loop quantities, which amounts to assuming that the geometry is homogeneous on the averaging scale. Initially, we work to quadratic order. We find g ˆ ·ρ =−  1 A  c  ¯ ) + 1 tc tc : ∇∇ψ(r ¯ ) tc : ∇ψ(r 2  c∈C  ¯ ¯ ) − 1G ˆ =g ˆ × ∇ψ(r 2 2 : ∇∇ψ(r ),  (F.1)  with ˆ2 = 1 G A  c c c  tt.  (F.2)  c∈C  ˆ vanishes on average and hence does not give a systematic correction. The tensor G To see this, consider the related tensor g ˆ2 = −(∇ × (rr)) 1 c c c = r r A c∈C  1 = A  c  (tc + r )(tc + r )  c∈C  ˆ2 + =G  1 A  c c  r r +  c∈C  1 A  c  r rc .  (F.3)  c∈C  In index notation1 , we find ˆ 2 · εˆ)ijk = (ˆ (ˆ g2 − G g · εˆT )ij (r )k + (ˆ g · εˆT )ik (r )j . 1  (F.4)  ˆij Bk . There is an insufficiency of notation here, since (ˆ g · εˆT )ik (r )j is not of the form A  136  However, discrete calculus implies g ˆ2  ijk  = − (∇ × (rr)) = −εij rk − εik rj = (ˆ g · εˆT )ij (r )k + (ˆ g · εˆT )ik (r )j .  (F.5)  ˆ = 0, and to quadratic order there is no renormalization of ρ; ρ = ρ. ¯ Hence G 2 In fact, this argument can be continued to all orders; only the tensor algebra ˆ for the part of A ˆ which is symmetric in its is more complicated. Let us write A n  last n indices. Then we define  g ˆn = − ∇ × (r n ) 1 c c n = (r ) A c∈C  1 = A  (tc + r )n  (F.6)  n (tc )k (r )n−k , k  (F.7)  c c∈C  so that 1 g ˆn = A n  n c k=1  c∈C  n  where the sum runs from k = 1 rather than k = 0 because correlations between r and other loop quantities, we have g ˆn  = n  1 A  n  n (tc )k k  c c∈C  k=1  = 0. Ignoring  r n−k n  n  =n g ˆ  c c∈C  · εˆT r n−1 + n  k=2  n k  1 A  c c∈C  (tc )k  (r)n−k .  (F.8)  n  However, with discrete calculus we observe that = −∇ × (r)n  g ˆn n  n  = −n (∇ × r)r n−1 n  = −n εˆ r n−1 . n  137  (F.9)  Since g ˆ  ˆ we find = δ, n k=2  n ˆ n−k Gk r = 0, k  (F.10)  n  with ˆk = G  1 A  c  (tc )k  .  (F.11)  c∈C  ˆ k = 0 for all k ≥ 2. Indeed, for n = 2 we find Inductively, (F.10) implies G ˆ2 = G ˆ 2 . This then implies 0 = G ˆ3 = G ˆ 3 , etc. So we find that ρ = ρ¯ to 0=G 2  3  all orders in n, up to neglected correlations.  138  Appendix G  Green’s Function for a Half-plane The solution to ∆∆ψ = 0 with σ ˆ = ∇ × ∇ × ψ and σ ˆ ·y by  y=0  = 2πA δ(x) is given  ψ(z, z¯) = 21 z¯f (z) + 21 g(z) + c.c., with f (z) = −A log z and dg/dz = A¯ log z. One can directly compute A¯ A¯ z + 2z 2 2z ¯ Az A ∂z¯2 ψ = 2 + = ∂z2 ψ 2¯ z 2¯ z A¯ A ∂z z¯ψ = − − . 2¯ z 2z ∂z2 ψ =  We use σxx = 2∂z z¯ψ − ∂z2 ψ − ∂z¯2 ψ σyy = 2∂z z¯ψ + ∂z2 ψ + ∂z¯2 ψ σxy = −i(∂z2 ψ − ∂z¯2 ψ) and the identities [z w] ¯ =w·z [z w] ¯ =w×z [z 2 ] = x2 − y 2 [z 2 ] = 2xy  139  (G.1)  with z = ( [z], [z]), etc. Then A·z z2 A ·z ∂z¯2 ψ = z¯2 A·z ∂z z¯ψ = − |z|2 ∂z2 ψ =  and 4x2 A · r |r|4 4y 2 A · r =− |r|4 4xyA · r =− |r|4  σxx = − σyy σxy  (G.2)  where we have switched to a full vector notation with r = (x, y). For the anisotropic Green’s function, we have ψ(z, w) = f (z) + g(w) with df (z) = A1 log z, dz  dg(w) = A2 log w. dw  We use the identities [Sad09] σxx = −f (z) + µ2 g (w) + c.c.  (G.3)  σyy = +f (z) + g (w) + c.c.  (G.4)  σxy = −if (z) − µg (w) + c.c.  (G.5)  and find A1 · r A3 · w + |r|2 |w|2 A1 · r A2 · w + =+ |r|2 |w|2 r × A1 A4 · w =+ − , |r|2 |w|2  σxx = −  (G.6)  σyy  (G.7)  σxy  with A3 = A2 µ2 and A4 = A2 µ.  140  (G.8)  Appendix H  Fabric Change due to Contact Opening Here we consider the renormalization of fabric due to contact opening. When the latter occurs, contact positions r c do not change, but loop positions do change, because of the definition z r = c∈C r c . We will use a basic property of the fabric tensor z Fˆ ≡  tc tc = c ∈C  rc rc − z r r . c ∈C  It follows that if loops and merge into by opening of their common contact ˆ c, then the change in NRC F is localized to the contributions from these loops. The change in NRC Fˆ is δ Fˆ− ≡ z Fˆ  − z Fˆ − z Fˆ tc tc −  =  tc tc − c ∈C  c ∈C  = −z r r  tc tc c ∈C  c c  − 2r r + z r r + z r r .  We have basic relations z  =z +z −2  and z r  rc +  = c ∈C \{c}  rc c ∈C \{c}  = z r + z r − 2r c , which imply z tc = z (r c − r ) = (z + z − 2)r c − (z r + z r − 2r c ) = z tc + z tc . 141  Hence δ Fˆ− = − z r + z r − 2r c r  − 2r c r c + z r r + z r r  = z r (tc − tc ) + z r (tc − tc ) − 2r c tc = z r + z r − 2r c tc − z r tc − z r tc = z (r c − tc )tc − z (r c − tc )tc − z (r c − tc )tc = −z tc tc + z tc tc + z tc tc . This form is suggestive, but it will be useful to eliminate tc completely, giving 1 z tc + z tc z tc + z tc + z tc tc + z tc tc z z z z z 1− z tc tc + 1 − z tc tc − tc tc + tc tc z z z  δ Fˆ− = − =  z −2 c c z −2 z z z tt + z tc tc − tc tc + tc tc z z z z z z c c z c c = tc − tc tc − tc − 2 t t −2 t t . z z z  =  142  Appendix I  Jacobians I In changing coordinates from f to ρ to ψ, Jacobians are introduced into the partition function. The relation of f to ρ can be understood by a classic subject in graph theory: cycle spaces and bond spaces [Bol79, BM76]. The contact network forms a graph, with grains as vertices and real contacts as edges. We consider f as a vector-valued function defined on the real contacts. The contacts are oriented, so each contact c = gg for two grains g and g . Formally, we have f c ec (c ),  f (c ) =  (I.1)  c  where ec (c ) is an abstract basis element for each contact. Each ec is a function such that ec (c ) = δcc . In this formalism, a change of coordinates consists of selecting NRC linear combinations of basis elements, as in linear algebra. Associated to the contact network there are certain natural bases, associated to loops, and cuts, defined below. To each loop we assign a basis element z such that z (c) = +1 if the contact c circulates anticlockwise around the loop, −1 if the contact circulates clockwise around the loop, and 0 otherwise. A cut p is a separation of the grains of the packing into two sets; call them Ap and Bp . We define the basis element up such that if the contact c goes from Ap to Bp , then up (c) = +1, if the contact goes from Bp to Ap , then up (c) = −1, and 0 otherwise. The number of linearly independent cuts is exactly N −1. Thus, using Euler’s formula NRC = NL +N −1, we can write f = ρ + η, where ρ and η are vector-valued functions defined on the loops, and the cuts, respectively. We have f (c ) =  η p up (c ).  ρ z (c ) +  (I.2)  p  The span of ρ is the cycle space of the graph, and the span of η is the bond space of the graph. These spaces are orthogonal, hence the coordinate change is 1 to 1. Moreover, the Jacobian is exactly ∂f 2 = NST , ∂(ρ, η) 143  (I.3)  where NST is the number of spanning trees in the contact network1 . This quantity has been well-studied for lattices, and is related to statistical mechanical problems in percolation [Wu01, SW00]. For large systems, it is exponential in the number of edges, NST = eqNRC . The constant q has been evaluated for various lattices. For the 4-8-8 (bathroom tile) and honeycomb planar lattices with z¯ = 3, we have q = 0.52 and q = 0.54, respectively, while for the square lattice with z¯ = 4 we have q = 0.58. All of these lattices can be constructed as contact networks with monodisperse grains, hence we expect that q for realistic contact networks is in this range. To apply these results to our partition function, we first note that ρ as defined above exactly corresponds to ρ as defined in the main text. Since loop forces automatically describe configurations in force balance, when we make a coordinate change from f to (ρ, η), the N vector constraints of force balance involve only the field η, except at the boundary of the packing. Solving the first N − 1 vector equations, we obtain η[ρ]. The final vector equation is a gauge fixing constraint on the ρ field. We have Z=  2 λNRC NST  ˆσ DρDη e−V α:ˆ Θ[f ]  2 λNRC NST  ˆσ Dρ e−V α:ˆ δG [ρ] δT B [ρ] ΘC [ρ],  T  =  (I.4)  T  where δG is the gauge fixing constraint. To introduce the field ψ, we first need to transform from the contact network to the Delaunay triangulation. To do so, we introduce the field ρ defined on the Delaunay triangles, with dρt δ(ρt − ρ  1=  (t)  ).  (I.5)  t  We insert this into Z, change the order of integration, write everything in terms of ρ defined on the triangles, and eliminate ρ defined on the loops by eliminating a δ function from each loop. We thus have remaining NT − NL vector equations equating ρ within a loop. By a graph identity NT − NL = NV C , these are exactly + − the NV C virtual contact constraints δ(f c ) = δ(ρt (c) −ρt (c) ) which say that there can be no contact forces on virtual contacts. Finally, we need to introduce new orthogonal fields ψ and φ and perform a similar transformation as above to eliminate the torque balance constraints. 1  This is almost exactly Corollary 12.4 in [BM76]. The only modification is that we have vector-valued functions, hence our result is squared.  144  However, since all fields are written in terms of the Delaunay triangulation, the resulting Jacobian will only depend on the slow geometry, not on the fast geometry T . Hence it contributes only an overall and irrelevant constant to the partition function. We thus obtain (6.8).  145  Appendix J  Jacobians II For each isostatic cluster, the integrals over the internal contacts give a Jacobian Dψ|C i δVi C .  Ji =  (J.1)  i = NVi C ((6.11)). From the fact that each cluster is isostatic, H i = 0, we have NRC Hence each real contact in the cluster can be paired with a virtual contact. This integral then involves two types of integrals: (i) integrals over virtual contacts c of c ), and (ii) integrals over real contacts c of δ(f c ) for a virtual contact c . For δ(fN T a large class of cluster geometries, contacts can be paired up in such a way that in (ii) the virtual contact c is adjacent to the real contact c. Let us first perform these integrals, assuming that they can all be done independently. Then for (i), we have c  dψ c δ( |A˜c| ψ c + · · · )  c dψ c δ(fN )=  = |A˜c |. +  (J.2)  −  Here 1/A˜c = 1/At (c) +1/At (c) and · · · refers to terms depending on ψ at adjacent contacts. We have also used the fact that | c | = 1 for real contact. In case (ii) we have dψ c δ(fTc ) = =  dψ c δ( | c |At , 2At  c  ×  c  c 1 | c |At ψ  + ···) (J.3)  since c and c form the sides of the triangle t. We have ignored the fact that when we sequentially perform the integrals, the earlier integrations affect the later ones; e.g. we solve ψ c in terms of ψ at its neighbours c , and then we perform a ψ c integral. This is generically affected by the original c integral. The key point is that the coupling has a random sign. We assume that these signs cancel, so that the Jacobian is the product of integrals of the above type. 146  The triangle areas are not broadly distributed so we replace each by their mean value A/NT = A/(2N ). Assembling results, we find i  J ≈  A 8N  147  NVi C  .  (J.4)  Appendix K  Mean-field Relation Between ρ and ρi. Within each isostatic cluster i, we make an approximation for the loop force in loop of the form δρM F = λ δρi ,  (K.1)  where ρi =  1 Ai  A ρ =− ∈Li  1 Ai  c  ψc,  (K.2)  c∈C ∂i  λ is a parameter. Here we show that, in the simplest mean-field approximation, λ = 1. First, the form of (K.1) should ensure that for any choice of the ψ c ’s, (K.2) holds with ρ → ρM F . This requires δρi = =  1 Ai 1 Ai  A δρM F ∈Li  A λ δρi ∈Li   =−  1 Ai  c∈C  1 Ai ∂i   Aλ  c  δψ c ,  ∈Li  hence 1 Ai  A λ = 1.  (K.3)  ∈Li  To determine the λ ’s, we note that when the cluster is a single triangle, we must have λ = Ai /A = 1. Conversely, when the cluster is large, we use the 148  c  c0 l  Figure K.1: Schematic of homotopy construction. Only grains on RC and a subset of C ∂i are shown. The green lines show paths from c ∈ C ∂i to c0 ∈ RC for M = 3. The dashed lines are the averages over j. c  exact discrete calculus relation ψ c − ψ c0 = − c0 dr × ρ to pair each boundary contact with one on the boundary of an loop , establishing a homotopy between the cluster boundary and the loop, allowing an estimate of the λ ’s. For each c ∈ C ∂i , we consider M paths from c to contacts c0 (c, j) ∈ RC , chosen so that M 1 c c = c z i /z (see Figure K.1). j ≡ M c∈C ∂i |c0 (c,j)=c j=1 c∈C ∂i |c0 (c,j)=c We have c  Ai δρi = −  δψ c  c∈C ∂i c c  =−  δψ c0 (c,j) −  =−  j c  zi z  dr × δρ c0 (c,j)  c∈C ∂i c  δψ c +  c  .  dr × δρ c0 (c,j)  c∈C ∂i  c ∈RC  (K.4)  j  The first term in (K.4) is simply (z i /z )A δρ . The second term is a discrete area integral over the cluster, outside of . Applying (K.1), c c c∈C ∂i  c  dr × δρM F c0 (c,j)  c  = j  c∈C ∂i  × δρi .  drλ c0 (c,j)  j  The simplest solutions will have, for 1 along the contour from c0 to c, λ 1 unc correlated with c c0 (c,j) dr. In a mean-field approximation, we therefore make a 149  replacement λ  1  ¯≡ →λ  assuming that the we let Ai δρi −  1 ’s  Ai  1 −A  A λ =1+ ∈Li , =  Ai  A (1 − λ ), −A  are sampled, on average, proportional to their area1 . Hence  zi ¯ A δρM F = λ z  c c c0 (c,j)  c∈C ∂i  ¯ =λ  c  × δρi  dr c  (r − r  j  c0 (c,j)  j)  × δρi  c∈C ∂i   ¯ =λ  c c  r −  c∈C ∂i  Using the geometric identity Ai δρi −  c∈RC  c rc    zi z  c  r c  × δρi  c ∈RC  = A εˆT , valid for any loop, we see that  i zi ¯ Ai − z A A λ δρi = λ z z  δρi ,  (K.5)  from which it follows that λ = 1.  1 In fact, loops will be sampled with a weight proportional to |r the same conclusion.  150  1  − r |−1 A 1 . This leads to  Appendix L  Contact-angle Anisotropy In the main text, we suppose that contact angles are isotropically distributed, and that the angles γ c vanish identically. Here we remove these assumptions. We will need to specify the joint distribution P(γ c , Ac , θc ) = P(γ c |Ac , θc ) P(Ac |θc ) P(θc ). We have performed simulations of nearly rigid disks in steady shear, and found that, with negligible error, the contact areas Ac are independent of γ c and θc , as previously observed [RB89, VRDEY09]. This implies that P(γ c |Ac , θc ) ≈ P(γ c |θc ) and P(Ac |θc ) ≈ P(Ac ). When crystallization is avoided, the distribution of γ c has 2 components: a δ-function peak at γ = 0, corresponding to grains which are locally in a triangular lattice configuration, and a diffuse peak, vanishing at γ − aγ sin(2θ) = ±π/2, with aγ an anisotropy parameter. As shown in Figure L.1, our simulations are relatively well-fit by P(γ c |θc ) = λδ(γ c ) + (1 − λ)  1 cosn (γ c − aγ sin(2θc )), Cn  (L.1)  for |γ c − aγ sin(2θc )| < π/2, with λ ≈ 0.1, n = 8, and Cn a constant1 . Deviations from the fit are visible near γ c = ±0.3. These are a signature of nascent crystallization. Because this is an artifact of monodispersity, which is absent in realistic packings, we do not refine (L.1) in this case. The anisotropy parameter aγ is small and can be written in terms of the anisotropy in the contact angle distribution, as discussed below. In particular, (L.1) should be considered only as valid to O(a2γ ). We have checked that in very strongly polydisperse disk packings, which look to the eye very different, (L.1) is robust, with the same values of the parameters λ and n. In polydisperse packings, the humps near γ c = ±0.3 are absent. Many studies have found that the contact angle distribution is well fit by its first two Fourier components, viz., [RDAR12, RR05] P(θc ) = 1  1 (1 + ac cos(2θc )) . 2π  The normalization constant Cn = π[(n − 1)(n − 3) · · · (1)]/[n(n − 2) · · · (2)].  151  (L.2)  −3  x 10  2.5  P(γc)  2  1.5  1  0.5  0 −1.5  −1  −0.5  0 γc  0.5  1  1.5  Figure L.1: Probability distribution of γ c from simulations of a simple shear flow (black), as compared with the fit (L.3) (red). The δ-peak extends above the displayed area. Because the Delaunay triangulation tiles the plane, the variables γ c and θc are not independent. In fact, by a geometrical identity, consistency of (L.1) and (L.2) to O(a2c ) requires that aγ (1 − λ) = ac . c =M ˆ into M11 +M22 = As in the main text, we decompose the matrix c c f¯N c c c c c . f¯N , M11 − M22 = cos(2θ )f¯N , and M12 + M21 = sin(2θ )f¯N We consider first the conditional expectation of normal forces at a given contact angle θ, denoted with a subscript γ|θ. Using (L.1), this is f¯N  γ|θ  =  2A NRC  P + T cos(2θ + 2ω) + T ac  Cn−2 sin(2θ) sin(2θ + 2ω) + O(a2c ) . Cn (L.3)  This expression is compared to data from numerical simulations in Figure L.2. The fit is extremely good, which is not surprising since the stress is a linear functional of the contact forces. It implies that 2A f¯N = P + 21 T ac cos(2ω) NRC A cos(2θc )f¯N = (T cos(2ω) + ac P ) NRC A sin(2θc )f¯N = − T sin(2ω). NRC 152  5  x 10 6  5  fN  4  3  2  1  0 0  0.5  1  1.5  2  2.5  3  3.5  θc  c from simulations of a simple shear flow Figure L.2: Probability distribution of fN (×), as compared with the theoretical expression (L.3) (o). Error bars show the magnitude of fluctuations across the packing.  As in the main text, we take coordinates so that   α + η 0  · Rν , α ˆ = R−ν ·  0 α−η where Rν is an anticlockwise rotation matrix with angle ν. We saw above that α ˆ must be positive definite; this requires α > |η| ≥ 0. By allowing η to take either sign, we can require 0 ≤ ν < π/2. In these coordinates, for a contact vector c = (cos(θ c ), sin(θ c )), we have αc = c · α ˆ · c = α + η cos(2θc + 2ν). Writing η/α = sin(2λ), with |λ| ≤ π/4, we find the equations of state 1 2A (1 − ac cos(2ν) tan(λ)) = P + 12 T ac cos(2ω) α cos(2λ) NF C 1 A − tan(λ) cos(2ν) − ac 12 cos(4ν) sec2 (λ) + sin2 (2ν) = (T cos(2ω) + ac P ) α cos(2λ) NF C 1 A tan(λ) sin(2ν) + ac − sec2 (λ) + cos(2ν) = − T sin(2ω). α cos(2λ) NF C  153  We notice that σ ˆ is no longer aligned with α, ˆ i.e., ω = ν. We also notice that stress anisotropy has coupled to geometrical anisotropy, even in the isotropic equation which relates P to α. It is likewise possible to find corrections to the fluctuations induced by geometrical anisotropy. We omit these, as the lengthy expressions offer little insight.  154  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items