# The influence of packing structure and interparticle forces on ultrasound transmission in granular media

See allHide authors and affiliations

Edited by David A. Weitz, Harvard University, Cambridge, MA, and approved June 1, 2020 (received for review March 7, 2020)

## Significance

Structure-property relations of granular materials are governed by the arrangement of particles and the chains of forces between them. These relations enable design of wave damping materials and nondestructive testing technologies. Wave transmission in granular materials has been extensively studied and demonstrates rich features: power-law velocity scaling, dispersion, and attenuation. However, the precise roles of particle arrangements and force chains on these features remain topics of continued research interest. Here, we employ X-ray measurements and analyses to show that velocity scaling and dispersion arise from both particle arrangements and force chains, while attenuation arises mainly from particle arrangements.

## Abstract

Ultrasound propagation through externally stressed, disordered granular materials was experimentally and numerically investigated. Experiments employed piezoelectric transducers to excite and detect longitudinal ultrasound waves of various frequencies traveling through randomly packed sapphire spheres subjected to uniaxial compression. The experiments featured in situ X-ray tomography and diffraction measurements of contact fabric, particle kinematics, average per-particle stress tensors, and interparticle forces. The experimentally measured packing configuration and inferred interparticle forces at different sample stresses were used to construct spring networks characterized by Hessian and damping matrices. The ultrasound responses of these network were simulated to investigate the origins of wave velocity, acoustic paths, dispersion, and attenuation. Results revealed that both packing structure and interparticle force heterogeneity played an important role in controlling wave velocity and dispersion, while packing structure alone quantitatively explained most of the observed wave attenuation. This research provides insight into time- and frequency-domain features of wave propagation in randomly packed granular materials, shedding light on the fundamental mechanisms controlling wave velocities, dispersion, and attenuation in such systems.

Stress wave propagation through granular media is important for detecting the magnitude of seismic events, locating oil and gas reservoirs, designing acoustic insulation and absorption materials, and designing materials for wave redirection (1⇓⇓⇓–5). Granular materials form disordered packing structures and feature heterogeneous interparticle forces controlled by nonlinear contact laws (6⇓⇓⇓–10). Ultrasound waves necessarily propagate through these interparticle contacts, making the prediction of wave behavior a highly nonlinear problem (11). Extensive research has been conducted using mathematical, numerical, and experimental tools to elucidate ultrasound wave behavior in granular materials in both the time and frequency domains. Some phenomena that have previously been studied include stress-dependent velocity scaling (8, 12), frequency filtering (8, 13, 14), wave dispersion (15⇓⇓–18), band gaps (8, 13, 19), rotational waves (17, 20⇓–22), wave focusing (23), and wave scattering (20, 21, 24, 25).

In the time domain, the sound velocity of waves traveling through granular media has been found to follow a power-law relationship with respect to the applied pressure with an exponent varying from 1/4 to 1/6 for one-dimensional (1D) chains (8, 13, 26) and two-dimensional (2D) (27) and three-dimensional (3D) ordered (15, 16) sphere packing. Here, the exponent value of 1/6 corresponds to theoretical predictions of Hertzian contact and has been confirmed in solitary wave propagation experiments on 1D chains of spherical particles (8). The pressure-induced transition of power-law exponent ranging from 1/4 to 1/6 has also been observed not only for randomly packed spheres (12, 28, 29) but also, for natural geomaterials, including dry and water-saturated soils and sands, under both preconsolidation and overconsolidation (28, 30⇓⇓–33). This transition has been qualitatively attributed to the deformation of interparticle asperities (28), force chain arrangements (27, 28, 33), stress heterogeneity (32, 34), packing reorganization (12, 16, 29), nonspherical contact geometries (27, 32), wave modes (such as compressional and shear wave), and wave amplitudes (12, 32, 33, 35). However, a study quantitatively linking packing structure disorder and contact force heterogeneity to wave velocity and velocity scaling in 3D disordered granular materials has not yet been performed.

