The role of force networks in granular materials

One of the most visually striking features of granular materials is their heterogeneous pattern of force transmission, commonly known as force chains. In this paper and its associated talk, we will review several experiments on two-dimensional photoelastic granular materials which highlight methods for making quantitative interparticle forces measurements, theoretical frameworks for interpreting such data, and the specific findings which result from these methods.


Introduction
Granular materials are inherently heterogeneous, and continuum models of properties such as the shear modulus and sound speed often fail to quantitatively capture their dynamics.This lack of understanding has serious consequences for the engineering of both geotechnical applications, and the industrial storage and handling of granular materials.One likely reason for these difficulties is that boundary stresses are transmitted via a chain-like network of strong forces, as shown in Fig. 1.
This paper presents several complementary methods of characterizing these systems: acoustic transmission at the particle-and network-scale [1,2], the force-moment ensemble [3], and network science [4,5].For each case, we provide a description of both the method itself, and the type of information that can be gained through that approach.Collectively, these measurements highlight the importance of understanding the role of the force chain network in controlling the bulk properties of the material.In this paper and associated talk, we will primarily focus on 2D experiments on photoelastic materials, since these readily allow for the quantitative determination of the particle-scale forces [6].The framework is general enough to be applied to 3D materials as similar data becomes available.

Photoelastic force measurements
The use of photoelastic materials, as shown in Fig. 1, allows for the quantification of the forces in a granular system by placing it between a pair of polarizing filters.In our experiments, we use darkfield photoelasticity: polarized light will only be be transmitted through the second polarizer if the stress on a particular grain has rotated the polarization of the light by an integer multiple of π/2.In the sample images shown here, dark particles are under e-mail: kdaniel@ncsu.edulow stress, and bright particles are under higher stress with an increasing pattern of dark/bright fringes as the stress increases.By using circular, rather than linear, polarized light, the resulting fringe pattern is independent of the relative angle between the applied stresses and the polarizer orientation.
These methods can be used either semi-quantitatively to measure the average pressure, or quantitatively to determine the vector force at each interparticle contact.Our experiments utilize either polyurethane or Vishay Photo-Stress particles, cut into centimeter-scale disks/ellipses using either a CNC mill or a water-jet.
By using a combination of unpolarized and polarized light (see Fig. 2), it is possible to obtain a snapshot of the entire state of the system.In the results presented below, we determine the positions of all particles from images taken in unpolarized light, and the pressure or contact forces from images taken in polarized, monochromatic light.In the most recent work described, this was done via LED lighting.In a single image, we obtain position-data in the red channel of a color camera, and the force-data in the green channel.An open-source code for making vector contact force measurements, along with a review of the underlying theory and practical advice for constructing a photoelastic setup, is described in an upcoming review article [7].

Acoustics
By using a combination of high-resolution, low-speed imaging (to measure the stress on each particle) and lowresolution, high-speed imaging (to measure the amplitude of sound waves), we can examine the relationship between the force chain network and the spatial pattern of sound propagation.To measure the local amplitude of the propagating sound wave, we calibrated changes in the photoelastic signal against measurements from a small number of particle-scale piezoelectric sensors known to have a linear response.To measure the pressure on individual particles, we located all particles within a high-resolution image and used a combination of the average intensity I, the squared average intensity gradient |∇I| 2 [8], and the position of photoelastic fringes in the vicinity of each contact to provide an empirical estimate of the normal component of the contact force.
As shown in Fig. 3, we observe that the amplitude A of the acoustic signal is on average largest within particles experiencing the largest forces, with an empirical relationship that is approximately a power law A ∝ P 0.35 .
For Hertzian contacts, the interparticle contact force is given by F ∝ δ β , where δ is the distance (measured from the edge) that each particle is compressed, and the exponent β depends on the particle geometry.The contact area between cylindrical particles is given by a ∝ F 1 2β , for a rectangular contact area.For our particles, we observe that the interparticle force obeys a Hertzian-like contact force law F ∝ δ 5/4 , a weak nonlinearity likely arising due to out-of-plane and asperity effects.For a Hertzian exponent of β = 5/4, the contact area between cylinders therefore grows as P 2/5 , consistent with the observed pressure dependence in Fig. 3.We therefore conclude that increases   3. Data collected at the particle-scale (small gray symbols), averaged to reveal a relationship between the pressure on a particle and the amplitude of sound propagation through that particle.The dashed line is ∝ P 0.35 , from a fit to all points.The large, dark circles are the average and standard error of this dataset, calculated for equally-populated bins.  in interparticle contact areas along force chains provides the main conduit for sound propagation.
Importantly, the force chain network does not remain a static structure during the propagation of sound.As shown in Fig. 4, the deformation of individual particles, for instance due to expansion via the Poisson ratio, can open and close individual contacts.This is exemplified by the particles labeled 1 and 2, which are clearly not touching in the static system (panel a, see elliptical region), yet transiently form a force chain during sound propagation (panel b, same elliptical region marked).For ν ≈ 1/2, the lateral expansion will be of the same order of magnitude as the compression δ; this is consistent with closing the gap shown.
Traditional discrete element methods, in which all particles in a simulation remain circular and "overlap" each other rather than deforming, are unlikely to be able to capture such effects.In fact, no preference for the sound waves to follow the force chains was seen in simulations that investigated the question [9], even though the effect would have been larger in that case (spherical particles with β = 3/2) than here.Future work will be necessary to resolve this discrepancy.

