Physical Review E 92, 022203 (2015)

Softening of Stressed Granular Packings with Resonant Sound Waves

C. J. Olson Reichhardt1, L.M. Lopatina1, X. Jia2, and P.A. Johnson1

1Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
2Institut Langevin, ESPCI ParisTech, CNRS UMR 7587, 1 rue Jussieu, 75005 Paris, France, EU

(Received 16 June 2014; revised manuscript received 25 February 2015; published 5 August 2015)

We perform numerical simulations of a two-dimensional bidisperse granular packing subjected to both a static confining pressure and a sinusoidal dynamic forcing applied by a wall on one edge of the packing. We measure the response experienced by a wall on the opposite edge of the packing and obtain the resonant frequency of the packing as the static or dynamic pressures are varied. Under increasing static pressure, the resonant frequency increases, indicating a velocity increase of elastic waves propagating through the packing. In contrast, when the dynamic amplitude is increased for fixed static pressure, the resonant frequency decreases, indicating a decrease in the wave velocity. This occurs both for compressional and for shear dynamic forcing, and is in agreement with experimental results. We find that the average contact number Zc at the resonant frequency decreases with increasing dynamic amplitude, indicating that the elastic softening of the packing is associated with a reduced number of grain-grain contacts through which the elastic waves can travel. We image the excitations created in the packing and show that there are localized disturbances or soft spots that become more prevalent with increasing dynamic amplitude. Our results are in agreement with experiments on glass bead packings and earth materials such as sandstone and granite, and may be relevant to the decrease in elastic wave velocities that has been observed to occur near fault zones after strong earthquakes, in surficial sediments during strong ground motion, and in structures during earthquake excitation.
I. INTRODUCTION
II. SIMULATION
III. RESULTS
IV. DISCUSSION
V. SUMMARY
References


I.  INTRODUCTION

Granular matter has very unusual properties and can exhibit liquidlike behavior by flowing under certain excitations, while it can have a solidlike resistance to shear for other excitations. The jamming phase diagram, originally proposed by Liu and Nagel [1], provides a concise description of the transition from jammed to unjammed states as a function of density, temperature, or loading. Granular matter can exhibit fragile properties in which the response depends on the loading history [2,3,4,5]. A number of studies have focused on the loading axis by applying a shear to the granular packing and studying the unjamming of the packing above a certain shear level [6,7,8,9,10,11,12,13,14]. Much work has also been performed on calculating the normal or soft modes of granular packings [15,16,17,18], with particular emphasis on the emergence of low frequency modes close to the jamming transition.
Most previous studies of granular matter under shear loading have considered primarily quasistatic shear involving plastic granular flow, applying a single direction of shear [6,7,8,9,10,11,12,13,14,19] or cyclic shear [20,21,22,23,24,25,26]. Relatively little work has been performed on the nonlinear elastic/plastic sound vibration regime of granular matter in the dense state where macroscopic plastic distortions of the grain positions do not occur [27]. Such dynamic shearing of dense granular packings is of particular interest in connection with prominent effects in surficial sediments from strong ground shaking from earthquakes [28,29] as well as the behavior of fault gouge material in response to earthquake forcing [30,31]. Gouge is a disordered granular matter that often exists along and within the fault core; it is produced by the long-term grinding of the tectonic plates against each other via a process known as communition [32]. It has been hypothesized to play a role in certain behaviors related to earthquake faults, such as a delayed triggering response in which a large distant earthquake can initiate an earthquake after a waiting time of days or months [33]. Moreover, large earthquakes have been observed to cause a long-lived, shaking induced depression of the elastic wave velocity in the mid to upper crust in localized areas, which slowly recovers over time [34,35] as well as in near surface sediments [36,37,38].
Experiments performed with glass bead packs [30,4] and on natural materials such as sandstone show a similar decrease in the elastic wave velocity under oscillatory or dynamic loading [40,3,39]. One common method for probing the softening of the elastic wave velocity is the use of nonlinear resonant ultrasound spectroscopy, which can measure the nonlinear elastic state of a rock or a glass bead pack [3]. The frequency of an applied wave of fixed amplitude A is swept or stepped across a resonant mode of the sample, and the resulting signal is measured on the opposite side of the sample [40]. In diverse materials including Berea sandstone, Lavoux limestone, or synthetic slate, the resonant frequency drops with increasing amplitude of the driving wave A [40,41,42,43], and this indicates a drop in the velocity at which an elastic wave pulse travels through the sample [46,44,4,45]. In granular matter, when the static confining pressure is increased, the elastic wave velocities increase [47,48,49,50]. Early work on elastic wave or sound propagation in a glass bead packing suggested that the detailed contact structure of grains within the packing play an important role in wave transmission [51,52], particularly in short-wavelength wave scattering [48]. For long-wavelength coherent waves, effective medium theory indicates a link between the coordination number (the average number of contacts per grain) and the elastic wave velocity [53,54,47,49]. Simulations and experiments with 3D packings indicated that the effective medium theory fails to account quantitatively for the shear elastic modulus when the affine approximation breaks down at low static pressures or high dynamical amplitudes [55,56,57]. Much is understood regarding grain behavior under shear [58,59]; however, despite a number of studies on sound wave propagation in two and three dimensional packings [60,61], a detailed microscopic understanding of the elastic wave velocity evolution with driving amplitude has not yet been obtained.
In this work, we study confined granular packings subjected to both a static pressure and to dynamic loading achieved by applying an ac compressional or shear loading to our model system. We show that the behavior of the elastic wave propagation matches what has been observed experimentally, and demonstrate that changes in the contact number of the grains are correlated with the elastic wave propagation changes. We illustrate the dynamical motion of the grains and discuss the implications of our work to dynamical triggering studies performed using earthquake catalogs.
Fig1.png
Figure 1: (Color online) Schematic of system showing the four confining walls. The bottom and side walls (dark grey) are fixed, while the top wall (light yellow) is subjected both to a static confining force Fsl=−p∧y (thick light red arrow) and a sinusoidal dynamic loading force (thin dark blue arrows) in either the compressional (center arrow) or shear (right arrow) direction.