In the frequency domain, granular materials have been found to feature dispersion: waves propagate at different velocities depending on their wavelength (or frequency) and the material’s particle size (24, 25). The acoustic signal measured through a disordered granular medium can be decomposed into coherent and incoherent parts. The coherent portion of the signal is the leading portion that represents the average of many measurements and provides a link to elastic properties and coordination number (16, 25, 36). The incoherent portion of the signal is the trailing portion that diminishes in magnitude when responses are averaged over an ensemble of the measured signals through the same packing configuration (16, 37). Theoretical calculations based on granular crystals that accounted for force balance and included Cosserat effects (17, 38) were able to provide reasonable predictions of wave dispersion for coherent waves in ordered granular materials. The differences between predictions and measurements were ascribed to nonlinearity of contacts and structural disorder. Unlike the propagation of plane waves in 1D chains and 2D and 3D granular crystals, the transmitted ultrasound signals for 3D disorder granular systems arise from a superposition of the signals traveling through different paths in this heterogeneous medium (23, 27, 39). Even small amounts of structural disorder can dramatically influence wave propagation, causing departures from theoretical force-dependent velocity scaling, dispersion, and attenuation (12, 16, 40, 41). Quantifying the influence of both structural disorder and force heterogeneity on dispersion may enhance the usefulness of using dispersion as a nondestructive probe of the state of a granular medium.

Both experiments and numerical simulations for ultrasound propagation have shown that coherent waves (including compressional, shear, and rotational waves) decay in amplitude with propagation distance. Frictional contact loss, viscoelastic/viscous contact dissipation, wave scattering, and material damping have been assumed to be responsible for wave attenuation (25, 35, 42, 43). Local contact damping has been incorporated in discrete element simulations (18, 35, 44⇓–46), enabling analysis of the energy dissipation during wave propagation. Additionally, a stiffness matrix incorporating all particle contacts of a granular packing has been used to build a lattice-based model to examine the wave attenuation and scattering for 2D granular systems (20, 21, 47). However, limited experimental data are available to highlight the quantitative roles of packing disorder and force heterogeneity on wave attenuation in 3D.

The primary objective of this paper is to quantitatively investigate the relative influence of packing structure disorder and interparticle force heterogeneity on ultrasound wave velocity, dispersion, and attenuation in a disordered granular material. To achieve this, we combined in situ X-ray computed tomography (XRCT), three-dimensional X-ray diffraction (3DXRD), and ultrasound wave measurements during uniaxial compression of a 3D disordered granular material in a unique loading frame (9, 10, 48). To evaluate the relative roles of packing structure and interparticle forces on wave behavior, we compared measurements of velocity, dispersion, and attenuation with simulations having contact networks controlled by measured interparticle forces and contact networks having uniform interparticle stiffness (equivalent to uniform forces). We further compared predictions with those made on isotropically loaded granular crystals to isolate the roles of packing structure disorder. Our results revealed that both packing structure and interparticle forces play important roles in controlling coherent wave velocities, velocity scaling, and dispersion, while packing structure alone can explain the majority of the observed wave attenuation.

## Results and Discussion

### Wave Velocity.

The acoustic velocities of the coherent portions of ultrasound waves transmitted through a random packing of sapphire spheres under increasing normal compression were computed using cross-correlation analyses (*Materials and Methods* and *SI Appendix*, section 1C) and are presented in Fig. 1*A*. The measured acoustic velocity across load steps was found to follow a power-law dependence on the average vertical sample stress,

To simulate wave propagation through the packing at each load step, systems of equations were constructed for a spring network representing the experimentally measured structural disorder and interparticle force network (called “heterogeneous”) and for a separate spring network with measured structural disorder but uniform stiffness (called “uniform”). An input signal was put into the contact network by forcing all particles contacting the top platen to oscillate following a single half-period of a sinusoidal with a fixed frequency of 0.4 MHz. The sum of all contact forces on the bottom platen was taken to be the received signal. The same cross-correlation analyses used in experiments were used to measure the wave velocity through the material in simulations. The uniform network features a reduced total velocity compared with the heterogeneous network, reflecting the influence of the horizontally aligned forces (used in finding an average force and thus, average stiffness), which are typically lower than the vertically aligned forces. The network with heterogeneous forces exhibits a similar velocity scaling to experimental measurements, suggesting that both structural disorder and force heterogeneity are needed to precisely predict the acoustic velocity.

