## Abstract

We investigate shock propagation in confined, frictionless granular media using discrete element simulations with an elastoplastic contact law. Depending on the level of confinement and loading, elastoplastic systems exhibit a weak or strong shock propagation response similar to an elastic Hertzian system although the details of the shock development differ markedly from the elastic case. Two modes of dynamic stress propagation are observed based on the shock intensity regime: weak shocks carry the stresses via the initial contact path while strong shocks form new contact networks behind the front. However, unlike for elastic shock propagation, there is an upper bound to the front velocity of strong shocks that depends on the maximum intergranular contact stiffness. Since elastoplastic contact is a dissipative process, results show that dissipation is enhanced with confining pressure in the weak shock regime.

## 1 Introduction

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*_{0}, 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.

## 2 Contact Law and Numerical Methodology

*α*= 1, 2) [24]. We assume that all spheres have the same constitutive properties: stiffness

*E*, Poisson’s ratio

*ν*, and yield stress

*σ*

_{y}. The loading force–displacement (

*F*

_{l}–

*δ*) relation adopted in this work has a smooth transition from an elastic Hertzian domain to a fully plastic flow domain, while the unloading force–displacement (

*F*

_{u}-

*δ*) relation is elastic:

*a*)

*b*)

*a*is the contact area,

*δ*

_{R}is the residual plastic displacement, $R\xaf$ is the deformed relative radius, and subscripts

*y*and

*p*, respectively, denote the value of the quantity of interest at the onset of yield and plastic flow. These quantities are calculated using the following equations:

*c*)

*E** and

*R** are the effective Young’s modulus and effective radius, respectively, and are given by

*E*= 115 GPa, Poisson’s ratio

*ν*= 0.3, and yield stress

*σ*

_{y}= 619 MPa. The Brake model’s specific parameters, the Meyer’s exponent

*n*, and contact pressure

*P*

_{0}are, to our knowledge, not tabulated for brass 260 and were extracted as

*n*= 2.042 and

*P*

_{0}= 1.017 GPa by curve fitting to the experimental data shown in Fig. 1 [23]. As reflected in that figure, the force–displacement relationship under loading is almost linear for higher contact forces. In this study, inter-granule contact is assumed to be frictionless and thus only involves normal interactions. The dissipated plastic energy associated with a given contact is the area under the loading and unloading curve and is given by

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^{5} 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.

^{−7}s was chosen to ensure stability and convergence of the numerical results. It is worth noting that this timestep value is an order of magnitude smaller than the transit time of a sound wave along the length of the granular domain. The first step in the simulation process consists of jamming the granules to generate a prescribed pressure in the granular ensemble. To that effect, granules are first randomly initialized with an Hertzian contact law such that the area fraction,

*A*

_{f}, of the total number of granules is close to the jamming transition point of 0.84 for 2D packings (Fig. 2(a)) [30]. Given the channel width and granular diameters, this packing involves 1375 granules. Next, the walls are moved inwards from this initial state to compress the granules and exceed the critical area fraction of 0.84, which causes the materials to jam and build-up an internal stress state,

*σ*

_{ij}, defined by

*F*

_{x}and

*F*

_{y}are the amplitudes of the horizontal and vertical force components associated with the pre-impact jamming phase,

*l*and

*w*, respectively, denote the length and width of the channel, and

*D*is the diameter of the larger granule (0.04 m). Artificial viscosity is applied at this stage to remove excess kinetic energy of the granules until the total kinetic energy equilibrates to 10

^{−8}J. The resulting stress state is effectively hydrostatic with a stress ratio between the vertical and horizontal directions between 0.95 and 0.98 due to the contacts being frictionless and the granules essentially behaving like a “granular liquid” [31]. The order in which the walls are moved does not affect the final confining pressure which is area fraction dependent [9] but may affect the underlying fabric of the force chain network. The confining pressure,

*P*

_{0}, can be computed as

*P*

_{0}= (1/2)|

*σ*

_{xx}+

*σ*

_{yy}|. In the next step, the contact law in the equilibrated system is switched from the elastic Hertzian contact law to the elastoplastic contact law. The change of contact law slightly perturbs the previously stabilized granular state, hence requiring the system to be re-equilibrated under artificial viscosity until the kinetic energy stabilizes back to approximately 10

^{−8}J. This state is considered as the starting jammed state as a strong force chain network forms to support and maintain the pressure when the artificial viscosity is removed (Fig. 2(b)).