II.  SIMULATION

We consider a two-dimensional (2D) packing of N=700 disks with Hertzian contact interactions [62,63]:
Fijgg=g⎡
⎣ 		1
2
 (Di + Dj) − rij⎤
⎦ 		3/2

 
^
r
 
ij 
(1)
where g=10 is the elastic constant, Di(j) is the diameter of particle i(j), rij=Ri−Rj, rij=|rij|, and ∧rij=(Ri−Rj)/rij. The value of g can be related to the Hertzian expression [62] as: g=(8 G √{Dg/2})/3(1−νg) where G is the grain shear modulus, Dg/2 is the radius of the smaller of the two grain sizes, and νg is the Poisson ratio of the material from which the grains are made. The two grains interact only when they are in contact with each other, for rij ≤ (Di+Dj)/2. To avoid crystallization of the packing, we use a bidisperse assembly of grains consisting of a 50:50 mixture of grains with a radius ratio of 1:1.4. We measure length in units of a0, the diameter of the smaller of the two sizes of grains. Following Gallas et al. [64] and similar works by Cundall-Strack [63] and Somfai et al. [60], we include viscous damping forces at the contacts to mimic the dissipation of sound waves in granular packings:
Fijn=−γn meff(rij ·vij)
^
r
 
ij 
 /|rij|
(2)
in the normal direction and
Fijt=−γs meff(tij ·vij)
^
t
 
ij 
 /|tij|
(3)
in the tangential direction. Here γn=0.1 and γs=0.1 are the dissipation coefficients, meff is the effective mass of the two grain system, vi(j) is the velocity of grain i (j), vij=vi−vj is the relative velocity, and
tij=⎛
⎜
⎜
⎝
−rijy
rijx
⎞
⎟
⎟
⎠ 		.
(4)
As mentioned in Ref. [64], the Coulomb friction proportional to the normal force is neglected in the present simple model of interaction; however, the tangential damping force in Eq. (3) may mimic to some degree the effect of static friction to halt the tangential relative motion. We employ a granular dynamics simulation technique to integrate the equations of motion for each particle, given by
Mi
⋅⋅
r
 
i 
 =
∑
j 
 Fijgg+ Fijn + Fijt + Fg
(5)
and
Ii
⋅⋅
ϕ
 
i 
 =
∑
j 
 δCij