We further examined the contact forces between particles contacting the bottom platen in simulations. We calculated a wave velocity for each contacting particle by dividing the height of the sample by the time between the input signal peak and the first peak in the particle–platen force (25). Distributions of these velocities for load steps 4, 6, 8, and 10 for a 0.4-MHz input signal are provided in Fig. 1 *B* and *C* for the heterogeneous and uniform networks, respectively. A significantly broader spread is observed for the network containing heterogeneous forces (Fig. 1*B*), which cause local variations in the acoustic velocities and acoustic paths. This spread signifies the influence of force heterogeneity on dispersion and loss of phase coherence, as described in the following sections.

To further understand the path of ultrasound waves propagating through the granular packing, we constructed a unit cell for each individual particle using the positions of its neighbors and the corresponding interparticle forces, extracted from X-ray measurements. We calculated the theoretical acoustic wave velocity in the vertical (z) direction through each unit cell. Based on this velocity, we estimated the sample’s theoretical wave velocity by 1) averaging all unit cells’ theoretical wave velocities using experimentally measured contact forces [labeled M(heterogeneous)] and uniform contact forces [labeled M(uniform)], 2) averaging only over unit cells corresponding to a percolating network of “strong” interparticle forces, and 3) identifying the fastest theoretical wave path or “chain” through the sample.

The percolating network of strong interparticle forces was identified as follows. First, a force threshold was specified, and all contacts with a force magnitude above that threshold were retained. Next, we searched for a continuous path through these contacts that connected the bottom platen to the top platen (ignoring lateral particle-boundary contacts). Finally, we incremented the force level by +0.01*B* for load steps 4, 6, 8, and 10) for percolation, and the corresponding path was identified as the percolating network of strong interparticle forces. The velocity through the percolating network was computed by averaging the velocity in the z direction across all unit cells whose central particles belong to this network (the numbers of particles in networks for load steps 4, 6, 8, and 10 are provided in Fig. 2*B*).

The fastest theoretical chain, called chain in Fig. 1, was obtained by comparing all possible acoustic paths through the sample from the top to bottom platen. One particle contacting the top platen was labeled as the first particle of this acoustic path. Any particles contacting this first particle with a z position lower than the first particle’s z position were treated as a viable path for wave transmission. This process was repeated until a path from the top platen to the bottom platen was constructed. All such paths were evaluated, and the time of flight for each path was computed by *B*.

Velocity predictions by all three methods (averaging over unit cells, on percolating strong force networks, through fastest chains) are shown in Fig. 1*D*, along with experimentally measured velocities. Velocities obtained for the fastest chain and the force percolation network, respectively, served as approximate upper and lower bounds for the experimentally measured velocities. The fastest chain prediction closely matched measured wave velocities, particularly at high loads. This finding indicates that the actual force magnitudes on the fastest chain are needed to explain the wave speed of the coherent pulse. In contrast, the predicted velocities averaged over all unit cells for both M(uniform) and M(heterogeneous) fall far below the measured wave velocities. The case of M(uniform) resulted in similar velocities to those of M(heterogeneous) and demonstrated smaller deviations. This finding indicates that averaging over all particles’ structural disorder and force heterogeneity when making velocity predictions is not appropriate, as such averaging significantly reduces expected wave speeds.

The three velocity predictions employ contact networks with different levels of average contact force and different levels of structural disorder. As illustrated in Fig. 2*B*, the percolating force network and fastest paths tend to develop along contacts with higher than average contact forces and particles with higher than average vertical compressive stress magnitudes. To further emphasize the role of contact forces in furnishing wave velocities, we evaluated the scaling of both the predicted and measured velocities with respect to the effective vertical contact force, *E*, the predicted and measured velocities with different methods generally follow power-law relationships with

