Shockwaves in Jammed Ductile Granular Media

In granular materials, the contact interactions between granules govern the mechanical response of the entire ensemble. The classical Hertzian contact theory, which describes elastic interactions between spherical particles, has paved the way to understand the unique wave propagation properties of granular media. Studies on 1D granular chains have shown that the inherent nonlinearity of the Hertzian force–displacement contact relationship makes uncompressed granular media behave as a sonic vacuum without a characteristic sound speed [1–3]. As a result, depending on the material properties and applied loading conditions, nonlinear wave propagation in granular media takes the form of singular solitary waves, solitary wave trains, or shockwaves.

Due to the nonlinear contact interactions among its constituents, wave propagation in 1D granular media can be tuned by simply confining the granular system through a lateral pressure or precompression. An individual contact without precompression following the Hertzian force–displacement relationship has zero stiffness at the origin, thereby making it impossible to propagate sound waves. However, by applying an initial confining pressure, P₀, to the 1D chain, the contact becomes linearized making it possible for small disturbances to propagate with a sound speed that scales as $P01/6$ [4,5]. As the amplitude of the disturbance is increased, wave propagation becomes increasingly nonlinear and eventually gives rise to strongly nonlinear solitary waves and strong shocks, in which the velocity of the shock front, v_s, scales with the maximum dynamic pressure, P, as P^1/6 [1,3,5].

In 2D and 3D packings of granular media, a further complexity arises where the application of an external pressure causes the granules to jam creating dense, heterogeneous force chains that support the externally applied loading [6–9]. Studies on elastic systems have shown duality in the scaling of the sound speed with the confining pressure. Under different experimental conditions, higher confining pressures show a $P01/6$ scaling of the sound speed similar to 1D Hertzian chains, whereas low confining pressures show a sound speed scaling of approximately $P01/4$ [10–14]. The crossover between these two regimes is configuration specific and is speculated to be a result of nonlinear effects such as dynamic changes in the coordination number, mixing of shear/compression wave modes, and contact asperities [11,13,15]. Although these results show a universal behavior irrespective of the underlying force chains, experiments probing even smaller, more localized length scales of the order of the force chains reveal that force chains act as conduits for stress propagation [16,17]. On the opposite end of the spectrum, strongly nonlinear wave propagation can occur for large impulses that produce dynamic forces much larger than the initial confining forces [18–20]. In those studies, the resulting shock front velocity was found to scale as P^1/6 similar to that observed in 1D chains.

A limitation of the Hertzian contact law for ductile granules is that it is only applicable for small amplitude loading conditions provided the contact surface has not been preconditioned (pre-yielding of the contact surface) [21]. If we consider non-preconditioned ductile granular spheres (granules in its most common form), the Hertzian contact law greatly overestimates the force–displacement relation beyond a relatively small critical displacement value as the contact region yields at low loads due to stress concentrations [22–24]. Over the years, several contact force–displacement laws have been proposed to capture elastoplastic contact for different material properties and have been used in discrete element method (DEM) simulations to study the dynamics of elastoplastic wave propagation [25,26]. While solitary waves in elastic systems propagate unattenuated, elastoplastic waves show severe attenuation due to plasticity [23,25,27]. Unlike their elastic counterparts, shockwaves in 1D uncompressed elastoplastic systems show fronts with smaller oscillations and have an upper bound for the wave velocity that depends on the maximum stiffness of the contact [25,26,28].

In this work, we use DEM simulations to study elastoplastic shockwave propagation in confined, random 2D granular packings made of two spherical grain sizes and make comparisons with corresponding elastic Hertzian granular systems. The primary goal of our simulations is to investigate the effect of confining pressure and impact velocity on plastic dissipation and front velocity. We first study the phenomenology of the shock by observing the shape of the front and force chain network. Next, we extract the front velocity from the simulations and predict the lower and upper bounds for 2D randomly packed granular systems by extending the approximations of these bounds for 1D chains. Finally, we show the dependence of plastic dissipation on the confining pressure and excitation velocity by investigating the dynamic pressure and changes to the contact network.

It should be noted that E_p increases monotonically with δ_max.

In our simulations, we study wave propagation along the length of a granular channel containing a bidisperse distribution of brass 260 granules with a 2:1 ratio of 0.02 m and 0.04 m diameters. These granules are confined in a 0.4 m × 4 m rectangular channel having smooth walls, which are approximated as granules with very large (10⁵ m) diameter, making these “wall granules” appear as flat elastoplastic surfaces while being compatible with the simulation framework detailed below. The material properties of the wall are assumed to be the same as those of the granules but yield a much stiffer contact due to the larger effective radius. The size of the simulation channel was chosen to capture the steady-state propagation of the front and doubling the width of the channel did not bring about a noticeable difference in the results.

