# On Phase Transition and Self-Organized Critical State in Granular Packings

## Abstract

We model two-dimensional systems of granular aggregates confined between two planes and demonstrate that at a critical grain volume fraction an abrupt rigidity transition occurs. This transition is observed both in static and shear tests. The grain volume fraction at which the transition occurs, $`\nu _c`$, decreases with increasing friction between the grains. Densely packed grains, with a volume fraction $`\nu >\nu _c`$, display an elastic-plastic rheology. Dilute packings, with $`\nu <\nu _c`$, display gas-like characteristics. It is shown that when volume fraction is allowed to change freely (using constant normal stress boundary condition), it evolves spontaneously to $`\nu _c`$ under a wide range of boundary conditions, exhibiting ’self-organized criticality’.

Granular media is a fundamental, yet not well understood, complex system with wide ranging applications to technological and natural systems. In recent years there has been much work on granular dynamics, with emphasis on how behavior of grain aggregates may resemble solids, liquids or gases . Already Reynolds in 1885  noted that loosely packed sands deform easily as fluids while dense packings resist shear as solids. The two phases are traditionally treated separately by ’kinetic-gas’ approaches for loosely packed grains  and elasto-plastic (often using associated plasticity ) theories for dense soils. In this Letter we numerically investigate the transition between solid and gas behaviors in granular aggregates, and show that the phase boundary between gas and solid is also an attracting point, to which systems naturally evolve. We numerically model grain aggregates using a version of the popular ’discrete element method’  which treats grains as inelastic disks with rotational and translational degrees of freedom. Two grains of radius $`R_i`$ and $`R_j`$ undergo an inelastic interaction when the distance separating them $`r_{ij}`$ is less than the sum of their radii. During the interaction the $`i`$th grain experiences a contact force that has both shear and normal components:

$$F_{ij}(t)=[k_n(R_i+R_jr_{ij})\gamma m(\dot{𝐫}_{\mathrm{𝐢𝐣}}\widehat{𝐧})]\widehat{𝐧}+[min(k_s\mathrm{\Delta }s,\mu (𝐅\widehat{𝐧}))]\widehat{𝐬}$$
(1)   
where $`\widehat{𝐧}=(𝐫_{ij}\widehat{𝐱},𝐫_{ij}\widehat{𝐲})/r_{ij}`$, $`\widehat{𝐬}=(𝐫_{ij}\widehat{𝐲},𝐫_{ij}\widehat{𝐱})/r_{ij}`$, are the unit vectors in the normal and tangential directions respectively. $`k_n,k_s`$ are the normal and shear elastic constants, $`m`$ is the grain mass, $`\gamma `$ is a damping coefficient ensuring inelasticity of the interaction, $`\mu `$ is the surface friction coefficient, and $`\mathrm{\Delta }s`$ is the shear displacement since the initial contact of two particles. Force is integrated through time to calculate grain position and velocity. For collisions governed by equation (1) energy loss is governed by a normal restitution coefficient $`e_n=\mathrm{exp}(\gamma t_{col}/2)`$, where $`t_{col}=\pi (2k_n/m\gamma ^2/4)^{1/2}`$ is the collision time, and by frictional work which depends on the amount of real slip (after the elastic limit is reached) and the frictional shear force, $`\mu F_n`$. Time is measured in units of undissipated elastic wave travel time $`t_0=\sqrt{m/k_n}`$, and distance in units of average disk diameter $`x_0=2\overline{R}`$. In simulations presented here $`k_n=1,\gamma =1,k_s=0.5`$ and $`m=1`$.