We note that the apparent Hertzian scaling of wave velocity with the effective vertical contact force on the fastest chain, *SI Appendix*, section 1A). The resulting fragments potentially introduce different contact laws and thus, different force-dependent velocity scaling. Other factors related to the contact law between particles, including rough contact (50), plastic deformation of particle contacts (8), and particle shape (27), can also contribute to the transition of the exponent value and were not investigated here.

The wave velocities predicted with the three methods described above were also compared in Fig. 1*E* with theoretical acoustic velocities of granular crystals to further evaluate the role of structural disorder. The granular crystals included simple cubic (SC; equivalent to monatomic chain), body-centered cubic (BCC), hexagonal close packed (HCP), and a 2D triangular lattice (Tri). Experimentally measured wave velocities for most load steps lie between those for SC and HCP structures and closely follow BCC predictions across load steps 4 to 11. This suggests that a BCC structure can be used to predict the scaling of ultrasound wave velocity with force at sufficient large force levels. However, the chain prediction more closely tracks the particular variations in the scaling of measured velocities with force, indicating that the particular details of local structure and force magnitudes on this chain can significantly improve predictions.

### Wave Dispersion.

To study wave dispersion, we computed spectrograms of received signals through the granular sample, received signals through pistons in contact, and received signals through simulated lattice networks. Corresponding spectrograms are shown in Fig. 3 *A*–*D* for the input signal of 0.8 MHz at load step 8. For the simulations (Fig. 3 *C* and *D*), experimental input signals (i.e., the Gaussian bursts shown in Fig. 7*A*) were scaled to limit the maximum particle displacement to <1 μm and injected into the heterogeneous and uniform networks. In both the experimental and numerical settings, wave reflection on the boundaries (e.g., platen–particle boundary at the top of the sample) cannot be avoided due to the sample size. Therefore, the spectrogram analysis is only performed for the first 5 μs (including three to five wave cycles) after the arrival of the coherent wave to minimize the effects of wave reflection, as in ref. 16. Contour lines at 20, 40, and 60% of the maximum signal amplitude in the first 5 μs of wave analysis are also shown in Fig. 3 *A*–*D*.

Theoretical dispersion predictions from averaging unit cell dispersion relations across the entire packing (with actual interparticle forces) and across only particles in the fastest chain are shown in Fig. 3 *E*–*H* as blue lines with square makers (the shaded region represents ±1 SD of arrival times) and brown lines with triangle markers, respectively. Results based on unit cells with uniform force are shown in green lines with diamond markers, the error bars of which (not shown) are slightly smaller than those for the unit cells with the actual interparticle forces. The dispersion relations based on network simulations are presented in pink lines with circle markers. Dispersion relations based on isotropically loaded SC and BCC granular crystals (having comparable coordination number with the studied sample) are provided for reference. For comparison, the numerical and analytical results were shifted to a common time axis with the experimental observations by adding the time of flight for the case of piston contact.

Compared with the input signal (similar to Fig. 3*A*), the spectrograms of received signals for experiments (Fig. 3*B*) and simulations (Fig. 3 *C* and *D*) demonstrated significant distortion indicative of dispersion, with delayed higher-frequency components. The earliest time, *B*, indicating frequency-dependent wave arrival time. These points are plotted for all input frequencies at load steps 4, 6, 8, and 10 in Fig. 3 *E*–*H*, with error bars representing SDs of

The measured dispersion varies significantly, not closely following any particular theoretical prediction. However, the trend of the curves, particularly for the 60% contours, more closely follows that from the average of unit cells in the entire packing with actual or uniform forces (green or blue curves) than, for instance, that of the fastest chain (brown curve). Specifically, the experimentally measured dispersion is stronger (i.e., arrival time is more delayed with frequency) than that predicted from the fastest chain. This suggests that the details of the entire packing, rather than the subset supporting the fastest wave, are needed to predict dispersion.