(6)
where Mi is the mass of grain i, Ii is the radius of gyration of grain i, ϕi is the angular degree of freedom of grain i, Fg is a gravitational force term, and δCij is the torque exerted on a grain through contact with other grains. We measure distances in units of a0, time in units of t0=√{a0/gg}, where gg is the gravitational acceleration constant, velocity in units of v0=√{gga0}, force in units of F0=m0gg, where m0 is the mass of the smaller of the two sizes of grains, elastic constants in units of k0=m0gg/a0, and stress in units of θ0=m0gg/a02.
The grains are confined within four walls in our simulation box as illustrated in Fig. 1. Wall interactions are modeled using image grains that are the reflection of a grain in contact with the wall to the other side of the wall. The bottom and side walls are held at fixed positions, while the top wall is a piston used to apply a static load Fls=−p∧y normal to the wall modulated by a dynamic load of the form Fld=A sinωt∧α, where α = y for compressional loading and α = x for shear loading. The position of the top wall is allowed to vary according to the sum of the total load Fload=Fls+Fld and the effective forces exerted on the image grains by the actual grains. A static pressure value of p=0.005 corresponds to a downward force on individual grains touching the top wall of 1.92×10−4. This pressure is transmitted throughout the packing and opposed by an effective force arising from the fixed bottom wall, so that individual grains move very little when the static pressure is modified. Similarly, a dynamic compressional amplitude of A=0.030 contributes an oscillating force of magnitude 1.15×10−4 to each grain touching the top wall, such that the motion of individual grains remains much smaller than a0. It is important to note that due to the confinement, once the grains have been prepared in the packing they are not able to rearrange their positions but can only make slight shifts relative to their neighbors, which do not change. We measure the net force exerted by the grains on the top wall, ft(t), and the bottom wall, fb(t), for fixed A while slowly stepping ω across a resonant frequency ω0. For each driving frequency ω, we collect data during a period of 50 drive cycles. We then compute the power spectrum S(ν)α=|∫fα(t)e−i2πνtdt|2, with α = t,b, of both ft(t) and fb(t), and obtain the response in the form of the relative or normalized amplitudes of the output to input signals at the driving frequency, η(A)=S(ν = ω/2π)b/S(ν = ω/2π)t. As a ratio, η is a dimensionless quantity that provides information about the relative spectral weight of the response; we use it to quantify changes in the resonant frequency of the system.
Fig2.png
Figure 2: (Color online) Results from the compressional dynamic simulation. (a) Scaled amplitude of detected response η vs driving frequency ω at A=0.025 for increasing static pressure p=0.0050, 0.0055, 0.0060, 0.0065, 0.0070, 0.0075, 0.0080, 0.0090, 0.0100, 0.0110, 0.0120, 0.0130, and 0.0140 (bottom to top), showing a shift of the resonant peak ω0 to higher frequencies with increasing p. (b) Resonant frequency ω0 vs static pressure p on a log-log scale, indicating an increase in the elastic wave velocity with increasing static pressure, for different values of the dynamic amplitude A=0.015, 0.020, 0.025, and 0.030, from top to bottom. Dashed lines are fits to ω0 ∝ pβ with β ≈ 0.35. (c) η vs ω at p=0.0050 for increasing dynamic amplitude A=0.010, 0.012, 0.014, 0.015, 0.016, 0.018, 0.020, 0.022, 0.024, 0.025, 0.026, 0.028, 0.030, 0.032, 0.034, 0.036, 0.038, and 0.040 (top to bottom), showing a shift of the resonant peak to lower frequencies with increasing A. (d) Resonant frequency ω0 vs dynamic amplitude A, on a log-log scale, indicating a decrease in the elastic wave velocity with increasing A, for different values of the static pressure p=0.005, 0.007, 0.009, and 0.011, from bottom to top. Dashed lines are fits to ω0 ∝ A−β with β ≈ 0.4.
To prepare our system, we remove the left wall of the sample, hold the top or piston wall in a fixed position, and fill the system with a granular gas. We add a gravitational force term Fg=migg∧x to each grain and allow the grains to settle into a dense packing. We then close the left wall and change the gravitational force to Fg=−migg∧y to force the grains toward the bottom wall of the packing; we then permit the piston or top wall to move and incrementally apply a static pressure to the piston, allowing the granular arrangement to settle to a state of no net motion between pressure increments. Once we have reached the desired static pressure level p, we add a sinusoidal term to the force exerted by the piston, resulting in a sinusoidal motion of the piston. We permit the system to oscillate for 20 cycles in order to eliminate any transient effects, and then measure the wall forces ft(t) and fb(t) during a period of 50 cycles. In a given run we perform a frequency sweep by holding the amplitude of the oscillation of the piston fixed but increasing the frequency of the oscillation to a new value after each set of 70 cycles. To change the static pressure or the magnitude of the dynamic forcing, we start with a fresh uncompacted sample in each case. This avoids a systematic increase in density that could otherwise occur after each frequency sweep. Our simulation measurement protocol is similar to that used in the experiment described in Ref. [30], where the resonance compressional P waves are observed.
Fig3.png
Figure 3: (Color online) Experimental glass bead pack results modified from Ref. [30]. (a) η vs f=ω/2π for increasing dynamical amplitude A = 10 mV, 70 mV, 130 mV, 190mV, 250 mV, 310 mV, and 370 mV, from top to bottom. (b) Normalized ∆ω0 vs A, in mV, for samples with increasing p from bottom to top. The elastic wave velocity decreases with increasing A in each case, but the overall magnitude of the decrease becomes smaller as p increases.

III.  RESULTS