In this Letter we investigate three related problems. The first problem is the general characteristics of static granular aggregates after compaction between two parallel plates. The second and third problems are the behavior of the same confined aggregates during shear using two different boundary conditions: constant distance between shearing walls, termed constant volume boundary conditions (CVBC), and constant normal stress applied to the wall, termed constant force boundary conditions (CFBC). Simulations are performed on square systems with $`n`$ disks. The top and bottom edges of the box are composed of grains glued together to form rigid rough walls of length $`l`$ (Figure 1). The box is periodic in the horizontal direction. Grain radii are randomly drawn from a Gaussian distribution that peaks at $`\overline{R}`$, with a standard deviation of $`0.5\overline{R}`$. Polydispersivity is introduced to discourage ordering effects. The system is initiated as tall loosely packed box, which is compacted vertically by normal stresses to a height $`l`$. After compaction the horizontal walls are allowed to move freely in the horizontal direction to relax forces. However, unrelaxed normal forces can be maintained on the walls in this case, since we do not allow the walls to move vertically, using CVBC. Global rearrangements during compaction and relaxation ensure (local) minimal energy configuration, as would occur during natural compaction.

After compaction and relaxation we measure properties of static configurations at different solid volume fractions, $`\nu =_{i=1}^n\pi R_i^2/l^2`$, ranging between 0.75 to 0.96. Measurements of three parameters, (a) the average number of grains touching (i.e. exerting a force on) a grain in the interior of the box, termed the coordination number $`Z`$, (b) The normal stress (normal force per unit length) $`N`$ operating on the upper and lower confining walls, and (c) the systems shear modulus, G, are presented in Figure 2. All measures show an abrupt change in behavior at a critical volume fraction $`\nu _c`$, which depends on the value of friction prescribed between the grains, but not on the system size. The coordination number is approximately zero for $`\nu <\nu _c`$. At $`\nu _c`$ the coordination number abruptly jumps and then increases as an empirical power-law, $`Z(\nu \nu _{c}^{}{}_{\mu }{}^{})^{\alpha _\mu }`$, where the subscript $`\mu `$ denotes different friction values used in simulations. Power-law fits (in solid curves) yields $`\nu _{c}^{}{}_{0}{}^{}=0.83\pm 0.01,\alpha _0=0.5\pm 0.1`$ and $`\nu _{c}^{}{}_{0.5}{}^{}=0.805\pm 0.01,\alpha _{0.5}=0.3\pm 0.05`$. Critical behavior with values of $`\nu _{c}^{}{}_{0}{}^{}=0.82\pm 0.02`$ were obtained for hard smooth disks , and visco-elastic 2D bubbles  for mono-dispersed and polydispersed systems, demonstrating that $`\nu _c`$ is fairly independent of the disk size distribution, and the interaction law in absence of friction. The difference in $`\nu _c`$ and $`\alpha `$ between frictional and smooth grains occurs because frictional grains tend to ’stick’, and thus cannot achieve the lower energy configuration of smooth disks.

The fact that $`N`$ and $`Z`$ (Figure 2a,b) are zero for $`\nu <\nu _c`$, indicates that below $`\nu _c`$ grains do not touch. At $`\nu _c`$ grains first touch and elastically repel each other, exerting normal stresses on the walls. For $`\nu >\nu _c`$, normal stress follows $`NZ(\nu \nu _c)`$ (shown in solid lines), a slightly modified form (accounting for the phase transition) of predictions from standard models of densely packed elastic disks . The observed transition is identified as a macroscopic rigidity transition in Figure 2c. The system’s elastic shear modulus $`G`$ is obtained by imposing a small homogeneous shear step strain of a magnitude $`ϵ`$ and measuring the resulting shear stress on the walls $`\sigma `$ ($`G=\sigma /ϵ`$ is independent of $`ϵ`$ for $`ϵ<10^4`$ of the grain radius, here we use $`ϵ=10^5`$). The procedure follows that outlined in . Figure 2c shows $`G`$ increasing from zero to a finite value, as the system passes through the phase transition at $`\nu =\nu _c`$. Physically, the rigidity should depend on average number of springs per disk, yet not on their compression, i.e. $`GZ`$ (solid curves). The rigidity-transition thus identified is a result of global geometrical constraints of grain packings, with the coordination number identified as the order parameter for the transition. It was previously shown that $`\nu =0.83\pm 0.02`$ marks the upper limit of compacity of disordered packings of smooth hard mono-sized and poly-sized disks , where long-range order in disk positions appearing for $`\nu >\nu _c`$.