To highlight the importance of force heterogeneity on dispersion, we made further dispersion predictions by averaging the power spectrum of all particles in the spring networks with heterogeneous and uniform forces after inputting an initial condition. The dispersion for the heterogeneous network is shown in Fig. 4.*A* and demonstrates a broader spread than that for the uniform network shown in Fig. 4*B*, emphasizing enhanced loss of phase coherence for the excited modes due to force heterogeneity.

### Wave Attenuation.

Attenuation of input signals was evaluated by the transmission ratio as a function of input frequency and load step, as shown in Fig. 5*A*. We define the transmission ratio as the ratio of maximum spectrogram amplitude of the received signal (in arbitrary units) through the granular packing to that of the received signal through the pistons in contact with one another (this latter amplitude is nearly frequency- and compression-level independent). The transmission ratio is below 0.3 at nearly all load levels for input frequencies 1.2 and 1.6 MHz. The transmission ratio rises above 0.5 and close to 1.0 at the highest load levels for input frequencies of 0.4 and 0.8 MHz. We suspect the observed attenuation arises mainly from two mechanisms: Rayleigh scattering (20, 21) and viscoelastic damping at particle contacts (44, 47, 51).

To evaluate attenuation due to Rayleigh wave scattering, we calculated the scattering attenuation coefficient (*SI Appendix*, section 4A), η, by fitting autocorrelation functions of the longitudinal velocity with **6** using the ODE45 solver in Matlab. Modal damping ratios, ζ, were computed by

To study the relative impact of interparticle dissipation and Rayleigh scattering on wave attenuation, we compared the measured wave attenuation with three types of simulations: 1) network simulations with inferred forces and *D*, we observed 1) increasing transmission ratios, thus decreasing attenuation, with increasing compression; 2) more attenuation with increasing frequency; 3) comparable transmission ratios for contact networks with heterogeneous and uniform forces; 4) a slight increase in wave attenuation due to interparticle dissipation; and 5) Rayleigh scattering as the main attenuation mechanism, capable of explaining much of the experimentally observed attenuation.

To further demonstrate these points, we present attenuation coefficients due to scattering in Fig. 5 *E* and *F* and the modal damping ratio due to interparticle dissipation in Fig. 5*G*. The attenuation coefficient due to scattering on a uniform network (Fig. 5*F*) was slightly lower than that of the heterogeneous one (Fig. 5*E*) for high frequencies, suggesting that disorder due to force heterogeneity only slightly enhances attenuation due to wave scattering. Examples of velocity correlation calculations for the scattering-only case (no contact damping) with different initial conditions are shown in Fig. 5 *B* and *C* to illustrate the dominant role of wave scattering on wave attenuation. On the other hand, for all frequencies between 0 and 3 MHz, ζ was found to be smaller than 0.07 with a restitution coefficient of 0.9 [within ranges tabulated for sapphire (53)], suggesting the secondary role of interparticle damping in attenuation, as first illustrated in Fig. 5*G*.

## Conclusions

In this study, ultrasound responses of externally stressed disordered granular materials were experimentally and numerically investigated. Gaussian burst ultrasound waves with different center frequencies were sent through a disordered 3D granular material compressed to different levels of sample stress. In situ XRCT and 3DXRD enabled the measurement of contact fabric, crystal orientation, average intraparticle stress tensors, and interparticle forces at the corresponding levels of sample stress. The packing configuration and inferred interparticle forces were used to construct a linear spring network. The ultrasound response of this network was simulated to quantitatively determine the relative importance of packing structure disorder and interparticle force heterogeneity in controlling wave velocity, acoustic paths, dispersion, and attenuation. The findings of this study are summarized as follows.

1) The experimentally measured velocity of the coherent wave was quantitatively predicted by the theoretical fastest velocity on any chain through the packing and scaled with interparticle forces measured on this chain according to Hertzian contact theory. We conclude that both the packing disorder and interparticle force heterogeneity are responsible for ultrasound wave velocities and velocity scaling in disordered granular media.

