# Shock Waves in Weakly Compressed Granular Media

###### Abstract

We experimentally probe nonlinear wave propagation in weakly compressed granular media, and observe a crossover from quasi-linear sound waves at low impact, to shock waves at high impact. We show that this crossover grows with the confining pressure , whereas the shock wave speed is independent of — two hallmarks of granular shocks predicted recently gomez . The shocks exhibit powerlaw attenuation, which we model with a logarithmic law implying that local dissipation is weak. We show that elastic and potential energy balance in the leading part of the shocks.

###### pacs:

Many disordered materials, including granular media gomez ; bobPRL ; Cheng ; bobfragile , foams gijsjam and emulsions emulsion ; brujic , lose their rigidity when their confining pressure is lowered. In almost all cases, the resulting unjamming transition goes hand in hand with the vanishing of one or both elastic moduli Bolton ; Durian ; makse ; OHernetal ; ellenbroek ; SS ; finsize ; review1 ; review2 ; footnote1 , and consequently, nonlinearities must dominate when such marginal systems are subjected to finite stresses gomez ; N&V ; ohern . For example, soft particles exhibit nonlinear rheology near jamming emulsion ; gijsflow , even when their local elastic and viscous interactions are linear olson ; brian , and marginally connected spring networks exhibit nonlinear elasticity near their critical points MW ; Fred . In these two cases, the vanishing of the elastic moduli is a collective phenomenon, closely connected to the isostatic character of the marginal state Bolton ; Durian ; makse ; OHernetal ; ellenbroek ; SS ; finsize ; review1 ; review2 ; footnote1 .

Here we experimentally probe a different scenario where local nonlinearities near unjamming lead to vanishing elastic moduli and nonlinear excitations: shock waves in granular media shocknote . Granular media have frictional interactions, and these prevent granular media to reach the isostatic limit. Consequently, there are no collective mechanisms leading to vanishing elastic moduli or nonlinearities review1 ; ellakcritical ; Kostya ; silke . Nevertheless, when granular media unjam as the pressure is lowered, the individual contacts weaken due to the nonlinear local Hertz contact law, which states that for elastic spheres, the contact forces scale with deformations as Johnson ; footnote2 . As a result, frictional granular media have a vanishing linear response at their unjamming point, and their elastic moduli and sound speed vanish as and , respectively makse ; ellaksoundwaves ; nesterenko ; weakening ; jia .

Recent simulations on frictionless Hertzian media by Gomez et al. suggest that sound waves give way to strongly nonlinear shock waves near the unjamming () point gomez . Three crucial questions remain open, as the numerical model of Gomez has no static friction, no dissipation and is in 2D. First, realistic granular media are frictional: do shock waves also arise for non-isostatic, frictional systems? Second, friction also leads to dissipation — do shock waves survive realistic levels of dissipation? Third, can such shock waves be excited in 3D experiments?

To answer these questions, we experimentally probe sound and shock waves by impacting a weakly compressed granular medium with a heavy mass, while measuring the propagation speed and front shape for a wide range of impact magnitudes (Fig. 1). We find a novel power law attenuation of the shock waves, which we can model with a simple model where local dissipation depends logarithmically on grain forces. Notwithstanding this weak dissipation, our shock waves exhibit the three main hallmarks of the numerically observed conservative shocks: (i) a crossover from linear waves to shock waves when the impact pressure exceeds the confining pressure ; (ii) independence of the shock speed on but powerlaw scaling with impact strength; (iii) a balance of kinetic and potential energies in the leading edge of the shock waves. The ease with which we can excite such granular shock waves suggests that they play an important role whenever loose granular media are strongly excited Thoroddsen ; detlef ; johnjet ; ChengSheet ; plow .