To investigate the transition in dynamic behavior of grains we conducted a set of simulation shearing the static 24x24 configurations, in couette flow. Here we show results of simulations with $`\mu =0.5`$. (Non-frictional grains and other values of friction show similar behavior with the transition occuring at $`\nu _{c}^{}{}_{\mu }{}^{}`$, and are thus not presented). The upper wall was moved at a constant velocity ($`v=10^3x_0/t_0`$) while $`l`$ was kept constant, maintaining CVBC. We observe different behavior as a function of the solid fraction: In configurations with $`\nu <\nu _c`$, momentum is transferred from the wall to interior grains mostly via short-lived collisions, observed in spiky fluctuations in stress measured on the wall and in coordination number, Fig. 3. Stresses are transmitted only within local clusters (as in Figure 1(t2)). The power spectra of the stress fluctuations time series approaches white noise demonstrating the uncorrelated nature of stress-transfer. For dense packings, with $`\nu >\nu _c`$, grains interact via long-lasting contacts. Global motion is characterized by elastic-plastic cycles: clusters of grains in contact accumulate recoverable elastic strain, but when stresses become too great, grains suddenly rearrange to relieve the stress. Continued shearing then begins accumulation of elastic strain on the new particle arrangement. This behavior leads to a stick-slip stress time series (Figure 3). Long clusters of these grains in contact form system-spanning ’stress chains’ that transmit forces from the boundaries into the interior (as in Figure 1(t1)). The stress fluctuation power spectra for $`\nu >\nu _c`$ follows $`f^\eta `$, where $`\eta 2`$, indicating long-time correlations and in agreement with experimental results  conducted in densely packed systems. Most interesting is the behavior at the transition point: When $`\nu =\nu _c`$ the system oscillates between a solid-like ’jammed’ state (Figure 1(t1)) and a gas-like behavior (Figure 1(t2)). At $`\nu =\nu _c`$, values of $`Z,\sigma ,N`$, fluctuate between values characteristic of the ’solid-phase’ and those characteristic of ’gas-phase’ (Figure 3). The transition density thus bridges the gap between liquid and solid by coexistence of the two phases. We note the relation to Kirkwood-Alder transition (KAT), a disorder-order transition of repulsive hard-spheres studied extensively in the context of collidol suspensions . The relation to KAT, as well as the transition characteristics, suggest that the granular rigidity transition is a first-order transition.

The relation between stress and strain rate (to be published elsewhere) is another property that changes across the phase boundary. When $`\nu <\nu _c`$ we measure $`(\sigma ,N)v^2`$, as expected from theory  and experiments . When $`\nu >\nu _c`$ measured average stresses are nearly independent of strain rate, as expected for elastic-plastic materials and as seen in experiments in densely packed granular systems . For $`\nu =\nu _c`$ stress-strain-rate curves resemble those of plastic materials, since high stresses in jammed states dominate time-averaged behavior.