The switch to the elastoplastic contact law inevitably yields contacts that carry loads greater than the yield force (Eq. (1*a*) 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.

Type of interaction | Yield force (N) |
---|---|

Small–small | 32.2 |

Small–big | 57.3 |

Big–big | 128.8 |

Small–wall | 128.8 |

Big–wall | 515.1 |

Type of interaction | Yield force (N) |
---|---|

Small–small | 32.2 |

Small–big | 57.3 |

Big–big | 128.8 |

Small–wall | 128.8 |

Big–wall | 515.1 |

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.

## 3 Results and Discussion

### 3.1 Structure of the Shock.

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.

### 3.2 Dependence on Confining Pressure and Piston Velocity.

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\u223cP01/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].

*u*

_{p}. This result is likely associated with the presence of a precursor traveling ahead of the front and creating a dynamic compression comparable to that observed for lower piston velocities. There also appears to be an upper bound for the front velocity as increasing the piston velocity more than two-fold from 70 m/s to 150 m/s increases the front speed by only 7%. This upper bound is due to the stiffness of the quasi-linear contact law for large intergranular forces making the granular fabric appear as if it were connected by quasi-linear springs. For 1D chains, the upper bound of the front velocity,

*V*

_{max,1D}, can be obtained from the long wavelength limit of the dispersion relation given by

*K*

_{max}is the maximum intergranular stiffness,

*D*is the granular diameter, and

*ρ*

_{avg}is the average density of the granules [26]. For 2D random packings, the maximal stiffness of a contact,

*c*, for a given granule,

*p*, must be averaged over all the specific angular orientations in the contact network. This averaging is performed using the granular stiffness tensor formulated by Liao and Chan [32]:

*V*is the volume encapsulating the granules,

*N*

_{c}is the number of contacts for granules

*p*,

*K*is the intergranular contact stiffness,

**is the normal vector connecting two contacts, and**

*n**l*is the spacing between two neighboring granules. Assuming macroscopic elasticity, a long wavelength approximation for the maximum front speed can be found using

*K*

_{max}, and

*ρ*

_{avg}is the average density of the granular assembly contained in

*V*. It should be noted that Eq. (4) can be derived from Eq. (6) by assuming a unit cell containing one contact and setting

**to be [1, 0]. As evident in Fig. 8, this upper bound agrees well with the numerical results for a piston velocity of 150 m/s. The practicality of using an elastoplastic contact law to describe wave propagation for impacts beyond 150 m/s is not feasible as the granules lose their spherical shape due to very large loads (>300 kN) for which the overlap between granules exceeds one-quarter of the granule diameter. Equation (6) can also be used to estimate the sound speed (**

*n**c*

_{0,predicted}) of a packing by replacing the maximal stiffness with the initial contact stiffness as a result of confinement [33]. Although the dependence of

*c*

_{0,predicted}on confining pressure (dashed curve in Fig. 8(a)) is similar to that obtained numerically, the analytical prediction overestimates the numerical results. This discrepancy has been observed elsewhere [13] when homogenization was performed to predict the sound speed in granular systems and is attributed to the localized modes associated with the initial force chain network, which cannot be captured in a long wavelength type analysis.

### 3.3 Energy Dissipation Characteristics.

*E*

_{p,tot}, and the work done by the piston,

*W*. Δ

*E*

_{p,tot}is computed by summing the change between the initial and final plastic energy dissipated (Eq. (2)) over all contacts. The work done by the piston at time

*t*is evaluated using

*F*

_{p}is the force exerted on the piston. Figure 9 shows the dependence of the Δ

*E*

_{p,tot}/

*W*ratio on the piston velocity for different values of the confining pressure. We observe that dissipation increases with confining pressure for weak shocks, but this trend gradually diminishes as the piston velocity increases and eventually disappears after the formation of strong shocks. Unlike the front velocity results, which showed a pressure dependence even in the strong shock regime, a clear transition is observed in Fig. 9 when

*u*

_{p}exceeds approximately 4 m/s. For reference, this transition speed is remarkably lower than the uniaxial shock speed in continuum brass, which is approximately 440 m/s [34]. This key difference is associated with the fact that the continuum must not only yield but also be in the fully plastic region, and this process requires very strong impacts.

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*_{0} = 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.

## 4 Conclusions

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.

## Acknowledgment

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

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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.