Athermal Statistical Ensembles
A key technique of statistical physics is to represent a system by an ensemble of allowed microstates, each occurring with a probability that depends on a small number of macroscopic state variables such as temperature, pressure, volume, and energy.Because granular materials exhibit reproducible ensembles of microscopic quantities, wellpredicted by a few macroscopic parameters, these methods are a promising route to describing the state of a granular system.Importantly, these methods seem to be successful even though (1) the thermodynamic temperature is too small (by a factor of 10 12 or so) to play a role in the state of a granular material, and (2) granular materials are both frictional and inelastic, ensuring that they are inherently non-equilibrium.
The Edwards ensemble and its later variations, including the force-moment ensemble, have gained significant attention as a promising approach to writing an equation of state for static granular materials [3].A common feature of these ensembles is that they introduce a temperaturelike variable of the form 1  T ≡ ∂S ∂E , where the entropy S is replaced by a configurational entropy and the temperature and energy (T, E) are replaced by a different pair of intensive (equilibrating) and extensive (additive) state variables.There has been considerable interest in the statistics of the locally-defined force moment tensor Σ, examining the extent to which it provides an extensive variable analogous to energy in equilibrium statistical mechanics.It is calculated by summing the outer products of the moment arm d and the force f between the particles (m) and their contacting neighbors (n): In two dimensions, it is convenient to decompose Σ into σ p = (σ 1 + σ 2 )/2 and the σ τ = (σ 1 − σ 2 )/2, where σ 1,2 are the principal stresses.The average hydrostatic pressure in the system is given by Γ N = Tr Σ/N, taken over clusters of N particles.Note that this formulation in terms of the force-moment tensor is related to the more familiar Cauchy stress tensor σ by the relation Σ ≡ V σ, where V is the total volume occupied by the particles [3].
In the force-moment ensemble, the corresponding temperature-like variable is known as the angoricity [10] where i = (p, τ) for either the pressure or shear components.The equivalent of the Boltzmann distribution would be given by Equilibration of the two components of the temperature-like variable angoricity, in experiments on a lowfriction subsystem (red ) placed within a high-friction bath (black •).The low-friction particles are monitored for location via markings visible in blacklight, and are returned to the subsystem if they stray out of it; no mechanical confinement is present.Values are measured as averages over clusters of N = 8 particles, to reduce statistical fluctuations [11].Pressure Γ is measured in units of [N/m].Figure adapted from [12].
where an exponential distribution is expected [11] instead of σ 2 i .We make measurements of the relative-angoricity of the system using the method of overlapping histograms.For a stress state σ (either pressure or shear), consider a measurement of P(σ) for the same system prepared at two different values of α, enumerated (1,2).The ratio of these histograms is conjectured to be given by [3] where (α, σ) can be either the pressure (p) or shear (τ) components.By taking logarithm of the ratio, and fitting the data to determine the slope, it is possible to measure changes in the angoricity.Since 1/α = 0 at σ = 0, this can also be used as an absolute measurement.
We have conducted laboratory tests of the validity of angoricity, examining the extent to which these ensemblebased formalisms provide well-defined state variables.For these experiments, we use a quasi 2D granular material levitated on an air cushion to eliminate basal friction.By analogy with the canonical ensemble in conventional statistical mechanics, we consider whether state variables equilibrate between a smaller subsystem and a larger bath.For a low-friction subsystem within a highfriction bath, subject to biaxial compression, we observed that the temperature-like variable compactivity (associated with the system volume in the original Edwards ensemble [13]) fails to equilibrate [12].This is a failure of the zeroth law of thermodynamics for this ensemble.However, for the force-moment ensemble, both α p and α τ are observed to equilibrate [12], as shown in Figure 5.
For a well-defined state variable, we would expect that the equation of state (1/α ∝ Γ) would be independent of loading history.To test this hypothesis, we performed experiments under three loading conditions -biaxial compression, uniaxial compression, and simple shear -as shown in Fig. 6.In the bottom half of the figure, we  report the volume and pressure measured over 8-particle clusters; these already show a strong history-dependence.In calculating the angoricities according to Eq. 2 and 3, we find that while σ p is exponentially-distributed, σ τ is Boltzmann-distributed (Eq. 3 and 4 have a σ 2 dependence instead of linear).
Subject to this modification, we observe (as expected) that both 1/α p and 1/α τ are linearly related to the internal pressure Γ.However, the proportionality constant for this equation of state is both volume-dependent (not shown) and history-dependent (see Fig. 7).Therefore, we find the stress-based ensemble to be volume-dependent and angoricity to be a variable of process.This complicates its use as a state variable, but highlights the importance of tracking both normal and tangential components of interparticle contact forces.