Fig4.png
Figure 4: (Color online) Results from the shear dynamic simulation. (a) η vs ω at A=2.0 for increasing static pressure p=0.3, 0.4, 0.5, 0.6, 0.7, and 0.8, from left maximum to right maximum, showing a shift in ω0 to higher frequencies with increasing p. (b) Resonant frequency ω0 vs static pressure p on a log-log scale, indicating an increase in the elastic wave velocity with increasing static pressure, for different values of the dynamic amplitude A=3, 5, 7, and 10, from top to bottom. Dashed lines are fits to ω0 ∝ pβ with β ≈ 0.25. (c) η vs ω at p=0.5 for increasing dynamic amplitude A=1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, and 10.0, from left maximum to right maximum, showing a shift of ω0 to lower frequencies with increasing A. (d) Resonant frequency ω0 vs dynamic amplitude A, on a log-log scale, indicating an increase in the elastic wave velocity with increasing A, for different values of the static pressure p=0.3, 0.4, 0.5, 0.6, 0.7, and 0.8, from top to bottom. Dashed lines are fits to ω0 ∝ A−β with β ≈ 0.1.
We first compare measurements of the dynamic response η in the compressional and shear oscillatory simulations and in experiments. In Fig. 2(a) we plot the normalized amplitude η as a function of driving frequency ω for fixed compressional dynamic amplitude A=0.025 and static pressures ranging from p=0.005 to p=0.0140. Here the resonant frequency ω0 shifts to higher values as the static pressure is increased. For the clamped boundary conditions we consider, the resonant frequency ω0 is related to the elastic wave velocity V by ω0=πV/L, where L is the sample thickness along the y-axis [3,30]. This indicates that the elastic wave velocity is increasing with increasing static pressure, in agreement with previous observations. By identifying the value of ω0 from each curve, we construct a plot of ω0 versus p shown in Fig. 2(b) for dynamic amplitudes ranging from A=0.015 to A=0.030. The resonant frequency increases with increasing static pressure roughly as a power law with slope β ≈ 0.35; however, there is an overall downward shift in the resonant frequency as the compressional dynamic loading A increases. In Fig. 2(c) we plot η versus ω for the compressed system at fixed static pressure p=0.0050 and dynamic amplitudes ranging from A=0.010 to A=0.040. Here, the peak value ω0 decreases in frequency with increasing dynamic amplitude A, indicating that the elastic wave velocity is decreasing with increased dynamic driving. This softening of the system with dynamic driving is more clearly shown in Fig. 2(d) where we plot ω0 versus A for values of p ranging from p=0.005 to p=0.011. The softening is very robust and appears for each value of p. For comparison, we illustrate in Fig. 3(a) the experimentally obtained values of η as a function of frequency for different dynamical amplitudes A. The resonant frequency decreases with increasing dynamic amplitude. This is more clearly shown in Fig. 3(b), where we plot ∆ω0, the shift in ω0 from a reference value, versus the dynamic amplitude A for different values of static pressure [30]. In each case, the resonant frequency decreases with increasing dynamic amplitude in agreement with the simulation results.
We find similar behavior for a system in which the top plate is dynamically sheared in the direction transverse to the applied static pressure. In Fig. 4(a) we illustrate representative η vs ω curves at A=2.0 and increasing static pressure p in the sheared system. Figure 4(b) shows a log-log plot ω0 versus p curves for values of A ranging from 3 to 10 in the same system. We observe a power law behavior ω0 ∝ pβ with β ≈ 0.25, a somewhat smaller exponent than in the dynamically compressed system. The resonant frequency increases with increasing p, indicating an increase in the elastic wave velocity with increasing static pressure. We note that significantly larger static pressures must be applied to the dynamically sheared system than to the dynamically compressed system in order to obtain a wave signal that propagates through the entire packing and is measurable on the bottom plate. For increasing dynamic amplitude, the resonant frequency decreases, as illustrated in Fig. 4(c) for p=0.5 and a range of values of A. The decrease is slower than linear, as shown in Fig. 4(d) where we plot ω0 versus A for different values of p in the dynamically sheared system. These simulation results are also in excellent agreement with our experimental observations on shear resonant modes [65].
Fig5.png
Figure 5: Results from the compressional dynamic simulation. The total magnitude of the frequency shift across our measured range of A, ∆ω = ω0(A=0.005)−ω0(A=0.030) vs p shows two regimes of frequency shift behavior. At low p, ∆ω0 increases with increasing static pressure, while for p > 0.005, ∆ω0 decreases with increasing p. The higher static pressure regime agrees with the experimental response.
In Fig. 3(b) we find that experimentally, the overall magnitude of the decrease in f0 with increasing dynamic amplitude, ∆ω0, becomes smaller when the static load p is increased. The same behavior occurs in the simulations, as shown in Figs. 2(d) and 4(d). For the compressional dynamic simulations, we find that if we decrease the static pressure p to very small values, ∆ω0 passes through a peak value and then begins to decrease with decreasing p instead of increasing. This is shown in Fig. 5, where we plot ∆ω0=ω0(A=0.005)−ω0(A=0.030) as a function of static pressure p in the compressional system. For p < 0.005, ∆ω0 increases with increasing static pressure, while for p > 0.005, ∆ω0 decreases with increasing static pressure. The higher p behavior agrees with the experimental results [30]. The behavior at low p < 0.005 will be further explored in a future work. In general, the elastic response of a granular packing under weakly confining pressure may seriously deviate from the predictions of effective medium approaches. For the remainder of this paper, we will focus on the higher pressure regime with p > 0.005 in the compressional dynamic simulation.
Fig6.png
Figure 6: (Color online) Results from the compressional dynamic simulation. (a) η vs ω for A=0.010, 0.015, 0.020, 0.025, and 0.030, from top to bottom, at fixed p=0.005. (b) Average contact number Zc in the packing vs ω for A=0.010, 0.015, 0.020, 0.025, and 0.030, from top to bottom, at fixed p=0.005. There is a pronounced dip in Zc that increases in magnitude with increasing A. (c) Value of 〈Zc〉 at ω = ω0 as a function of dynamic amplitude A for p=0.005, 0.007, 0.009, and 0.011, from bottom to top.
We next compare the response of the system with the average coordination number 〈Zc〉 = N−1∑Zi of the packing, where Zi is the number of particles in direct contact with particle i. In Fig. 6(a) we plot the normalized amplitude versus driving frequency in a compressional system with p=0.005 for dynamic amplitudes ranging from A=0.010 to A=0.030. As before, we observe that the resonant frequency ω0 decreases with increased A. In Fig. 6(b) we show the corresponding 〈Zc〉 versus driving frequency. For larger driving amplitudes A > 0.015, near the resonance frequency ω ≈ 7×10−6 there is a dip in 〈Zc〉 which increases in magnitude with increasing A, indicating that the packing is becoming looser as the dynamic amplitude increases. The decrease of 〈Zc(ω0)〉, the value of 〈Zc〉 at the resonant frequency, is shown in Fig. 6(c) as a function of A. There is a slight increase in 〈Zc(ω0)〉 as the static pressure increases, but there is a clear decrease in 〈Zc(ω0)〉 with increasing A. As suggested by the effective medium theory, the elastic wave velocity is proportional to the coordination number [53,54,47,49]. Thus the reduction of the number of the contacts in the packing at a large driving (A > 0.015) is the main physical reason for the decrease in the elastic wave velocity with increasing dynamical amplitude in the granular packing. Higher static pressure forces more grains into direct contact. In contrast, larger amplitudes of dynamical forcing A tend to break contacts in the packing. At small A, we note that the dip in 〈Zc〉 occurs below the resonant frequency. This suggests that the elastic softening in this regime is likely determined both by the contact stiffness decrease and the contact number reduction [4]; weak contacts between grains may be important for low driving amplitudes but are destroyed at higher driving amplitudes.
Fig7.png
Figure 7: (Color online) Contour plot showing the regions of the sample undergoing the largest amount of motion in the dynamically compressed system at p=0.005 at the resonant frequency. Colors (dark to light) indicate the magnitude of the motion at each spatial location; the colorbar scale has been multiplied by 100. (a) A=0.010. (b) A=0.020. (c) A=0.025. (d) A=0.030.
Fig8.png
Figure 8: (Color online) Images of grain positions in the dynamically compressed system at p=0.005. Individual grains are colored according to their average net displacement (dark blue: smallest motion; light blue: largest motion); in each case, the net displacement is much smaller than the grain radius. Grains that are outlined in light red have zc < 4. Lines indicate the force contact network; force contacts with the walls are not shown. Heavy lines are strongly stressed bonds that experience forces that are greater than one standard deviation above the mean force in the packing. Heavy light red lines indicate strongly stressed bonds that pass through a grain with Zi < 4. Left column: Low frequency state with ω = 5.0×10−6. Center column: Resonant frequency state at ω = ω0. Right column: High frequency state with ω = 1.475×10−5. (a,b,c) A=0.010, where ω0=1.2×10−5. (d,e,f) A=0.020, where ω0=8.75×10−6. (g,h,i) A=0.025, where ω0=8.5×10−6. (j,k,l) A=0.030, where ω0=7.5 ×10−6.
Finally, we find that the motion of the grains under excitation takes two forms. The bulk of the packing responds collectively, with a large section of the packing moving coherently in response to the dynamic forcing at and near resonance. We also observe isolated soft spots or rattler areas where an individual grain has a much higher amplitude of motion than the grains that surround it. These soft spots tend to contribute additional damping to the propagating elastic wave signature. To identify the soft spots in the compressional dynamic simulation, we compute δri = max(ri(t)−ri(0)), which is the maximum displacement of an individual particle from its average equilibrium position ri(0) for a given driving frequency, amplitude, and static pressure. In Fig. 7 we show contour plots of the value of δri at the resonant frequency for packings with p=0.005 and A=0.010 to 0.030. The position in the packing is indicated on the x and y axes, while the coloring indicates the value of δri. The number and density of soft spots, indicated by local maxima in δri, increases with increasing driving amplitude, with a single spot in Fig. 7(a), two in Fig. 7(b), three in Fig. 7(c), and more than four in Fig. 7(d). This proliferation of soft spots contributes to the drop in ω0 with increasing driving amplitude. The soft spots can also be imaged as in Fig. 8, where we illustrate the grain positions along with the force chain network in the packing, and color each grain according to its relative value of δri. The soft spots in Fig. 7 coincide with rattler particles in the center column of Fig. 8, where we show the grain configurations at the resonant frequency ω0. For comparison, in the left column of Fig. 8 we show the system response at a low frequency of ω = 5.0×10−6, and in the right column we show the response at a high frequency of ω = 1.475×10−5. Away from resonance, there is significantly less motion throughout the packing. The highest stresses in the packing become more heterogeneously distributed at resonance as the driving amplitude increases.

