Citation
Ellenbroek, W. G. (2007, June 14). Response of granular media near the jamming transition. Leiden Institute of Physics, Institute-Lorentz for Theoretical Physics, Faculty of Science, Leiden University. Retrieved from https://hdl.handle.net/1887/12083
Version: Corrected Publisher’s Version
License: Licence agreement concerning inclusion of doctoral thesis in the Institutional Repository of the University of Leiden
Downloaded from: https://hdl.handle.net/1887/12083
Note: To cite this publication please use the final published version (if applicable).
near the Jamming Transition
W. G. Ellenbroek
particles that have an increased (decreased) force between them. The back cover shows a similar network for a global shear deformation.
near the Jamming Transition
P R O E F S C H R I F T
ter verkrijging van
de graad van Doctor aan de Universiteit Leiden,
op gezag van Rector Magnificus prof. mr. P.F. van der Heijden,
volgens besluit van het College voor Promoties
te verdedigen op donderdag 14 juni 2007
klokke 13.45 uur
door
Wouter Gerhardus Ellenbroek
geboren te Leiden
in 1979
Promotor: Prof. dr. ir. W. van Saarloos Co-Promotor: dr. M. L. van Hecke
Referent: Prof. dr. J. M. J. van Leeuwen Overige leden: Prof. dr. G. T. Barkema
Prof. dr. H. J. F. Knops (Radboud Universiteit Nijmegen) dr. D. van der Meer (Universiteit Twente)
Prof. dr. B. Nienhuis (Universiteit van Amsterdam) Prof. dr. J. M. van Ruitenbeek
Prof. dr. H. Schiessel
Casimir PhD Series, Delft-Leiden, 2007-05 ISBN 978-90-8593-029-7
Dit werk maakt deel uit van het onderzoekprogramma van de Stichting voor Fundamenteel Onderzoek der Materie (FOM), die financieel wordt gesteund door de Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO).
iv
They were pretty to look at and marbles Was always my favourite game
We played all the summer days In the stony alleyways
In the playground after class We would trade the coloured glass More valuable than diamonds More magical than diamonds
Steve Hogarth
1 Introduction 1
1.1 Jamming . . . 1
1.2 Granular media . . . 3
1.2.1 Static granular media . . . 5
1.3 Jamming in granular media . . . 5
1.4 This thesis . . . 8
1.4.1 Linear response . . . 8
1.4.2 Force network ensemble . . . 9
1.4.3 Organization . . . 10
2 Linear response of frictionless granular packings 13 2.1 Introduction . . . 14
2.2 Numerical procedure . . . 15
2.2.1 Packing generation . . . 15
2.2.2 Linear response . . . 16
2.3 Macroscopics and continuum elasticity . . . 19
2.3.1 Point response . . . 19
2.3.2 Bulk deformations . . . 22
2.4 Microscopics: local floppiness . . . 24
2.4.1 The displacement angle distribution . . . 25
2.4.2 Divergence of local displacements . . . 28
2.4.3 The scaling of ∆z . . . . 30
2.5 Critical length scale and finite size scaling . . . 32
2.5.1 Observation of ℓ∗in inflation response . . . 34
2.5.2 Inhomogeneity of global response . . . 36
2.5.3 Finite size effects . . . 37
2.6 Discussion . . . 39
2.6.1 Elastic moduli . . . 39
2.6.2 P (α) and steric hindrance . . . . 41
2.A The Hertzian interaction . . . 42
2.B 2D elasticity and point response . . . 43
2.B.1 Response to a point force . . . 44
3 Critical and non-critical jamming of frictional grains 47
3.1 Introduction . . . 48
3.2 Jamming of frictional particles . . . 48
3.3 Numerical frictional packings . . . 50
3.3.1 Packing generation . . . 50
3.3.2 Variation of z . . . . 51
3.4 Density of states . . . 52
3.4.1 Scaling of ω∗ . . . 54
3.5 Elastic moduli . . . 54
3.6 Discussion . . . 56
3.A The frictional dynamical matrix . . . 58
4 Theory of the Force Network Ensemble 63 4.1 Introduction . . . 64
4.1.1 Definition of the force network ensemble . . . 67
4.2 Regular packings: balls in a snooker triangle . . . 68
4.2.1 Three balls . . . 69
4.2.2 Numerical analysis of larger systems . . . 72
4.2.3 P (f) and phase space geometry . . . . 74
4.3 General formulation for arbitrary packing geometry . . . 76
4.3.1 Mathematical definition of the ensemble . . . 77
4.3.2 Calculation of P (f) . . . . 79
4.4 Exact results for small crystalline packings . . . 81
4.4.1 2D triangular packings with periodic boundaries . . . 81
4.4.2 3D FCC packing with periodic boundaries . . . 83
4.4.3 6 balls in a snooker-triangle . . . 84
4.5 Beyond packings . . . 85
4.5.1 Random matrices . . . 85
4.5.2 Perturbing a physical packing matrix . . . 89
4.6 Discussion . . . 90
4.A The pressure constraint . . . 95
4.B The Gaussian random matrix P (f) . . . . 97
5 Sheared force-networks: anisotropies, yielding and geometry 101 5.1 Introduction . . . 102
5.2 Imposed-shear force network ensemble . . . 105
5.3 Angle resolved force distributions . . . 105
5.4 Analytical bounds on the shear stress . . . 107
5.5 High-dimensional ensemble . . . 108
5.6 Discussion . . . 111
6 Bounds on the shear load of cohesionless granular media 113
6.1 Introduction . . . 114
6.2 From contact forces to macroscopic stress . . . 117
6.3 Frictionless packings . . . 119
6.3.1 Realistic force anisotropy . . . 119
6.3.2 The limit N→ ∞ and tensile stresses . . . 120
6.4 Frictional packings: τmax(µ) . . . 122
6.4.1 The optimization problem . . . 122
6.4.2 Analytical solution . . . 124
6.5 Fabric anisotropy . . . 126
6.6 Discussion . . . 127
6.A Steric exclusion and the coordination number . . . 129
6.B Frictionless case at arbitrary order . . . 130
6.C The anisotropic frictionless optimization . . . 131
Bibliography 133
Summary (in Dutch) 139
Publications 143
Curriculum Vitae (in Dutch) 145
Acknowledgments (in Dutch) 147
Introduction
1.1 Jamming
Two of the dictionary explanations of being jammed are “packed or squeezed tightly into a space” and “blocked through crowding”. Anyone who has ever packed a suitcase knows that these two meanings are closely related: The more stuff you squeeze in there, the less likely it is that the contents will move around through the suitcase after it has been closed. The same thing happens in traffic every morning.
The situation that particles become stuck in such a way that they start to re- sist stress is what defines jamming [1]. This means they cannot be made to flow by applying a small stress. In physics, the most well-known situation in which this happens is ordinary freezing: In a liquid, molecules move about more or less freely, while in a solid they are nicely ordered in a crystal lattice, and move only small distances due to thermal vibrations. A liquid flows easily, while a solid does not. However, the term jamming has been reserved to describe situations in which both the flowing phase and the jammed phase are disordered. The classical example of such a situation is a glass: there is no way of telling a glass from a liquid by just looking at the structure — the distinction is purely based on the fact that a glass does not flow. The commonly used experimental definition of the glass transition is the point where the viscosity of the fluid becomes immea- surably large, i.e. ν = 1014Pa s. Obviously, this definition contains a somewhat arbitrary choice of “immeasurably large”. This is true in general for jamming: to determine whether something can flow or not, one has to apply a shear stress and see if the material starts to flow, or if the stress relaxes away. Of course this definition is limited by the time the experimentalist has (or is willing to take) to
Figure 1.1: The jamming diagram, as proposed by Liu and Nagel (figure taken with permission from Ref. [1]).
wait and see what happens after applying the stress. Play around with the toy known as “silly putty” and you realize how important this is: throw it on the floor and it will bounce like rubber, leave it on the table for a while and it will spread out like a fluid.
There are many more examples of disordered materials that display such a jamming transition. Shaving foam really flows like a fluid when it comes out of the can, but when it hits your hand it “freezes” in a fraction of a second. If you gently vibrate your hand, the blob of foam will vibrate along more or less elastically and keep its shape. However, it can be easily made to flow again, by rubbing it onto your skin. To make mayonnaise, one has to create tiny droplets of oil in a mixture of egg-yolk and lemon juice. At the beginning there are only a few such droplets and the mixture is still fluid, but after adding more and more oil at some point one ends up with an emulsion that appears solid: the droplets are packed so closely together that they cannot move anymore. Finally, take a walk on the beach. You do not sink (at least not very deeply), because the sand appears solid. Then try to build a steep hill from dry sand. It does not work, because the sand keeps flowing down whenever you try to make the hill steeper than some critical angle.
These everyday examples may not sound very surprising, but in any case they lead one to wonder to what extent the jamming of all these different systems is similar. What can we learn from one by investigating the other? Is there a general way to characterize jamming in these systems, just like we have characterized normal freezing by looking at crystalline ordering? To facilitate the discussion about these questions, Liu and Nagel introduced a jamming diagram [1], which
summarizes the jamming behavior of disordered materials in terms of three pa- rameters that typically control it: density ρ, temperature T , and load Σ. The load Σ is usually thought of as a shear loading, like in the case of rubbing shaving foam against skin. This diagram, displayed in figure 1.1, is a kind of “phase dia- gram”, in the sense that it described the state of a system in terms of its control parameters. It is important to note that it is not a phase diagram in the ordi- nary thermodynamic sense, because the systems at hand are not necessarily in thermodynamic equilibrium. Often, these systems have a significant history de- pendence: It is not so easy to start an avalanche, but once it is going, it is pretty hard to stop it.
Let us end this section with a few small remarks about the axes of the dia- gram, which are important for the context of this thesis.
The density axis is often related to a pressure that keeps the system together:
This is air pressure for the shaving foam, osmotic pressure for emulsions, and mechanical pressure from gravity in the case of beach sand. In most of our numerical simulations, we will use a prescribed value of the confining pressure to set the density of the system.
Unjamming by shear loading is also known as yielding. The question how much shear stress a material can sustain before it starts to flow is obviously im- portant for various fields, such as soil mechanics and civil engineering, but also for industry (e.g. food and cosmetics).
For a given system, not all axes of the jamming diagram need to be equally relevant. For example, if the particles are so large that they do not really feel the thermal fluctuations (as in sand), the temperature axis plays no role. However, there are other (non-equilibrium) ways to “excite” such systems, using shakers or air flow to drive the particles, and many efforts are being made to describe the resulting behavior in terms of effective temperatures. A self-consistent description of these systems using effective temperatures would contribute greatly to the unificationidea behind the jamming diagram.
1.2 Granular media
Granular media are materials that consist of many separate particles. These
“grains” are so large that they have no thermal motion: Their thermal ener- gies are much smaller than the gravitational or elastic energies involved in their behavior. Real-life examples are sand, rice, sugar, and so on, but physicists of course often prefer to study simpler systems such as spherical glass beads, or even quasi two-dimensional systems of cylinders between plates. Because the particles are so large, they have many internal degrees of freedom, which gener- ally cause them to behave very inelastically in collisions: When sugar is poured in a jar, the grains almost immediately come to a rest. The grains have no time to find their configuration of lowest energy, and therefore the system is always out of equilibrium. In principle there can be attractive forces between the particles,
Figure 1.2: Left: The force network in a disordered granular packing, obtained using photoelastic particles between crossed polarizers. Particles that carry a large load light up more brightly. Clearly the load is distributed very inhomogeneously. Figure courtesy of R.P. Behringer. Right: The force network in a numerical granular packing, obtained by the procedure described in section 2.2.1. The thicknesses of the lines represent the amount of force carried by the contact.
due to water bridges (wet granular media), Van der Waals forces (fine powders) or electrostatic charges, but in this thesis we consider only dry granular media, in which grains interact only repulsively.
The presentation of granular media as materials that can behave as “an un- usual gas”, “an unusual liquid”, and “an unusual solid” has by now become clas- sical [2]. Indeed, grains that are shaken around in a closed container at low density, colliding only occasionally, behave in some sense like a gas. It is not a normal gas, though, because energy is lost in every collision, which leads to clusteringof grains and to the fact that the system will come to rest pretty quickly after one stops shaking it. At higher densities, granular media sometimes flow like a liquid, and sometimes form static packings that resist stress like a solid.
The transition between these two is of course what we introduced in the pre- vious section as the jamming transition. Again, these states are pretty unusual compared to ordinary liquids and solids [2, 3]. In fast granular flows, the non- trivial behavior is usually attributed to dissipation (through inelastic collisions, or friction between the grains), or the effect of the air that is always present be- tween the particles. Slow granular flows are no less interesting, especially from the point of view of the jamming transition: A bucket full of glass beads behaves like a solid in the sense that it can support a heavy object. The surprise comes when the material is sheared very slowly deep beneath the surface: There is no observable motion of the grains at the surface, but the ability to support a load disappears: a steel ball placed on top now immediately sinks in [4]. In other words: the yield stress has become zero. This is clearly another example of a sit- uation where a system can go from solid-like to liquid-like without any apparent
structural changes.
1.2.1 Static granular media
Static granular media are interesting solids. When grains are confined under a static pressure, they all have to exert forces onto each other, and they do so in a rather inhomogeneous manner. A small fraction of the particles carries a large part of the load, as is seen from the force networks in figure 1.2. Naturally, these large forces line up a bit — because the total force on each grain is zero, a large force from one side implies a large force from the other side. This lin- ing up has been called the emergence of force chains, sometimes implying that these are long-range correlated force-bearing structures. However, in isotropi- cally compressed samples no definition of “force chains” has led to correlation lengths longer than a few grains.
The key quantity to characterize the inhomogeneity of the contact forces is the probability density of the magnitude of the contact forces, P (f). Consistent with the observation that some grains hardly carry any load, this distribution has a finite value for f → 0. In jammed systems it usually has a peak around the average valuehfi [5]. The statistics of the large forces has been a subject of debate: early experiments clearly found an exponential tail for P (f) [6, 7], as do some simple models like the q-model [8], but many numerical simulations found P (f )that bend slightly downward on a log-lin plot, indicating a slightly faster than exponential decay [5, 9–11]. The distinction is usually difficult to make, because the large forces are rare and hence it is hard to get good statistics for f >5hfi.
To obtain a continuum theory for stresses in granular media, one has to coarse grain the force network to obtain a stress field. It has been shown that a consis- tent elastic theory for granular media can be written down by using the proper definitions for stress and strain tensor [12]. Other models for stress propagation have been based on wave-like equations [13,14], usually focusing on force chains as the main actor. Recently it was proposed that the latter work better in small systems, while elasticity does a better job in large systems [15]. However, it has remained unclear how exactly the length scale that separates the two behaviors depends on the parameters of the system. For a more elaborate discussion on the nature of stress propagation in granular media see Ref. [16] and references therein.
1.3 Jamming in granular media
By studying jamming in granular media, we can expect, first of all, to learn some- thing new about granular media. However, we can even hope to learn things about jamming in general, because the jamming transition in granular model
(a)
jammed
unjammed
1/φ J
T
Σ
jammed
unjammed
(b)
potentialenergy
Ri+ Rj distance
i j
0 0
Figure 1.3: (a) The jamming diagram for dry granular media, with the transition point at zero temperature and zero applied stress indicated as “point J” [17]. (b) A typical interaction potential for dry granular media: particles i and j only interact when their centers are closer than the sum of their radii Ri+ Rj.
systems turns out to be very clean and well-defined, and understanding a com- plicated system often starts with understanding a conceptually simpler system that displays similar behavior.
Granular media are ideal systems to study the jamming transition. The grains do not interact when they do not touch, so that the number of contacts in the sys- tem is well-defined. When they touch, the interaction is repulsive. Therefore, the jamming transition at zero temperature and zero applied stress becomes sharply defined [17]. Consider a collection of frictionless spheres in a box. If the density is below a critical value φc, the particles fit in the box without touching each other, so that they do not interact, and the system is unjammed. Above φc, the grains cannot avoid touching, and start repelling one another (figure 1.3b). A pressure builds up and the system is jammed. Exactly at φc, the spheres just start to touch, but are not pressed into each other yet. For large systems φcconverges to a well-defined value that does not depend on the initial conditions of the pack- ing procedure. This makes the jamming transition on the density axis, marked
“point J” in the jamming diagram in figure 1.3, sharply defined.
Because at point J the particles have to touch but are not yet deformed, we can determine the average number of neighbors z that a grain must have. It is exactly the minimum number that is needed to be mechanically stable. If we have N spheres1 in d dimensions, there are dN coordinates that describe the positions of the particles. There are zN/2 contacts in the system, each of which gives an equation constraining the coordinates, in the sense that
|rj− ri| = Rj+ Ri for i and j in contact,
where ri and Ri are the position and radius of particle i. For this system of
1Here, we do not include rattlers, particles that happen to have no contacts at all, in the counting of N.
equations to have at least one solution, we must have z ≤ 2d. On the other hand, in the jammed system particles are in mechanical equilibrium: On each sphere there are dN equations of force balance, constraining the values of the contact forces. For frictionless contacts these forces are scalars, so there are zN/2 contact force variables. For this system of equations, we need z ≥ 2d to have a solution. At the jamming transition, where φ = φcand the pressure p = 0, both inequalities must hold and we get
z =2d . (1.3.1)
This value is called the isostatic contact number ziso, because it represents the minimum number of contacts needed to have mechanical equilibrium. Jammed packings, which have a nonzero confining pressure p, have z > ziso and are called hyperstatic. In this analysis we have assumed that the equations relating the positions are linearly independent. If this is not the case, packings can be hyperstatic even at zero pressure. This happens for example in a hexagonal lattice of particles in two dimensions, which has z = 6 while ziso=4.
The fact that packings at point J are isostatic has the consequence that break- ing a single contact is enough to make the system unstable: Each broken contact generates a zero-energy mode called floppy mode. This is easy to see: for particles with a radial potential energy to move freely, the distance between each pair of interacting particles must be unchanged. When z < ziso, there are less than dN of these pairs, which allows to choose the dN components of the displacement field (a vector uifor each particle) such that these distances indeed stay fixed. In jammed packings these floppy modes are not present, but it has been shown that their proximity is important for the properties of packings close to point J [18].
The numerical model systems studied by O’Hern et al. have an interaction po- tential between the grains that goes as the overlap δij = Ri+ Rj− |rj− ri| to some power. The two common choices are the three-dimensional Hertzian interaction V ∼ δ5/2 [19] and the harmonic interaction V ∼ δ2. However, when studying the properties of the packings as a function of distance to the jamming transi- tion, one obtains many scaling relations between them which are independent of the choice of interaction. The two most interesting ones, for which no simple explanation exists, are the scaling of the excess contact number ∆z ≡ z − ziso
and the scaling of the ratio of elastic constants G/K:
∆z ∼ (φ − φc)1/2 (1.3.2)
G/K ∼ (φ − φc)1/2, (1.3.3)
where G and K denote the shear modulus and bulk modulus, respectively. These relations deviate from what one would expect if the displacements of the parti- cles are affine, i.e. such that their coordinates change in a way which smoothly follows the global deformation. For a random set of points in a box, one would naively expect that changes in distance between the points would go linear in the change in density, and therefore ∆z ∼ (φ − φc). Similarly, assuming affine
displacements, the shear modulus and bulk modulus should be of the same or- der, and G/K should be independent of φ. These models do not work because the deformations are not affine, but a clear characterization of the nonaffinity is lacking, and therefore the microscopic origin of the scaling relations is unclear.
The fact that the observed scaling relations are power laws reminds us of critical phenomena. This has attracted a lot of attention, because if point J really is analogous to a critical point2, it may have an influence on a large part of the jamming diagram around it. An important aspect of this discussion is the existence of a length scale that diverges upon approaching the transition. Such a length scale has been identified in the vibrational properties of these packings, first through a crossover frequency in the density of states [20, 21], and later by looking at the spectral properties of the eigenmodes [22]. This length scale diverges as
ℓ∗∼ 1
∆z ,
and has been suggested to be the length beyond which one can treat the packing as an elastic continuum [21]. We will show in chapter 2 of this thesis that this is indeed the case. Therefore, it is definitely also of interest to those in the gran- ular community who worry about the nature of stress transmission in granular media [16].
1.4 This thesis
In this thesis we study the jamming behavior of granular media, using two differ- ent approaches: linear response and the force network ensemble. Both are based on the idea to simplify the problems, which makes them more tractable and often allows to disentangle various effects which would happen all at the same time if the full problem is studied. Let us first look at the basic properties of these tools, and then explain the organization of the thesis.
1.4.1 Linear response
The changes in a granular medium that result from a very small perturbation are given by its linear response. This means all interactions between grains are lin- earized around a starting configuration, so that effects of e.g. breaking or creation of contacts are ignored. The main advantage with respect to full-scale numeri- cal simulations (Molecular Dynamics) is that it is numerically cheap, because for each calculation only one matrix equation needs to be solved. Furthermore, by comparing the results to those obtained using other methods, one can determine whether a particular phenomenon can or cannot be caused by nonlinear effects.
2At least, to a certain extent — remember that the jamming transition is not a thermodynamic transition.
The linearization of the interparticle interactions leads to a matrixM, known in condensed matter physics as the dynamical matrix, which contains all infor- mation about the geometry of the contact network, the stiffnesses of the springs that (after the linearization) describe the contact interaction, and the effect of pre-stress3 in the system. It is a square matrix with dimension dN for friction- less spheres (see chapter 2) and d(d + 1)N/2 for frictional (rough) spheres (see chapter 3). It has two main uses, which I will introduce shortly here, and which will be described in more detail when they are used in chapters 2 and 3.
The first one is a many-particle equivalent of Hooke’s law:
fext=Mu , (1.4.1)
where fextis an externally applied force, and u contains the displacements of all particles. Solving this equation tells us how all particles move as a result of the external force. Then, using the interaction potential, we obtain the changes in the contact forces from these displacements, which in turn allow to determine the changes in the stress tensor.
The second use of the dynamical matrix is a many-particle equivalent of the equation of motion of a harmonic oscillator:
m¨u +Mu = 0 , (1.4.2)
where m¨uis a vector that contains the product of mass and acceleration for each particle. From this equation one obtains the eigenfrequencies and eigenmodes of the packing.
Even though both types of calculation use the same linearization, we will use the term linear response equation only for equation (1.4.1).
1.4.2 Force network ensemble
The force network ensemble is a recently introduced tool to study the statistics of forces in granular media [23]. The first application was to calculate the force distribution P (f) for several packings of discs. The resulting P (f) had all the fea- tures that are typically observed in experiments and numerics, and subsequently the model has received a lot of attention [24–27].
The motivation behind the ensemble is that packings of relatively hard parti- cles show a separation of length scales: the deformations corresponding to typ- ical forces in the packing are orders of magnitude smaller than the sizes of the particles. Exploiting this scale separation, the positions of the particles are kept fixed, and one analyzes the ensemble of force networks that satisfy mechanical equilibrium on each grain. On the one hand, this constitutes a considerable in- crease in physical relevance with respect to models for force statistics that do not incorporate the local conditions of mechanical equilibrium. On the other hand,
3I.e. the forces that the particles exert on each other in the starting configuration.
by fixing the positions of the particles the ensemble can be used to disentangle complex phenomena in granular media, such as yielding, in which the forces and the contact network are evolving at the same time.
As was explained in section 1.3, jammed packings are hyperstatic, which means they have more contacts per particle on average than the isostatic value ziso. Therefore, the conditions of force balance are not sufficient to determine all contact forces, and we get a ∆zN/2-dimensional solution space of force net- works that satisfy force balance on each grain. The force balance conditions read simply
X
j(i)
fji=0 ,
where the sum is over all neighbors j of particle i and fjiis the force that j exerts on i. They are usually written in matrix form as
A ~f = ~b ,
where ~f is a vector of zN/2 components representing the contact forces, and
~b is a vector of dN zeroes. It is important to note that in dry granular media all forces are repulsive, which is incorporated by demanding all components of f~ to be positive. The solution space can then be bounded by adding additional equations toA that fix the pressure in the system, or even each component of the stress tensor separately. This stress constraint has a nonzero right-hand side, so that the general solution is
f = ~~ f0+
∆zN/2
X
k=1
αk~v(k), (1.4.3)
where ~v(k) is the set of solutions toA ~f = ~0. The coefficients αk are restricted by the condition that all components of ~f need to be positive. If the contact network thatA corresponds to was made by compressing real repulsive particles, it is often convenient to choose the particular solution ~f0 to be the actual force network of this real system.
Finally, to turn this “solution space” into an “ensemble”, we have to define a measure. There is no a priori obvious measure for this problem, so a flat measure is being used [23]. In the spirit of the Edwards ensemble, all balanced force networks are taken to be equally probable. Because of the fixed contact geometry, however, the force network ensemble is much more restricted than the Edwards ensemble, which considers “all blocked states” [28]. Using the force network ensemble with this flat measure, realistic P (f) have been obtained, both on hexagonally ordered packings and on disordered packings [23].
1.4.3 Organization
We study the jamming transition from various angles. Throughout this thesis, a crucial role is played by the distance to isostaticity, ∆z≡ z − ziso.
In chapter 2, we study the jamming transition in frictionless granular media, the basics of which were summarized in section 1.3, to uncover the microscopic origin of the nontrivial scaling relations between the confining pressure, elastic moduli, excess contact number, and the length scale ℓ∗. We determine the linear response of these packings to local forces, applied to a single particle, and to global deformations. We characterize the nonaffinity of the response by a local measure of the correlations between neighboring particles, the displacement an- gle distribution, and show that the nonaffinity has a strong signature of floppy modes. This leads to the observation that particles in contact tend to slide past each other when the system is deformed, by an amount that diverges as the jamming transition is approached. We discuss how these sliding motions are re- sponsible for the creation of new contacts when the system is compressed, and how they should be limited by finite size effects. Furthermore, we show that granular media, when stress fields are averaged over many packings, are well- described by linear elasticity theory. In a single packing, however, the response is dominated by the local disorder of the packing, and one needs to go to length scales larger than ℓ∗to see a continuum-like response.
To show that it is really the distance to isostaticity that controls the scaling re- lations, we study jamming in frictional packings in chapter 3. Frictional packings are not automatically isostatic at jamming: for moderate values of the micro- scopic friction coefficient, ∆z goes to a constant as the pressure goes to zero. We determine the crossover frequency ω∗ in the density of states, and the ratio of elastic constants G/K, for a range of pressures and friction coefficients. Both of these quantities are known to scale with ∆z in frictionless systems. We find that they are proportional to ∆z in these frictional systems as well, and that they do not scale with, for example, the confining pressure.
In chapter 4, we study the mathematical structure of the force network en- semble, by applying it to simple hexagonally ordered packings. This chapter is not so much related to jamming, but more to statistics of forces in granular me- dia. There is therefore no other role for ∆z than to determine the dimension of the force space. We develop a general method to solve for the force statistics of a packing using the ensemble, which we apply to several small regular packings.
Larger packings are treated numerically, using a simulated annealing procedure.
Already for small packings, we find the typical features of the force distribution P (f ): a finite value at f = 0, and a small peak close to f =hfi, and a tail that is in between a Gaussian and an exponential. To determine how special the ma- tricesA have to be to obtain this sort of force statistics, we apply the ensemble sampling to several types of random matrices and perturbed packing matrices, and find that the typical features of real P (f) quickly disappear.
In chapter 5, we study force networks under shear stress, using the force network ensemble on disordered contact geometries for various values of ∆z.
The maximum shear stress that a packing can sustain is important for yielding of granular materials. The contact geometries that we use as input for the ensemble are isotropic, so that we can disentangle the role of force anisotropies from that
of contact (also called fabric) anisotropies. We study the anisotropy of the force network, starting from the angle-resolved force statistics Pφ(f ). We relate the maximum possible value of the shear stress to the size of the force space, and find that this maximum is proportional to ∆z: strongly hyperstatic contact networks have a much larger force space than those with z close to ziso, and therefore they can accommodate more shear stress. Strongly hyperstatic packings develop a force anisotropy which gets very close to a bound that we derive analytically.
In chapter 6 we derive this bound in more detail, and generalize it to the more general case of frictional particles and anisotropic contact networks.
Linear response of frictionless
granular packings
abstract
We study the microscopic origin of the scaling behavior in frictionless granular media above the jamming transition by analyzing their linear response. The response to local forcing is non-self-averaging and fluctuates over a length scale that diverges at the jamming transition. The response to global deformations becomes increasingly nonaffine near the jamming transition. This is due to the proximity of floppy modes, the influence of which we characterize by the local relative displacements, using the displacement angle distribution. We show how the local response governs the anomalous scaling of elastic constants and contact number.
2.1 Introduction
As we described in the general introduction, many disordered materials, such as foams, colloidal suspensions, granular media and glasses, display a jamming transition between solid-like and fluid-like states that is typically controlled by density, temperature and the amount of shear stress in the system [1, 29]. Gran- ular media have been found to be ideal model systems to study jamming at zero temperature, because they consist of large, athermal particles that easily lose kinetic energy in collisions [17]. In collections of frictionless spherical particles that interact through repulsive contacts, the jamming transition as a function of the packing fraction φ is sharply defined: In the limit of large system sizes, a well- defined packing fraction φc(marked as “point J” in the jamming diagram) sep- arates the regime of jammed packings that resist shear-loading from the regime of loose particles.
A crucial aspect of packings of frictionless spheres is the fact that they are marginally stable (isostatic) at point J: the system has just enough contacts to be mechanically stable (z = ziso = 2d in d dimensions). From there, any bro- ken contact will give rise to a zero-energy displacement mode known as a floppy mode [30], as we explained in section 1.3. The vibrational properties of jammed packings (away from point J) have been shown to be strongly influenced by these floppy modes: there is an excess of low-frequency vibrational modes, which shows up in the density of states as a plateau that extends down to a crossover frequency ω∗∼ ∆z [20, 21], where ∆z = z − ziso.
The mechanical and geometric properties of jammed packings close to point J were found to exhibit power law scaling as a function of φ− φc[17]. These scal- ings suggest a critical nature of the jamming transition, so that properties of the critical point J may control a large part of the jamming diagram around it [17].
This is of interest to several communities: Firstly, because understanding what happens at point J may improve our understanding of the mechanical properties granular media themselves. Secondly, we may learn new things about jamming in general by first studying this rather clean model system.
Understanding the microscopic origin of the scaling relations near the jam- ming transition has remained a great challenge. Under the assumption that these packings deform affinely, some of them have simple explanations. Other relations have proven to be harder to understand, in particular the scaling of the excess contact number ∆z and the ratio of elastic constants G/K. The puzzle becomes even bigger if one realizes that, close to point J, deformations are not at all affine, so that even the “simple” relations may have to be revisited.
In this chapter, we investigate the microscopic origin of the scaling of the elastic bulk modulus K, shear modulus G, and excess contact number ∆z, by calculating the linear response to external perturbations. We show that these are
strongly related to the nonaffine deformations, and we introduce the displacement angle distributionto quantify the correlations between neighboring particles. The nonaffinity is such that particles predominantly slide past each other, indicating a growing influence of floppy modes as the jamming transition is approached.
Furthermore, we show that the response of granular packings to a local point force can be described by linear elasticity. The correspondence is very good when an average over many packings is considered; in a single packing the response is dominated by packing disorder up to a length scale ℓ∗. This length scale, which diverges upon approaching point J, was first introduced to explain the excess of low-frequency vibrational modes [21], and observed in the spatial structure of these modes [22]. We show that it sets the scale above which a granular medium can be treated as a continuum, both from the disorder in the changes in the contact forces in linear response to a localized perturbation, and from an analysis of the coarse graining length needed to obtain a smooth stress tensor in linear response to a global deformation.
The chapter is organized as follows: We start with an introduction of our numerical tools, explaining how we prepare the packings and how we calcu- late their linear response. In section 2.3 we show that the linear response to a point force can consistently be described in terms of linear elasticity. Section 2.4 deals with the microscopic details of the response, introducing the displacement angle distribution, and relating the elastic moduli to the typical relative displace- ments of neighboring grains. There we also discuss how these displacements are responsible for the scaling of ∆z, by considering the response to a global com- pression. In section 2.5 we present evidence for the diverging length scale ℓ∗and discuss possible implications for finite size effects. We conclude with a discussion on the scaling of the elastic moduli and on an interpretation of P (α) in terms of steric hindrance, including implications for experiments on disordered media close to the jamming transition.
2.2 Numerical procedure
All numerical results presented in this paper concern the linear response of a given starting configuration to a small perturbation. It is therefore explicitly a two-step problem. First, we prepare a packing at a given density or confin- ing pressure, using a Molecular Dynamics simulation (section 2.2.1). Then, we analyze its properties using calculations of the linear response to perturbations (section 2.2.2).
2.2.1 Packing generation
The starting configuration that will later serve as input for the linear response calculations, is prepared using a molecular dynamics simulation with discs in
two dimensions. The discs interact via the three-dimensional Hertzian potential,
Vij = (
ǫij
1−rσijij5/2
rij≤ σij
0 rij≥ σij
, (2.2.1)
where i, j label the particles, σij is the sum of the particle radii Ri+ Rj, and rij
is the distance between the particle centers. The energy scale ǫijdepends on the radii and effective Young moduli of the particles, see Appendix 2.A. We will refer to the quantity between brackets as the dimensionless overlap, δij=1− rij/σij.
Using the 3D Hertzian potential makes the system a quasi-2D packing of discs with round edges. We use zero gravity to have a homogeneous pressure, and the radii of the discs are drawn from a flat distribution between 0.4 < R < 0.6, thus creating a polydispersity of±20% around the average particle diameter d = 1 (our unit of length) to avoid crystallization. The simulation starts from a loose granular gas in a square box with periodic boundaries, which is compressed to a target pressure p by changing the radii of the particles while they are moving around. The radii are multiplied by a common scale factor rs, which evolves in time via the damped equation ¨rs=−4ω0˙rs− ω20[p(t, rs)/p− 1]rs, where ω0∼ 6 · 10−2, p(t, rs)is the instantaneous value of the pressure and p the target pressure.
This ensures a very gentle equilibration of the packings. Energy is dissipated through a drag force and inelastic collisions. The simulation stops when the accelerations of all grains have dropped below a threshold which is 10−10in our reduced units. This way we generate packings in mechanical equilibrium for pressures ranging from p = 10−6to p = 3· 10−2, in units of the modified Young modulus of the individual grains [19]. For one particular type of calculation we use packings which are only periodic in the x-direction, and have hard walls on top and bottom. These are generated by compressing the system using a hard piston, as described in Ref. [31].
Because there is no gravity in our simulations, at the end of the simulation there will usually be particles without neighbors, or due to numerical precision effects, particles with fewer neighbors than is needed to make them rigidly con- nected to the rest of the packing. These rattlers or floaters do not contribute to the rigid backbone of the packing. They are removed from the system when determining the contact number z or calculating its linear response to small per- turbations.
2.2.2 Linear response
We calculate the response of a packing to a mechanical perturbation, which can be either an external force or a deformation of the periodic box. The response in general involves displacements of all particles in the packing, which we analyze in the linear regime.
Linear response equation
The total potential energy of the system is a sum over all pairs of particles in contact,
E =X
hi,ji
Vij(rij), (2.2.2)
where we made the dependence on the relative positions of the particles rij = rj− riexplicit. The change in rij due to displacements uiof the particles is
∆rij =q
(rij+ uk,ij)2+ u2⊥,ij− rij, (2.2.3) where ukis the relative displacement along the line connecting the centers of the particles and u⊥is the relative displacement perpendicular to that. In the linear regime the change in energy is expanded in these relative particle displacements as
∆E = 1 2
X
hi,ji
kij
u2k,ij− fij kijriju2⊥,ij
. (2.2.4)
Here fij = −dVij/drij and kij = d2Vij/drij2 so that both the initial force f and the stiffness k are nonnegative for our potential (2.2.1). There are no terms linear in the displacements because the starting configuration is in mechanical equilibrium, which makes them sum to zero. The u2k-term represents the sim- ple stretching of the linearized spring. The u2⊥-term is only there if there is a pre-stress or initial force and corresponds to a lowering of the energy due to a transverse displacement of the particles.
For any potential which has the form V ∼ δαij, the prefactor of the second term in equation (2.2.4) can be written as δij/(α− 1). The closer to point J, the smaller these dimensionless overlaps, so this prefactor will be small close to point J. However, this does not allow us to ignore this term, because, as we will see, the typical u⊥are going to be much larger than the typical uk, in the limit of approaching the jamming transition.
Writing the change in energy in the independent variables of the problem, the displacements uiof the particles, defines the dynamical matrixM:
∆E =1
2Mij,αβui,αuj,β, (2.2.5)
whereM is a dN × dN matrix, the indices α, β label the coordinate axes and we use the summation convention. The dynamical matrix contains all information on the elastic properties of the system. It can be diagonalized to study the vibra- tional properties [17, 32], and it can be used to study the response of the system to an external force (defined on each particle) [33]
Mij,αβuj,β= fi,αext. (2.2.6)
The dynamical matrixM is very sparse for large systems, because the only non- zero elements are those for which i and j are in contact and those for which i = j. Therefore we can compute the response efficiently by using the Conju- gate Gradient Method [34]. This method has been used before by O’Hern et al.
to study the quasistatic response of granular systems to a global shear or com- pression [17]. The difference between our numerical approach and theirs is that we are considering the linear regime where changes in the contact topology are ignored. In this regime calculations are numerically cheaper, allowing to obtain averages over large ensembles, while the main results for the scaling behavior are recovered (see section 2.3.2).
There has of course also been a lot of work on the response of granular ma- terials using full Molecular Dynamics (also called Discrete Element Methods in this field, to distinguish it from continuum modeling) [15, 35–37]. Such simu- lations probe the actual dynamics of the system, where grains have inertia and contacts are formed and broken all the time. Whereas these methods are prob- ably essential when simulating granular flows, the basic physics of quasistatic deformations of jammed granular media is captured well within the linear re- sponse. In addition, one of the key quantities we use to analyze the response to bulk deformations, the relative sliding of neighboring particles, is extracted straightforwardly from the linear response, but requires a more subtle treatment when used on MD data.
Forces and stresses
Solving the linear response equation (2.2.6) yields the displacements ui of all particles. From this we calculate the local relative displacements u⊥,ij, uk,ij, and the change in contact force ∆fij using
uk,ij = cos φijuij,x+sin φijuij,y (2.2.7) u⊥,ij = cos φijuij,y− sin φijuij,x (2.2.8)
∆fij = −kuk,ij , (2.2.9)
where φij is the angle the bond vector rij makes with the x-axis. The change in the stress tensor is obtained using a local coarse graining [12]:
∆σαβ(r) =X
hi,ji
∆fij,αrij,β
Z1 0
ds Φ r− rj+ srij , (2.2.10)
where we use a Gaussian coarse graining function Φ with a standard deviation ξCG, which in our case is of the order of one particle diameter. Taking Φ = 1/V in this expression (where V is the volume of the system) gives the well- known expression for the average of the change of the stress tensor over the entire system.
2.3 Macroscopics and continuum elasticity
There has been a long-standing debate about the use of continuum elasticity to describe the behavior of granular materials under small deformations. The views on stress propagation in granular media can roughly be divided in two classes:
wave-like approaches, often focusing on force chains and their propagation [13, 38–41], and approaches focusing on continuum elasticity [12, 15, 37, 41–44]1.
In this section we show that the changes in the stress field of a granular packing in response to a directional point force are, after averaging over an en- semble of a few hundred packings, well described by linear continuum elasticity.
The elastic moduli used in the fitting process are equal to those obtained from a global shear or compression of a granular packing. We do not consider spa- tial variations in the elastic moduli, assuming the elasticity tensor is reasonably homogeneous. All elastic moduli in this section are therefore spatially averaged moduli.
2.3.1 Point response
We calculate the linear response of packings of N = 10 000 particles to a force in the y-direction. The packings are periodic in the x-direction and have hard walls on top and bottom to carry the load. The confining pressure used to create the packings ranges from p = 10−6to p = 10−2, corresponding to contact numbers ranging from z = 4.05 to z = 5.41. We have 10 different packings at each pressure, and use each one about 20 times by applying a point force to different particles, all of which are closer than 0.1 to the line that marks the center of the system. Examples of the resulting response networks, depicting the changes in the contact forces, are shown in figure 2.2.
The continuum solution that we are comparing the granular point response to, is obtained from the static Navier-Cauchy equation [45]
G∆u + K∇(∇ · u) = −fext. (2.3.1) Here K and G denote bulk and shear modulus and we use a two-dimensional formulation of elasticity theory. This equation is the direct equivalent of the lin- ear response equation (2.2.6), because solving it yields displacements in terms of external forces. We solve equation (2.3.1) for fextequal to a unit point force in the y-direction, see Appendix 2.B for details. The resulting stress field σαβ(r) only depends on the ratio K/G of elastic constants, since the overall scale drops out when relating an imposed force to a resulting stress. Therefore, there is only one free parameter when fitting the stress field of the granular system, equa- tion (2.2.10), to the continuum expression. As described in Appendix 2.B, we use the the effective Poisson ratio ν = (K− G)/(K + G). To fit the overall scale
1See also Ref. [16] and references therein.
σxx
σxx σxx
σxx
σxy
σxy σxy σxy
σyy
σyy σyy σyy
(a)
(b)
(c)
(d)
Figure 2.1:Stress fields for the response to a point force. The greyscale images represent the ensemble averaged granular stress field; the thick contours denote the fitted contin- uum stress field. For the σyy thin contours for the granular stress field are included as well. The components of the stress tensor are plotted for (a) p = 10−2, (b) p = 10−4, (c) p =10−6. (d) Same as (c), but when leaving the pre-stress term (∼ u2⊥)in the energy expansion (2.2.4) out of the calculation.
p =10−2 p =10−4 p =10−6
Figure 2.2:Force response networks for a point loading with pressure as indicated. Blue (red) lines indicate positive (negative) changes in contact force, the thickness indicating the amount. The particles themselves are not drawn. The figures show only the part of the packing close to the forced grain, and contain about 1800 of the 10 000 particles.
of the moduli, we fit the displacement field. In particular, we fit the the average y-displacement of all particles in a strip at height y,
Hgran(y):=ui,y
yi≈y , (2.3.2)
to its continuum counterpart:
H(y):= 1 Lx
ZLx/2
−Lx/2
uy(x, y)dx = Ly− 2|y|
4Lx(K + G), (2.3.3) which is obtained by taking fext = δ(x)δ(y)ˆyand integrating equation (2.3.1) over x, leaving a standard Green’s function problem with the boundary condition that uy vanishes at the top and bottom wall, i.e. at y =±Ly/2.
Figure 2.1 displays the ensemble averaged stress fields ∆σαβ(r) (equa- tion (2.2.10)) and the fitted continuum stress fields for p = 10−2, p = 10−4, and p = 10−6. To illustrate the effect of the pre-stress term in the energy expan- sion (2.2.4), we also show the stress fields for a different calculation at p = 10−6 where we explicitly left out the u2⊥-term. This term is strictly negative in a sys- tem of purely repulsive particles, and therefore tends to destabilize the system and lower the energy associated with deformations. The destabilization is seen from comparing figure 2.1c with figure 2.1d: the spatial fluctuations from the continuum elasticity stress fields are much smaller if we leave out the pre-stress term. The lowering of the energy can also be shown by calculating the elastic moduli without the pre-stress term: the resulting shear modulus is higher than the real shear modulus. We will come back to the relative importance of the two terms in the energy expansion in section 2.4.
The fits of the average vertical displacements Hgran(y) are shown in fig- ure 2.3; the resulting elastic moduli are presented in figure 2.4.
y/Ly
H(y)·(K+G)
0.0 0.2 0.4 0.6
0.0 0.2 0.4
−0.2
−0.4
Figure 2.3:The fitting procedure used to obtain K +G: The average vertical displacement of particles at height y, equation (2.3.2), is fitted to the triangle-shaped function H(y), given by equation (2.3.3). The functions are rescaled by the fitted K + G and vertically offset for clarity.
2.3.2 Bulk deformations
A more conventional way to extract elastic constants is to apply a compression or shear deformation to the entire system. In packings with periodic boundaries this is done by imposing a relative displacement on all bonds that cross the boundary of the periodic box. For example, to impose a globally uniform compressional strain ǫxx = ǫyy = ǫ, all terms in the energy expansion of equation (2.2.4) that represent a bond that crosses the periodic x-boundary are changed such that each occurrence of uj− uiis replaced by uj− ui± ǫLxˆx, where the±-sign is given by the sign of rij,x. The y-boundary is treated analogously. Writing the linear response equation (2.2.6) from this energy expression for the deformed system, the terms proportional to ǫ end up on the right hand side of the equation and act like an effective fext. The response of the system to this shape or volume change of the periodic box is again calculated by solving equation (2.2.6) for this effective external force. The elastic moduli then follow from the resulting change in stress tensor according to equation (2.2.10) with the trivial coarse graining function Φ = 1/V . The bulk modulus is extracted from a uniform compressional strain ǫxx= ǫyy= ǫthrough
K = σαα
2ǫαα
= σxx+ σyy
4ǫ , (2.3.4)
and the shear modulus from a uniform shear strain ǫxy= ǫthrough G = σxy
2ǫxy
= σxy
2ǫ . (2.3.5)
10−5 10−4 10−3 10−2 10−1 1
10−710−610−510−410−310−210−1
p ∆z
K,G G/K
0.01 0.1 1
0.01 0.1 1 10
Figure 2.4: (a) Bulk modulus K and shear modulus G as a function of pressure. The squares (K) and diamonds (G) are obtained using the bulk response described in this section. The crosses (K) and plus signs (G) follow from the fits of the point response stresses and displacements from the previous section. The error bars on the bulk response data indicate the intervals spanned by the median 50% of the data; the actual standard error of the mean is much smaller than the symbol size. The error bars on the point response data are error estimates from the fitting procedure explained in section 2.3.1.
(b) The ratio of elastic moduli G/K scales approximately as ∆z.
The shear deformation is applied in the form of a pure shear, i.e. by having a displacement in the y-direction imposed on all bonds that cross the x-boundary andvice versa.
Using this procedure, we calculate the elastic moduli for packings of N = 1000 particles at pressures ranging from p = 5· 10−6to p = 3· 10−2, and average over 100 different packings at each pressure. The result is shown in figure 2.4, together with the moduli obtained in the previous section, from fitting the point response. All data points except one fall on the same scaling curves
K ∼ p0.36±0.03 (2.3.6)
G ∼ p0.70±0.08, (2.3.7)
which shows that we can consistently describe the response in terms of con- tinuum elasticity. To our knowledge, this is the first time the connection be- tween the localized response and independently obtained elastic moduli has been shown explicitly. The deviating point is the lowest pressure data point from the point response. This is due to the fitting procedure: we fit the Poisson ratio ν from the stress field ∆σij(r). As p→ 0, we observe ν → 1, its maximum value in two dimensions. In this limit, the fitting of the stress tensor becomes less accu- rate for numerical reasons, which is made even worse by the fact that the data become more noisy at low p. In the bulk response we see a similar effect: the fluctuations around the average values are much larger at lower pressure.
A simple estimate of the elastic moduli of disordered solids relies on the as- sumption that the displacements of the particles in the packing are affine. This
