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 ma- terials 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 homo- geneous 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 iso- staticity 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 con- tact 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 Ele- ment 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 Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Symbols - Chapters 1-4 . . . . . . . . . . . . . . . . . . . . . . x List of Symbols - Chapters 5-6 . . . . . . . . . . . . . . . . . . . . . . xiii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Outline of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 Discrete Calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1 Intrinsic formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Extrinsic formulation . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Discussion and extensions . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.1 Physical interpretation . . . . . . . . . . . . . . . . . . . . 20 2.3.2 Polydisperse packings . . . . . . . . . . . . . . . . . . . . . 23 2.3.3 Force law . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 Isostaticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1 Geometric formulation . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Discrete calculus equations . . . . . . . . . . . . . . . . . . . . . . 36 3.3 Final equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 iv 3.3.1 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3.2 Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.3 Homogeneous and isotropic fabric . . . . . . . . . . . . . . 43 3.3.4 Homogeneous and anisotropic fabric . . . . . . . . . . . . . 45 3.3.5 Inhomogeneous and anisotropic fabric . . . . . . . . . . . . 47 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4 Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1 Kinematic variables . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Constitutive laws . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3 Geometric compatibility . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3.1 Floppy modes . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.3.2 General deformations . . . . . . . . . . . . . . . . . . . . . 60 4.4 Fabric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.5 Summary of Chapters 2-4 . . . . . . . . . . . . . . . . . . . . . . . 65 4.6 Homogeneous solutions . . . . . . . . . . . . . . . . . . . . . . . . 69 4.6.1 Simple shear flow . . . . . . . . . . . . . . . . . . . . . . . 70 4.6.2 Pure shear flow . . . . . . . . . . . . . . . . . . . . . . . . 71 4.6.3 Isotropic expansion . . . . . . . . . . . . . . . . . . . . . . 72 4.6.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5 Hyperstaticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.1 Phase space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2 Measure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.2.1 Edwards’ approach . . . . . . . . . . . . . . . . . . . . . . 76 5.2.2 Force Network Ensemble . . . . . . . . . . . . . . . . . . . 78 5.3 Properties of the enlarged Force Network Ensemble . . . . . . . . 79 5.3.1 Contact potential . . . . . . . . . . . . . . . . . . . . . . . 79 5.3.2 Variational principles . . . . . . . . . . . . . . . . . . . . . 80 5.3.3 Microcanonical v.s. Canonical Ensembles . . . . . . . . . . 81 5.3.4 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.4 A note on percolation . . . . . . . . . . . . . . . . . . . . . . . . . 83 6 Mean-field Theory of Hyperstatic Packings . . . . . . . . . . . . 85 6.1 Gauge fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.2 Isostaticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.3 Isostatic cluster partition function . . . . . . . . . . . . . . . . . . 92 6.3.1 Mean-field approximations . . . . . . . . . . . . . . . . . . 93 6.3.2 Computing the partition function . . . . . . . . . . . . . . 95 6.4 Global partition function . . . . . . . . . . . . . . . . . . . . . . . 97 v 6.4.1 Mean-field contact forces . . . . . . . . . . . . . . . . . . . 98 6.4.2 Free-energy minimization . . . . . . . . . . . . . . . . . . . 99 6.4.3 Saddle-point equations . . . . . . . . . . . . . . . . . . . . 100 6.5 Free energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.5.1 Equation of state . . . . . . . . . . . . . . . . . . . . . . . 103 6.5.2 Fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 . . . . . . . . . . . . . . . . . . . . . . . 58 vii List of Figures 1.1 Sketch of 2D packing of rigid grains, showing the position of a contact, rc, and the position of a grain, rg. The hatched grains are rattlers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 Delaunay triangulation and Voronoi tesselation . . . . . . . . . . . 16 2.2 Maxwell-Cremona Diagram . . . . . . . . . . . . . . . . . . . . . . 21 3.1 Contact-network tessellation and its dual . . . . . . . . . . . . . . 27 3.2 Notation convention for `c, tc`, and s ` g vectors. We also have s c = r` ′ − r` = tc` − tc`′ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 Geometric interpretation of ψ. . . . . . . . . . . . . . . . . . . . . 32 3.4 Averaging cells Ω`m. . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.5 Simple-shear simulation snapshot. . . . . . . . . . . . . . . . . . . 34 3.6 Convergence of fabric tensor ĝ to δ̂. . . . . . . . . . . . . . . . . . 35 3.7 Convergence of fabric tensor F̂ . . . . . . . . . . . . . . . . . . . . . 41 3.8 Green’s function to a normal forcing on an isotropic geometry . . . 45 3.9 Green’s function to a normal forcing on an anisotropic geometry . 48 3.10 Green’s function to a shear forcing on an isotropic geometry . . . . 48 3.11 Green’s function to a shear forcing on an anisotropic geometry . . 49 4.1 Mean-field phase diagram for rigid, quasistatic, 2D granular materials. 67 6.1 Simple-shear simulation snapshot. . . . . . . . . . . . . . . . . . . 85 6.2 Isostatic cluster decomposition . . . . . . . . . . . . . . . . . . . . 91 6.3 Simulation results: √ P vs. z and ∆z vs φ. . . . . . . . . . . . . . 107 6.4 Proportion of rattlers . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.5 Prediction of 1/α . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.6 Ratio of α predictions . . . . . . . . . . . . . . . . . . . . . . . . . 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 L.2 Probability distribution of f cN . . . . . . . . . . . . . . . . . . . . . 153 ix List of Symbols - Chapters 1-4 Symbol Definition Equation Reference Ac sc × `c (3.18) A Maxwell-Cremona area (2.19) Cg Set of contacts belonging to grain g (1.5) d Dimension of space (1.7) D Mean linear dimension of grains (1.2) D̂ 12(∇w + (∇w)T ) F g Body force exerted on grain g (1.5) F̂ ` Fabric tensor on loop ` (3.32) F̂ ′ Fabric tensor (3.39) f cg Contact force exerted on grain g at contact c (1.5) f cN Normal component of contact force at c (4.8) f cT Tangential component of contact force at c (4.8) Ĝc Contact normal fabric tensor at c (4.12) g Gravity vector (4.55) ĝ` Fabric tensor (3.11) H Generalized hyperstaticity (2.5) H ′ Hyperstaticity (2.5) I Inertial number (1.2) Lc Set of loops adjacent to contact c (2.9) `c Contact vector Figure 3 `cg Contact vector Figure 2.1 M̂± δ̂ ± 1µ ε̂ (4.17) NL Number of loops (2.6) NRC Number of real (force-bearing) contacts (2.5) NRG Number of real (force-bearing) grains (2.5) NS Number of sliding contacts (1.7) ns Number of sliding contacts per grain (1.7) P Pressure (2.28) qcg Weight (4.22) Continued on next page x Table 1 – continued from previous page Symbol Definition Equation Reference <[z] Real part of z (3.48) rc Position of contact c (1.5) rg Centre-of-mass of grain g (1.5) r` Position of loop ` (2.11),(3.9) sc Cross-contact vector Figure 3 stg Vector pointing across triangle t (2.13) tcg r c − rg (4.1) tc` r c − r` (3.3) tr Tensor trace (1.12) u Continuum velocity field (4.45) V Volume of the packing (1.8) V c Relative velocity of g′(c) with respect to g(c) at c (4.2) vc Weight (3.5) v Continuum grain velocity (4.24) vg Grain linear velocity (4.1) w Continuum contact velocity (4.24) wcg Contact velocity of g at c (4.1) z̄ Contact number (1.7) z̄L Loop contact number (2.7) ziso Isostatic contact number (1.7) z` Loop contact number (3.8) Γ̂ Continuum Cosserat strain rate tensor (4.43) Γ̂c Cosserat strain rate tensor at c (4.3) ∆z̄ z̄ − ziso ε̂ Permutation tensor (Levi-Civita symbol) (1.11) κ Stiffness number (1.1) µ Coulomb friction coefficient (1.6) ρ Continuum loop force (3.27) ρ`g Loop force on loop `, oriented to act on grain g (2.9), (3.1) ρt Loop force on triangle t (2.14), (2.15) σ̂ Stress tensor (1.8) φ Area fraction (4.53) ϕ Continuum torque potential (3.27) ϕ`g 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 Table 1 – continued from previous page Symbol Definition Equation Reference ψc Airy stress function on contact c (2.16) ψc` Airy stress function on contact c in loop ` (3.7) ξ √ trF̂ (3.40) ω Continuum grain angular velocity (4.24) ωg Grain angular velocity (4.1) 〈·〉 Coarse-grained average (3.13) (∇ϕ)c Discrete gradient of ϕ` at c (3.18) (∇v)c Discrete gradient of vg at c (4.31) (∆ϕ)` Discrete Laplacian of ϕ` at ` (3.17) (∆v)g Discrete Laplacian of vg at g (4.30) xii List of Symbols - Chapters 5-6 Symbol Definition Equation Reference Â Angoricity (5.4) C 〈(P − 〈P 〉)(T − 〈T 〉)〉 D Multiple integration (5.7) F [P] Free-energy functional (5.11) f̄ Mean-field contact force f̄ cN −`c · f̄ c, normal component of f̄ c H Generalized hyperstaticity (5.1) H ′ Hyperstaticity (5.2) J i Cluster Jacobian (6.16) `c Contact vector Figure 3 `0(c) Loop adjacent to c, exterior to cluster (6.20) NFC Number of contacts on force chains (6.12) NGC Number of potential contacts NRC Number of real (force-bearing) contacts (5.2) NRG Number of real (force-bearing) grains (5.2) NS Number of sliding contacts (5.1) NST Number of spanning trees (6.8) NV C Number of virtual contacts (6.11) ns Number of sliding contacts per grain (5.1) P Pressure (6.44) q log(NST )/NRC (6.8) S[P] Entropy functional (5.10) T Shear part of σ̄ (6.44) T Topology of contact network (5.7) X Compactivity (5.4) Z Partition function (5.7) Zi Isostatic cluster partition function (6.14) ZMF Mean-field partition function (6.13) z̄ Contact number (5.2) Continued on next page xiii Table 2 – continued from previous page Symbol Definition Equation Reference ziso Isostatic contact number, d+ 1 (1.7) α Isotropic part of α̂ (6.46) αc α̂ : `c`c (6.26) α̂ Inverse of angoricity (5.4) β 〈αc(f̄ cN )0〉 (6.36) γc Weight (6.27) ∆z̄ z̄ − ziso δ[f |V C ] Virtual contact constraint δ-function (6.8) δf f − f̄ δρ ρ− ρ̄ δψ ψ − ψ̄ η Shear part of α̂ (6.46) λ Contact potential (5.5) ν Principal axis angle of α̂ (6.46) ρ̄ Mean-field loop force ρi Cluster loop force (6.22) ρ` Loop force on loop ` (6.7) ρt 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] κ ≡ ( E P )2/3 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 = UL/ν, 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 τ = µeP 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)2Iα with α ≈ −1 and C0 a dimensionless constant [MiD01]. We find Re = UD ν ≈ C0 I 2+α/2 µ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 func- tions δ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 1We 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 computa- tional work, mostly performed in engineering departments, has led to sophisti- cated 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 gen- eration of physicists began performing controlled experiments on granular mate- rials [MOB96, VHBC99, GHL+01, GLBH01, GRCB03]. In particular, several experiments measured the normal stress under a sand pile, with the surpris- ing 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 2We 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]. rg rc Figure 1.1: Sketch of 2D packing of rigid grains, showing the position of a contact, rc, and the position of a grain, rg. 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: F g = − ∑ c∈Cg f cg , 0 = ∑ c∈Cg (rc − rg)× f cg . (1.5) 4 Here f cg is the contact force exerted on grain g at contact c, r c is the position of c, rg 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, respec- tively. For frictional grains we have dNRC −NS degrees-of-freedom in the contact forces. The number of linearly-independent constraints must not exceed the num- ber 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̄ ≥ 1dns + 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. 3z̄ is also called the mean coordination number. 4Friction can also be modelled by geometrical asperities on otherwise spherical particles, in- teracting with normal, repulsive forces, without Coulomb friction. Such packings are found to be exactly isostatic, obviating the need to consider sliding contacts [POS12]. 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 ∑ g∈G ∑ c∈Cg (rc − rg)f cg , (1.8) 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 coarse- graining 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 fric- tionless5 packings of spheres become exactly isostatic in the rigid limit κ → ∞ [Mou98, Rou00, MTDT02, OSLN03]. Numerical studies have shown that fric- tional 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 5For 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 nu- merical simulations of frictionless ellipsoids violate the rigidity inequality anal- ogous 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 = ncnc · f c at each contact, with nc the contact normal. In the spherical limit, contact normals are collinear with con- tact 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 de- pendence 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 count- ing is modified to account for these dependencies, then these packings can be con- sidered 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 illus- trative 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 sta- bility. When only statics are considered, this definition is tautological; it is only well-defined when dynamics are considered. Operationally, one usually defines rat- tlers 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, ∆δ δ = αLD∆T D/κ , (1.10) where αL is the coefficient of linear thermal expansion; αL ≈ 8 × 10−6K−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, 6For frictionless spheres, one can make sharper statements. Hilbert showed that in d dimen- sions, 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 gran- ular 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, tri- angles, 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 Cg 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. Tensors are denoted with a hat. The identity tensor is written δ̂. Outer (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 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 indices; e.g., the kmth component of Â : ∇(B̂ (∇· Ĉ)) is ∑i,j,lAij∂i(Bjk(∂lClm)). The 3D cross-product is defined with the Levi-Civita symbol. We take the free index to be the first, e.g., (Â×B̂)ijm = ∑ k,lAikεjklBlm. To conform to convention 9 we make an exception to the above rules for the tensor curl in 3D, letting (∇ × B̂)ij = ∑ k εikl∂kBjl. We make repeated use of the 2D cross product (Â × B̂)il = ∑ j,k AijεjkBkl, with ε̂ the 2D Levi-Civita symbol, ε12 = −ε21 = 1, ε11 = ε22 = 0. We frequently use the properties ε̂T = −ε̂ and ε̂ · ε̂T = δ̂. The 2D curl is (∇×u)ij = ∑ k εik∂kuj , i.e., the 2D curl only differs from the gradient by multiplication by ε̂: ∇× = ε̂ ·∇. In both 2D and 3D, the condition that a 2-tensor Â be symmetric can be written Â = ÂT ⇐⇒ ε̂ : Â = 0. (1.11) Note that the trace of a 2-tensor can be written tr(Â) = δ̂ : Â. (1.12) Contours are always oriented anticlockwise, and we use +/− to indicate geo- metric 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 amor- phous 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 tri- angulation. 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: 0 = ∑ c∈Cg f cg , 0 = ∑ c∈Cg (rc − rg)× f cg . (2.1) As above, f cg is the contact force exerted on grain g at contact c, r c is the position of c, and rg 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 V ∑ g∈G ∑ c∈Cg (rc − rg)f cg , (2.2) where V is the volume of the packing. In what follows we will decompose this into contributions from each grain g, σ̂g ≡ 1/V g ∑c∈Cg(rc − rg)f cg , where V g 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∂kjl∂lψ. In component form, σ̂ = ( ∂yyψ −∂xyψ −∂xyψ ∂xxψ ) . In three dimensions (3D), (2.4) also holds, with ψ replaced by a tensor ψ̂, known 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 1We 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 contacts 2, 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) 2 NRG, (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.6) 2These 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/z̄L = 1− 2/z̄ + 1/NRC , or 1 z̄ + 1 z̄L = 1 2 + 1 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) 2 NRG − d, (2.8) indicates that phase space is spanned by a vector field defined on the loops, pro- viding dNL degrees of freedom, subject to torque balance, providing 1 2d(d−1)NRG constraints, and a single extra vector constraint. To exhibit this formulation ex- plicitly, 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 f cg = ∑ `∈Lc ρ`g, (2.9) 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. 4Euler’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. 5In 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) 2 NL − d(d− 1) 2 NRC − d(d+ 1) 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 rc × f cg = ∑ `∈Lc ( ϕ`g + r ` × ρ`g ) , (2.11) 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 in- ferred from ρ` must equal those determined from ϕ`, which requires 0 = ∑ `∈Lc ( ϕ`g + ( r` − rc)× ρ`g) . (2.12) 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 topo- logical divergence and curl operators. The former, defined for vector fields on the contacts, is (d∗f)g ≡ ∑c∈Cg f cg , while the latter, defined for pseudovector and pseudoscalar fields on the loops, is (d∗ρ)cg ≡ ∑ `∈Lc ρ ` g. These operators sat- isfy 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é’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 accom- plished 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. 18 19 20 21 22 23 24 25 21.5 22 22.5 23 23.5 24 24.5 25 25.5 26 26.5 t `cg c stg g 18 19 20 21 22 23 24 25 21.5 22 22.5 23 23.5 24 24.5 25 25.5 26 26.5 g scg 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 6Also known as the Dirichlet tessellation 7This 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 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 dis- crete calculus on these graphs in which (1.5) appear as differential relations9. By Poincaré’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: ∫ g dA ∇ · σ̂ ≡ Ag (∇ · σ̂)g ≡ − ∑ c∈Cg scg × σ̂c ≡ − ∮ ∂g dr × σ̂. 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)`cf 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. 8In 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. 9Elements 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∫ g dA ∇× ρ ≡ Ag (∇× ρ)g ≡ − ∑ t∈T g stg ρ t ≡ − ∮ ∂g dr ρ. (2.13) Noting that 2stg = ` c+ g − `c − g , we see that σ̂ g = (∇× ρ)g if f cg = ρ 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 Cg, 0 = f cg = ρ 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]∫ Ω dA σ̂ ≡ Aσ̂ = − ∑ t∈BT stρt ≡ − ∮ ∂Ω 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 `, N `V G − N `V C + N `T = 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∫ g dA ∇ · ρ ≡ Ag (∇ · ρ)g ≡ − ∑ t∈T g stg × ρt ≡ − ∮ ∂g dr × ρ. Comparing this equation with (2.13), and using σ̂g = (∇× ρ)g, we see by inspec- tion that σ̂g = (σ̂g)T is equivalent to (∇ ·ρ)g = 0. This motivates a search for ψc such that ρt = ∇× ψc. We define∫ t dA ∇× ψ ≡ At (∇× ψ)t ≡ − ∑ c∈Ct `ct ψ c ≡ − ∮ ∂t dr ψ. 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, 3N `T = 2N ` V C + N ` RC . 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 · rc + 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 ρtn − ρt0 = n∑ i=1 f ci = n∑ i=1 sci × σ̂ci ≡ ∫ tn t0 dr × σ̂, (2.15) where the forces are exerted from the left side of the contour to the right side. Similarly, ψcn − ψc0 = − n−1∑ i=0 sti × ρti ≡ − ∫ cn c0 dr × ρ. (2.16) Force and torque balance ensure that the line integrals in (2.15) and (2.16), re- spectively, 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` − rc)× ρ`, where z`r` ≡∑c∈RC` rc. Comparing with (3.3) we see that ϕ` = ψ̄` and the consistency constraints in the intrinsic formulation are equivalent to ψc+` = ψ c −`. In the next chapter, we will exploit this observation. With the Voronoi and Delaunay constructions, it is possible to build an ex- trinsic formulation in 3D as well. This is presented in Appendix B. 2.3 Discussion and extensions 2.3.1 Physical interpretation As gauge variables, ρ and ψ do not admit a unique physical interpretation. How- ever, in 2D, the inversion formulae suggest a macroscopic interpretation of po- tential 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 f cg can be decomposed as f cg = ( 1 A`+ + 1 A`− ) `cgψ c + ∑ c′∈Cc 1 A`(c′) `c ′ ψc ′ , (2.17) where Cc 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 force- and torque-balance. c ρ̄g ρ̄gext ρ` − ρ` + 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 = 1 2(ρ `− − ρ̄g)×(ρ`+−ρ̄g). For interior contacts, adding the contributions from each adjacent grain, we have ac ≡ acg+ + acg− = 12f cg− × (ρ̄g + − ρ̄g−). (2.18) The total area of the tiling is A = ∑ c∈RC ac − ∑ c∈EC acgext , (2.19) 21 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 = σ̂ × rg, 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∑ c∈RC ac = 12(ε̂ · σ̂ · ε̂) : ∑ c∈RC f cg` c g (2.20) = −12(ε̂ · σ̂ · ε̂) : Aσ̂T (2.21) = Adet(σ̂), (2.22) so that A = Adet(σ̂) + 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 ε̂ ·rg, where P = 1/2 tr(σ̂) is the pressure, then we find ac = −det(σ̂)/(2P ) f cg− · `cg− . In this representation, repulsion of the contact forces is equivalent to positivity of the areas. Moreover, ∑ c∈RC a c = Adet(σ̂) 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 par- ticular 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 conve- nient 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 σ̂ = ∇× ρ. How- ever, ψ 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, inter- acting without cohesion. In fact, these assumptions are inessential. We considered hard disks and spheres, and framed this chapter as a discus- sion 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 arbi- trary 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) The same relations hold in 3D with ρ→ ρ̂ and ψ → ψ̂. To obtain physical force 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 limit, σ̂ will satisfy σ̂ = ∇×∇× ψ, with ψ → ψ̂ in 3D. From this expression we see that the pressure P ≡ 1 d tr(σ̂) = { 1 2∇2ψ in 2D 1 3∇2tr(ψ̂)− 13∇ · (∇ · ψ̂) in 3D. (2.28) 24 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 micro- scopic 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 geom- etry. 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 bound- ary 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 solu- tion 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 f cg = ∑ `∈Lc ρ`g, (3.1) rc × f cg = ∑ `∈Lc ( ϕ`g + r ` × ρ`g ) . (3.2) 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 0 = ∑ `∈Lc ( ϕ`g − tc` × ρ`g ) , (3.3) 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). ` `c g sc 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 anti- clockwise around loops `. In its dual, the vectors scg circulate anticlockwise around grains g. Summing around a grain g, we find 0 = ∑ c∈Cg rc × ∑ `∈Lc ρ`g, (3.4) 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 vc, we find 0 = ∑ c∈C` vc (∑ `′∈Lc ϕ` ′ g − tc`′ × ρ` ′ g ) . (3.5) We will specify vc 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. c g `′` g′ `c s`g tc` Figure 3.2: Notation convention for `c, tc`, and s ` g vectors. We also have s c = r` ′ − r` = tc` − tc`′ . In the plane, we can orient all loops anticlockwise; (3.3) becomes 0 = ϕ` ′ − ϕ` − tc`′ × ρ` ′ + tc` × ρ`. (3.6) 3.1 Geometric formulation Let us first establish the relationship between ρ and ϕ. We introduce auxiliary variables ψc` ≡ ϕ` − tc` × ρ`, (3.7) 1This is closely related to MacLane’s planarity criterion, a result in graph theory [Bol79, Mac37, O’N73] 2In 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` = rc−r`, the definition (3.7) says that for each loop, we create a facet of a surface which passes through (rc, ψ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∈C` ψc` (3.8) provided r` = 1 z` ∑ c∈C` rc, (3.9) 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 ρ` = − ∑ c∈C` `ctc` · ε̂ −1 ·∑ c∈C` `cψc` = (ĝ`)−1 · −1 A` ∮ ∂` dr ψ = (ĝ`)−1 · (∇× ψ)` (3.10) where ĝ` = 1 A` ∑ c∈C` `ctc` · ε̂ = 1 A` ∑ c∈C` `crc · ε̂ (3.11) 29 is a fabric tensor, the first of several that we will encounter. To understand the behaviour of ĝ`, we use the remarkable geometric identity A`δ̂ = ∫ ` dA ∇r = ε̂ · ∮ ∂` dr r = ε̂ · ∑ c∈C` `c 12(r g′ + rg), (3.12) 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 rc = 1 2(r g′+rg) for each c, then ε̂ · ĝ` · ε̂T = δ̂ and hence ĝ` = δ̂. Defining grain centres by centre-of-mass, the condition that contacts lie midway between adjacent grain centres is always satisfied by monodisperse disks. Hence ĝ` − δ̂ 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 ĝ` 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 m th 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 simu- lations 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 ĝ` converges to δ̂ as averaging size increases, as anticipated above. 3Details of the numerics are discussed in Appendix D. 30 Note that with discrete calculus, ĝ` can be written as ĝ` = ∮ ∂` dr r · ε̂ = (∇× r)` · ε̂T = ε̂ · (∇r)` · ε̂T . (3.14) Since ∇r = δ̂ in the continuum, convergence of 〈ĝ`〉 to δ̂ is consistent with this ex- pression; indeed, it could have been anticipated without reference to the powerful identity (3.12). 31 30 40 50 60 70 −10 0 10 20 30 40 50 60 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 (rc, ψc). Each hole in the surface corresponds to a grain. 32 −5 0 5 10 15 20 25 30 35 30 32 34 36 38 40 42 44 46 48 50 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 0 10 20 30 40 50 60 0 10 20 30 40 50 Figure 3.5: Frame of a highly polydisperse packing undergoing simple shear. Line thickness is proportional to the normal force. 34 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 m g ij Figure 3.6: Convergence of fabric tensor ĝ to δ̂. From simulations of highly polydisperse disks, we compute averages of the fabric tensor 〈ĝ`〉 over neighbour- hoods 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 ĝ (gxx (red,x), gxy (green,o), gyx (blue,+), gyy (black,)). We see that, on average, ĝ 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 0 = ∑ `∈Lg s`g × ρ` ≡ ∮ ∂g dr × ρ ≡ −Ag(∇ · ρ)g, (3.15) where s`g = r c+ − rc− circulates around the grain (Figure 3). This already has a natural continuum meaning as (∇ · ρ)g = 0. The loop equations become 0 = ∑ c∈C` vc ( ϕ` ′ − ϕ` ) − ∑ c∈C` vc ( tc`′ × ρ` ′ − tc` × ρ` ) , (3.16) which will take some work to be put into a continuum form. We begin with the ϕ term. There is a special choice of vc that makes A`(∆ϕ)` ≡ ∑ c∈C` vc ( ϕ` ′ − ϕ` ) 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 ∆ϕ ≡ ∑ `∈L A`(∆ϕ)` = ∑ c∈∂Ω vc ( ϕ` ′ − ϕ` ) , 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 vc so that (3.17) holds identically. A natural definition of ∇ϕ on the contacts is (∇ϕ)c = 1 Ac ε̂ · `c ( ϕ` ′ − ϕ` ) , (3.18) 4For 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 theo- rem5. Indeed, for a contour of adjacent loops (`0, `1, . . . , `n), separated by contacts (c1, c2, . . . , cn), we find ϕ`n − ϕ`0 = n∑ i=1 sci · (∇ϕ)ci ≡ ∫ `n `0 dr · (∇ϕ). (3.19) With this definition,∮ ∂Ω dr ×∇ϕ ≡ ∑ c∈∂C `c × (∇ϕ)c = − ∑ c∈∂C |`c|2 Ac ( ϕ` ′ − ϕ` ) . To enforce (3.17), we should take vc = |`c|2/Ac. Similarly, we can write∑ c∈C` vc ( tc`′ × ρ` ′ − tc` × ρ` ) = ∑ c∈C` vc ( rc × (ρ`′ − ρ`)− (r`′ × ρ`′ − r` × ρ`)) = ∑ c∈C` `c × ((∇ρ)c × rc + (∇(r × ρ))c) = ∮ ∂` dr × ((∇ρ)c × rc + (∇(r × ρ))c) = −A`∇ · ((∇ρ)c × rc + (∇(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) and (3.22). Indeed, at each contact, |f cT | ≤ µf cN can be written as the pair of inequalities (δ̂ ± 1µ ε̂) : ncf c ≤ 0, (3.24) 5Note 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 ĝ` 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 interpreta- tion, 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 depo- sition 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 grain- scale 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, 〈ĝ`〉 = δ̂. (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 imple- mented. However, it is analytically intractable. Hence to proceed theoretically, we will need to make approximations on some averages. Recalling that ρ` = (ĝ`)−1 · (∇×ψ)`, this is already in a discrete calculus form ρ = 〈(ĝ`)−1 · (∇× ψ)`〉. (3.29) Since fluctuations of ĝ` are small, it is natural to consider a mean-field approxi- mation ρ = 〈ĝ`〉−1 · 〈(∇× ψ)`〉 = ∇× ψ. Likewise, since ϕ` = 1 z` ∑ 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 〈ϕ`〉 = ψ̄(r`) in terms of this smooth surface. But positivity of pressure implies that the 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 ψ̄(r`), and hence ϕ 6= ψ in the continuum. Discrete calculus allows us to obtain self-consistent corrections to the mean- field assumption ϕ = ψ̄. Since the coarse-grained surface ψ̄(r) is assumed smooth, we can Taylor expand ψ̄(rc) around ψ̄(r`), and compute ϕ` = 1 z` ∑ c∈C` ψ̄(r c) exactly, introducing fabric tensors which characterize the local geometry. Here we will fit ψ̄(r) to a quadratic polynomial; this is physically equivalent to ho- mogeneously distributing stress across a loop. Later we will discuss the effect of higher-order terms6 . We have ψ̄(r) = ψ̄(r`) + (r − r`) · ∇ψ̄(r`) + 12(r − r`)2 : ∇∇ψ̄(r`), (3.30) 6For a typical loop, a quadratic polynomial is the lowest order polynomial that fits through all ψc. To see this, consider an mth degree polynomial ψ̄(r) = ∑m n=0 ân :: n· · · :: (r − r0)n, in which (r − r0)n is the n-fold outer product of r − r0 and :: n· · · :: is shorthand for n contractions. Each coefficient ân is a symmetric n th-order tensor, with n+1 degrees-of-freedom, so the polynomial has j(m) = 1 2 (m+1)(m+2) degrees-of-freedom. Within a loop, the requirement that ψ̄ goes through all (rc, ψc) gives z` constraints. The identity 1 2 = 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. 39 and hence ϕ` = ψ̄(r`) + 1 2z` ∑ c∈C` tc`t c ` : ∇∇ψ̄(r`) = ψ̄(r`) + 12 F̂ ` : ∇∇ψ̄(r`), (3.31) defining a fabric tensor F̂ ` = 1 z` ∑ c∈C` tc`t c `. (3.32) Note that F̂ can be rewritten F̂ ` = 1 z` ∑ c∈C` tc`r c = 1 z` ∑ c∈C` rcrc − z`r`r` . (3.33) 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. We may likewise compare ρ̄(r`) with ρ`. In Appendix F we show that up to 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 ρ = ∇× ψ ϕ = ψ + 12 F̂ : ∇∇ψ. (3.35) Hence 0 = ∇ · ρ is identically satisfied. However, the loop equation gives the leading-order stress-geometry equation 0 = ∆ϕ−∆ψ = 12∆ ( F̂ : ∇∇ψ ) (3.36) 40 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, 12∆ϕ, must equal the pressure in the grains, 1 2∆ψ; 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 ∑ `∈Ω ∑ c∈C` tc`t c `. (3.37) 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 from the packing displayed in Figure 3.5. We notice that F̂ is close to δ̂, but has a systematic deviation which expresses fabric anisotropy. We also notice that fluctuations are much larger than in ĝ, since there is no exact cancellation on averaging. 0 1 2 3 4 −0.1 0 0.1 0.2 0.3 0.4 m F i j 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 12〈|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̂ ′ = 1 NRC ∑ c∈RC `c`c. Since tc`/|tc`| = ε̂ · `c/|`c| up to fluctuations, over a sufficiently large area we have F̂ ∝ ε̂ · F̂ ′ · ε̂T , (3.39) with the constant of proportionality reflecting different normalizations. This im- plies 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 ψ ∼ L2P, F̂ ∼ ξ2, ∇ ∼ 1 L . (3.40) Consider the Taylor expansion ψ̄(rc) = ψ̄(r`) + ∑ n≥1 1 n! (tc`) n :: n· · · :: ∇nψ̄(r`), (3.41) with (tc`) n the n-fold outer product of (tc`), and :: n· · · :: shorthand for n contrac- tions. 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 δnL2P . In the macroscopic limit δ → 0, we need only keep leading order terms in δ. Noting that ∆ψ ∼ P, ∆(F̂ : ∇∇ψ) ∼ δ2P, (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 σ̂1 = g · rδ̂. The Coulomb condition P ≥ 0 applies to the total stress, hence if σ̂ = σ̂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 The simplest geometry is an isotropic and homogeneous fabric, F̂ ∝ δ̂. In this 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 solu- tions 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. 7See 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̄ = x− iy, ∂z = 1 2(∂x − i∂y), ∂z̄ = 12(∂x + 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)] = 12(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 = 2piA δ(x), (3.49) where the vector A indicates the direction of forcing, and the factor of 2pi is there for later convenience. This amounts to 4 scalar equations: on y = 0, x 6= 0, we have σ̂ · y = 0, and the traction on an infinitesimal curve around the origin is 2piA. Since σ̂ = ( ∂yyψ −∂xyψ −∂xyψ ∂xxψ ) these are equivalent to ∂2xψ ∣∣ y=0 = +2piAyδ(x), ∂x∂yψ|y=0 = −2piAxδ(x). (3.50) For these boundary conditions, the potentials take the form f(z) = −A log z, dg(z) dz = Ā log z, (3.51) where A = Ax + iAy, Ā = 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 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 2 4 6 8 10 12 14 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 2 4 6 8 10 12 14 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 −10 −5 0 5 10 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0 5 10 15 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∂ 2 x∆ψ + Fy∂ 2 y∆ψ = 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 elas- ticity [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 perpendic- ular 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 1− a ( −b± √ 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 differ- ences 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 re- lated 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) dz = A1 log z, dg(w) dw = A2 logw, (3.58) with A1 = Ax + µAy pi(1 + iµ) , A2 = −Ax + iAy pi(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. (3.61) 47 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0 5 10 15 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0 5 10 15 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 −10 −5 0 5 10 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 −20 −10 0 10 20 Figure 3.9: Stresses resulting from a normal point forcing at the origin of a semi- infinite half plane, with anisotropic fabric. (σyy (top left), σxx (top right), σxy (bottom left), P (bottom right)) −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 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 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 −5 0 5 10 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 −15 −10 −5 0 5 10 15 20 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0 5 10 15 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 −20 −10 0 10 20 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 pi/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 corre- lation lengths decrease away from isostaticity, and hence mean-field assumptions should become better away from isostaticity. Of course, it is possible that corre- lations 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〉 ≈ A piNL = 2Ag pi(z̄ − 2)φ(1− x0) , (3.62) where have used NLpiξ 2 ≈ A and A = AgNRG/((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 vari- ables are subject to geometric compatibility conditions, expressing the grain-scale relationships between grain movement, grain rotation, and contact movement. Af- ter 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 vg 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 rc, on grains g′(c) and g(c). Under deformation, the contact point c on grain g moves at a rate wcg = v g − ωgε̂ · tcg, (4.1) where tcg = r c − rg. The relative velocity at a contact is V c = wcg′ −wcg = vg ′ − ωg′ ε̂ · tcg′ − vg + ω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 vg = 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 de- formation is a combination of rolling and rotation of the contact normal nc. 51 When V cN = n c · V c < 0, grains g and g′ are moving closer together; when V cN > 0, the grains are separating. Finally, sliding deformation is described by V cT = −nc × V c 6= 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 Ac ε̂ · scV c, (4.3) which has units of strain rate1. Following Kruyt, we will call Γ̂ the Cosserat strain 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 g0 dr · Γ̂ = − n∑ i=1 `ci × 1 Ac sciV ci = n∑ i=1 V ci , (4.4) which should be compared with∫ `n `0 dr × σ̂ = − n∑ i=1 sc ′ i × 1 Ac ′ i `c ′ if c ′ i = n∑ i=1 f c ′ i , (4.5) for a contour of loops (`0, `1, . . . , `n) separated by contacts (c ′ 1, c ′ 2, . . . , c ′ n). This shows that ε̂T · Γ̂ is to V as σ̂ is to f . When an open contour is taken which shrinks to a point in the continuum, (4.4) becomes n · Γ̂ = ncV , where the unit vector n is directed from g0 to gn, 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 Γ̂ ∼ −δ̂. The significance of Γ̂ is further illustrated by energy considerations. Under the action of external forces F g and torques τ g, the rate at which work is done in a quasistatic deformation is dW dt = ∑ g∈Ω (F g · vg + τ gωg) = − ∑ g∈Ω ∑ c∈Cg ( f cg · vg + (tcg × f cg )ωg ) = ∑ c∈Ω f c · V c (4.6) = − ∑ c∈Ω 1 Ac `cf c : ( ε̂ · scV 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, Rou00, Kru03, Sat04]. It indicates that Γ̂ describes strain which costs energy. 4.2 Constitutive laws Because Γ̂ is directly related to deformation, it is the kinematic variable that is di- rectly related to stress. The relation between Γ̂ and σ̂ 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 V cN + ε̂ · nc V cT (4.8) f c = −nc f cN − ε̂ · nc f cT , (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 either the contact bears a repulsive force f cN > 0 and the normal relative velocity vanishes, V cN = 0, or the normal relative velocity is repulsive, V c N > 0, and the normal force vanishes, f cN = 0. This condition is called the Signorini condition and can be written [RBR96] f cN ≥ 0, V cN ≥ 0, f cNV cN = 0. (4.10) Similarly, Coulomb friction requires that if |f cT | < µf cN , then no sliding can occur, V cT = 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 |f cT | ≤ µf cN , f cTV cT ≥ 0, (|f cT | − µf cN)V cT = 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 f cN = 0 or V c N = 0, and similarly for the tangential components. Nevertheless, 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 = ncnc, (4.12) so that f cN = Ac nc · `c Ĝ c : σ̂c V cN = Ac sc × nc Ĝ c : Γ̂c (4.13) f cT = Ac nc · `c (Ĝ c · ε̂T ) : σ̂c V cT = Ac sc × nc (Ĝ c · ε̂T ) : Γ̂c. (4.14) For convex grains, we have `c · nc > 0, and we expect that sc × nc > 0 so that f cN ≥ 0 ⇐⇒ Ĝc : σ̂c ≥ 0 (4.15) V cN ≥ 0 ⇐⇒ Ĝc : Γ̂c ≥ 0. (4.16) Let M̂± = δ̂ ± 1µ ε̂. Then we also have f cN ± 1µf cT ≥ 0 ⇐⇒ ( Ĝc · M̂± ) : σ̂c ≥ 0. (4.17) 54 Note that these inequalities imply f cN ≥ 0. Since (σ̂c)T · Γ̂c = −f cV c/Ac, the equalities of the Signorini and Coulomb conditions can be written 0 = ( (σ̂c)T · Γ̂c) : Ĝc, 0 = ( (σ̂c)T · Γ̂c) : (M̂+ · Ĝc · ε̂)((σ̂c)T · Γ̂c) : (M̂− · Ĝc · ε̂). The second is nonlinear because of the inherent nonlinearity of |f cT |. Under the simplest mean-field assumptions the above conditions become Ĝ : σ̂ ≥ 0 Ĝ : Γ̂ ≥ 0( Ĝ · M̂± ) : σ̂ ≥ 0, (4.18) and 0 = ( σ̂ · Γ̂) : Ĝ, (4.19) 0 = ( σ̂ · Γ̂) : (M̂+ · Ĝ · ε̂)(σ̂ · Γ̂) : (M̂− · Ĝ · ε̂). These forms are not unique: for example, the discrete equation f cNV c N = 0 could also give a continuum equation (Ĝ : σ̂)(Ĝ : Γ̂) = 0, using (4.13). This equa- tion differs from (4.19) in the mean-field assumptions implicit in writing it. We comment further on this below. If σ̂, Γ̂, and Ĝ vary smoothly in the continuum, then the domain will be 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 0 = ( σ̂ · Γ̂) : (M̂± · Ĝ · ε̂) (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 NRC ∑ c∈RC ncnc = 1 ξ ε̂T · F̂ · ε̂. (4.21) 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 The fields v, ω, w, and Γ̂ are related through discrete equations (4.1) and (4.2); 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 vg and ωg can be recovered from wcg. Consider ∑ c∈Cg qcgw c g = v g (∑ c∈Cg qcg ) − ωgε̂ · ∑ c∈Cg qcgt c g. For convex particles, we can always choose positive coefficients qcg so that the vector ∑ c∈Cg q c gt c g = 0. Then we have vg = 1 qg ∑ c∈Cg qcgw c g, (4.22) with qg = ∑ c∈Cg q c g, so that v is a weighted grain average of w. Similarly,∑ c∈Cg scgw c g = −ωg ∑ c∈Cg scg ε̂ · tcg, so that ωg = 1 2Ag ∑ c∈Cg scg ·wcg = 1 2Ag ∮ ∂g dr ·w = −12tr (∇×w)g, (4.23) where Ag = 12 ∑ c∈Cg t c g × scg is an area associated to a grain. Under the simplest mean-field closures on ω = 〈ωg〉, w = 〈wc〉, and v = 〈vg〉, these become ω = −12tr(∇×w) v = w. (4.24) Note that the relation between ω and w can be rewritten as ωε̂ = −12tr(∇×w)ε̂ = −12 ε̂ε̂ : (∇w)T = 12(∇w − (∇w)T ), (4.25) 56 hence grain rotation is equivalent to macroscopic vorticity in the mean-field ap- proximation2. 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 K = d(d+ 1) 2 NRG − ( d(NRC −NS) + (d− 1)NS ) = d 2 NRG ( ziso + 1 dns − z̄ ) , (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 iso- static 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 infinitesi- mal deformation. This property of isostaticity has been discussed extensively by Combe and Roux [CR00, Rou00, RC02]. 2Here we have used the identity ijkl = δikδjl − δilδjk. 57 Now let us find continuum equations for floppy modes. Initially we will con- sider 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 = wcg′ −wcg, the floppy mode condition V c = 0 is equivalent to continuity of w across a contact. There is a clear analogy between wcg 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 ρ` wcg ψ c ` wcg′ = w c g ψ c `′ = ψ c ` (∇× ω)` = 0 (∇ · ρ)g = 0 (∆v)g = · · · (∆ϕ)` = · · · Table 4.1: Deformation-stress duality Summing V c = 0 around a loop, we find 0 = ε̂ · ∑ g∈G` ( tc + g − tc − g ) ωg = ε̂ · ∮ ∂` dr ω = ε̂T ·A` (∇× ω)` , (4.28) so 0 = ε̂T · (∇× ω)` = (∇ω)`. Summing around grains, with coefficients zc, we find 0 = ∑ c∈Cg zc ( vg ′ − vg)− ε̂ ·∑ c∈Cg zc ( tcg′ω g′ − tcgωg ) . (4.29) 58 As before, we write Ag(∆v)g ≡ − ∮ ∂g dr × (∇v) = ∑ c∈Cg sc × (∇v)c, (4.30) with sc oriented clockwise around the grain, and (∇v)c ≡ − 1 Ac ε̂ · sc(vg′ − vg). (4.31) Green’s theorem indicates that we should take zc = |sc|2/Ac. We have∑ c∈Cg zc ( tcg′ω g′ − tcgωg ) = ∑ c∈Cg rcsc × (∇ω)c − ∑ c∈Cg sc × (∇(rω))c = Ag (∇ · ((∇ω)r))g −Ag (∇ · ∇(ω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 homoge- neous 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 D̂ ≡ 12(∇w + (∇w)T ) = +12(∇u+ (∇u)T ) (4.37) ωε̂ = 12(∇w − (∇w)T ) = −12(∇u− (∇u)T ), (4.38) hence the strain rate D̂ described by u is the same as that described by w, but 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 degrees- of-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, topologi- cally, 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 trian- gles 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 3Mathematically, 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∈C` V c = ∑ c∈C` `c × ( − 1 Ac scV c ) = − ∮ ∂` dr × (ε̂ · Γ̂) = A`(∇ · (ε̂ · Γ̂))`, (4.40) while summing around a grain with coefficients zc = |sc|2/Ac, we find∑ c∈Cg zcV c = ∑ c∈Cg sc × ( − 1 Ac ε̂ · scV c ) = − ∮ ∂g dr × Γ̂ = Ag(∇ · Γ̂)g. (4.41) The terms involving Γ̂ modify (4.34) and (4.33) to (∇ · (ε̂ · Γ̂))` = (∇ω)` (∇ · Γ̂)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 ω = −12tr(∇×w), we find 0 = −∇ · (ε̂ · Γ̂)+∇ω (4.43) 0 = ∇ · (− Γ̂ +∇w − ωε̂). (4.44) Using (4.25) we find ∇w − ωε̂ = 12(∇w + (∇w)T ) = D̂, hence 0 = ∇ · (− Γ̂ + D̂). We see that the grain equation constrains the strain while the loop equation con- strains the vorticity. The latter can be written 0 = ∇ · (−ε̂ · Γ̂ + 12(∇w · ε̂T )T ), which implies Γ̂ = 12(∇w)T − 12∇u (4.45) 61 for a velocity field u. The grain equation becomes 0 = ∇ · (− Γ̂ + 12∇w + 12(∇w)T ) = ∇ · (− Γ̂ + Γ̂T + 12(∇u)T + 12(∇w)T ) = −∇ · (Γ̂− Γ̂T )+ 12∇(∇ · u) + 12∇(∇ ·w) = −∇ · (Γ̂− Γ̂T )+ 12∇(∇ · u) +∇trΓ̂ + 12∇(∇ · u), so that u must solve ∇(∇ · u) = ∇ · (Γ̂− Γ̂T )−∇tr(Γ̂), (4.46) which is the analog of the equation which, for true floppy modes, resembled linear elasticity. Here we see no resemblance for general Γ̂. By similar manipulations, it can be shown that this equation is equivalent to ∆u = ∆w, (4.47) which may be more convenient if w is known. The 6 degrees-of-freedom in Γ̂ and 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 dt F̂ ` = 1 z` d dt ∑ c∈C` tc`t c ` = 1 z` d dt ∑ c∈C` tc`r c = 1 z` d dt ∑ c∈C` rcrc − z`r`r` = 1 z` ∑ c∈C` ( drc dt rc + rc drc dt ) − dr ` dt r` − r`dr ` dt , 62 where we have used the fact that z` is constant during the continuous part of deformation. We need an expression for drc/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, wcg′ and w c g, corresponding to the movement of each grain participating in the contact. The true velocity of the contact c, w̄c, depends on the local curvature of the grains. We can proceed as in the previous Chapter by Taylor expanding the continuum field w̄(rc) about the loop centre. It will be sufficient to expand to linear order. Hence we let drc dt = w̄(rc) = w(r`) + tc` · ∇w(r`) and dr`/dt = w(r`), which implies d dt F̂ ` = (∇w(r`))T · 1 z` ∑ c∈C` tc`r c + 1 z` ∑ c∈C` rctc` · ∇w(r`) = (∇w(r`))T · F̂ ` + F̂ ` · ∇w(r`). 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̂ dt = (∇w)T · F̂ + F̂ · (∇w). 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 NRCF̂ is δF̂− = zz′ z′′ (tc` − tc`′)(tc` − tc`′)− 2z z′′ tc`t c ` − 2z′ z′′ tc`′t c `′ , (4.48) where z = z`, z′ = z`′ , and z′′ = z`′′ . We can obtain a mean-field form by making replacements z, z′ → z̄L, tc` − tc`′ → 2tc` tc`t c ` → F̂− tc`′t c `′ → 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̂− ≈ 2z̄LF̂− (4.49) Under the same assumptions, contact creation will lead to a change in NRCF̂ of −2z̄LF̂+, where F̂+ is the mean fabric which is created by contact closing. These processes modify the fabric evolution equation to d dt (NRCF̂ ) = NRC(∇w)T · F̂ +NRCF̂ · (∇w) + 2z̄L ( J−F̂− − J+F̂+ ) , (4.50) where J± are the rates of contact opening and closing, which must satisfy J+ − J− = dNRC dt . (4.51) The quantities J± and F̂± should be related to the strain rate D̂ by consideration of 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 F̂± = 12ξ 2 δ̂ ∓ 1√ | det D̂| D̂ , (4.52) where ξ2 = 〈|tc`|2〉 = tr(F̂ ). The rates of contact opening and closing are deter- mined by the proximity of neighbouring grains and the strength of velocity fluctu- ations; 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 J± = NRGI± z̄, φ, tr(D̂)√ |det D̂| √| det D̂|, (4.53) 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 = AgN/φ to see that dφ dt = − φ A dA 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 com- plete 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 F̂ and Ĝ, the contact number z̄, the area fraction φ, the contact velocity w, and 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) the stress-geometry equation 0 = ∆ ( F̂ : (2P δ̂ − σ̂)), (4.57) the frictional rigid-contact constitutive laws 0 = ( σ̂ · Γ̂) : Ĝ (4.58) 0 = ( σ̂ · Γ̂) : (M̂± · Ĝ · ε̂), (4.59) the equations of geometric compatibility 0 = ∇ · (ε̂ · Γ̂)−∇ω (4.60) 0 = ∇ · (− Γ̂ +∇w − ωε̂), (4.61) geometry evolution equations dF̂ dt = dF̂ dt [F̂ ,w, z̄], (4.62) dz̄ dt = dz̄ dt [w, z̄, F̂ , φ] (4.63) dφ dt = −φ∇ ·w (4.64) and miscellaneous constitutive laws Ĝ = 1 ξ ε̂T · F̂ · ε̂ (4.65) ω = −12tr(∇×w) (4.66) ξ = √ tr(F̂ ) (4.67) M̂± = δ̂ ± 1µ ε̂. (4.68) 65 Admissible solutions are also constrained by inequalities Ĝ : Γ̂ ≥ 0( Ĝ · M̂± ) : σ̂ ≥ 0. (4.69) These equations can be simplified by writing σ̂ and Γ̂ in terms of potentials ψ and u as σ̂ = (g · r)δ̂ +∇×∇× ψ (4.70) Γ̂ = 12(∇w)T − 12∇u. (4.71) Written in this way, σ̂ and Γ̂ identically solve g = ∇ · σ̂ (4.72) 0 = ε̂ : σ̂ (4.73) 0 = ∇ · (Γ̂− 12(∇w)T − 12∇w) (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 dz̄/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 σ̂ 6= 0, but in the quasistatic limit, we have σ̂ = 0. Thus, when ∆z̄ < 0 the rigid-contact constitutive laws are trivially satisfied, independent of Γ̂. Note also that the Coulomb inequality is saturated. 66 In this Chapter, we counted generalized floppy modes, which are deformations that preserve the steric exclusions. For rigid grains, these are the only deforma- tions possible. We found that the equations of geometric compatibility can be solved if ∆z̄ ≤ 0. If this is satisfied, Γ̂ can be written in terms of w and u, 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 forbidden with rigid grains. Thus we expect Γ̂ = 0 and w = 0, so that the contact 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 me- chanical equilibrium and geometric compatibility; indeed, these are closely re- lated 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, defor- mation 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]. z̄ ns I ziso solid (Γ̂ = 0) gas (σ̂ = Γ̂ = 0) flu id (∆ z̄ = 0) 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. In the fluid phase, both Γ̂ and σ̂ 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. In the solid phase, ∆z̄ > 0, the packing is static and the fabric is fixed: Γ̂ = 0, 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 boundary- value problem violates the Coulomb inequality, then the packing fails along the line ( Ĝ ·M̂± ) : σ̂ = 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 contacts pose no constraints in this case, ns = 0 and Γ̂ = 0. In this phase 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 inequality, then a jamming transition occurs along the curve Ĝ : Γ̂ = 0. In 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 con- tact 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 equa- tions needs to be solved. It is convenient to use Pauli matrices ς̂1 ≡ 0 1 1 0 , ς̂2 ≡ 0 −i i 0 , ς̂3 ≡ 1 0 0 −1 , which satisfy ς̂a · ς̂b = δabδ̂ + iabc ς̂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 (4.76) σ̂ = P (δ̂ + τaς̂a) (4.77) Ĝ = 12 δ̂ + gaς̂a, (4.78) with implied sums over a. Symmetry of σ̂ and Ĝ implies that τ2 = g2 = 0. Substituting these relations into the contact constitutive laws (4.19) and (4.20) gives, after some algebra, 0 = ΓV ( 1 + 2τaga ) + 2iγ2 ( τ1g3 − τ3g1 ) + γ1 ( 2g1 + τ1 ) + γ3 ( 2g3 + τ3 ) (4.79) and 2ΓV ( τ1g3 − τ3g1 ) + γ1 ( 2g3 + τ3 ) + iγ2 ( 1− 2τ1g1 − 2τ3g3 )− γ3(2g1 + τ1) = ± 2 µ [ΓV + γ1τ1 + γ3τ3] . (4.80) Using the identity ε̂ · Â · ε̂ = ÂT − δ̂ tr(Â) and (4.65) we can check that F̂ : (2P δ̂ − σ̂) = F̂ : (ε̂ · σ̂ · ε̂T ) = (ε̂T · F̂ · ε̂) : σ̂ = ξĜ : σ̂, hence the stress-geometry equation is equivalent to 0 = ∆ ( ξĜ : σ̂ ) = 2∆ ( Pξ(12 + g1τ1 + g3τ3) ) . (4.81) 69 We can also write the Signorini and Coulomb inequalities in these variables. They are Ĝ : σ̂ = 2P ( 1 2 + g1τ1 + g3τ3 ) ≥ 0 (4.82) Ĝ : Γ̂ = 2 ( 1 2ΓV + g1γ1 + g3γ3 ) ≥ 0 (4.83)( Ĝ · M̂± ) : σ̂ = 2P ( 1 2 + τ1(g1 ± 1µg3) + τ3(g3 ∓ 1µg1) ) ≥ 0. (4.84) We now want to write the equalities as partial differential equations in u and w. To do so, from (4.45) we simply compute 2ΓV = tr(Γ̂) = 1 2∇ · (w − u) (4.85) 2γ1 = ς̂ T 1 : Γ̂ = 1 2∂y(wx − ux) + 12∂x(wy − uy) (4.86) 2γ2 = ς̂ T 2 : Γ̂ = i 2∂y(wx + ux)− i2∂x(wy + uy) (4.87) 2γ3 = ς̂ T 3 : Γ̂ = 1 2∂x(wx − ux)− 12∂y(wy − uy), (4.88) 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 σ̂, Ĝ, ∇w and ∇u constant. We will assume that Ĝ is fixed and look for solutions as Ĝ is varied. Note that for homogeneous stress and fabric, the stress-geometry equation is identically satisfied. 4.6.1 Simple shear flow A natural flow problem is homogeneous shear flow w = (w0y, 0). It is easy to see that for any constant Γ̂, we can find a function u = (ax+by, cx+dy) that satisfies ∆(w − u) = 0 and Γ̂ = 12(∇w)T − 12u. 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 Γ̂. Boundary conditions on slip velocity V c can be applied to Γ̂ by considering n · Γ̂ = ncV , which relates the mean relative velocity V across contacts oriented with contact normal n to Γ̂, with nc the contact density in the direction n. Under 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 = 1 2V0(1− β), and iγ2 = 12V0(1 + β). In steady simple shear flow, 70 it is observed that the fabric has principal axes at pi/4 to the shear direction, so we let g3 = 0 [OK74]. The solution to the constitutive laws is then τ1 = 4(β2 + 1)µg1 −µ+ 2µβ − µβ2 + 4g21µ+ 8µg21β + 4µβ2g21 + 4g1 − 4β2g1 , (4.89) τ3 = − µβ 2 − 4β2g1 − µ− 4g21µ+ 8βg1 + 4µβ2g21 − 4g1 −µ+ 2µβ − µβ2 + 4g21µ+ 8µg21β + 4µβ2g21 + 4g1 − 4β2g1 . (4.90) Symmetry requires that the principal stress direction is at pi/4 to the shear direc- tion, 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 = (vxx,−vyy) [AR10, Kru10]. As above, we can always satisfy the geomet- ric compatibility equations relating w to Γ̂, so that the latter must be constrained separately. By reflection symmetry of pure shear flow, Γ̂ must take the form Γ̂ = V0 1 0 0 −β , (4.92) where V0 and β are constants. Again, by symmetry, we must have g1 = τ1 = 0. The constitutive laws then become 0 = (1− β)(1 + 2τ3g3) + (1 + β)(2g3 + τ3) (4.93) 0 = (1− β) + (1 + β)τ3. (4.94) 71 These imply τ3 = (β−1)/(β+1) and 0 = g3(1−τ23 ). The latter equation indicates that either the fabric is isotropic, g3 = 0, or one of the normal stresses vanishes, τ23 = 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 = τ1g1 + τ2g3, τ1g3 = τ3g1 ± 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 µ2 ≥ 0, 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, and becausew can always be found which is compatible with the given Γ̂, they offer 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. As discussed above, the relation between the geometric variables Γ̂, w, F̂ , Ĝ, z̄, and ns is incomplete. Further investigation into the contact opening and 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 = NRGz̄/2 force-bearing contacts. Force and torque balance on each grain constrain contact forces to lie in a space of dimension H ′ = dNRC − d(d+ 1) 2 NRG = d 2 NRG(z̄ − ziso), (5.1) with ziso = d + 1. Sliding constraints further reduce the dimension of available contact force solutions to H = dNRC −NS − d(d+ 1) 2 NRG = d 2 NRG(z̄ − 1dns − ziso), (5.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 1dns+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 fluctua- tions 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 dissi- pates energy through friction and inelastic collisions and forms a static packing in metastable mechanical equilibrium. Hence the dynamics can be viewed as a pro- cess 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 ∼ µPD3 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 1Near 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. 2In 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 Vc 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 = {f c1} and f2 = {f c2}. 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, Vc is convex, hence each f c lies in these convex sets. The field f thus lies in the intersection of all the Vc, 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 probabil- ity 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 increas- ing 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 in- duced. 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 hy- pothesis 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. 3Of course, β = 1/(kBT ) in terms of Boltzmann’s constant kB and temperature T . 76 Hence the microcanonical Edwards probability density is Pmc(r,ω,f) = 1 Zmc δ ( V − V [r,ω])δ(σ̂ − σ̂[r,f ])Θ[r,ω,f ], (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 〈A〉 = ∫ Dr ∫ Dω ∫ Df A[r,ω,f ] Pmc(r,ω,f), where Dr is shorthand for integration over all rg. Inserting 1 = ∫ dσ̄δ(σ̄ − σ̂) and 1 = ∫ dV̄ δ(V̄ − V ) into this expression, we find 〈A〉 = 1 Zmc ∫ dσ̄ ∫ dV̄ A(V̄ , σ̄) eS(V̄ ,σ̄), with S(V̄ , σ̄) = log [∫ Dr ∫ Dω ∫ Df δ(σ̄ − σ̂[r,f ])δ(V̄ − V [r,ω])Θ[r,ω,f ] ] the microcanonical entropy. As in usual Hamiltonian mechanics, the microcanon- ical 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 macro- scopic observables V and σ̂ are replaced by fixed temperature-like variables com- pactivity X = (∂S/∂V )−1 and angoricity Â = V (∂S/∂σ̂)−1 [Edw05]. Note that angoricity is a tensor. Given the validity of a microcanonical ensemble, the jus- tification for the canonical ensemble is standard and discussed in, for example, [BDD06, HC09]. The canonical Edwards probability density is Pc(r,ω,f) = 1 Zc e−V [r,ω]/X−V α̂:σ̂[r,f ] Θ[r,ω,f ], (5.4) where α̂ = Â−1. 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 mechan- ically 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 simula- tions [ABOS12, GBO06, GBOS09]. These authors found that microstates were accessed with a frequency that was far from uniform, suggesting that (5.3) is in- correct: 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 justi- fication, 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 4However, 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 Z λNRCe−V α̂:σ̂[f ] Θ[T ,f ], (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, Θ[T ,f ] = ∏ g [ δ (∑ c∈Cg f cg ) δ (∑ c∈Cg rc × f cg )]∏ c Θ ( f cN − 1µ |f cT | ) ≡ δFB[f ] δTB[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 Z = ∑ T λNRC ∫ Df e−V α̂:σ̂[f ] Θ[T ,f ], (5.7) 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 Properties of the enlarged Force Network Ensemble 5.3.1 Contact potential It is clear from the expression for Z that λ is necessary for dimensional homogene- ity. By defining an elementary phase space volume, it plays a role very similar to ~ in classical statistical mechanics. Indeed, we can write √ λ 3NRG Z = ∑ T ∫ D( √ λf) e −V α̂√ λ :σ̂[ √ λf ] Θ[ √ λf ], (5.8) where all forces in the right-hand side now appear in the dimensionless combina- tion √ λ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(α̂, λ) = 1√ λ 3NRG Z(α̂/ √ λ, 1). (5.9) 5.3.2 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 restricted to O5, for which 〈σ̂〉P = σ̄ and 〈NRC〉P = N̄RC , we have 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 V ∂F ∂α̂ , α̂ = 1 V ∂S ∂〈σ̂〉 〈NRC〉 = −λ∂F ∂λ , − log λ = ∂S ∂〈NRC〉 (5.14) where S = S[Peq], and 〈 〉 denotes an ensemble average taken with respect to Peq. 5Meaning 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- sembles. Inserting 1 = ∫ dσ̄δ(σ̂ − σ̄) and 1 = ∑N̄RC δN̄RC ,NRC(T ) into (5.7), we find Z = ∑ N̄RC λN̄RC ∫ dσ̄e−V α̂:σ̄ ∑ T δN̄RC ,NRC(T ) ∫ Dfδ(σ̂[f ]− σ̄)Θ[f ] = ∑ N̄RC λN̄RC ∫ dσ̄e−V α̂:σ̄eSmc(σ̄,N̄RC), (5.15) where Smc(σ̄, N̄RC) = log [∑ T δN̄RC ,NRC(T ) ∫ Dfδ(σ̂[f ]− σ̄)Θ[f ] ] (5.16) is the entropy of the micro-canonical eFNE corresponding to fixed stress σ̄ and fixed contact number N̄RC . In the thermodynamic limit, (5.15) is dominated by its saddle-point contribution (σ̂0, N0), satisfying α̂ = 1 V ∂Smc ∂σ̄ ∣∣∣∣ σ̄=σ̂0,N̄RC=N0 (5.17) and − log λ = ∂Smc ∂N̄RC ∣∣∣∣ σ̄=σ̂0,N̄RC=N0 (5.18) In this limit, F = − logZ = −N0 log λ+ α̂ : V σ̂0 − Smc(σ̂0, N0) +O(logS′′mc), (5.19) hence σ̂0 = 1 V ∂F ∂α̂ = 〈σ̂〉, N0 = − ∂F ∂ log λ = 〈NRC〉 (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 space on the shell σ̂[f ] = σ̄, NRC(T ) = N̄RC . Hence, computing Z is equivalent to knowing the geometry of O(T ) as σ̄ and N̄RC are varied. Before presenting a 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 SFNE(σ̄) = log [∫ Dfδ(σ̂[f ]− σ̄)Θ[f ] ] , (5.22) the entropy of the Force Network Ensemble on a strictly fixed contact network. Making a change of coordinates in (5.22) to σ̄ = kσ̂′, it is easy to show that SFNE(kσ̂ ′) = (12NRGd∆z − d2) log(k) + SFNE(σ̂′), (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 ZFNE(α, η̂) = α −NRGd∆z/2 ∫ dσ̄e−(δ̂+η̂):V σ̄eSFNE(σ̄), (5.25) hence 〈P 〉 = 1 d Tr〈σ̂〉 = 1 dV ( ∂FFNE ∂α + 1 α ∑ i ∂FFNE ∂ηii ) . (5.26) The tensor ∂FFNE/∂η̂ is function of η̂. On an isotropic geometry, we can only form tensors from δ̂ and η̂, hence its most general form is ∂FFNE/∂η̂ = f(I1, I2, . . . , Id) δ̂ + g(I1, I2, . . . , Id) η̂, (5.27) where In is the n th invariant of the tensor η̂, and f and g are unknown functions. Under reflection symmetry r → −r, we have F → F , In → (−1)nIn, 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 ∂FFNE/∂η̂|η̂=0 = 0. 82 Using F = − logZ 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 ma- terial discussed in this Chapter derive from the work of Edwards, and are not closely connected to models commonly appearing in the statistical mechanics lit- erature. 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 pre- cisely, links have been drawn between granular materials and rigidity percolation, in which one considers a regular lattice of springs, diluted at random with a prob- ability 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 net- work 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 clus- ters 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 0 10 20 30 40 50 60 0 10 20 30 40 50 Figure 6.1: Frame of a polydisperse packing undergoing simple shear. Line 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 mi- croscopic force configurations consistent with macroscopic constraints. As dis- cussed in the previous Chapter, we consider a canonical distribution Peq (f) = 1 Z λNRCe−Aα̂:σ̂[f ] Θ[f ], (6.1) where volume V has been replaced by area A in two dimensions, and Θ[f ] = ∏ g [ δ (∑ c∈Cg f cg ) δ (∑ c∈Cg rc × f cg )]∏ c Θ ( f cN − 1µ |f cT | ) ≡ δFB[f ] δTB[f ] ΘC [f ]. (6.2) 1In 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 Z = ∑ T λNRC ∫ Df e−Aα̂:σ̂[f ] Θ[f ], (6.3) 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 = − logZ. Hence our main goal is compu- tation 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 in- teractions 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 f cg = ρ 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∈Ct `ct ψ c, (6.6) 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∈RC` `c` ψ c, (6.7) 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 |V C ] ≡ ∏ c∈V C δ (f 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 Z = ∑ T λNRCJ ∫ Dψ e−Aα̂:σ̂ δGδ[f |V C ]ΘC . (6.8) 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 factor which depends only on the slow geometry, it is exactly N2ST , the square of the number of spanning trees in the contact network. The latter is exponential in the system size, NST = e qNRC , 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 at each vertex. The Boltzmann factor e−Aα̂:σ̂ weights these surfaces by their cur- 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 (0,f c). Hence the Coulomb inequality |f cT | ≤ µf cN is a constraint on the curvature 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 theo- ries, 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 includ- ing 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 af- fecting 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 2After 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 f cN ≥ 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 prepara- tion 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 degrees- of-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 = N iC − 2N iV C = 3 + 2N iRC − 3N i − zi, (6.11) where the equality is a consequence of Euler’s relation. Here, N iRC and N i V C are the number of real and virtual contacts strictly inside the cluster, zi is the number of contacts around the cluster boundary, N i is the number of non-rattling grains inside the cluster, and N iC = N i RC + N i V C . The crucial difference between (6.11) and (6.9) is in the boundary terms. For a single loop, H i = 3 − zi ≤ 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 perco- lating 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 5 10 15 20 25 30 35 40 45 50 10 15 20 25 30 35 40 45 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 NFC = N∆z +N∂C , (6.12) a consequence of assuming H ′i = 0 for each cluster. For systems large enough that √ N∆z 1, there will be many, large isostatic 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 Z ≈ ZMF = ∑ T λNRCe2qNRC 〈 1 ZFC ∏ i Zi 〉 C , (6.13) where 〈·〉C denotes an average over cluster decompositions, Zi is the partition func- tion for cluster i, and ZFC will correct for overcounting, as discussed below. Here, and in the following, the subscript MF 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 clus- ter boundary, RCi the real contacts in the interior of i, V Ci the virtual contacts in the interior of i, and Ci = RCi ∪ V Ci. For each isostatic cluster, we have a partition function Zi = ∫ Dψ|C∂i ∫ Dψ|Ci e−α̂:(Aσ̂) i ΘiCδ i V C . (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 de- termined by those on the boundary. These interior integrations disappear, giving a Jacobian J i = ∫ Dψ|Ci δiV C . (6.15) As shown in Appendix J J i can be estimated for a class of geometries, with the result J i ≈ ( A 8N )N iV C . (6.16) Hence we have the cluster partition function Zi = J i ∫ Dψ|C∂i e−α̂:(Aσ̂) i ΘiC . (6.17) 92 Note that Zi is holographic: the properties of the cluster are determined by the values of the gauge field ψ on its boundary. In the geometric interpretation, Zi 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 Zi, 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 To compute Zi, we need to express all quantities appearing in ΘiC and σ̂ i in terms of the values of ψ on the boundary C∂i. The stress tensor for the cluster is (Aσ̂)i = − ∑ c∈Ci `cf c = ∑ c∈Ci Acσ̂c = ∫ i dA σ̂ = ∫ i dA ∇× ρ = − ∫ ∂i dr ρ = − ∑ c∈C∂i `cρ`0(c) (6.18) 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 ρ `0(c) MF = ρ̄ `0(c) + 1 ĀL δψc`c, (6.19) where, for simplicity, we assume that the external loop has the mean loop area ĀL = A/NL. Hence, as a consequence of ∇ · σ̂ = 0, the stress tensor is easily written in terms of boundary ψ. We find (Aσ̂)iMF = (Aσ̄) i − ∑ c∈C∂i 1 ĀL δψc`c`c, (6.20) with (Aσ̄)i = −∑c∈C∂i `cρ̄`0(c) = −∑c∈Ci `cf̄ c. 93 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 ∑ `∈Li A`ρ` = − 1 Ai ∑ c∈C∂i `cψc, (6.22) 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 δρ`MF = λ `δρ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 f cMF = f̄ c − δρi + 1 ĀL δψc`c. (6.24) Equivalently, we can write δf cMF = ∑ c′∈C∂i Gcc ′ `c ′ δψc ′ , with Gcc ′ = 1/Ai + { 1/ĀL if c ′ = c 0 if c′ 6= c and δf = f − f̄ . For interior contacts c ∈ RCi, δf cMF = 0. Hence each interior contact in i has its contact force set entirely by its mean value f̄ c, and the fluctu- ations 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 ĀL δψc−`c− − δρi ] + 12` c+ × [ 1 ĀL δψc+`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. It is easily checked that the Coulomb condition Θc(f cMF ) ≥ 0 is equivalent to δψc ≤ ψcM ≡ ĀL ( f̄ cN + ` c · δρi) . (6.25) 6.3.2 Computing the partition function Extracting the mean-stress contribution Z̄i = exp(−α̂ : (Aσ̄)i), we find eSi ≡ Zi/Z̄i = J i ∏ c∈C∂i ∫ dψce 1 ĀL αcδψc ΘiC , (6.26) 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 ĀL )2 ∫ dρi (2pi)2 dkie Ai ĀL (iki+νi)·(ρi−ρi[ψ]) , 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 Zi, changing the order of integration, and performing a trivial change of variables to δψ and δρ, we find eSi = J ′i ∫ dki d(δρi) (2pi)2 e Ai ĀL (iki+νi)·δρi ∏ c∈C∂i ∫ ψcM −∞ d(δψc)e 1 ĀL (αc+γc)δψc , with γc = (iki + νi) · `c (6.27) and J ′i = J i(Ai/Ā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 ∫ dki d(δρi) (2pi)2 e Ai ĀL (iki+νi)·δρi ∏ c∈C∂i 1 αc + γc e 1 ĀL (αc+γc)ψcM = J ′′i ∫ dki ∏ c∈C∂i 1 αc + γc e(α c+γc)f̄cN ∫ d(δρi) (2pi)2 ea i·δρi , (6.28) with ai = Ai ĀL (iki + νi) + ∑ c∈C∂i (αc + γc)`c and J ′′i = J ′i(ĀL)z i . Convergence of the δρi integral requires that Re[ai] = 0, and hence 0 = Ai ĀL νi + ∑ c∈C∂i (αc + νi · `c)`c. (6.29) The δρi integral then becomes δ(Im[ai]), so that 0 = Ai ĀL ki + ki · ∑ c∈C∂i `c`c = ki · M̂ . (6.30) The matrix M̂ is non-singular and approximately equal to (Ai/ĀL+ 1 2z i)δ̂. Hence the only solution for ki is the zero wavenumber ki = 0 3. We find eSi = J ′′′i ∏ c∈C∂i 1 αc + νi · `c e (αc+νi·`c)f̄cN , (6.31) 3This is a consequence of the neglect of inter-cluster coupling. Future work should explore the effect of the latter on k. 96 with νi evaluated with (6.29) and J ′′′i = J ′′i/det(M̂) ≈ J ′′i/(Ai/ĀL+ 12zi)2. Note that νi is a linear functional of the αc’s, which represents an additive renormal- ization of αc from coupling with the cluster. Solving (6.29), νi = − δ̂ + ĀL Ai ∑ c∈C∂i `c`c −1 · ĀL Ai ∑ c∈C∂i αc`c. 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 6= 0 only from coupling between stress anisotropy and geometrical anisotropy. In this work we consider the physically relevant case of large clusters, zi 1, for which νi is negligible independent of α̂. Indeed, Ai ∼ (zi)2, hence νi ∼ 1/zi 1. Neglecting νi, we have our result eSi = J i ∏ c∈C∂i ĀL αc eα cf̄cN . (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 ZFC = (ĀL) NFC [ ∏ c∈FC 1 αc e−α̂:A cσ̄ceα cf̄cN ] . 6.4 Global partition function Hence we obtain the global mean-field partition function ZMF = ∑ T λNRCe2qNRC ( A 8N )NV C 〈( A NL )NFC e−α̂:Aσ̄ ∏ c∈FC 1 αc eα cf̄cN 〉 C . 4Our 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 basic variation of f̄ needed to sustain an applied stress σ̄ can be obtained by solving σ̂ = ∇ × ρ = ∇ × ∇ × ψ for a constant σ̂ = σ̄ in the continuum. The result is ψ̄(r) = −12r × σ̄ × 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 the applied stress. The discrete mean fields are then obtained by ψ̄c ≡ ψ̄(rc) and ρ̄` ≡ ρ̄(r`), for each real contact c and loop `, where r` is the centroid of the loop. The basic mean-field variation of contact forces is then f̄ c0 = ρ̄ `− − ρ̄`+ = −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̄ c0 = (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 −NFC/(NRC −NFC) on RC - FC, and varies smoothly in between, 5To satisfy the Coulomb constraint, one could project f̄ inside the Coulomb cone. We leave this for future work. 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 − NFC)/NFC . 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 + k A ∑ c∈FC `cf̄ c0 − k A NFC NRC −NFC ∑ c∈RC−FC `cf̄ c0 = σ̄0 ( 1− kNFC NRC + k NFC NRC −NFC NRC −NFC NRC ) = σ̄0, where we have used self-averaging to replace, e.g., −∑c∈FC `cf̄ c0 by Aσ̄0NFC/N . 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 pre- scription: 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 this minimization is that while σ̄ in Z comes from all contacts, the f̄ cN term comes only from contacts on force chains. Hence f̄ cN = (1 +k)(f̄ c N )0 for contacts on force chains, and ZMF = ∑ T λNRCe2qNRC ( A 8N )−NRC 〈( A NL )NFC e−α̂:Aσ̄0 ∏ c∈FC 1 αc e(1+k)α c(f̄cN )0 〉 C , where we have ignored an irrelevant constant term. For each cluster decomposi- tion, 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〉, (6.35) 99 and β ≡ 〈αc(f̄ cN )0〉 = α̂ : 〈`c`c(f̄ cN )0〉 ≥ 0. (6.36) Assuming that cluster decompositions are self-averaging, the average over decom- positions gives ZMF = ∑ T λNRCe2qNRC ( A 8N )−NRC e−α̂:Aσ̄0eβ(1+k)NFC ( A ζNL )NFC . The k-dependence is via a term exp(kβ). Since β ≥ 0, exp(β) ≥ 1. For any topol- ogy T , this minimizes the free energy (maximizes Z) when k takes its maximal value (NRC −NFC)/NFC . 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: ∑ T = NGC∑ 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 write ZMF = e −α̂:Aσ̄0 ∑ NRC eS0 and seek the maximum value of S0 over NRC . 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λe2qeβ A ) + (2NRC − 3N) log ( A ζNL ) , (6.37) where we have used the identity (6.12) to eliminateNFC , and ignored sub-extensive boundary terms, irrelevant in the thermodynamic limit. We seek ∂S0/∂NRC = 0. 100 Applying Stirling’s approximation logN ! = N log(N/e) + 12 logN + O(1/N) as N →∞, it is easy to see that log ( N K ) = N logN −K logK − (N −K) log(N −K) +O(logN). The extremal condition is 1 = NGC −NRC NRC ( 8Nλe2qeβ A )( A ζNL )2 e−(2NRC−3N)/NL = Z − z̄ z̄ 32Aλe2qeβ Nζ2(z̄ − 2)2 e −2∆z/(z̄−2) = g(z̄;Z) 32Aλe2qeβ Nζ2 , (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 h̃(λ, α̂) = 32A N e2qeβλ ζ2 = 8pi φ e2qeβλ ζ2 , where φ is area fraction. Here we have used A = AgN/φ with Ag = pi/4 the average area per grain. The final expression for the partition function is ZMF = e −α̂:Aσ̄eS0 , 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(ZMF ) = α̂ : Aσ̄ − S0 = α̂ : Aσ̄ −NGC logNGC +NRC logNRC + (NGC −NRC) log(NGC −NRC) −NRC log(8Nλeβe2q/A)−NRG∆z log(A/(ζNL)). (6.39) In this expression, σ̄, q, A, NRG, and NGC are constants; NRC = NRC(λ, α̂), as specified implicitly by (6.38), ζ = ζ(α̂) = exp(〈logαc〉), and β = β(α̂) = 〈αc(f̄ cN )0〉; and z̄ = 2NRC/N , NL = NRC − N + 1, αc = `c · α̂ · `c. From F 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〉 = ∫ dz̄ ∫ dσ̂ λNRCe−Aα̂:σ̂ Ω(σ̂, z̄)Q(σ̂, z̄), where Ω = eS is the density of states. We can write Ω as Ω = eS = eS1eS2eS3 , where S1 = NGC logNGC −NRC logNRC − (NGC −NRC) log(NGC −NRC), S2 = NRC log(8Ne β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(8Ne 2q/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 allowable normal force (f cN = 0). In a model with a finite coefficient of friction, this would be generalized to β = 〈αc(f̄ cN − 1µ |f̄ cT |)〉, with the same interpretation. 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 The mean-field equation of state is 〈σ̂〉 = σ̄. Using (5.13), σ̄ = 1 A ∂F ∂α̂ = σ̄ − NRC A ∂β ∂α̂ + N∆z A ∂ log ζ ∂α̂ , since the terms from ∂NRC/∂α̂ vanish by virtue of ∂S0/∂NRC = 0. This simplifies to 2∆z z̄ 〈 `c`c `c · α̂ · `c 〉 = 〈 `c`cf̄ cN 〉 , (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`cf̄ cN 〉 = β. (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(λ, α̂) = 8pie2q φ λ ζ2 , where φ is area fraction. The contact potential equation implies that at isostaticity, the geometrical contact number is Ziso = ziso + 3φ 8pie2q ζ2 λ . (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 σ̄ = R−ω · P + T 0 0 P − T ·Rω (6.44) = P + T cos(2ω) −T sin(2ω) −T sin(2ω) P − T cos(2ω) , (6.45) 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 ≤ ω < pi/2. The vectors sc connect the centres of neighbouring loops, and are, on average, perpendicular to `c. We let sc/|sc| = −`c×Rγc , 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 = 12s c× `c = 12 |sc| cos(γc). Writing contact vectors as `c = (cos(θc), sin(θc)), we have f̄ cN/|sc| = −`c ×Rγ c × σ̄ · `c = cos(γc)P + cos(2θc + 2ω − γc)T . For simplicity, here we will consider only macroscopically isotropic geometries for which the distribution of contact angles is P(θc) = 1/(2pi); 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 to decompose the matrix 〈`c`cf̄ cN 〉 = M̂ into M11 + M22 = 〈f̄ cN 〉, M11 −M22 = 〈cos(2θc)f̄ cN 〉, and M12 +M21 = 〈sin(2θc)f̄ cN 〉. Then with the above prescriptions 〈f̄N 〉 = 2A NRC P 〈cos(2θc)f̄N 〉 = AT NRC cos(2ω) 〈sin(2θc)f̄N 〉 = − AT NRC sin(2ω), which gives the right-hand-side of the tensorial equation of state. To compute the left-hand-side, we likewise take coordinates so that α̂ = R−ν · α+ η 0 0 α− η ·Rν , (6.46) 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 ≤ ν < pi/2. In these coordinates, we have αc = `c · α̂ · `c = α+ η cos(2θc + 2ν). Writing η/α = sin(2ξ) (6.47) 104 with |ξ| ≤ pi/4, we find the isotropic-geometry equations of state 1 α cos(2ξ) = 2AP ∆zNRG −tan(ξ) cos(2ν) α cos(2ξ) = AT ∆zNRG cos(2ω) −tan(ξ) sin(2ν) α cos(2ξ) = AT ∆zNRG sin(2ω). These can be rewritten more simply as P = NRG∆z 2A 1 α cos(2ξ) , (6.48) T/P = −2 tan(ξ), (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 〈mv2〉 = kBT . 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 A2 ∂2F ∂α̂∂α̂ = NFC A2 〈 `c`c`c`c (`c · α̂ · `c)2 〉 . (6.51) 105 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 P and T . Using α̂ : σ̄ = 2αP + 2ηT cos(2ω − 2ν), one can verify directly from (6.3) that (δP )2 ≡ 〈(P − 〈P 〉)2〉 = −1 4A2 ∂2F ∂α2 , (δT )2 ≡ 〈(T − 〈T 〉)2〉 = −1 4A2 ∂2F ∂η2 , C ≡ 〈(P − 〈P 〉)(T − 〈T 〉)〉 = −1 4A2 ∂2F ∂α∂η . After some algebra we find (δP )2 = 1 α2 cos3(2ξ) NRG 4A2 [ ∆z − 18g(z̄) z̄2g′(z̄) cos2(ξ) ] (6.52) (δT )2 = 1 α2 cos3(2ξ) NRG 4A2 [ ∆z ( 1− 2 sin 2(ξ) tan2(2ξ) ) − 18g(z̄) z̄2g′(z̄) cos(2ξ) cos2(ξ) ] (6.53) C = − 1 α2 cos3(2ξ) NRG 4A2 [ ∆z sin(2ξ) + 18g(z̄) z̄2g′(z̄) tan(ξ) cos(2ξ) ] (6.54) In each of these covariances, the term with ∆z comes directly from stress fluctu- ations 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 N2RG ∂2F ∂(log λ)2 = 4 N2RG ∂NRC ∂(log λ) = − 2 NRG g(z̄) 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. 3 3.5 4 4.5 0 500 1000 1500 2000 2500 3000 z P1 /2 ×106 0.8 0.82 0.84 0.86 0 0.2 0.4 0.6 0.8 1 1.2 φ (∆ z)2 Figure 6.3: Results of numerical simulations showing the scaling of √ P with z and (∆z)2 with φ. 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 force- bearing grain, and hence NGC = NRC+x0N , 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 sys- tem. Any such subsystem must be sufficiently large for the canonical distribution to hold. Numerical simulations of frictionless packings indicate that the angoricity 6Note 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 4 6 8 10 x 106 0 0.05 0.1 0.15 0.2 P x 0 Figure 6.4: Proportion of rattlers as a function of P , as obtained from numerical simulations. 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 x 107 ∆ z 1/ α 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 x 107 ∆ z 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. 0 0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 14 ∆ z α F/α E 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 pack- ing 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 do- main [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 piξ2I and 2piξI . If there are NI isostatic clusters, we find NIpiξ 2 I = A and NI2piξI = NFC , which implies ξI = 2A NFC ∼ 2A NRG 1 ∆z , where we have used NFC = 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 equilib- rium, it is natural to suppose that the isostatic clusters are statistically inde- pendent. 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 7This precise identification of the minimal changes between configurations in mechanical equi- librium 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 car- ries 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 sim- ulations. 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 elas- 112 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 granu- lar materials, which shed some light on the development of fabric, and the jamming and unjamming transitions between solid and gaseous phases. Owing to incom- plete understanding of contact opening and contact closing processes, the defor- mation 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 fluc- tuations 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 ap- proximations 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 sim- ulations 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 con- tinuum 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ément. 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 Shat- tuck. 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öhlich. Diseases of triangulated random surface models, and possible cures. Nuclear Physics B, 257:433–449 [AKV+04] E. Almaas, B. Kovács, T. Vicsek, Z. N. Oltvai, and A.-L. Barabási. 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ël 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éma and Farhang Radjäı. Stress-strain behavior and ge- ometrical properties of packings of elongated particles. Physical Re- view 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 pseudo- elastic 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 re- view letters, 96(12):120601 [BKLS00] A. Barrat, J. Kurchan, V. Loreto, and M. Sellitto. Edwards’ mea- sures 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, vol- ume 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 dimen- sional quantum gravity. Nuclear Physics B, 438(1):320–344, 1995. [CBD12] C Coulais, RP Behringer, and O Dauchot. Dynamics of the con- tacts reveals widom lines for jamming. EPL (Europhysics Letters), 100(4):44005, 2012. [CCL97] F. Calvetti, G. Combe, and J. Lanier. Experimental micromechan- ical analysis of a 2d granular material: relation between structure evolution and loading path. Mechanics of Cohesive-frictional Mate- rials, 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. Physical Review Letters, 109(20):205501–, 11 2012. [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 frot- tement de leurs parties et A la roideur des cordages. Bachelier, Paris, 1821. [CR00] Gaël Combe and Jean-Noël 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öter, 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. Jam- ming, 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 be- haviour. 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ël Forterre and Olivier Pouliquen. Flows of dense granular media. Annual Review of Fluid Mechanics, 40(1):1–24, 2008. [GBO06] Guo-Jie Gao, Jerzy B lawzdziewicz, 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 Shat- tuck. Experimental demonstration of nonuniform frequency distri- butions of granular packings. Physical Review E, 80(6):061304, 2009. [GF82] B.J. Gellatly and J.L. Finney. Characterisation of models of multi- component 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ément, and S. Luding. Footprints in sand: The re- sponse 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ément, and R. P. Behringer. Green’s function measurements of force transmission in 2d granular materi- als. Physica D: Nonlinear Phenomena, 182(3–4):274–303, 8 2003. 117 [GRH+90] Etienne Guyon, Stéphane 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 imagi- nation. 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 prin- ciples. Journal of the Mechanics and Physics of Solids, 11(5):357– 372, 9 1963. [Hil98] Rodney Hill. The mathematical theory of plasticity, volume 11. Ox- ford university press, 1998. [HOC07] Silke Henkes, Corey S. O’Hern, and Bulbul Chakraborty. Entropy and temperature of a static granular assembly: An Ab Initio ap- proach. 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, liq- uids, 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 ma- terials. Comptes Rendus Mécanique, 338(10–11):596–603, 2010/11// 2010. [Kuh99] M.R. Kuhn. Structured deformation in granular materials. Mechan- ics of materials, 31(6):407–429, 1999. [Kur00] J. Kurchan. Emergence of macroscopic temperatures in systems that are not thermodynamical microscopically: towards a thermo- dynamical 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 jam- ming/rigidity transition of a granular glass. EPL (Europhysics Let- ters), 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. Funda- menta 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 measure- ments 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 Exper- iment, 2011(07):L07002, 2011. [MD95] C. Moukarzel and P. M. Duxbury. Stressed backbone and elas- ticity 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ër, and R. De- lannay. 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 ellipse- shaped 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. Cam- bridge, 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 compres- sional deformation of sand. Soils and Foundations, 12(2):1–18, 1972. [OK74] Masanobu Oda and Junichi Konishi. Microscopic deformation mech- anism of granular material in simple shear. Soils and Foundations, 14(4):25–38, 1974. [OKNN82] Masanobu Oda, Junichi Konishi, and Siavouche Nemat-Nasser. Ex- perimental micromechanical evaluation of strength of granular ma- terials: 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ák Somfai, and Bernard Nienhuis. Scale invari- ance 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 temperature- like variables in jammed granular subsystems. arXiv preprint arXiv:1207.7349, 2012. [PEKK12] C. Pérez-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écile Appert-Rolland, François Chevoir, Philippe Gondret, Sylvain Lassarre, Jean-Patrick Lebacque, and Michael Schreckenberg. Stability and Jamming Transition in Hard Gran- ular 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éphane Roux. Nonsmooth- ness, indeterminacy, and friction in two-dimensional arrays of rigid particles. Physical Review E, 54(1):861–873, 07 1996. [RC02] Jean-Noël Roux and Gaël 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éma, and S. Roux. Fabric evolution and accessible geometrical states in granular materials. Granular Matter, 14(2):259–264, 2012. [REYR06] Vincent Richefeu, Moulay Säıd El Youssoufi, and Farhang Radjäı. 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 granu- lar materials. Journal of non-crystalline solids, 352(42):4505–4508 0022–3093, 2006. [Rou00] Jean-Noël Roux. Geometric origin of mechanical properties of gran- ular materials. Phys. Rev. E, 61:6802–6836, Jun 2000. [RR05] Farhang Radjäı and Stéphane 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 mate- rials. 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 quanti- ties for granular assemblies. Int. J. of Sol. and Struc., 41(21):5775 – 5791, 2004. [SEmcG+02] Leonardo E. Silbert, Deniz Ertaş, Gary S. Grest, Thomas C. Halsey, and Dov Levine. Geometry of frictionless and frictional sphere pack- ings. 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 differ- ence 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. [SRC+01] D. Serero, G. Reydellet, P. Claudin, É. Clément, and D. Levine. Stress response function of a granular layer: Quantitative compar- ison between experiments and isotropic elasticity. The European Physical Journal E: Soft Matter and Biological Physics, 6(2):169– 179, 2001. [SVC02] D. Segré, 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 fric- tional 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 en- sembles 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ás Unger, János Kertész, and Dietrich E. Wolf. Force indetermi- nacy 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 distribu- tions under sandpiles. Phys. Rev. E, 60:R5040–R5043, 1999. [VRDEY09] C. Voivret, F. Radjai, J.Y. Delenne, and MS El Youssoufi. Multi- scale 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 Grin- spun. 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 amor- phous 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 8 9 10 11 12 13 14 12 12.5 13 13.5 14 14.5 15 15.5 16 16.5 17 g stg ℓc + ℓc − ℓc 0 Figure A.1: Geometry used in torque balance computation. 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 in Figure A.1, we note that stg is parallel to ` c0(g,t), so that Ag(∇ · ρ)g = ∑ t∈T g stg × ρt = ∑ t∈T g stg × −1 At ( −`c−ψc− + `c+ψc+ + `c0ψc0 ) = ∑ t∈T g 1 2 (`c + − `c−)× −1 At ( −`c−ψc− + `c+ψc+ ) = ∑ t∈T g 1 At ( Atψc − −Atψc+ ) = 0 on summation around the grain. 127 Appendix B Discrete Calculus in 3D c g stgc `cg 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:∫ g dV ∇ · σ̂ ≡ V g (∇ · σ̂)g ≡ ∑ c∈Cg Ac ˜̀cg · σ̂c ≡ ∫ ∂g dS · σ̂, 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 = −`cf 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 · ˜̀cf c = −f cg , 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∫ g dV ∇× ρ̂ ≡ V g (∇× ρ̂)g ≡ ∑ t∈T g dStg × (ρ̂t)T ≡ ∫ ∂g dS × ρ̂T . The natural way to define dStg is to decompose the contact faces on either side of t into triangles of area Atc± and set dS t g = A t c+ ˜̀c+ g + A t c− ˜̀c− g ≡ ∑ c∈Ctg A t c ˜̀c g. It is easily seen from Figure B.1 that Atc ˜̀c g = 1 2(r t − rc) × stgc, where rt is the 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 = −12 ∑ t∈T g ∑ c∈Ctg ( (rc − rt)× stgc )× (ρ̂t)T = −12 ∑ t∈T g ∑ c∈Ctg ( (rc − rg)× stgc )× (ρ̂t)T = −14 ∑ c∈Cg ∑ t∈T c ( `cg × stgc )× (ρ̂t)T , where we have used the identity ∑ c∈Ctg s t gc = 0 to exchange r t with rg. Using (1.8), we now see that σ̂g = ∇× ρ̂t if `cgf cg = 12 ∑ t∈T c ( `cg × stgc )× (ρ̂t)T . This is equivalent to f cg = −12 ∑ t∈T c ρ̂t · stgc ≡ −12 ∮ (∂c)g ρ̂ · dr, (B.1) and 0 = ∑ t∈T c s t gcρ̂ t · `cg. The latter equation is satisfied identically if we let ρ̂t = − 2|st|2ρtgc stgc. This allows (B.1) to be rewritten f cg = ∑ t∈T c ρ t gc, 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 ρ t gc → ρtgc + (∆ρ)tgc with (∆ρ)tgc = Bv + − Bv− , for any vector field Bv 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]∫ Ω dV σ̂ ≡ V σ̂ = ∑ c∈BC Aclc × (ρt)T ≡ ∮ ∂Ω dS × ρT , where the l’s are oriented outwards. We now seek ψ̂c such that ρ̂t = ∇× ψ̂c. Using Stokes’ theorem, we set∫ t (∇× ψ̂) · dS ≡ At ( ∇× ψ̂ )t · s̃t ≡ ∑ c∈Ct ψ̂c · `c ≡ ∫ ∂t ψ̂ · dr, where s̃ = s/|s| and the contour is oriented anticlockwise around s̃t. A natural choice to guarantee ρ̂t = ∇ × ψ̂c is to set ψ̂c = 14 ˜̀c ˜̀cψc; then we have Atρtgc = −18 |st| ∑ c∈Ct ` c gtψ 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 NgV −NgT +NgC = 2. Each vertex is the confluence of three edges of the polyhedron, so that 3NgV = 2N g T . 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 Bv 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 = rc · B, which lead to (∆ρ)tgc = −stgc×B. HereB is a constant, but it may be derived from a fluctuating 130 field on the triangles, viz. B = ∑ t∈T cB t. This gives 3NT unknowns constrained by 3NC equations. Second, we can have ∆ψ c = |`|−1 ˜̀c · (∆ψg+ − ∆ψg−), with ∆ψg = rgB, where again B is a constant. In this case it can be derived from B = ∑ c∈Cg 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. Writing Ĥ ≡ ∇∇ψ, convexity of ψ is equivalent to positive definiteness of Ĥ. The latter is equivalent to the statement that Ĥ has positive eigenvalues. Writing λi and ui for the eigenvalues and eigenvectors of σ̂, we have tr(σ̂) = ∑ i λi. From (C.1), Ĥ ·uj = ∑ i 6=j λiuj , so uj is an eigenvector of Ĥ with eigenvalue ∑ i 6=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 f cN = kNδ c df cT /dt = kTV c T , 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 |f cT | ≤ µf cN , such that contacts begin sliding when this inequality is saturated. The code di- rectly solves Newton’s laws mg d2rg dt2 = ∑ c∈Cg f cg Ig dωg dt = ∑ c∈Cg rc × f cg , where mg and Ig 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 quan- tities 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 ĝψ = δ̂ + (∇ψ)(∇ψ), (E.1) which measures distances along the surface, and the shape tensor2 I = 1√ 1 + |∇ψ|2∇∇ψ, (E.2) which measures curvature. The Gaussian curvature is K = det I det ĝψ = det(∇∇ψ) +O(ψ4) (E.3) Using ∇∇ψ = ε̂ · σ̂ · ε̂T , we have K = det(σ̂) +O(ψ4), (E.4) while the mean curvature is H = 1 det ĝψ tr ( ĝψ · ε̂ · I · ε̂T ) = 1 1 +O(ψ2) tr ( σ̂ +O(ψ2) ) = 2P +O(ψ2). (E.5) 1Also called the first fundamental form. 2Also called the second fundamental form. 134 The celebrated Gauss-Bonnet theorem states that∫ Ω dA K = ∫ ∂Ω kg ds+ 2pi, (E.6) where Ω is the surface and kg is the geodesic curvature of the boundary. The physical interpretation of this theorem is that A = ∫ Ω dA det(σ̂) (E.7) depends only on boundary quantities. In fact, A is exactly the area of the Maxwell- Cremona reciprocal tiling discussed in Chapter 2; compare equation (2.20). 135 Appendix F Renormalization of ρ. Here we compare ρ̄(r`) with ρ`. Throughout, we will ignore correlations between 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 ĝ` · ρ` = − 1 A` ∑ c∈C` `c ( tc` : ∇ψ̄(r`) + 12tc`tc` : ∇∇ψ̄(r`) ) = ĝ` ×∇ψ̄(r`)− 12Ĝ`2 : ∇∇ψ̄(r`), (F.1) with Ĝ`2 = 1 A` ∑ c∈C` `ctc`t c `. (F.2) The tensor Ĝ vanishes on average and hence does not give a systematic correction. To see this, consider the related tensor ĝ`2 = −(∇× (rr))` = 1 A` ∑ c∈C` `crcrc = 1 A` ∑ c∈C` `c(tc` + r `)(tc` + r `) = Ĝ`2 + 1 A` ∑ c∈C` `crcr` + 1 A` ∑ c∈C` `cr`rc. (F.3) In index notation1, we find (ĝ`2 − Ĝ`2 · ε̂)ijk = (ĝ` · ε̂T )ij(r`)k + (ĝ` · ε̂T )ik(r`)j . (F.4) 1There is an insufficiency of notation here, since (ĝ` · ε̂T )ik(r`)j is not of the form ÂijBk. 136 However, discrete calculus implies 〈ĝ`2〉ijk = −〈(∇× (rr))`〉 = −εijrk − εikrj = 〈(ĝ` · ε̂T )ij〉〈(r`)k〉+ 〈(ĝ` · ε̂T )ik〉〈(r`)j〉. (F.5) Hence 〈Ĝ`2〉 = 0, and to quadratic order there is no renormalization of ρ; ρ = ρ̄. In fact, this argument can be continued to all orders; only the tensor algebra is more complicated. Let us write Â ⌋ n for the part of Â which is symmetric in its last n indices. Then we define ĝ`n = − (∇× (rn))` = 1 A` ∑ c∈C` `c(rc)n = 1 A` ∑ c∈C` `c(tc` + r `)n (F.6) so that ĝ`n ⌋ n = 1 A` ∑ c∈C` `c n∑ k=1 ( n k ) (tc`) k(r`)n−k ⌋ n , (F.7) where the sum runs from k = 1 rather than k = 0 because ∑ c∈C` ` c = 0. Ignoring correlations between r` and other loop quantities, we have 〈ĝ`n ⌋ n 〉 = 〈 1 A` ∑ c∈C` `c n∑ k=1 ( n k ) (tc`) k 〉 rn−k ⌋ n = n 〈 ĝ` 〉 · ε̂Trn−1 ⌋ n + n∑ k=2 ( n k )〈 1 A` ∑ c∈C` `c(tc`) k 〉 (r)n−k ⌋ n . (F.8) However, with discrete calculus we observe that 〈ĝ`n ⌋ n 〉 = −∇× (r)n⌋ n = −n (∇× r)rn−1⌋ n = −n ε̂ rn−1⌋ n . (F.9) 137 Since 〈ĝ`〉 = δ̂, we find n∑ k=2 ( n k ) Ĝkr n−k ⌋ n = 0, (F.10) with Ĝk = 〈 1 A` ∑ c∈C` `c(tc`) k 〉 . (F.11) Inductively, (F.10) implies Ĝk = 0 for all k ≥ 2. Indeed, for n = 2 we find 0 = Ĝ2 ⌋ 2 = Ĝ2. This then implies 0 = Ĝ3 ⌋ 3 = Ĝ3, etc. So we find that ρ = ρ̄ to 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∣∣ y=0 = 2piA δ(x) is given by ψ(z, z̄) = 12 z̄f(z) + 1 2g(z) + c.c., (G.1) with f(z) = −A log z and dg/dz = Ā log z. One can directly compute ∂2zψ = Az̄ 2z2 + Ā 2z ∂2z̄ψ = Āz 2z̄2 + A 2z̄ = ∂2zψ ∂zz̄ψ = − Ā 2z̄ − A 2z . We use σxx = 2∂zz̄ψ − ∂2zψ − ∂2z̄ψ σyy = 2∂zz̄ψ + ∂ 2 zψ + ∂ 2 z̄ψ σxy = −i(∂2zψ − ∂2z̄ψ) and the identities <[zw̄] = w · z =[zw̄] = w × z <[z2] = x2 − y2 =[z2] = 2xy 139 with z = (<[z],=[z]), etc. Then ∂2zψ = A · z z2 ∂2z̄ψ = A · z z̄2 ∂zz̄ψ = −A · z|z|2 and σxx = −4x 2A · r |r|4 σyy = −4y 2A · r |r|4 σxy = −4xyA · r|r|4 (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) dz = A1 log z, dg(w) dw = A2 logw. We use the identities [Sad09] σxx = −f ′′(z) + µ2g′′(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 σxx = −A1 · r|r|2 + A3 ·w |w|2 (G.6) σyy = + A1 · r |r|2 + A2 ·w |w|2 (G.7) σxy = + r ×A1 |r|2 − A4 ·w |w|2 , (G.8) with A3 = A2µ 2 and A4 = A2µ. 140 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 rc 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̂ ` ≡ ∑ c′∈C` tc ′ ` t c′ ` = ∑ c′∈C` rc ′ rc ′ − z`r`r`. It follows that if loops ` and `′ merge into `′′ by opening of their common contact c, then the change in NRCF̂ is localized to the contributions from these loops. The change in NRCF̂ is δF̂− ≡ z`′′F̂ `′′ − z`F̂ ` − z`′F̂ `′ = ∑ c′∈C`′′ tc ′ `′′t c′ `′′ − ∑ c′∈C` tc ′ ` t c′ ` − ∑ c′∈C`′ tc ′ `′t c′ `′ = −z`′′r`′′r`′′ − 2rcrc + z`r`r` + z`′r`′r`′ . We have basic relations z` ′′ = z` + z` ′ − 2 and z` ′′ r` ′′ = ∑ c′∈C`\{c} rc ′ + ∑ c′∈C`′\{c} rc ′ = z`r` + z` ′ r` ′ − 2rc, which imply z` ′′ tc`′′ = z `′′(rc − r`′′) = (z` + z` ′ − 2)rc − (z`r` + z`′r`′ − 2rc) = z`tc` + z `′tc`′ . 141 Hence δF̂− = − ( z`r` + z` ′ r` ′ − 2rc)r`′′ − 2rcrc + z`r`r` + z`′r`′r`′ = z`r`(tc`′′ − tc`) + z` ′ r` ′ (tc`′′ − tc`′)− 2rctc`′′ = ( z`r` + z` ′ r` ′ − 2rc)tc`′′ − z`r`tc` − z`′r`′tc`′ = z` ′′ (rc − tc`′′)tc`′′ − z`(rc − tc`)tc` − z` ′ (rc − tc`′)tc`′ = −z`′′tc`′′tc`′′ + z`tc`tc` + z` ′ tc`′t c `′ . This form is suggestive, but it will be useful to eliminate tc`′′ completely, giving δF̂− = − 1 z`′′ ( z`tc` + z `′tc`′ )( z`tc` + z `′tc`′ ) + z`tc`t c ` + z `′tc`′t c `′ = ( 1− z ` z`′′ ) z`tc`t c ` + ( 1− z `′ z`′′ ) z` ′ tc`′t c `′ − z`z` ′ z`′′ ( tc`′t c ` + t c `t c `′ ) = z` ′ − 2 z`′′ z`tc`t c ` + z` − 2 z`′′ z` ′ tc`′t c `′ − z`z` ′ z`′′ ( tc`′t c ` + t c `t c `′ ) = z`z` ′ z`′′ ( tc` − tc`′ )( tc` − tc`′ )− 2 z` z`′′ tc`t c ` − 2 z` ′ z`′′ tc`′t c `′ . 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′) = ∑ c f c ec(c ′), (I.1) 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′) = ∑ ` ρ` z`(c ′) + ∑ p ηp up(c ′). (I.2) 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∂(ρ,η) ∣∣∣∣ = N2ST , (I.3) 143 where NST is the number of spanning trees in the contact network 1. 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 = e qNRC . 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 = ∑ T λNRCN2ST ∫ DρDη e−V α̂:σ̂ Θ[f ] = ∑ T λNRCN2ST ∫ Dρ e−V α̂:σ̂ δG[ρ] δTB[ρ] ΘC [ρ], (I.4) 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 1 = ∏ t ∫ dρt δ(ρt − ρ`(t)). (I.5) 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. 1This 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 J i = ∫ Dψ|Ci δiV C . (J.1) From the fact that each cluster is isostatic, H i = 0, we have N iRC = N i V C ((6.11)). 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 δ(f cN ), and (ii) integrals over real contacts c of δ(f c′ T ) for a virtual contact c ′. For 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 ∫ dψcδ(f cN ) = ∫ dψcδ( |` c| Ãc ψc + · · · ) = |Ãc|. (J.2) Here 1/Ã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δ(f c ′ T ) = ∫ dψcδ(`c × `c′ 1|`c|Atψc + · · · ) = |`c|At 2At , (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 J i ≈ ( A 8N )N iV C . (J.4) 147 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 δρ`MF = λ `δρi, (K.1) where ρi = 1 Ai ∑ `∈Li A`ρ` = − 1 Ai ∑ c∈C∂i `cψc, (K.2) λ` 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 ρ` → ρ`MF . This requires δρi = 1 Ai ∑ `∈Li A`δρ`MF = 1 Ai ∑ `∈Li A`λ`δρi = − 1 Ai ∑ c∈C∂i 1 Ai ∑ `∈Li A`λ` `cδψc, hence 1 Ai ∑ `∈Li A`λ` = 1. (K.3) 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 lc0 c 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. exact discrete calculus relation ψc − ψc0 = − ∫ cc0 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 〈∑c∈C∂i|c0(c,j)=c′ `c〉j ≡ 1M ∑Mj=1∑c∈C∂i|c0(c,j)=c′ `c = `c′zi/z` (see Figure K.1). We have Aiδρi = − ∑ c∈C∂i `cδψc = − ∑ c∈C∂i `c 〈 δψc0(c,j) − ∫ c c0(c,j) dr × δρ 〉 j = −z i z` ∑ c′∈RC` `c ′ δψc ′ + ∑ c∈C∂i `c 〈∫ c c0(c,j) dr × δρ 〉 j . (K.4) The first term in (K.4) is simply (zi/z`)A`δρ`. The second term is a discrete area integral over the cluster, outside of `. Applying (K.1),∑ c∈C∂i `c 〈∫ c c0(c,j) dr × δρMF 〉 j = ∑ c∈C∂i `c 〈∫ c c0(c,j) drλ 〉 j × δρi. The simplest solutions will have, for `1 along the contour from c0 to c, λ `1 un- correlated with `c ∫ c c0(c,j) dr. In a mean-field approximation, we therefore make a 149 replacement λ`1 → λ̄ ≡ 1 Ai −A` ∑ `′∈Li,`′ 6=` A` ′ λ` ′ = 1 + A` Ai −A` (1− λ `), assuming that the `1’s are sampled, on average, proportional to their area 1. Hence we let Aiδρi − z i z` A`δρ`MF = λ̄ ∑ c∈C∂i `c 〈∫ c c0(c,j) dr 〉 j × δρi = λ̄ ∑ c∈C∂i `c(rc − 〈rc0(c,j)〉j)× δρi = λ̄ ∑ c∈C∂i `crc − z i z` ∑ c′∈RC` `c ′ rc ′ × δρi Using the geometric identity ∑ c∈RC` ` crc = A`ε̂T , valid for any loop, we see that Aiδρi − z i z` A`λ`δρi = λ̄ ( Ai − z i z` A` ) δρi, (K.5) from which it follows that λ` = 1. 1In fact, loops will be sampled with a weight proportional to |r`1 − r`|−1A`1 . This leads to the same conclusion. 150 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θ) = ±pi/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 Cn cosn(γc − aγ sin(2θc)), (L.1) for |γc − aγ sin(2θc)| < pi/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 crystal- lization. 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 2pi (1 + ac cos(2θ c)) . (L.2) 1The normalization constant Cn = pi[(n− 1)(n− 3) · · · (1)]/[n(n− 2) · · · (2)]. 151 −1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 x 10−3 γc P( γc ) 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. As in the main text, we decompose the matrix 〈`c`cf̄ cN 〉 = M̂ into M11+M22 = 〈f̄ cN 〉, M11 −M22 = 〈cos(2θc)f̄ cN 〉, and M12 +M21 = 〈sin(2θc)f̄ cN 〉. 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ω) + Tac Cn−2 Cn sin(2θ) sin(2θ + 2ω) +O(a2c) ) . (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 〈f̄N 〉 = 2A NRC ( P + 12Tac cos(2ω) ) 〈cos(2θc)f̄N 〉 = A NRC (T cos(2ω) + acP ) 〈sin(2θc)f̄N 〉 = − A NRC T sin(2ω). 152 0 0.5 1 1.5 2 2.5 3 3.5 0 1 2 3 4 5 6 x 105 θc f N Figure L.2: Probability distribution of f cN from simulations of a simple shear flow (×), 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 α̂ = R−ν · α+ η 0 0 α− η ·Rν , 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 ≤ ν < pi/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 |λ| ≤ pi/4, we find the equations of state 1 α cos(2λ) (1− ac cos(2ν) tan(λ)) = 2A NFC ( P + 12Tac cos(2ω) ) − 1 α cos(2λ) ( tan(λ) cos(2ν)− ac ( 1 2 cos(4ν) sec 2(λ) + sin2(2ν) )) = A NFC (T cos(2ω) + acP ) 1 α cos(2λ) ( tan(λ) sin(2ν) + ac (− sec2(λ) + cos(2ν))) = − A NFC T sin(2ω). 153 We notice that σ̂ is no longer aligned with α̂, i.e., ω 6= ν. 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 geo- metrical anisotropy. We omit these, as the lengthy expressions offer little insight. 154
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Continuum limits of granular systems
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Continuum limits of granular systems DeGiuli, Eric 2012
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Continuum limits of granular systems |
Creator |
DeGiuli, Eric |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | 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 homo- geneous 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 iso- staticity 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 con- tact forces ubiquitously observed in experiments, arise from an entropic instability which favours localization of contact forces. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-05-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0071936 |
URI | http://hdl.handle.net/2429/44393 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2013-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2013_fall_degiuli_eric.pdf [ 3.1MB ]
- Metadata
- 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
- 24-1.0071936-fulltext.txt
- Citation
- 24-1.0071936.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0071936/manifest