2) Experimentally measured velocities were found to be dependent on frequency, indicating significant dispersion. No dispersion predictions closely matched those that were experimentally measured. Simulations with uniform and heterogeneous forces suggested that packing disorder significantly increases dispersion, with force heterogeneity playing a secondary role primarily discernible in dispersion relation plots. We conclude that packing disorder is primarily responsible for dispersion, with force heterogeneity playing a secondary role.

3) Additionally, the measured force heterogeneity and corresponding simulation results suggest that wave attenuation was primarily due to wave scattering. Force heterogeneity enhanced loss of phase coherence but played a negligible role in overall wave attenuation. Therefore, we conclude that packing disorder is primarily responsible for attenuation of longitudinal ultrasound waves.

The findings of this research provide insight into the microscopic features associated with particular wave behaviors and may support the development of materials for wave manipulation and methods for nondestructive testing.

## Materials and Methods

### Experimental Setup.

The granular material we studied was a sample of 621 single-crystal sapphire spheres (Bird Precision) with radii ranging from 94 to 103 μm and uniform roughness below 15 nm (for rms roughness of the particle surface). The granular material was prepared by pouring particles into an aluminum cylinder of 1.5-mm inner diameter, 3.00-mm outer diameter, and 8-mm height. The sample developed locally crystallized regions but was largely disordered. The sample was then placed in a custom-built uniaxial load frame described next and shown in Fig. 6*A*. The load frame features a linear actuator with encoder (Haydon Kerk Size 34 Stepper Motor), a broadband piezoelectric transducer and receiver (Olympus V133-RM), each in an aluminum encasement, stainless steel compression platens, and a load cell (Futek LCM200), all mounted to an aluminum plate and a polycarbonate tube. To drive and monitor the load frame components, we employed a National Instruments (NI) system, which contained a computer (PXIe-1082 Express Chassis), a motion controller (PXI-7342 and UMI-7764), an arbitrary waveform generator (PXIe-5423), and an oscilloscope (60 MHz; PXIe-5105). The system is shown in its operational configuration outside the X-ray hutch in Fig. 6*B*. The load frame was mounted on the translation and rotation stage present at the Advanced Photon Source (APS) beamline 1-ID-E (Fig. 6*C*) throughout the experiment, and XRCT and 3DXRD measurements were performed between increments of sample strain. Fig. 6*D* illustrates the experimental geometry schematically.

### Mechanical Loading.

The granular sample was loaded mechanically in small strain increments, ranging from 0.005 to 0.01, by microstepping the linear actuator, which drove the stainless steel platen located above the granular material downward in the −*z* direction. The loading rate was limited to approximately 1 μm/s to ensure that the loading remained quasistatic (strain rate of

### Ultrasound Transmission.

After each 3DXRD measurement, an ultrasound signal was transmitted through the sample using the piezoelectric transducer above the sample. The signal was received using the piezoelectric receiver below the sample. Four types of signals were transmitted, each of which was transmitted and received five times. Each signal was generated by sending a time-varying voltage to the piezoelectric transducer: a sinusoid convolved with a Gaussian, yielding “Gaussian bursts” (16) with center frequencies of 0.40, 0.80, 1.20, and 1.60 MHz. The center frequencies were selected to fall below the cutoff frequencies of a 1D chain of 200-μm diameter sapphire spheres under axial compression of 0.50 N, which were calculated to be 3.27 MHz (*SI Appendix*, section 5). Features of input signals (in the time domain and frequency domain) are provided in Fig. 7 *A* and *B*. Note that all frequencies used in this study are angular frequencies unless otherwise specified.

Wave velocities for the coherent portion of transmitted waves were calculated by combining cross-correlation analysis of the average received signals and the input signal and cross-correlation analysis of the signals received themselves (*SI Appendix*, section 1C). The resulting cross-correlation coefficients at center frequencies of 0.4 and 0.8 MHz typically featured a sharp rise at the arrival of the coherent signal. However, this sharp rise could not be captured in some signals with center frequencies of 1.2 and 1.6 MHz that demonstrated significant attenuation. At these frequencies, we therefore inspected the received signals, and we selected a peak in the cross-correlation that corresponded to a cross-correlation below 0.75 but was consistent with arrival times at previous load levels. To compute the wave speed through the granular sample, we subtracted the arrival time of the coherent wave obtained from two pistons in contact,

### Inference of Interparticle Forces.

By combining XRCT and 3DXRD datasets, all interparticle and particle-boundary force vectors were calculated using a modified version of the minimization procedure described in refs. 9, 10, 54, and 55. The minimization procedure involved an effort to satisfy equilibrium (Eqs. **2** and **3**) and stress–force relations (Eq. **4**) for each particle:**4**) and **2** and **3**) are *SI Appendix*, section 2. We denote the overall number of contacts in the system *B*. These renderings employ linear scaling of force vector lengths, widths, and opacities with force vector magnitude.