Setup — Our setup consist of a large metal container ( cm) filled with glass beads (diameter mm) and covered by an 2 kg aluminum top plate. To excite shock waves, we impact a freely sliding cylindrical piston (A), (diameter 10 cm, length 18 cm, mass 3.7 kg), which makes contact with the grains through a circular hole in the side of the container (B), with a mass (C, 1.2 kg) suspended from a pendulum (mass 2 kg, length 1.3 m). We detect wave propagation throughout the material via pressure sensors and accelerometers burried in the granular medium; the sensors are attached to steel wires to ensure their correct positioning and orientation. We thus probe the local excess pressure (our pressure sensors have no DC response) and acceleration at distances and from the impact zone (Fig. 1). See Supplemental Material at [URL will be inserted by publisher] for experimental details.

Phenomenology — In Fig. 2 we show typical time traces of the local pressures and velocities detected at locations cm and cm, for a strong impact and low confining pressure. The waves take the form of fronts with a clearly identifiable leading edge where the pressures and particle velocities peak, followed by a long tail with a complex structure. Here we will focus on the nonlinear regime, where the propagation is set by the amplitude. We characterize our waves by the peak pressures, and , and peak particle velocities, and , at and . We determine the front speed from the time of travel , where is the interval between reaching its 50% value at and .

Before embarking on a systematic study of the propagation of these fronts, we briefly discuss the typical scales shown in Fig. 2. First, the typical grain velocities are of the order of cm/s, the impact speeds of the plunger are of the order of m/s whereas is of the order of hundreds of m/s. Second, the grain displacements at at the time of the pressure peak are of the order of 10 m, consistent with purely elastic deformations (See Supplemental Material at [URL will be inserted by publisher] for details), whereas the total motion of the plunger into the sand for such impacts is of the order of a mm.

The picture that emerges is that upon impact, a rapid front is formed, which, for the cases when is larger than the sound speed, we will call a shock wave. We note that we do not observe the breaking up of our fronts in several distinct solitons — in the simulations of gomez , the stability of the shocks was attributed to the disorder of the granular packings. The shock being nonlinear, its speed is set by its amplitude, and the leading edge of the shock remains leading, as it outruns whatever happens in the tail. The interactions in the leading edge are dominated by elastic, Hertzian interactions. However, long after the shock wave has outrun our system, the plunger is still slowly penetrating the sand bed, leading to a very long tail, where the vast majority of rearrangements and dissipative events take place.

Propagation speed — To probe the nature of the waves excited in this system, we varied the pressure at the depth of the sensors from 0.9 to 6 kPa and have determined the propagation speed for a wide range of impact strengths. To compare our data to the theoretical predictions of gomez , we plot our data in a so-called Hugoniot plot. Due to attenuation (see Fig. 2), and as a measure of the impact strength we use their geometric mean, .

Fig. 3 shows as a function of and , and that our data is well fit by the form , which is very close to the more involved analytical formula for shock waves presented in gomez .

The main three features of our data are (i) For strong impacts, becomes independent of , and scales consistently as — the dependence on is a hallmark of nonlinear waves, and the exponent is the one predicted for Hertzian shock waves. (ii) For weak impacts, becomes independent of the impact strength — the hallmark of linear waves — but increases with pressure due to the nonlinear local interaction law. (iii) The crossover from the linear to nonlinear regime is expected to arise when , because when , the Hertzian interaction can be linearized, whereas for , linearization fails and shocks arise gomez . Indeed we find that grows with .

All these features are in good qualitative and quantitative agreement with the earlier numerical findings of Gomez et. al: the waves we excite for large impacts are indeed shock waves.

Attenuation — The peak pressure diminishes whilst the shock propagates through the material. In Fig. 4 we compare the peak pressures and at cm and cm for the same experiments as shown in Hugoniot plot. Strikingly, the attenuation varies significantly with impact strength: for the weakest impacts, , while for larger impacts the relation between and is consistent with power law scaling:

(1) | |||||

(2) |

where , , the characteristic pressure Pa and .