IV.  DISCUSSION

Different regimes of the elastic wave velocity c (compressional or shear) through 2D or 3D bead packings as a function of applied static pressure pext have been observed previously in experiment, with c ∝ p1/6 at high pressures but c ∝ p1/4 at low pressure [53,47,66,67,68,69,70,71]. The results shown in Figs. 2 and 4 are more consistent with the low pressure regime. For the case of monodisperse disk packings, this effect was treated analytically in Ref. [72]. This scaling behavior might be correlated to the change in static pressure with the average contact number in the packing [47,73] as has been confirmed numerically [49,60,74] and in experiments [24].
Regarding the magnitude of the change in the wave velocity c with dynamical amplitude, it is generally larger for low static pressure than for high static pressure [44,30], in agreement with our results in Figs. 2 and 4. In Ref. [57], this behavior was suggested to result from significant rearrangements of the contact network, resulting in a change in the average contact number but without significant motion of particles or a significant change in the packing density. Indeed, as we observe here, small shifts in the positions of individual grains can modify the local contact number enough to change the effective velocity of elastic waves in the system. We find a reduction in the average contact number when the amplitude of the dynamical forcing is increased, consistent with the experimental results. This resembles the "acoustic fluidization" effect (initially introduced for describing frictional weakening [75]) that has been observed in which the elastic wave velocity can soften under large wave amplitudes even when significant contact reorganization and sliding do not occur [4,76]. If the wave amplitude were large enough to drop the average contact number below the jamming threshold in a significant portion of the sample, a "sonic vacuum" state could occur in which transmission of elastic waves would become impossible [77]. Above these amplitudes, the entire packing fluidizes, as in Refs. [78,79].
Granular packings often exhibit heterogeneous responses due to their highly disordered internal contact structure. In a 2D idealized granular packing, based on the response of a single grain to a sinusoidal driving frequency, localized normal modes at high frequency were predicted to occur, likely due to the interference between scattered plane waves [80]. Evidence for localized soft spots has been observed in Hertzian packings where the velocity distribution functions for the motion of individual particles have fat tails, indicating strongly non-Gaussian behavior [81]. These soft spots found at relatively low frequencies have been connected particularly with highly nonlinear responses such as glass-like behavior and non-affine displacement fields in granular packings [15,16,17,18,82,83,84].
Numerous studies have employed granular packings as a surrogate for the complex behavior occurring along fault zones in Earth [32,30,85,31,86]. Our results may suggest that the decrease in velocity observed along and near a fault after a large earthquake is analogous to a change in the granular packing to a state with a reduced number of contacts, even if no significant rearrangements of the grain positions have occurred. These contacts could gradually reconnect over time, in analogy with the slow recovery of the sound velocity that has been observed in the earth. Indeed, in the earth, the fault blocks surrounding a fault zone contain fractures at many scales. These are analogous to the grain contacts in our simulation and laboratory experiments. As wave amplitudes increase, slip is mobilized along the fractures resulting in a bulk modulus softening of the rock. This behavior is followed by slow dynamics where contacts in fractures are re-established, as demonstrated in laboratory experiments [30,4].