The switch to the elastoplastic contact law inevitably yields contacts that carry loads greater than the yield force (Eq. (1a) with δ = δ_y). Contact interactions involving smaller granules show more yielding as they have a lower yield force due to a small effective radius. Analysis of the plastic contacts prior to the impact event indicates that all the contact forces lie within the loading curve of the contact law irrespective of yielding. As shown in Fig. 3, increasing the pressure results in a progressive increase in the fraction of yielded contacts with nearly 5% of the contacts being yielded for a 25 kPa confining pressure and nearly 90% by 350 kPa. For completeness, we list in Table 1 the yield force values for the different contacts involved in the simulated granular medium.

The final step in our simulation procedure involves moving the left wall inwards like a piston with constant velocity, u_p, to generate a shock. The other three walls are stationary and do not interact with the piston or one another. The simulation is terminated before the shock reaches the opposite end of the channel to prevent reflection.

The rightward motion of the piston generates a pressure front that travels from left to right into the granular medium, with the granules located behind the front moving at approximately the velocity of the piston. Visualizing the propagation of the front by zooming into a portion of the channel after jamming but before impact (Fig. 4(a)) and for piston velocities of 0.8 m/s (Fig. 4(b)) and 70 m/s (Fig. 4(c)) reveals how the impact amplitude, i.e., the piston velocity, controls the dynamic creation of the force network. In Fig. 4(a), we observe that a significant force chain network develops in the jammed material with almost every contact point being active. In Fig. 4(b), this pre-existing contact network carries the dynamic loading with force magnitudes higher than the static pre-load as they provide a transfer path for the contact forces. As a result, the void locations and shapes also remain mostly intact. In contrast, Fig. 4(c) shows significant formation of new contacts accompanied by a decrease in void area fraction behind the propagating front. Since both Fig. 4(b) and Fig. 4(c) are taken at the same time instant, the front speed appears to be an increasing function of the piston speed.

Since the granules are discrete, uniquely identifying the wave front and extracting the wave speed are challenging. Adopting the approach of Gómez et al. [19], we obtain the front profile by dividing the length of the channel into 80 uniformly spaced bins and averaging the horizontal velocity, v_x, of the granules in each bin. A spatially varying velocity profile is then obtained by fitting a cubic spline through the bin-averaged granule velocity data. Figure 5 shows the evolution of the front at similar spatial positions for elastic and elastoplastic shocks for piston velocities of 0.8 m/s (top figures) and 70 m/s (bottom figures) and for a confining pressure of 212 kPa. For the lower value of u_p, we observe that the elastic front (Fig. 5(a)) is sharper and less smooth than the elastoplastic front (Fig. 5(b)). In addition, in the elastic case (Fig. 5(a)), the horizontal velocity of the granules reaches a plateau, albeit noisy, at about the imposed piston velocity, u_p. In the plastic case, at least for this length of channel, no plateau in granule velocity is achieved as a result of plastic dissipation.

Comparing the strong and weak shock profiles, we observe that the strong elastic (Fig. 5(c)) and elastoplastic (Fig. 5(d)) fronts are sharper than their weak counterparts with the elastoplastic front being smoother and increasingly wider. Here the shock width refers to the length over which the profile extends from the point it first goes above 0 m/s to u_p. This result implies that the “base” of the front travels faster similar to an elastic precursor (annotated in Fig. 5(d)) in homogeneous material. In uncompressed 1D elastoplastic chains, the spreading of a shock front is not observed even though loading is dominated primarily by the linear part of the force–displacement curve, which is responsible for wave dispersion [26]. This 1D result is due to the fact that loading at any contact must first “travel” up the nonlinear Hertzian part of the contact law before reaching the linear plastic contact domain. However, the formation of a contact network in the jammed state initializes contacts to have a linearized stiffness depending on the amount of load they carry. The presence of a confining pressure therefore tends to accentuate features of dispersion, including the presence of a precursor, since it eliminates the intrinsically nonlinear response associated with the unconfined state.

The observation of a precursor is also apparent by extracting the front velocity at different positions along the front profile. The position-dependent front velocity can be quantified by tracking the position, x, and time, t, at which the velocity of each bin exceeds a certain threshold. Figure 6 presents the x–t curves obtained for five velocity threshold values ranging from 0.005 to 35 m/s for the shock displayed in Fig. 5(d). The linear nature of these curves points to a constant front velocity associated with each velocity threshold, with the lower-threshold velocities yielding higher front velocities confirming the existence of an elastic precursor.

