Dynamic induced softening in frictional granular material investigated by DEM simulation

A granular system composed of frictional glass beads is simulated using the Discrete Element Method. The inter-grain forces are based on the Hertz contact law in the normal direction with frictional tangential force. The damping due to collision is also accounted for. Systems are loaded at various stresses and their quasi-static elastic moduli are characterized. Each system is subjected to an extensive dynamic testing protocol by measuring the resonant response to a broad range of AC drive amplitudes and frequencies via a set of diagnostic strains. The system, linear at small AC drive amplitudes has resonance frequencies that shift downward (i.e., modulus softening) with increased AC drive amplitude. Detailed testing shows that the slipping contact ratio does not contribute significantly to this dynamic modulus softening, but the coordination number is strongly correlated to this reduction. This suggests that the softening arises from the extended structural change via break and remake of contacts during the rearrangement of bead positions driven by the AC amplitude.

💡 Research Summary

This paper presents a comprehensive Discrete Element Method (DEM) investigation of dynamic softening in a three‑dimensional granular assembly composed of frictional glass beads. The authors first construct a realistic particle model using the Hertz normal contact law together with a Mindlin‑type tangential force that incorporates Coulomb friction (μₛ = 0.22). Viscous damping is introduced both at the contact level (γₙ, γₜ) and globally (γₐ = 10⁻⁷ kg s⁻¹) to stabilize the numerical integration. The particle size distribution is bimodal (300 µm and 400 µm radii, 60 %/40 % by mass) and the material properties correspond to typical soda‑lime glass (Y = 65 GPa, ν = 0.25, G = 43.3 GPa, loss factor β = 0.0163). Approximately 4100 particles are confined in a rectangular cell (10 cm × 25 cm × 10 cm) with periodic boundaries in the horizontal directions and rigid walls in the vertical direction.

Quasi‑static characterization
The packing is first compacted by moving the bottom wall at high speed, then a static vertical stress σ is imposed (10 kPa – 1.44 MPa). Small cyclic strains (Δε ≈ 10⁻⁶) are applied in uniaxial compression and simple shear to extract the bulk modulus K and shear modulus G. The results show K ∝ σ^0.379, close to the Hertzian prediction K ∝ σ^1/3, while G ∝ σ^0.233, a somewhat weaker scaling than reported in earlier frictionless or low‑friction simulations. The coordination number Z_c varies from 4.38 to 4.95 and the packing fraction from 0.60 to 0.62, consistent with isotropic packings of frictional spheres.

Dynamic testing protocol
Dynamic response is probed by imposing a sinusoidal vertical displacement on the top wall: A_drive sin(ωt). For each driving frequency (6 × 10³ rad s⁻¹ – 6 × 10⁵ rad s⁻¹, step ≈ 1.9 × 10² rad s⁻¹) the system is allowed to evolve for 100 cycles; the last 60 cycles are averaged to obtain a steady‑state strain field in a probe layer located away from the walls (0.2 l_y – 0.4 l_y). The drive amplitude A_drive is varied from 5 × 10⁻⁸ m to 6 × 10⁻⁵ m, corresponding to drive strains ε_drive = A_drive/h_ly ranging from 4 × 10⁻⁶ to 5 × 10⁻³. Local strain ε_local is extracted by subtracting a moving‑average baseline and fitting the residual oscillation to a sinusoid, yielding amplitude and phase for each ω.

Results – resonance and softening
At the smallest ε_drive the resonance frequencies are constant and higher‑order harmonics appear, reproducing the linear regime observed experimentally in granular media. As ε_drive exceeds a threshold (different for each confining stress), the primary resonance peaks shift to lower frequencies and broaden, indicating a reduction of the effective elastic modulus (dynamic softening). Two low‑order modes are identifiable at low σ (≈ 10 kPa) and high σ (≈ 695 kPa); at high σ the second mode merges with the first as the drive amplitude grows, leaving a single dominant resonance.

Mechanistic insight – role of contacts
To pinpoint the origin of softening, the authors monitor two micro‑structural indicators in the probe layer: (i) the slipping‑contact ratio (SCR), i.e., the fraction of contacts where the tangential force reaches the Coulomb limit, and (ii) the average coordination number Z_c. SCR remains essentially unchanged across the full range of drive amplitudes, demonstrating that increased sliding does not drive the observed modulus reduction. In contrast, Z_c systematically decreases with increasing ε_drive, and the magnitude of the frequency shift correlates almost linearly with the loss of contacts. This strongly suggests that the softening originates from a restructuring of the contact network—contacts break and reform as particles undergo small, sub‑particle‑size rearrangements induced by the acoustic drive.

Discussion and implications
The findings challenge the conventional effective medium theory (EMT) assumption of a static, homogeneous contact network. Even in the presence of friction, the dominant non‑linear effect is not frictional dissipation but the dynamic evolution of the fabric. The observed scaling of K and G with σ, together with the coordination‑number‑driven softening, aligns with recent theoretical work on marginally jammed systems where the contact network controls elastic response near unjamming. The work also bridges observations from laboratory acoustic experiments (e.g., resonance softening, harmonic generation) with particle‑scale mechanisms, providing a clear pathway to interpret non‑linear acoustic signatures in geophysical contexts such as fault gouge or engineered granular dampers.

Conclusions
Through systematic DEM simulations, the paper demonstrates that (1) granular packings of frictional glass beads exhibit a clear dynamic softening of their elastic moduli when driven at sufficiently large acoustic amplitudes; (2) this softening is not caused by an increase in sliding contacts but is tightly linked to a reduction in the average coordination number, i.e., a restructuring of the contact network; (3) the effect is more pronounced at higher confining stresses where the contact network is initially more robust. These insights advance our understanding of non‑linear wave propagation in granular media and highlight the importance of micro‑structural evolution in governing macroscopic dynamic properties.