Network Science
The two experiments described above approach granular problems from the standpoint that the particle-scale statistics control the bulk behaviors.However, the image in Fig. 1 clearly shows that the force chains create a secondary, mesoscale structure within the bulk.It remains an open question how best to extract length and force scales in a consistent, quantitative framework.One approach to quantifying these structures is to consider a network representation in which the nodes (all particles) are connected by weighted edges with values corresponding to the measured contact forces.We construct the weighted force network W from a list of all inter-particle forces, obtained from photoelastic analysis.If particle m and n are in contact, then W mn = f mn / f , where f mn is the normal force (or, alternatively, the tangential force) between them and f is the average force in the whole network.
Such network representations -commonly used throughout the social, biological, and physical sciencesprovide a mathematical framework which allows for the comparison of features at different spatial scales.A large catalog of statistical measures have been developed by the network science community [5], allowing for the quantification of particles, force-chains, or the whole system, as desired by the particular application.Here, we focus on developing measures which correspond to the force network.
We have observed [2] that domains of force-chains, located via automated community-detection [4], are more highly associated with sound propagation than are other smaller/larger scale measures.It is important, however, to have a way to determine the optimum length scale on which to define the communities.Community-detection algorithms work by optimizing the partition of the network into sub-units (communities) based on whether a node is more-strongly connected to members of its community than to others (the null model) [4,15,16].To identify communities in a network, we use an open-source code [17] to maximize the modularity Q, defined as where node (particle) m is assigned to community g m , node (particle) n is assigned to community g n , the Kronecker delta δ(g m , g n ) = 1 if and only if g m = g m , the quantity γ is a resolution parameter (discussed below), and P mn (the null model, discussed below) is the expected weight of the edge that connects node m and node n.These methods require a well-justified choice of both the null model, and the resolution of comparison to that null model (for both the spatial and temporal scales).Force networks in granular materials are an unusual class of network in that they are spatially-embedded: a grain can only be in contact with its geographical neighbors.Therefore, we have modified existing communitydetection tools to include this constraint directly, and better-develop existing techniques so that they provide physically-meaningful results.The first modification we make is to improve the null model used in the communityoptimization process (Eq.5).Instead of using the more common Newman-Girvan null model, in which any particle could be connected to any other particle, we instead define a geographical null model [18] in which only a particle's geographical neighbors are possible edges (contacts).
Second, the choice of resolution parameter γ sets whether the optimization process probes small vs. large scales.For a single network (e.g Fig. 1), Eq. 5 provides the optimization.In order to follow the dynamics of such networks over time or strain, it is necessary to use a multilayer network which is a composite of many such networks.In this case, there are two resolution parameters: one for the spatial scale, and the other for the temporal scale.We have developed two methods to chose a spatial resolution parameter which best-captures the force chain network, illustrated in Fig. 8ab.The first method seeks a value for γ for which there is the lowest correlation between the physical (Euclidean) distance between two nodes and the 'hop' distance (measured only along network edges); this is illustrated in Fig. 8a [18].The second method seeks a value for γ which maximizes the ratio between the volume of the convex hull (see Fig. 8b) and the volume of the particles within each network [19,20].In either case, the optimum value corresponds to choosing γ ≈ 1: particles with the mean force are equally likely to be included/excluded from a network.Importantly, the optimization process does not impose a hard threshold.
Community detection techniques also work on on multilayer networks; we have adapted these methods to be used to understand granular dynamics by representing each layer as a step in time [20].This is illustrated schematically in Fig. 8c, with the interlayer coupling constant ω controlling whether two adjacent layers are considered independent (ω = 0) or coupled (ω > 0).For ω = 1, the interlayer coupling strength is matched to the mean edge weight (interparticle force).In that limit, we find that there is no flexibility in the community detection: a single set of communities is frozen throughout the dynamics.For ω = 0, the flexibility is maximized: there is no correlation in community assignments for adjacent steps.We find that choosing a value of ω ≈ 0.1 to 0.2 gives approximately half the maximum flexibility and is a reasonable choice [20].
For each community detected, we can calculate the average size, strength (mean interparticle force), and sparsity (via the complex hull).This provides a rich set of characterizations through which we can quantify how the force chain network depends on particle properties, strain history, etc.These analyses can be performed on either the normal force network or tangential force network.One such comparison is shown in Fig. 9, in which we have analyzed the same dataset shown in Fig. 5 (from [12]).
As a test-case, we examine whether community detection is able to distinguish the subsystem of low-friction particles from the bath of higher-friction particles (see §4) when tangential forces are taken into account.We observe that the increasing community strength, measured from normal forces, is not significantly different between the subsystem and bath.However, the tangential network displays a clear distinction after a few steps, suggesting that the multilayer network model and community detection are sensitive to differences in friction.Importantly, the direction of this sensitivity agrees with the expectation from local force information alone (the modular structures found in the low friction bath are characterized by lower tangential intracommunity force), but the difference is more pronounced.