We will now determine a simple local model that captures this remarkable powerlaw attenuation. We ignore disorder and imagine the signal to propagate over a linear chain of beads, where each bead attenuates the signal by a factor , where the local attenuation depends on the local pressure . After some algebra, we find that we can capture the power law relation between and when has the following logarithmic form (See Fig. 5a):

(3) | |||||

(4) |

where is a material constant — See Supplemental Material at [URL will be inserted by publisher] for details.

A surprising consequence of this model is that it predicts that the exponent should depend on the distance between and . To test this, we have performed several sets of experiments where we vary this distance. Fig. 5b shows that the log slope relating and indeed increases for larger propagation distance — these trends are captured by our model without additional fit parameters.

Local Elastic Motion and Energy Balance — As Fig. 5a illustrates, even for the strongest impacts we can access, the attenuation per particle is not more than 10% — it is likely smaller, as our model underestimates the number of contacts per unit length in a disordered media. Hence, we believe that the interactions in the leading edge of the shock wave are predominantly elastic. Additional evidence for this comes from a direct comparison between the contact forces and particle displacements; these are consistent with purely elastic (Hertzian) deformations (See Supplemental Material at [URL will be inserted by publisher] for details), showing that sliding and rearrangements do not play a dominating role in the leading edge of the shock.

This motivates us to probe the balance between elastic and potential energies in our shock waves. For the conservative shocks studied in the numerical simulations gomez , the kinetic and potential energies balance in the shock regime, whereas away from the shock regime, the kinetic energy tends to zero, but the potential energy saturates at its lower bound .

To probe this balance in our experimental data, we present scatter plots of the peak velocities and vs the peak pressures and in Fig. 6. The kinetic energy simply scales where is the maximum local velocity, whereas the potential energy for Hertzian contacts scales , where and are the local maximum deformations and forces — note that the contact force should contain both the DC component and AC component . Our data shows a strikingly good qualitative agreement with the numerical data by Gomez et al. (Fig. 3d of gomez ): it is fully consistent with a power law in the strongly nonlinear regime notesaturate , and a pressure dependent deviation from this law in the weakly nonlinear regime. Moreover, an estimate of the kinetic and potential energies per particle based on the data shown in Fig. 6, shows that they are of similar order, just as in the simulations gomez . We finally note that this balance is equally good at and , even though the energies in are smaller than in due to attenuation. All this suggest that the attenuation does not significantly upset the balance between kinetic and potential energies.

Discussion— By varying the impact strength, we can tune our waves from pressure dependent sound waves at low impact, to pressure independent nonlinear shock waves at higher impact, similar to what was predicted for dissipationless Hertzian particles gomez . We find that propagating from one grain to the next, a small amount of energy is dissipated, which leads to novel powerlaw scaling of the attenuation of the shocks amplitude over several decades. Is this attenuation caused by scattering or dissipation? While we cannot rule out scattering, we note that the overall decrease of the pressure and velocity profiles shown in Fig. 2 favor a dissipative picture; in addition, the width of the leading edge does not change very much under propagation. Notwithstanding this weak dissipation, the magnitude of displacements in the shock, and the balance of kinetic and potential energy, show that the physics in the leading edge of the shock is dominated by elastic interactions.

The shock waves we observe here are therefore qualitatively different from the slow granular densification waves known as “plowing” plow ; cornstarch . Plowing is associated with densification through highly dissipative rearrangements, and such densification fronts propagate with velocities far below the sound speed plow . Similarly slow and dissipative events may occur in the tail of our waves, but the point is that the leading edge of our shock waves propagate faster than the sound speed. They are also qualitatively different from the weakly nonlinear waves observed under continuous driving weakening ; jia .

Finally we note that there are two pressure scales that set the physics of granular shock waves in experiments— the external pressure that sets the crossover from linear to shock waves, and the characteristic pressure above which attenuation sets in. In our experiment, , but it is conceivable that for more elastic particles, or in microgravity, one can reach , in which case, virtually dissipation free granular shock waves could be produced.