V.  SUMMARY

We characterize the evolution of the internal characteristics of bidisperse two-dimensional granular packings under large amplitude dynamic forcing and varied confining pressure. We find that the resonant frequency or fundamental mode of the frequency decreases with increasing dynamic amplitude at constant static pressure, in agreement with laboratory and field experiments. For fixed dynamic amplitude, the resonant frequency increases with increasing confining pressure, also in agreement with experiment. We show that the average contact number Zc of the packing decreases both at resonance and for increasing dynamic amplitude. We characterize the heterogeneity of the packing response by measuring the vibration displacement of each grain, and find regions of high and low displacements. Our approach provides insight into the elastic nonlinear nature of unconsolidated materials such as granular packings, as well as consolidated materials such as sandstone.
This work was supported by Institutional Support (LDRD) at Los Alamos National Laboratory. This work was carried out under the auspices of the NNSA of the U.S. DoE at LANL under Contract No. DE-AC52-06NA25396.

References

[1]
A.J. Liu and S.R. Nagel, Nature 396, 21 (1998).
[2]
M.E. Cates, J.P. Wittmer, J.-P. Bouchaud, and P. Claudin, Phys. Rev. Lett. 81, 1841 (1998).
[3]
A. Ostrovsky and P.A. Johnson, Riv. Nuovo Cimiento 24, 1 (2001).
[4]
X. Jia, Th. Brunet, and J. Laurent, Phys. Rev. E 84, 020301(R) (2011).
[5]
C. Reichhardt and C.J. Olson Reichhardt, Soft Matter 10, 2932 (2014).
[6]
D. Howell, R.P. Behringer, and C. Veje, Phys. Rev. Lett. 82, 5241 (1999).
[7]
D.M. Mueth, G.F. Debregeas, G.S. Karczmar, P.J. Eng, S.R. Nagel, and H.M. Nagel, Nature 406, 385 (2000).
[8]
O. Dauchot, G. Marty, and G. Biroli, Phys. Rev. Lett. 95, 265701 (2005).
[9]
P. Olsson and S. Teitel, Phys. Rev. Lett. 99, 178001 (2007).
[10]
A. Tordesillas, Phil. Mag. 87, 4987 (2007).
[11]
C. Heussinger and J.-L. Barrat, Phys. Rev. Lett. 102, 218303 (2009).
[12]
P. Olsson, Phys. Rev. E 81, 040301(R) (2010).
[13]
D. Bi, J. Zhang, B. Chakraborty, and R.P. Behringer, Nature 480, 355 (2011).
[14]
P. Olsson and S. Teitel, Phys. Rev. E 83, 030302(R) (2011).
[15]
L.E. Silbert, A.J. Liu, and S.R. Nagel, Phys. Rev. Lett. 95, 098301 (2005).
[16]
M. Wyart, L.E. Silbert, S.R. Nagel, and T.A. Witten, Phys. Rev. E 72, 051306 (2005).
[17]
L.E. Silbert, A.J. Liu, and S.R. Nagel, Phys. Rev. E 79, 021308 (2009).
[18]
N. Xu, V. Vitelli, A.J. Liu, and S.R. Nagel, EPL 90, 56001 (2010).
[19]
J.T. Jenkins and M.W. Richman, Phys. Fluids 28, 3485 (1985).
[20]
Y.M. Bashir and J.D. Goddard, J. Rheol. 35, 849 (1991).
[21]
S. Luding, M. Nicolas, and O. Pouliquen, in Compaction of Soils, Granulates and Powders, edited by D. Kolymbas and W. Fellin (Balkema, Rotterdam, 2000), p. 241.
[22]
M. Nicolas, P. Duru, and O. Pouliquen, Eur. Phys. J. E 3, 309 (2000).
[23]
J. Tejchman and E. Bauer, Granular Matter 5, 201 (2004).
[24]
J. Zhang, T.S. Majmudar, A. Tordesillas, and R.P. Behringer, Granular Matter 12, 159 (2010).
[25]
V. Magnanimo and S. Luding, Granular Matter 13, 225 (2011).
[26]
J. Ren, J.A. Dijksman, and R.P. Behringer, Phys. Rev. Lett. 110, 018302 (2013).
[27]
L.R. Aarons, J. Sun, and S. Sundaresan, Ind. Eng. Chem. Res. 49, 5153 (2010).
[28]
E.H. Field, Y. Zeng, P.A. Johnson, and I.A. Beresnev, J. Geophys. Res. 103, 26869 (1998).
[29]
I.A. Beresnev, G.M. Atkinson, P.A. Johnson, and E.H. Field, Bull. Seism. Soc. Am. 88, 1402 (1998).
[30]
P.A. Johnson and X. Jia, Nature 437, 871 (2005).
[31]
P.A. Johnson, H. Savage, M. Knuth, J. Gomberg, and C. Marone, Nature 451, 57 (2008).
[32]
S.L. Karner and C. Marone, J. Geophys. Res. 106, 19319 (2001).
[33]
J. Gomberg and P. Johnson, Nature 437, 7060 (2005).
[34]
F. Brenguier, M. Campillo, C. Hadziioannou, N.M. Shapiro, R.M. Nadeau, E. Larose, Science, 321 1478 (2008).
[35]
A.A. Delorey, K. Chao, K. Obara, and P.A. Johnson, unpublished.
[36]
D.P. Schaff and G.C. Beroza, J. Geophys. Res. 109, B10302 (2004).
[37]
J.L. Rubinstein, N. Uchida, and G.C. Beroza, J. Geophys. Res. 112, B05315 (2007).
[38]
J.L. Rubinstein, Bull. Seism. Soc. Am. 101, 275 (2011).
[39]
J.A. TenCate, D. Pasqualini, S. Habib, K. Heitmann, D. Higdon, and P.A. Johnson, Phys. Rev. Lett. 93, 065501 (2004).
[40]
R.A. Guyer and P.A. Johnson, Physics Today 52(4), 30 (1999).
[41]
E. Smith and J.A. TenCate, Geophys. Res. Lett. 27, 1985 (2000).
[42]
P. Johnson and A. Sutin, J. Acoust. Soc. Am. 117, 124 (2005).
[43]
J.A. TenCate, Pure Appl. Geophys. 168, 2211 (2011).
[44]
E.I. Mashinskii, J. Geophys. Eng. 1, 295 (2004).
[45]
M.W. Knuth, H.J. Tobin, and C. Marone, Granular Matter 15, 499 (2013).
[46]
P.A. Johnson, B. Zinszner, and P.N.J. Rasolofosaon, J. Geophys. Res. 101, 11553 (1996).
[47]
J.D. Goddard, Proc. R. Soc. Lond. A 430, 105 (1990).
[48]
X. Jia, C. Caroli, and B. Velicky, Phys. Rev. Lett. 82, 1863 (1999).
[49]
H.A. Makse, N. Gland, D.L. Johnson, and L. Schwartz, Phys. Rev. E 70, 061302 (2004).
[50]
M.A. Zimmer, M. Prasad, G. Mavko, and A. Nur, Geophys. 72, E1 (2007).
[51]
C.-h. Liu and S.R. Nagel, Phys. Rev. Lett. 68, 2301 (1992).
[52]
C.-h. Liu and S.R. Nagel, Phys. Rev. B 48, 15646 (1993).
[53]
J. Duffy and R.D. Mindlin, J. Appl. Mech. 24, 585 (1957).
[54]
P.J. Digby, J. Appl. Mech. 48, 803 (1981).
[55]
H.A. Makse, N. Gland, D.L. Johnson, and L.M. Schwartz, Phys. Rev. Lett. 83, 5070 (1999).
[56]
I. Agnolin and J.-N. Roux, Phys. Rev. E 76, 061304 (2007).
[57]
S. van den Wildenberg, M. van Hecke, and X. Jia, EPL 101, 14004 (2013).
[58]
B. Ferdowsi, M. Griffa, R.A. Guyer, P.A. Johnson, C. Marone, and J. Carmeliet, Geophys. Res. Lett. 40, 4194 (2013).
[59]
B. Ferdowsi, M. Griffa, R.A. Guyer, P.A. Johnson, C. Marone, and J. Carmeliet, Phys. Rev. E 89, 042204 (2014).
[60]
E. Somfai, J.-N. Roux, J.H. Snoeijer, M. van Hecke, and W. van Saarloos, Phys. Rev. E 72, 021301 (2005), and references therein.
[61]
Y. Khidas and X. Jia, Phys. Rev. E 81, 021303 (2010).
[62]
K.L. Johnson, Contact mechanics. (Cambridge University Press, Cambridge, 1985).
[63]
P.A. Cundall and O.D.L. Strack, Geotechnique 29, 47 (1979).
[64]
J.A.C. Gallas, H.J. Herrmann, and S. Sokolowski, Phys. Rev. Lett. 69, 1371 (1992).
[65]
J. Laurent, Ph.D. thesis, Université Paris-Est Marne-La-Vallée, July 2011, in French.
[66]
B. Gilles and C. Coste, in Powders and Grains 2001, edited by Y. Kishino (Balkema, Lisse, 2001), p. 113.
[67]
X. Jia and P. Mills, in Powders and Grains 2001, edited by Y. Kishino (Balkema, Lisse, 2001), p. 105.
[68]
B. Velicky and C. Caroli, Phys. Rev. E 65, 021307 (2002).
[69]
B. Gilles and C. Coste, Phys. Rev. Lett. 90, 174302 (2003).
[70]
C. Coste and B. Gilles, Phys. Rev. E 77, 021302 (2008).
[71]
T. Brunet, X. Jia, and P.A. Johnson, Geophys. Res. Lett. 35, L19308 (2008).
[72]
S.R. Pride and J.G. Berryman, Acta Mech. 205, 185 (2009).
[73]
S. Henkes and B. Chakraborty, Phys. Rev. Lett. 95, 198002 (2005).
[74]
T.S. Majmudar, M. Sperl, S. Luding, and R.P. Behringer, Phys. Rev. Lett. 98, 058001 (2007).
[75]
H.J. Melosh, Nature (London) 379, 601 (1996).
[76]
D. Espindola, B. Galaz, and F. Melo, Phys. Rev. Lett. 109, 158301 (2012).
[77]
L.R. Gomez, A.M. Turner, and V. Vitelli, Phys. Rev. E 86, 041302 (2012).
[78]
S. Luding, H.J. Herrmann, and A. Blumen, Phys. Rev. E 50, 3100 (1994).
[79]
S. Luding, Phys. Rev. E 52, 4442 (1995).
[80]
M. Leibig, Phys. Rev. E 49, 1647 (1994).
[81]
E.T. Owens and K.E. Daniels, Soft Matter 9, 1214 (2013).
[82]
M. Tsamados, A. Tanguy, C. Goldenberg, and J.-L. Barrat, Phys. Rev. E 80, 026112 (2009).
[83]
M.L. Manning and A.J. Liu, Phys. Rev. Lett. 107, 108302 (2011).
[84]
K. Chen, M.L. Manning, P.J. Yunker, W.G. Ellenbroek, Z. Zhang, A.J. Liu, and A.G. Yodh, Phys. Rev. Lett. 107, 108301 (2011).
[85]
K.E. Daniels and N.W. Hayman, J. Geophys. Res. 113, B11411 (2008).
[86]
P.A. Johnson, B. Carpenter, M. Knuth, B.M. Kaproth, P.-Y. Le Bas, E.G. Daub, and C. Marone, J. Geophys. Res. 117, B04310 (2012).


File translated from TEX by TTHgold, version 4.00.
 Back to Home