By extracting the front velocity using the values from the lowest threshold in Fig. 6, we first verify our system by comparing the results of the elastic (Hertzian) system with theoretically expected scalings. The weak and strong shock propagation regimes can readily be observed in the dependence of the shock front velocity on the confining pressure obtained for different piston velocities presented in Fig. 7(a). For weak shocks obtained at low piston velocities, the front velocity scales with the confining pressure as $vs∼P01/6$⁠, as shown by the solid line in Fig. 7(a). As mentioned earlier, this front velocity is effectively the speed of sound of the confined granular system for very small piston velocities.

As the piston velocity is increased, a transition to the strong shock regime is visible when the front velocity loses dependency on the confining pressure. Figure 7(b) shows another representation of the data in Fig. 7(a) by viewing the front velocity as a function of piston velocity u_p. That figure indicates that this regime transition takes place for u_p > 4 m/s and that the front velocity emerges to scale with the theoretically expected scaling of $up1/5$ [19].

To shed further insight on the dissipation characteristics in the two shock regimes, we characterize the pressure behind the front by computing the traction acting on the piston. As shown in Fig. 10, the dynamic pressure is comparable to the confining pressure in the weak shock regime. As expected, the trend of dynamic pressure with increasing piston velocity mimics the independence of the dissipated plastic energy to piston work ratio on the confining pressure observed in Fig. 9, and all the data points coalesce in the strong shock regime when the dynamic pressure is an order of magnitude higher than the confining pressure.

The increase in dynamic pressure is partly due to the increase in the intergranular contact forces and partly due to the formation of new contacts, N_nc. Figure 11(a), which tracks the evolution of N_nc for different piston velocity values given the same initial confining pressure (P₀ = 102 kPa), shows that N_nc is effectively 0 in the weak shock regime and that there is a steady increase in N_nc with time for strong shocks. This result implies that in the weak shock regime, the dynamic stresses are carried by the initial force chain network as seen in Fig. 4. Figure 11(b), which presents the evolution of N_nc normalized by the length traveled by the front, shows that, for the range of confining pressures in this study, the formation of new contacts is dependent only on the piston velocity.

Therefore, the increase in dissipation with confining pressure in the weak shock regime is due to individual contacts within the initial force chain network being capable of dissipating more energy with increasing pre-load. This is enabled by the fact that a strongly pre-loaded contact dissipates more energy for an incremental increase in force due to the monotonically increasing relationship between plastic dissipation and force. However, in the strong shock regime, the dynamic forces far supersede the confining forces, thus making the effect of confinement negligible.

In this paper, we have numerically studied elastoplastic shock propagation in frictionless, 2D jammed granular media with the aid of DEM simulations based on an elastoplastic contact law between granules. Similar to their elastic counterparts, elastoplastic shock propagation shows weak and strong shock regimes with the transition between them dependent on the piston velocity and the level of confinement. In the weak shock regime, the shock propagation speed shows a dependance on the initial confining pressure, whereas the shock speed is mostly affected by the piston velocity in the strong shock regime.

Based on the granular properties and conditions considered in this study, strong shocks can be created in granular systems with substantially (about two orders of magnitude) smaller piston velocities than in a homogeneous medium made of the same material. In the weak shock regime, dynamic stresses are carried by force chains formed due to confinement, while a strong shock is characterized by a considerable rearrangement of granules and the formation of new contacts. When compared with shockwaves in elastic granular media, elastoplastic shockwaves have smoother, less sharp profiles, with an elastic precursor traveling ahead of the front. Elastoplastic shockwaves also propagate considerably slower than elastic shockwaves and the shock velocity is bounded by two linearized approximations of the stiffness of the contact network. The lower bound (effectively the sound speed of the system) and upper bound are functions of the volume averaged intergranular contact stiffness due to confinement and maximum contact stiffness, respectively.

Energy dissipation is enhanced with confinement in the weak shock regime as strongly pre-loaded contacts dissipate more energy. This finding can extend the already known capability of a granular medium to dissipate more energy per volume than a continuum medium made out of its constituent material [28]. However, in the strong shock regime, the dynamic stress behind the shock front supersede that associated with the initial force chain network, making dissipation independent of the confining pressure.

This material is based upon work supported by the National Science Foundation under grant number NSF CMMI 17-61243.

There are no conflicts of interest.

The data and information that support the findings of this article are freely available. The authors attest that all data for this study are included in the paper.