###### Acknowledgements.

We acknowledge discussions with J. Burton, L. Gomez, H. Jaeger, X. Jia, and V. Vitelli, and funding from FOM, SHELL and NWO.## References

- (1) L.R. Gomez, A.M. Turner, M. van Hecke, and V. Vitelli, Phys. Rev. Lett. 108, 058001 (2012)
- (2) T. S. Majmudar, M. Sperl, S. Luding and R. P. Behringer, Phys. Rev. Lett. 98, 058001 (2007).
- (3) X. Cheng, Phys. Rev. E 81, 031301 (2010).
- (4) D. Bi, J. Zhang, B. Chakraborty and R. P. Behringer, Nature 480 355 (2011).
- (5) G. Katgert and M. van Hecke, EPL 92 34002 (2010).
- (6) K.N. Nordstrom, E. Verneuil, P.E. Arratia, A. Basu, Z. Zhang, A.G. Yodh, J.P. Gollub and D.J. Durian, Phys. Rev. Lett. 105, 175701 (2010).
- (7) I. Jorjadze, L. L. Pontani, and J. Brujic, Phys. Rev. Lett. 110, 048302 (2013).
- (8) F. Bolton and D. Weaire, Phys. Rev. Lett. 65, 3449 (1990).
- (9) D. J. Durian, Phys. Rev. Lett. 75, 4780 (1995).
- (10) H. A. Makse, N. Gland, D. L. Johnson and L. M. Schwartz Phys. Rev. Lett. 83 5070 (1999); D. L. Johnson, H. A Makse, N. Gland and L. Schwartz, Phys. B 279, 134 (2000); H. A. Makse, N. Gland, D. L. Johnson and L. Schwartz, Phys. Rev. E 70 061302 (2004).
- (11) C. S. O’Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
- (12) W. G. Ellenbroek, E. Somfai, M. van Hecke and W. van Saarloos, Phys. Rev. Lett. 97 258001 (2006).
- (13) S. Dagois-Bohy, B.P. Tighe, J. Simon, S. Henkes, and M. van Hecke Phys. Rev. Lett. 109, 095703 (2012).
- (14) C. P. Goodrich, A. J. Liu and S. R. Nagel, Phys. Rev. Lett. 109, 095704 (2012).
- (15) M. van Hecke, J. Phys. Cond. Matt. 22 033101 (2010)
- (16) A. J. Liu and S. R. Nagel, Ann. Rev. Cond. Matt. Phys. 1, 347 (2010).
- (17) The only counter example of a system retaining finite elastic moduli near jamming we know of are frictional harmonic particles, which are not isostatic review1 .
- (18) V. Vitelli and M. van Hecke, Nature 480, 325 (2011)
- (19) C. F. Schreck, T. Bertrand, C. S. O’Hern, and M. D. Shattuck, Phys. Rev. Lett. 107, 078301 (2011).
- (20) G. Katgert, M. E. Möbius and M. van Hecke, Phys. Rev. Lett. 101, 058301 (2008); G. Katgert, A. Latka, M. E. Möbius and M. van Hecke, Phys. Rev. E 79 066318 (2009).
- (21) P. Olsson and S. Teitel, Phys. Rev. Lett. 99, 178001 (2007).
- (22) B.P. Tighe, E. Woldhuis, J.J.C. Remmers, W. van Saarloos, and M. van Hecke, Phys. Rev. Lett. 105, 088303 (2010).
- (23) M. Wyart, H. Liang, A. Kabla and L. Mahadevan, Phys. Rev. Lett. 101, 215501 (2008).
- (24) C. P. Broedersz, X. Mao, T. C. Lubensky and F. C. MacKintosh, Nat. Phys. 7, 983 (2011).
- (25) We use the term shock waves here in a loose sense, meaning strongly nonlinear wave.
- (26) E. Somfai, M. van Hecke, W. G. Ellenbroek, K. Shundyak and W. van Saarloos, Phys. Rev. E. (R) 75, 020301 (2007)
- (27) K. Shundyak, M. van Hecke and W. van Saarloos Phys. Rev. E. (R) 75, 010301 (2007).
- (28) S. Henkes, M. van Hecke and W. van Saarloos, EPL 90, 14003 (2010).
- (29) K. L. Johnson, Contact Mechanics (Cambridge University Press, 1987).
- (30) For purely 2D Hertzian contacts, with , although 2D discs in experiments have somewhat larger exponents bobPRL .
- (31) E. Somfai, J.-N. Roux, J. H. Snoeijer, M. van Hecke, and W. van Saarloos, Phys. Rev. E 72, 021301 (2005).
- (32) V. F. Nesterenko, J. Appl. Mech. Tech. Phys. 24, 733 (1984); ibid., Dynamics of Heterogeneous Materials (Springer-Verlag, New York, 2001).
- (33) X. Jia, Th. Brunet and J. Laurent, Phys. Rev. E 84 020301(R) (2011)
- (34) S. van den Wildenberg, M. van Hecke and X. Jia, EPL 101 14004 (2013).
- (35) S. T. Thoroddsen and A. Q. Shen, Phys. Fluids 13, 4 (2001).
- (36) G. Caballero, R. Bergmann, D. van der Meer, A. Prosperetti and D. Lohse, Phys. Rev. Lett. 99 018001 (2007).
- (37) J. R. Royer, E. I. Corwin, A. Flior, M.-L. Cordero, M. L. Rivers, P. J. Eng and H. M. Jaeger, Nat. Phys. 1 164 (2005).
- (38) X. Cheng, G. Varas, D. Citron, H. M. Jaeger and S. R. Nagel, Phys. Rev. Lett. 99, 188001 (2007).
- (39) J. R. Royer, B. Conyers, E. I. Corwin, P. J. Eng, and H. M. Jaeger, EPL 93, 28008 (2011).
- (40) For the highest impacts, our accelerometers saturate, making an accurate estimate of impossible — such points have been omitted in Fig. 4, which explains the somewhat smaller range of in comparison to Fig. 3.
- (41) S. R. Waitukaitis and H. M. Jaeger, Nature 487, 205 (2012).

Supplementary Material

Experimental Details — To detect the propagation of waves throughout our granular material, we employ several Brüel & Kjaer sensors. First, at each location , we use pairs of accelerometers (Deltatron 4508 B002, 1 V/g and 4508 B001, 10 mV/g, resonance at 26 kHz, dimensions 1x1x1.6 cm, weight 4.8 g). Together these cover the wide range of accelerations that we encounter (from less than 0.1 m/s to larger than 7000 m/s), although they saturate for some of the very strong impacts, which we conclude to lead to accelerations exceeding 700 g! The accelerometers allow us, by integration, to measure the local particles velocities and displacements.

We also use force sensors (Type 8230, 112.5 mv/N at and 82301-001, 22.05mV/N at , resonance at 75 kHz, cylindrical shape 19 mm long, circular contact area 1.33 cm, weight 30 g). These do not saturate for our experiments and we focus on the signals from these sensors in most of our analysis. The sensors are much stiffer (2kN/m) than the typical stiffnesses of the Hertzian contacts we encounter, and are immersed in the grains, thus acting as pressure sensors. They have no DC response, so that they only are sensitive to the deviation from the confining pressure. To convert force to pressure, we divide by their circular contact area.

In our experiments, we vary the pressure at the depth of the sensors from 0.9 to 6 kPa by varying the load on the topplate. For each , we perform several series of experiments, each consisting of hundreds of impact experiments. We vary the impact strength (characterized by the peak pressure) over four decades: lightly tapping the resting pendulum generates the weakest impacts we can detect, while swinging the pendulum with full force from its maximum height generates the strongest impacts — only at the lowest confining pressures and at the strongest impacts, the piston penetrated the packing significantly, limiting the number of experiments in this regime. Before each series of experiments the container was freshly filled and then tapped to stabilize the packing, while after each series the container was emptied and the correct placing of the detectors was verified. No systematic dependence on filling procedure or ”age” of the sample was detected. For each impact, the signals from all six sensors were recorded at 250 kHz, and data sets ranging from 300 point before to 2000 points after impact were stored. To accurately determine the timing where the pressures reach their 50% peak magnitude, we perform linear fits to the 20% - 80% leading slope of the force signals.

Wave profiles — Our sensors allow us to extract the time evolution of the local grain pressure, acceleration, velocity and displacement. Fig. 7 shows representative examples of these profiles, both for a weak and a strong impact, qualitatively illustrating important aspects of the phenomenology of the waves we observe.

These signals allow us to extract typical scales characterizing the wave profiles. The first striking observation is that even for strong impacts, the displacements are rather small, and consistent with Hertzian deformations. Assuming that approximately 10 particles make contact with the sensor, we find that for the string shock wave shown in Fig. 7e-h, the peak force at is of order 0.6 N (and is reached at time ms). Using Hertz law, which states that the contact force varies with particle deformation as Johnson , and taking Gpa (typical value for glass), we estimate that m. The location is 12 particle diameters 24 contact zones away, so the cumulative motion at , assuming a fairly flat pressure profile, is predicted to be of order m, which is close to the measured displacement at ( m, see Fig. 7h). Such estimates are also reasonably good for other nonlinear waves, which suggests that in the leading edge of the wave, the deformations are predominantly elastic, and grain sliding and rearrangements have not arisen yet (these do arise in the long tail of the wave weakening ; jia ). The elastic nature of the contacts also explains why the pressure and velocity signals look qualitatively similar — local forces are proportional to the pressure, while local deformations of the grains are given by the gradient of , which, assuming that the pulse moves with constant speed , is proportional to the local grain velocity .

Model for Attenuation — In our simple model for attenuation, we ignore disorder and thus imagine the signal to propagate over a linear chain of beads, where each bead attenuates the signal by a factor — see Fig. 8a. We take as our length scale (), define , and in Eq. (5-9) below drop the tildes for simplicity, so that

(5) |

and after taking the continuum limit we obtain:

(6) |

To infer from the observed relation between and (Eq. 1-2) is far from trivial, as follows from integrating Eq. (6) over a finite distance .

To determine we note that solutions of Eq. (6), , can all be derived from a single mastercurve : , where the shift is adjusted to satisfy the initial condition of Eq. (6) — see Fig. 8a. Hence, we obtain the following relation between , and :

(7) |

Focussing on the power law relation between and for , we need to determine that satisfies: , which suggest that should be a double log function — so as to turn the exponent into the additive term . After some algebra, we obtain that is a double exponential function:

(8) | |||||

(9) |

where we stress that asymptotes to , indicating the absence of dissipation when (Fig. 8b).

This solution satisfies Eq. 6 when has a logarithmic form, which in dimensional quantities is given by Eqs. (3) and (4).

Potential and Elastic Energies — The peak kinetic energy per particle measured at equals the , where for average grains ( mm), kg. The peak potential energy per contact equals , where is the maximum deformation and the contact force, which we assume to be Hertzian: , where . Hence, . Assuming that 10 particles are in contact with the force sensor of area m, and taking into account that the total pressure is the sum of and , we can estimate the individual contact forces as .

Now assuming that the peak potential energy of contacts equals the peak kinetic energy, we arrive at

(10) | |||

(11) |

The straight lines in Fig. 6 correspond to , which is a reasonable number for frictional particles in 3D, which, in the limit of hard frictional spheres have .