Conclusions
We have used photoelastic contact force measurements to examine the network of forces within granular materials.We observe that the resulting force networks control the patterns of sound propagation, via the increase in contact area at the strongest contacts.Temperature-like variables associated with the force-moment ensemble obey some key thermodynamic-like features (e.g.obey the zeroth law) but their equations of state are invalidated through the presence of history-dependence.This suggests the presence of internal degrees of freedom -likely related to friction -which store the memory of how a granular packing was assembled.The application of network science tools within a granular context opens new horizons for investigating mesoscale phenomena within a granular material.In future work, the development of network techniques which allow for vector-weighed edges (to simultaneously account for both normal and tangential forces) will provide additional insight.Finally, while this paper focuses on measurements within 2D systems, due to the relative ease of making particle-scale force measurements, recent advances in 3D systems will soon make it possible to extend the techniques.

Figure 1 .
Figure 1.Sample image of force chains in a quasi-2D granular material composed of photoelastic disks.Dark particles carry little force, and bright particles carry more.

Figure 2 .
Figure 2. The use of multiple light fields -(a) unpolarized light and (b) polarized, monochromatic light -allows for the determination of all particle locations and all vector contact forces.Figure adapted from [3].

Figure
Figure3.Data collected at the particle-scale (small gray symbols), averaged to reveal a relationship between the pressure on a particle and the amplitude of sound propagation through that particle.The dashed line is ∝ P 0.35 , from a fit to all points.The large, dark circles are the average and standard error of this dataset, calculated for equally-populated bins.Figure adapted from[1] Figure3.Data collected at the particle-scale (small gray symbols), averaged to reveal a relationship between the pressure on a particle and the amplitude of sound propagation through that particle.The dashed line is ∝ P 0.35 , from a fit to all points.The large, dark circles are the average and standard error of this dataset, calculated for equally-populated bins.Figure adapted from[1]

Figure 4 .
Figure 4. (a) Image of force chains, subject to acoustic driving, and (b) magnitude of difference in two images, taken during the transmission of an acoustic pulse.The small ellipse marks the formation of a new contact, which repeatedly opens and closes during each wave cycle while a transiently strong force chain forms to its right.Figure reprinted from [1].
Figure 5.Equilibration of the two components of the temperature-like variable angoricity, in experiments on a lowfriction subsystem (red ) placed within a high-friction bath (black •).The low-friction particles are monitored for location via markings visible in blacklight, and are returned to the subsystem if they stray out of it; no mechanical confinement is present.Values are measured as averages over clusters of N = 8 particles, to reduce statistical fluctuations[11].Pressure Γ is measured in units of [N/m].Figure adapted from[12].

Figure 6 .
Figure 6.Left: Sample images of force chains until three loading conditions: biaxial, uniaxial, and shear.Right: Corresponding two-dimensional histograms of the volume V [AU] and pressure Γ [N/m], again measured over 8-particle clusters.Figure credit: [14] Figure 6.Left: Sample images of force chains until three loading conditions: biaxial, uniaxial, and shear.Right: Corresponding two-dimensional histograms of the volume V [AU] and pressure Γ [N/m], again measured over 8-particle clusters.Figure credit: [14]

Figure 8 .
Figure 8.(a) Schematic plot of force chain communities with different degrees of correlation between hop distance (along the edges of a network) and physical distance (as the crow flies).(b) Illustration of force chain communities with a low or high ratio between their convex hull area and the physical area of the particles.(c) Schematic image of a multilayer network comprising multiple strain steps.Adapted from [20, 21].

Figure 9 .
Figure 9. Curves showing the average community (a,d) size, (b,e) strength, and (c,e) sparsity of the subsystem (red curves) vs. bath (black curves) for both the normal (a,b,c) and tangential (d,e,f) force networks.Each strain step corresponds to a change in packing fraction of 0.009.Figure reprinted from [20].