We finally demonstrate the role of the critical solid fraction when the system is allowed to evolve to it’s own prefered solid fraction. Soil mechanists have long known of a ’critical density’. If a granular aggregate is over-compacted and sheared under CFBC, it will deform while shearing and expand to this ’critical density’. If it is initially under-consolidated it will compact while shearing until it reaches the same ’critical density’. We perform simulations, similar to these experimental conditions, by applying a constant normal stresses to the rigid boundaries (CFBC), while shearing the upper wall at velocity $`v`$. The systems expand or contract, depending on the initial porosity, and finally reach a steady state where variables ($`\nu ,\sigma ,Z`$) fluctuate around a constant value. Steady state solid fractions, $`<\nu >`$, are shown in Figure 4(a) as function of $`N`$ and $`v`$, for simulations with $`\mu =0.5`$. For a wide range of applied normal stresses (range widens with decreasing velocity) systems attain the critical volume fraction for frictional grains, i.e. $`<\nu >\nu _{c}^{}{}_{0.5}{}^{}=0.805`$. For smooth disks, simulations converge on $`<\nu >\nu _{c}^{}{}_{0}{}^{}=0.835`$, in a similar manner. Though $`<\nu >`$ is fairly independent of $`N`$ in the ’critical regime’, the average coordination number $`<Z>`$ increases (Figure 4(b1)(b2)). In those systems having $`<\nu >\nu _c`$ system-spanning stress-chains coexist with unstressed ’gas islands’. With increasing $`N`$, the number and size of gas islands decreases, and the connectivity of chains increases, and thus $`<Z>`$ increases. Deviations from the critical state occur in two cases: 1) At high enough normal stresses, ’solidification’ occurs: $`<\nu >`$ increases above $`\nu _c`$, islands virtually vanish and chains become fully connected ($`N=10^2`$ in Figure 4a). 2) When the ratio of inertial to normal forces is high inertial forces cause decompaction, $`<\nu >`$ becomes smaller than $`\nu _c`$ (simulations performed at higher velocities and lower normal stresses: Figure 4a, $`v=10^3`$ and $`N10^5`$), with islands growing to divide stress chains, resulting in ’liquification’. We also observed the gas to solid transition in the power spectra of stress fluctuations, where power follows $`f^\eta `$ with $`\eta `$ increasing continuously from $`\eta 0`$ for $`N=10^7`$ to $`\eta 2`$ for $`N10^2`$. In these CFBC simulations the system, although at the critical state, does not flip between gas and solid phases as in CVBC, since ’jamming’ episodes may be avoided by slight dilation (producing fluctuations around mean porosity and height). Instead of temporal coexistence seen in CVBC the critical state in CFBC is marked by spatial coexistence of two phases: gas-islands and stress-bearing chains.

Why is the transition between gas and solid also the ’critical-density’ to which the system is attracted? any finite normal stress acting on grains would cause grain compaction and contact, and thus loss of ’gas-like’ properties, approaching $`\nu _c`$ from below. Solid fractions of $`\nu >\nu _c`$ are characterized by finite elastic deformation ($`>1\%`$) of grains, which require extremely high normal stresses for stiff natural granular material as rocks. Thus $`\nu _c`$ is the rigid limit, where stresses are accommodated by efficient load bearing structures. Experimentally, a coexistence of liquid and solid regimes has been observed  in shearing granular systems, and oscillations between ’jamming’ and ’flowing’ states occur spontaneously in a variety of systems, from hoppers to natural and experimental land-slides , suggesting proximity of these systems to the phase boundary. An attracting phase boundary may also explain ’fragile materials’, a recent term used to describe ’jammed’ states which may be unjammed by small fluctuations .

To summarize, we identified the rigidity transition for granular media in $`2D`$ at $`\nu _c=0.800.84`$. Two-dimensional results can be mapped to corresponding 3D volume fractions using a relation for inter-particle voids, $`\nu _{3D}=4\nu _{2D}^{3/2}/3\pi ^{1/2}`$ , predicting that in 3D $`\nu _c=0.540.58`$ with the lower number representing frictional grains. Experiments  confirm that immense stiffening of rapidly shearing grain aggregates occurs at $`\nu _c0.54`$. The critical solid-fraction is a phase-boundary between gas and solid regimes of behavior, where the two phases were observed to coexist in space and/or time. The critical state is self-organizing, and corresponds to the ’critical density’ known in soil mechanics: Sheared aggregates tend to this critical volume fraction under a wide range of normal loads. We used normal stresses ranging between $`10^210^7`$, which using characteristic young modulus for Earth materials (e.g. rock), translates to $`N=110^5`$ KPa, (corresponding for soils to burial depth of $`110^5m`$). Based on this we suggest that many natural granular deformation processes will occur at a solid fraction which constitutes a phase boundary, and neither gas nor elasto-plastic descriptions will fully capture these systems behavior. However, such criticality in natural systems is easily missed, since spatial and temporal averages of stresses tend to be dominated by stress chains, and thus have the mark of solid deformation. It is clear that much more theoretical and experimental work needs to be done, including investigating much larger systems, and examining the critical state more closely.