### Dynamical Responses of Contact Network.

With the experimentally measured contact network and interparticle forces described above, we constructed a spring network to further study the system’s dynamical response at each level of applied sample stress in the experiment (20, 21). We assumed that the contact fabric remained unchanged by wave propagation. We further assumed that, during ultrasound wave propagation, particles oscillated with respect to their equilibrium positions (identified experimentally). Finally, we calculated a linear contact stiffness at each interparticle contact by linearizing a Hertzian contact law at the measured interparticle force.

The governing equations of wave propagation in the granular sample are given by*SI Appendix*, section 3 and were notably constructed using the experimentally measured packing structure and in some calculations, the experimentally measured interparticle forces. Eq. **6** was numerically integrated using the ordinary differential equation solver ODE45 in Matlab, for given initial conditions.

We focused numerical simulations on studying longitudinal wave propagation along the vertical (*z*) direction of the sample to mirror input and output signals controllable and measurable in experiments. To investigate the relative influence of packing structure and interparticle force heterogeneity on observed wave behaviors, we simulated wave propagation in two ways. The first way used both the experimentally measured contact networks and interparticle forces in the construction of mass, damping, and Hessian matrices. The second way used the experimentally measured contact network but assumed a uniform contact stiffness, derived from the average interparticle force in the packing, at all interparticle contacts. Although the latter system would be out of equilibrium if its interparticle contact forces were derived from the corresponding interparticle stiffness values, it provided a convenient method of decoupling the influence of force heterogeneity from that of structural disorder. Further theoretical calculations with granular crystals then provided a comparison of these results with scenarios without structural disorder.

### Data Availability.

All data are described in *SI Appendix*, section 6. Particle data are provided as Datasets S1–S5, and all data, including particle and waveform data, can be accessed through www.zenodo.org (60).

## Acknowledgments

C.Z. and R.C.H. acknowledge support from The Johns Hopkins University’s Whiting School of Engineering. E.B.H. and R.C.H. acknowledge support for the design and construction of the load frame from Lawrence Livermore National Laboratory’s Laboratory Directed Research and Development Program Grant 17-LW-009. We thank Dr. J. Almer, Dr. P. Kenesei, and Dr. J. S. Park for help in executing the experiments. We acknowledge the APS for synchrotron beamtime under proposal GUP-59127. Use of the APS, an Office of Science User Facility operated for the US Department of Energy (DOE) Office of Science by Argonne National Laboratory, was supported by US DOE Contract DE-AC02-06CH11357.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: rhurley6{at}jhu.edu.

Author contributions: C.Z., E.B.H., and R.C.H. designed research; C.Z., E.B.H., and R.C.H. performed research; C.Z. and R.C.H. analyzed data; E.B.H. provided extensive feedback on manuscript; and C.Z. and R.C.H. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

Data deposition: All raw data have been uploaded through the repository Zenodo, www.zenodo.org (10.5281/zenodo.3785083).

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2004356117/-/DCSupplemental.

- Copyright © 2020 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- I. L. Vér,
- L. L. Beranek

- I. L. Vér

- ↵
- ↵
- ↵
- G. Michlmayr,
- D. Cohen,
- D. Or

- ↵
- G. Gantzounis,
- M. Serra-Garcia,
- K. Homma,
- J. Mendoza,
- C. Daraio

- ↵
- ↵
- ↵
- V. Nesterenko

- ↵
- R. Hurley,
- S. Hall,
- J. Andrade,
- J. Wright

- ↵
- C. Zhai,
- E. Herbold,
- S. Hall,
- R. Hurley

- ↵
- ↵
- M. Hiraiwa,
- S. Wallen,
- N. Boechler

- ↵
- E. Herbold,
- J. Kim,
- V. Nesterenko,
- S. Wang,
- C. Daraio

- ↵
- B. P. Lawney,
- S. Luding

- ↵
- ↵
- C. Coste,
- B. Gilles

- ↵
- A. Merkel,
- S. Luding

- ↵
- H. Cheng,
- S. Luding,
- K. Saitoh,
- V. Magnanimo

- ↵
- E. Wang et al.

- ↵
- S. Gelin,
- H. Tanaka,
- A. Lemaître

- ↵
- K. Saitoh,
- R. K. Shrivastava,
- S. Luding

- ↵
- P. Poorsolhjouy,
- A. Misra

- ↵
- A. Leonard,
- L. Ponson,
- C. Daraio

- ↵
- ↵
- V. Langlois,
- X. Jia

- ↵
- C. Coste,
- B. Gilles

- ↵
- E. T. Owens,
- K. E. Daniels

- ↵
- J. Goddard

- ↵
- H. A. Makse,
- N. Gland,
- D. L. Johnson,
- L. Schwartz

- ↵
- B. O. Hardin,
- F. Richart Jr

- ↵
- B. O. Hardin,
- W. L. Black

- ↵
- M. A. Zimmer,
- M. Prasad,
- G. Mavko,
- A. Nur

- ↵
- A. Gheibi,
- A. Hedayat

- ↵
- J. Santamarina,
- G. Cascante

- ↵
- J. O’Donovan et al.

- ↵
- S. Van den Wildenberg,
- M. van Hecke,
- X. Jia

- ↵
- E. Somfai,
- J. N. Roux,
- J. H. Snoeijer,
- M. Van Hecke,
- W. Van Saarloos

- ↵
- ↵
- ↵
- R. K. Shrivastava,
- S. Luding

- ↵
- R. F. Waymel,
- E. Wang,
- A. Awasthi,
- P. H. Geubelle,
- J. Lambros

- ↵
- ↵
- E. B. Herbold,
- V. F. Nesterenko

- ↵
- D. O. Potyondy,
- P. Cundall

- ↵
- A. Di Renzo,
- F. P. Di Maio

- ↵
- A. Di Renzo,
- F. P. Di Maio

- ↵
- D. L. Johnson,
- Y. Hu,
- H. Makse

- ↵
- R. Hurley et al.

- ↵
- ↵
- K. L. Johnson,
- K. L. Johnson

- ↵
- K. J. Hanley,
- C. O’Sullivan

- ↵
- ↵
- ↵
- ↵
- R. Hurley,
- J. Lind,
- D. Pagan,
- M. Akin,
- E. Herbold

- ↵
- R. Hurley,
- E. Marteau,
- G. Ravichandran,
- J. E. Andrade

- ↵
- R. Hurley,
- K. Lim,
- G. Ravichandran,
- J. Andrade

- ↵
- E. Oberg,
- F. D. Jones,
- H. L. Horton,
- H. H. Ryffel,
- J. H. Geronimo

- ↵
- C. Sandeep,
- K. Senetakis

- ↵
- C. Zhai,
- E. B. Herbold,
- R. C. Hurley

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics