Efficacy of simple continuum models for diverse granular intrusions

Abstract

Granular intrusion is commonly observed in natural and human-made settings. Unlike typical solids and fluids, granular media can simultaneously display fluid-like and solid-like characteristics in a variety of intrusion scenarios. This multi-phase behavior increases the difficulty of accurately modeling these and other yielding (or flowable) materials. Micro-scale modeling methods, such as DEM (Discrete Element Method), capture this behavior by modeling the media at the grain scale, but there is often interest in the macro-scale characterizations of such systems. We examine the efficacy of a macro-scale continuum approach in modeling and understanding the physics of various macroscopic phenomena in a variety of granular intrusion cases using two basic frictional yielding constitutive models. We compare predicted granular force response and material flow to experimental data in four quasi-2D intrusion cases: (1) depth-dependent force response in horizontal submerged-intruder motion; (2) separation-dependent drag variation in parallel-plate vertical-intrusion; (3) initial-density-dependent drag fluctuations in free surface plowing, and (4) flow zone development during vertical plate intrusions in under-compacted granular media. Our continuum modeling approach captures the flow process and drag forces while providing key meso- and macro-scopic insights. The modeling results are then compared to experimental data. Our study highlights how continuum modeling approaches provide an alternative for efficient modeling as well as a conceptual understanding of various granular intrusion phenomena.

---

1 Introduction

Granular intrusion is common in various natural and human-made situations.^1–3 In many cases, the flow response to intrusion leads to a multi-phase behavior of the granular media, where the media shifts among solid-like, fluid-like, and gas-like states.^4,5 This behavior depends on material properties and local variables such as pressure and granular packing fraction.⁶ The space and time dependence of these fields, as well as the multi-phase behavior of the media, makes granular flow response challenging and computationally expensive to model.

Discrete Element Methods (DEM) are reliable and established methods for modeling granular media.^7,8 In DEM, the granular medium is modeled at the particle level, with the frictional, elastic contact forces between adjacent and interacting particles obtained using standard contact models like Hertzian and Hookean contacts,^9,10 generally with an assumption of no plastic yielding of the individual grains. The elastic and external forces are integrated via Newton's laws of motion to calculate particle velocities and positions. In the past few decades, DEM has played a major role in the advancement of the field of granular media.^11,12 However, the precision of DEM comes at a computational cost. As an example, a 0.1 mm particle diameter granular system with a 0.5 m per side cubic domain would require evaluating ∼9 × 10¹¹ DOF per time step (at a packing fraction of 0.64 with 6 degrees of freedom, DOF per particle). While current large-scale granular DEM studies remain in the range of ∼10⁷ particles (∼10⁸ DOF),^13,14 efforts are being made to extend this range. Most recently, some have reached the ∼10⁹ DOF range with the use of supercomputers and GPUs.¹⁵

In recent years, there has been a growing interest in modeling increasingly large granular systems, such as those relevant in ballistic impacts, wheeled locomotion, and agricultural plowing.^1–3,16,17 This has led to a push for computationally faster alternatives to DEM, such as using continuum descriptions of granular media and further reduced models such as Resistive Force Theory.^18–20 At the fundamental level, cohesionless granular media can be described as a plastic material with a frictional yield criterion and no cohesion strength (hence no ability to support tension). But, a comprehensive continuum model for fully describing a granular material requires incorporating a variety of behaviors shown by the media arising from the constituent grain shape, size, stiffness, distribution, and inter-particle interaction properties,⁵ as well as effects due to the surrounding air (pressure and flow-drag) on the granular media.²¹ Thus, often scenario-specific simplifications are made during constitutive model development. Conversely, even for a simple granular material, a constitutive framework encapsulating all granular phenomena would have to merge many different pieces together, e.g. dissipative kinetic theory^22,23 in dilute regimes; strengthening based on packing fraction²⁴ and internal state variables (like fabric and hysteresis) evolution;²⁵ rate-dependence of the strength (such as μ(I) rheology²⁶); non-trivial yield surface shapes in stress space;²⁷ and particle size effects^28,29 in the dense flow regime.

A variety of computational methods have been developed and utilized recently to aid in simulating continuum models for large deformation processes, such as the Material Point Method (MPM),^18,30–32 the Particle Finite Element Method (PFEM),³³ and Smoothed Particle Hydrodynamics (SPH).³⁴ In this regard, a recent study on high-speed locomotion of wheels²⁰ in grains, a process which involves complex trans-phase characteristics of the soil and non-trivial grain motions, revealed that a continuum treatment implemented with MPM captures the essential behaviors and agrees well with experimental data. In the present work, we highlight the capability of the continuum approach in modeling and explaining diverse granular intrusion phenomena using a very minimal set of constitutive ingredients. The phenomena we focus on herein have been independently studied in the literature, and our goal is to demonstrate that they can be quantitatively modeled and unified under a family of basic constitutive assumptions. We first introduce two constitutive representations of non-cohesive granular media: a non-dilatant plasticity model (NDPM), and a dilatant plasticity model (DPM). The DPM model permits dilatancy during dense flow while the NDPM model assumes dense flow is a constant-volume process. Thus, the NDPM is suitable for modeling steady-state granular behaviors, while the DPM is a more suitable model for transient flow processes. The latter is also a more elaborate model that converges to the former model under certain limits. Both of these models have their respective advantages depending on the scenarios they simulate. We do not include micro-inertial μ(I) effects²⁵ (effects of grain-level inertia on material properties) in these models; thus both of these constitutive models are rate-insensitive.

We demonstrate the utility of each model by considering granular flow and force response in four fundamental intrusion cases which have been studied in the literature: (1) depth-dependent force response in horizontal submerged intruder motion; (2) separation-dependent drag variation in parallel plate vertical intrusion; (3) initial density-dependent drag fluctuations in free surface plowing; and (4) flow zone development during vertical plate intrusions in under-compacted granular media (Fig. 1). We use NDPM in the first two cases (which focus on a steady response), and DPM in the latter two cases (where transient effects are the focus). Our continuum modeling approach captures the dynamics of such granular flows and gives a deeper macroscopic insight into each of these cases. Thus, our study highlights the efficacy of continuum modeling for predictive purposes, as well as in developing macroscopic conceptual understanding of diverse granular intrusion phenomena.

Fig. 1 Schematics of four test cases considered in this study: (A) horizontal dragging of a submerged cylinder at different depths, (B) multiple-plate vertical granular intrusion, (C) free-surface plowing in under- and over-compacted granular media, and (D) vertical intrusion in under- and over-compacted granular media. Cases (A) and (B) are be modeled with the non-dilatant plasticity model (NDPM), and cases (C) and (D) are modeled with the dilatant plasticity model (DPM) to model transient effects.

2 Material and methodology

We use two constitutive models to represent granular volumes in this study, namely: a non-dilatant plasticity model (NDPM), and a dilatant plasticity model (DPM). We first explain the two models in detail and then explain the numerical method we use to implement the schemes.

Non-dilatant plasticity model (NDPM)

This model is taken from Dunatunga and Kamrin.^30,35 The model assumes a rate-insensitive, Drucker–Prager yield criterion, incompressible plastic shear flow (with no dilatancy), and cohesionless response in extension whereby the material becomes stress-free below a ‘critical density’, ρ_d. The constitutive flow equations representing the above behaviors are shown below. These simultaneous constraints describe the material's separation behavior, shear yield condition, and tensorial co-directionality, respectively:

(ρ − ρ_d)P = 0 and P ≥ 0 and ρ ≤ ρ_d, (1)

[small gamma, Greek, dot above](τ − μ_sP) = 0 and [small gamma, Greek, dot above] ≥ 0 and τ ≤ μ_sP, (2)

for i, j = 1, 2, 3, where, ρ is the local granular density, ρ_d is the critical density of close-packed grains, μ_s is the internal friction coefficient, P = −σ_ii/3 is the local hydrostatic pressure, is the equivalent shear rate. is the equivalent shear stress, is the deviatoric part of Cauchy stress tensor, D_ij = (∂_iv_j + ∂_jv_i)/2 is the flow rate tensor.

For simplicity, the equations above are expressed in rigid-plastic form, where the flow rate is approximately all due to plastic flow. However, we also include a small elastic component which ensures that below the yield criterion the grains have a well-defined solid constitutive behavior. The model evolves the flow by solving the momentum balance equations, ∂_jσ_ij + ρg_i = ρv̇_i. We also assume a constant friction coefficient μ_f between the granular continuum and solid-body surfaces. The internal friction value μ_s is often measured from the tangent of the angle of repose for the granular media or from a direct shear test.³⁶ We refer to this model as the ‘Non-Dilatant Plasticity Model’ (NDPM). We provide the input material properties used in various cases in the relevant sections. A basic implementation of this continuum modeling approach can be downloaded for Matlab.³⁷

Dilatant plasticity model (DPM)

While NDPM assumes a single critical granular density cutoff for the entire granular volume in the system, this is often not the case in real granular systems. Granular media can support stress in a range of densities including close and loose packings.⁵ The second model we present addresses this limitation of the NDPM model and allows for density evolution in the media.⁶ This model takes inspiration from the family of critical state models.^6,24 We use a simple rigid-particle version of the Critical State model in which grains are not modeled to crush under hydrostatic loads; rather, they only have a hydrostatic elastic response in this case. When pressure is non-zero, the material is deemed to be in a ‘dense state’ and the density evolution occurs through shearing-based Reynolds dilation/contraction,³⁸ representing particle rearrangements at the granular level. In a dense state, the packing fraction ϕ is equal to an evolving state variable ϕ_d. Every shearing action in dense media drives the local value of ϕ_d towards a limiting steady-state critical packing fraction, ϕ_c. Furthermore, when the packing fraction varies, the material's friction coefficient also changes; material with low packing fraction has low consolidation and thus has less strength than more consolidated media.

Thus, in terms of constitutive relations, the simultaneous constraints describing the material's separation behavior and shear yield condition (eqn (1) and (2), respectively) remain the same with the two exceptions that the close-packed density, ρ_d, and internal friction value, μ_s, are non-constant local state variables (which vary with time and space), as follows:

ρ_d = ρ_gϕ_d, μ_s = μ_c + (ϕ_d − ϕ_c)χ (4)

Here, ρ_g represents the solid grain density, χ represents a dimensionless scaling constant, and μ_c represents the critical-state internal friction value reached by material at steady state (when ϕ_d = ϕ_c). When pressure is positive (the packing is dense), the evolution of ϕ_d is modeled as:where, β = (1/3)(ϕ_d − ϕ_c)χ is a local dilatancy variable. In the absence of confining pressure, the material can leave the dense state (ϕ < ϕ_d), but reconsolidates at the current value of ϕ_d unless the material opens up beyond a global lower limit, ϕ^min_d, which defines a minimum density possible for loaded media, below which no connected material states exist. We set this value to 0.45 in all simulations. Thus, the material can exist in a dense state only for ϕ ≥ 0.45 and the local state variable ϕ_d gets reset to ϕ^min_d(= 0.45) whenever ϕ drops below this value.

In addition to the material's separation behavior and shear yield condition (eqn (1) and (2)), the material flow rule has a deviatoric dependence as before (eqn (3)) but due to dilation now also has a corresponding volumetric component. The volumetric strain-rate is given by the dilatancy variable, β. The plastic flow rate tensor, thus, has the form:

Thus, the combination of eqn (1), (2), (4)–(6) collectively represent the Dilatant Plasticity Model (DPM). Note that DPM reduces to NDPM in the limit of χ → 0. Regarding its usage in cases 3 and 4, χ is calibrated for each material qualitatively based on experimental observations.

While DPM is a more elaborate representation of granular media compared to NDPM, NDPM is easier to implement and is computationally more efficient to solve. In this study, we aim to highlight the advantages and sufficiency of using each of these constitutive forms on a case-by-case basis.

For the numerical implementation, we use a continuum simulation approach based on the Material Point Method (MPM),³⁹ a derivative of the fluid-implicit-particle (FLIP) method for implementing a continuum description of granular media. The MPM implementation uses a combination of Lagrangian material points, which contain the state of the continuum, and a background computational mesh that is used to solve the equations of motion. Since the state is saved on the material points, the mesh is reset at the beginning of each computational step without any information loss, and thus the method allows for large deformations in the media without the error associated with large mesh deformations. A schematic representation of an explicit time-integration MPM-step is shown in Fig. 2. We use the MPM implementation described in Dunatunga and Kamrin³⁰ to implement the different constitutive equations representing granular volumes assuming 2D plane-strain motion. We choose the spatial grid resolution (Δx) in each system on a case-by-case basis, such that all the geometric features are well represented and all system outputs converge. Also, to ensure numerical stability, we choose a time step (Δt) that satisfies the CFL condition based on minimum element size, maximum Young's modulus, and minimum medium density.

Fig. 2 Sample explicit time integration step in MPM: a sample explicit time-integration step used in an MPM implementation. The Lagrangian tracers representing material points (solid red circles) carry the state of material over time and space. The background grid (green squares) assists in integrating the motion on the simulation domain and is reset each step. More details can be found in Dunatunga and Kamrin.³⁰

3 Results and discussion

3.1 Drag and lift on submerged cylinder dragging

Drag and lift forces on submerged objects in granular media are relevant in processes such as mixing, mining, soil-buried pipelines, and animal locomotion.^40–42 Over the past two decades, several studies^43–46 have explored the mechanisms and the variation of these forces by considering the horizontal dragging of submerged, rigid cylinders (at different depths) as a test case. We specifically consider the work of Guillard et al.,⁴⁴ who observed that the horizontal drag force on cylindrical objects moving horizontally increases with increasing depth while the vertical lift force plateaus in such motions for depths greater than a O(1) factor of the cylindrical diameter. The observations consider a depth range deep enough so that the cylinders are completely immersed in the granular bed even at the minimum depth. At very low depths near the free surface, the dragging motion of the cylinders results in an accumulation of the media in front of the intruders. This accumulation can augment the depth-dependence of drag forces and thus can change the force trends mentioned above. Thus, we do not consider such a near free-surface depth range in this study (as in Ding et al.⁴⁶) and focus on depths ≥2.5D.

The experimental drag and lift results in Guillard et al.⁴⁴ were obtained by evaluating torques and vertical lifts on horizontal cylinders inside large granular beds rotated about the vertical axis. After a half rotation or so about the vertical axis, the authors observed a reduction in the force response as the cylinder begins re-interacting with the disturbed material as it is rotated multiple times. Here, we use the quasi-steady data obtained from the first half rotation of the cylinder, before it interacts with the disturbed media. Presumably, density variations play a smaller role in this portion, and we choose to compare the behavior observed with that of the NDPM model.

A schematic representation of the cylinder drag case is provided in Fig. 1A. For the NDPM implementation, we use a grain density ρ_g of 2520 kg m⁻³, and a critical packing fraction of 0.60. We calibrate the internal friction value (μ_s) for glass beads to 0.48 to accurately model the initial slope of depth vs drag graph from the experiments. This is in line with the expected range of μ_s for the glass beads which lies in the 0.39–0.55 range (corresponding to an angle of repose range ∼22–29°, depending on the surface roughness).⁴⁷ The calibration absorbs possible inconsistencies between our 2D simulations with 3D experiments and also incorporates the variations due to indirect measurements of experimental drag and lift from rotating cylinder experiments.⁴⁴

The cylinder is modeled as an elastic body with a high elastic modulus, thereby acting as an approximate rigid body. The media-cylinder interface friction coefficient (μ_f) is set to 0.35. The cylinder diameter, D is 4 mm, and the presumed out-of-plane length, w is 1 m (as the simulations are 2D plane-strain). We use a 0.20 m × 0.16 m granular bed, and a 5 × 10⁻⁴ m spatial resolution (Δx) for simulating these cases.

Fig. 3A and B show the comparison of Guillard's experimental data to our continuum results. The linearly increasing nature of drag vs. depth graphs, along with the plateauing of lift vs. depth graphs are well captured. The visualizations of plastic strains and strain rates (Fig. 3C and D) qualitatively agree with those reported in previous numerical and experimental studies.^43,44 Notably, the variation of the equivalent plastic strain shows a strong semblance with the work of Wu et al.⁴⁵ Additionally, we also plot the variation of the local hydrostatic pressure at four different depths in Fig. 3E. The pressure distribution shows an asymmetry about the intruder, consistent with the existence of both drag and lift forces, and as the intruder's depth increases the vertical asymmetry diminishes, which is consistent with the plateauing of the lift force and continued growth of drag force.

Fig. 3 Case 1: experiments vs. simulations for drag and lift on submerged cylinders (using NDPM model): comparison of experimental data⁴⁴ (in blue) and the continuum simulation results (in red) for variation of (A) drag, and (B) lift on a submerged cylinder moving in granular media at depth |z|. The drag and lift forces are non-dimensionalized with characteristic force π(D/2)²ρ_dgw where D is the cylinder diameter, ρ_d is the effective medium density, g is the gravity, and w is out-of-plane cylinder width. Both experiments and continuum results use D = 4 mm, ρ_d = 1512 kg m⁻³, and g = 9.8 m s⁻². Similarly, the depth is non-dimensionalized with cylinder diameter, D. Variation of macroscopic state variables: (C) equivalent plastic strain, (D) equivalent plastic strain rate, and (E) local hydrostatic pressure at four different depths, |z| = [10,30,51,71] mm. The gray circle indicates the cylinder position and the drag direction is from left to right. See Movie S1 (ESI†) for visualizing material flow over time for cylinder drag cases considered in (C) and (D).

Similar to past DEM studies by Guillard et al.⁴⁴ and experimental studies by Wu et al.,⁴⁵Fig. 3C and D show high localization of the flow field for the depths greater than the lift turnover depth. A further increase in depth after this point results in a minimal-to-no change in the material flow profile. The drag and lift force behaviors are also related to this flow localization behavior. While the drag forces continue to increase linearly with depth, the lift forces become saturated near the depth the flow localizes, indicating a possible correlation between the flow and force. Fig. 3E indicates a reducing asymmetry of the pressure variation around the cylinder as depth increases, consistent with the notion that drag forces continue to increase but lift forces plateau. A deeper analysis of this variation (i.e. the distribution and the relative magnitude of pressure) to understand the mechanics of the force distributions on different faces of the intruders is reserved for future study.

Thus, the case of dragging submerged cylinders provides an example where NDPM captures the observed intrusion behavior.

3.2 Vertical drag in two-plate granular intrusions

Multi-body intrusions offer another commonly encountered scenario in real-life situations. Several researchers^48–50 have examined the dynamics of granular intrusions involving multiple intruding bodies. The case we study specifically takes inspiration from the work of Swapnil et al.⁵⁰ who examined the variation of vertical drag forces on a pair of parallel rigid plates during vertical downward intrusions as a function of the separation between them. They observed a peak in the average vertical drag (measured as work per unit depth) as the separation between the plates was increased. A similar ‘cooperative effect’ was observed by Merceron et al.⁴⁹ during the upward motion of parallel intruders in the granular media. When the separation between the intruders was below a critical value, jammed media between the intruders was observed. We use NDPM to investigate the phenomena observed by Swapnil et al.⁵⁰ and focus on the regime where intruders and separations are large enough that the role of particle size effects is minimal. We also conduct experiments to verify claims that follow from theoretical analysis of the NDPM model.

A schematic representation of the case is given in Fig. 1B. The case is studied in two parts. In the first part, we use generic material properties (provided next) to determine the behavior of the continuum model in the single-plate intrusion case. And in the second part, a direct quantitative comparison of simulation results with the experiments is done in the two-plate intrusion case (details later). We use a set of generic NDPM fitting properties from the earlier case, i.e. a grain density, ρ_g of 2520 kg m⁻³, a critical packing fraction of 0.6, but choose the internal friction coefficient to be μ_s = 0.4, which is more common for glass beads.⁴⁷ The plates were modeled as stiff elastic bodies with vertical displacement of control points assigned; thus, the plates act as quasi-rigid objects with a common fixed downward velocity. The intrusion velocity was set to 0.1 m s⁻¹ in all simulations, and the media-plate surface friction was set to 0.35. We use a 1.2 m × 0.6 m granular bed, and a 1.25 × 10⁻³ m spatial resolution (Δx) for simulating these cases.

Understanding single plate intrusions. We begin by first analyzing vertical drag forces in single plate intrusions. Fig. 4A shows the drag forces on vertical plate intrusions under NDPM for various plate widths at a constant intrusion velocity of 0.1 m s⁻¹. Note that the same model and implementation was used previously by Dunatunga and Kamrin³⁵ for modeling high-speed impacts of a circular intruder in granular media, where it was shown to match the flow and force data of Clark and Behringer.⁵¹ Before analyzing the simulation results, we perform a scaling analysis of the problem assuming the NDPM equations hold, in order to predict the dependence of drag on various system parameters according to NDPM. In this case, we expect the resulting vertical drag force, F^D on a vertical intruder of width L and out-of-plane width w at a depth of z to depend on the various system parameters, i.e. intruder dimensions (w and L), depth z (distance between the intruder bottom and the free-surface), close-packed media density ρ_d, gravity g, the media's internal friction μ_s, the media–intruder interface friction μ_f, and the velocity of intrusion v. Using scaling analysis, with base units of length as L, time as , and mass as ρ_dL³, we obtain:where, f̂ represents some unknown function. Recall that we are scaling based on NDPM, so certain properties such as the particle diameter do not appear in the above relation.

Fig. 4 Case 2: vertical intrusion of single plates (using NDPM model): (A) variation of drag forces with intrusion depth for plates of various widths, L [0.05, 0.10, 0.15, 0.20, 0.25] m during single-plate vertical granular intrusion. The arrow shows the direction of increasing L. The simulations are plane-strain and glass bead properties are used for simulating granular media (ρ_g = 2520 kg m⁻³, ϕ = 0.6, μ_s = 0.4). Each data series on the left graph was time-averaged over a 0.5 ms window to remove high-frequency force fluctuations. The vertical black dotted line in (A) shows the depth after which this averaging window includes sufficient data. Dependence of (B) initial force peaks, and (C) average drag forces on plate widths in single plate intrusions (blue squares). The corresponding quadratic and linear fits for (B) and (C) are shown as red dotted lines. The force peaks are shown with arrows in the left graph, and the time window used for calculating the average forces is the highlighted blue region (force averaging window) on the left graph.

The dependence of F^D on z can be divided into two regimes. (1) At low depths z ≪ L, the variable z/L is negligible and can be ignored. And (2), at larger depths, the drag forces F^D are known to show a linear dependence on depth (after an initial jump in the vertical drag near free surfaces^52–55). In both of these regimes, we set the dependence of F^D on w to be linear assuming the plane-strain nature of the intrusions. We also focus on intrusions that are sufficiently slow such that F^D is independent of velocity v (low-velocity regimes) as seen in the work of Swapnil et al.⁵⁰

With the above assumptions, in the z ≪ L regime, we obtain:

And in the latter, moderate-depth regime with F^D ∝ z, we obtain:Thus, from the scaling analysis, we expect the drag force per unit out-of-plane width F (= F^D/w) (unit Nm⁻¹) to have a quadratic dependence on plate width, L, near the free surface, and a linear dependence on plate width at larger depths. The scaling analysis does not provide information on the exact form of the variation of F^D in z in the region connecting the two regimes. However, we do expect the variation to be non-monotonic in z since in the vanishing depth limit the force scales as L² and in the deeper limit it scales with zL.

Using the simulation's force output, in Fig. 4B and C we plot the observed relationships between F and L at, respectively, a smaller depth (where initial peaks occur) and in a deeper regime, where linear force variations are predicted. The forces in Fig. 4C are averaged over a depth window (z_i − z_f), i.e.. The expected trends from the dimensional analysis are apparent. Besides our own simulation results, the linear dependence of drag force on intruder area, once deep enough, is a well studied relation.^52–54 We reiterate that our simulations are all in the quasistatic regime. Thus, the velocity contributions are negligible in comparison to static force contributions in our study. Thus, the initial peak in the force response is not related to inertial drag as seen in faster intrusions.^53,56

Physically, the initial force maxima (F^peak) in the force vs. displacement graphs of Fig. 4A correspond to the force requirements for initiating media flow in the system. The F^peak force contributions enable the beginning of flow by initiating the shearing of the media, which exists at some finite strength due to finite pressure under the plate. From slip-line theory for granular intrusions,⁵⁷ the flow develops a wedge-shaped no-shear-zone below the intruder, which has an acute angle of π/2 − tan⁻¹(μ) (where μ represents the material internal friction). Thus, the requirement of shearing media of a finite strength along the wedge-shaped shear zone is responsible for these force contributions. We can also obtain an intuition for the initial L² dependence of force by considering the conventional limit analysis for indentation of materials with yield stress Y, common in manufacturing processes such as drawing and blanking.⁵⁸ In such cases, indenter force per unit out-of-plane width, F, varies as F ∝ Y × L for L the indenter plate length. In frictional media, the strength is pressure-sensitive. With Y = μ × P and supposing P grows linearly along the edge of the no-shear-zone wedge, the mean strength along the edge grows ∝ L. Substitution of Y ∝ L in the limit analysis formula then gives the observed quadratic dependence of the drag force on L. We do not further investigate these forces but identify that the existence of small regions of under-compacted granular media near free surfaces is expected to suppress the growth of such forces in many cases (see Fig. 8D). Our simulations indicate that the drag forces’ dependence of F on L enter linear regimes at depths z ∼ O(10⁻¹)L. These forces are unique from the ‘added mass’ effects⁵⁹ and other macro-inertial effects^52,53 in granular impacts that are common in high-speed intrusion (and vary ∝ v²) since our intrusion velocities are small.

Two plate intrusions. The above scaling analysis gives important insights for the case of multibody intrusions. In view of eqn (9), we expect that two plates intruding far from each other and at moderate depth, will experience the same net vertical, drag F^D, as that experienced by a single plate with an equivalent surface area—the linearity in L (eqn (9)) indicates that this equality should hold for any ratio of areas between the two plates. Thus, any combination of plate widths should experience the same force both when they are infinitely far apart and when they have no separation. This analysis does not provide any information on the force variations when the plates are close but are at a finite distance from each other. However, physical intuition suggests that the presence of plates in the vicinity of each other will restrict the material flow, which should increase the drag on each plate. Given that the plates experience equal net forces at infinite and zero separation between them, we expect there to be a value of separation at which the force response is maximal or minimal (unless the force response is constant). The work of Swapnil et al.⁵⁰ explored this variation and found there is a peak in the force response at a low value of separation between the plates. Note that near the free surface (z ≪ L), we do expect drag force at zero separation to be higher than at infinite separation due to the quadratic dependence of the force on plate size in this regime—the ratio of forces for plate lengths of L₁ and L₂ would be F₀/F_∞ = (L₁ + L₂)²/(L₁² + L₂²) which is always greater than 1.

Fig. 5A and D show the variation of force for different combinations of plate widths and plate separations. Continuum modeling shows the existence of force peaks for both equal plate cases (Fig. 5A) and unequal plate cases (Fig. 5D). Our observations are in accord with similar experiments and DEM simulations by Swapnil et al.'s⁵⁰ for equal plates. As the continuum modeling successfully captures the behavior, the detailed material states in these simulations can help identify the macro-mechanical origins of the phenomena. We visualize the material flow by plotting snapshots of the plastic strain in a few of these cases. The plastic strain fields before, at, and after the force peak in Fig. 5B and E, suggest a macro-mechanical picture. We observe higher granular flow interaction between the two plates as the separation between the plates is decreased. For a single plate intruder, any neighboring flow restriction is expected to make it more difficult to push material during the intrusion. Thus, decreasing plate separation results in increasing drag on each plate. Once the plates are sufficiently close, a large wedge-shaped rigid zone forms spanning the plates, causing the two plates to act as a large single plate. Thus, any further reduction in the separation does not result in additional flow restriction. Instead, it leads to a reduction in the effective area of the merged plate systems, and thereby the drag forces decrease upon a further decrease in separation. It is interesting to note that the material flow profiles when the two plates are together or far-separated are similar to the classical plasticity solution of Prandtl⁶⁰ for yielding of metals upon indentation under a plate, characterized by a single rigid wedge of media under the indenter and flow emanating from both edges of the wedge. However, when the plates are separated but very close, the flow profile looks similar to the indentation plasticity solution of Hill,⁶¹ characterized by two smaller wedges under the plate and flow emanating only from the two outermost edges. We also plot the variation of local hydrostatic pressure in the media for the two cases in Fig. 5C and F. The ratio of local pressure magnitudes with equivalent hydrostatic pressure, ρ_dgz, is large (∼40) in accord with the observations of Brezinski et al.⁵⁴ To understand the behavior of the flow solutions, we re-derive the scaling relations for combinations of plates of different lengths L₁ and L₂ (L₁ ≤ L₂) assuming NDPM, similar to the single plate case. Using base units of length as L_avg = 1/2(L₁ + L₂), time as , and mass as ρ_dL_avg³, we obtain:

Fig. 5 Case 2: vertical intrusion of two parallel plates (using NDPM model): (A) variation of vertical drag (normalized) with plate separation, s (normalised to plate width L) for equally sized plates (L = 0.10 m each). (B) The corresponding equivalent plastic strain fields and (C) local hydrostatic pressure fields at different separations. (D) Variation of vertical drag with plate separation, s (normalised with average plate width L_avg = (L₁ + L₂)/2) for unequally sized plates  (L₁ = 0.05 m and L₂ = 0.15 m). (E) The corresponding equivalent plastic strain fields and (F) local hydrostatic pressure fields at different separations. The forces F^D in (A) and (D) are averaged over a depth range (z_i − z_f) of 0.04–0.08 m for all the cases and an average depth, z̄ = 0.06 m, w = 1 m, g = 9.8 m s⁻², ρ_d = 1512 kg m⁻³ (= ρ_g × ϕ) is used. The simulations are plane-strain and glass bead properties are used for simulating granular media (ρ_g = 2520 kg m⁻³, ϕ = 0.6, μ_s = 0.4). See Movie S2 (ESI†) for visualizing material flow over time for the two cases considered in (B) and (E).

With similar assumptions as before in the moderate-depth regime (F ∝ z),

F^D = ρ_dgzL_avgwf̂(μ_s,μ_f,s/L_avg,L₁/L₂).

Defining a non-dimensional force, F̃ = F^D/ρ_dgzL_avgw, we get

F̃ = f̂(μ_s,μ_f,s/L_avg,L₁/L₂) (10)

and for equal plates, L₁ = L₂, L_avg = L and we obtain:

F̃ = f̂(μ_s,μ_f,s/L). (11)

Therefore, if the NDPM model suffices to describe the physics of two-plate intrusion, we expect the above relation to describe a master curve that collapses data for plate-pairs of various L intruding into the same media.

Based on eqn (10), for F ∝ z, we obtain following relations for peak separation (s_p) and corresponding peak force (F_p) values:

s_p = L_avgf̂₁(L₁/L₂,μ_s,μ_f), (12)

and

F_p = ρ_dgzL_avgf̂₂(L₁/L₂,μ_s,μ_f). (13)

For equal plates (L₁ = L₂ = L = L_avg) we obtain

s_p = Lf̂₁(μ_s,μ_f), (14)

and

F_p = ρ_dgzwLf̂₂(μ_s,μ_f). (15)

Note that F_p has units of force per out-of-plane length. Also note that, under these relations, the only material properties that influence the peak force and separation are the friction coefficient(s) and density; the sole length scale comes from the plate itself. This makes physical sense when the smallest in-plane feature (min(w,L,s)) is sufficiently large compared to the grain diameter. Swapnil et al.⁵⁰ explored the dependence of force-peak separation s_p/d with intruder size, L/d, when the intruder size is close to the grain diameter, d (L/d range 1–6). Interestingly, their data agree with our proposed linear dependence s_p ∝ L at intruder sizes as low as 3–4 grain diameters.

We also verify these drag force variations and the peak separation scaling relation with new vertical intrusion experiments and compare them to calibrated continuum simulations (see Fig. 6 for the details). A DENSO VS087 robot arm intruded an apparatus that held pairs of steel plates at various separations into a bed of loosely consolidated Quikrete Pool Filter Sand with grain density ρ_g = 2520 kg m⁻³, and effective close-packed density ρ_d = 1512 kg m⁻³. The internal friction value of this medium was μ_s = 0.72 based on angle of repose measurements from experimental tilting tests, with the media globally fluidized for 15 seconds to an initial packing fraction of ϕ ≈ 0.58 for all trials. All intrusions were performed at 11 mm s⁻¹, where speed dependence of the force response is negligible. Over 128 trials, the net resistive forces on the pair of intruding plates were measured using an ATI Mini40 force transducer. We observed a satisfactory match between the experiments and the simulations. The fact that the experimental data for different values of L collapse onto a single dimensionless master curve supports the robustness of the scaling relation implied by NDPM in eqn (11). Additionally, we observe that although both the Fig. 5(A) and 6 show peaks in normalised force responses of the media at low separations, the shapes of the graphs are not ‘identical’. This variation is expected because eqn (11) indicates that the graph between normalized drag (F̃) and normalized separation (s/L) depends on material friction properties (μ_s and μ_f) and the two cases use different internal friction (μ_s) values.

Fig. 6 Case 2: experimental verification of peak force phenomenon in two-plate intrusions: the comparision of experimental data (dotted lines) and calibrated continuum simulations (solid line) for two (equal) plate intrusion experiments. The experiments use plates of width (L) 15 mm (blue data with ● marker), 20 mm (orange data with * marker), 26 mm (yellow data with ■ marker), 29 mm (violet data with Δ marker), and the continuum model uses a plate of width 10 mm (green data with ★ marker). Quikrete Pool Filter Sand with grain density ρ_g = 2520 kg m⁻³ was used (inset shows microscopic view). The effective critical density ρ_d and angle of repose of the sand was found to be 1512 kg m⁻³ and 36° ± 1 resp. All of the paired plates have a 1:5 horizontal aspect ratio. The continuum results correspond to plane-strain two-plate intrusions with drag forces F^D averaged over a depth range of 0.06–0.08 m (an average depth, z̄ = 0.07 m), with vertical axis normalized by w = 1 m, g = 9.8 m s⁻², ρ_d = 1512 kg m⁻³ (= ρ_g × ϕ). The material properties are calibrated based on experimental data with ρ_g = 2520 kg m⁻³, ϕ = 0.6, μ_s = 0.72 (= tan θ_repose).

We also briefly explore, in our simulations, the effect of changing plate width ratios on the peak separation distance (s_p), and the peak separation force (F_p) in our study (see Fig. 5D). We find that F_p monotonically decreases from a maximum value to F_∞ as the plate ratio decreases from 1 to 0 (note that plate-ratio L₁/L₂ is always ≤1 as L₁ ≤ L₂ by definition). The normalized peak separation distance (s_p/L_avg) increases with decreasing plate-ratios (L₁/L₂). The peak separation relations are also expected to be a function of μ_s from scaling analysis (eqn (13)). Thus, an in-depth shape and material property dependence characterization of this phenomenon is relegated to future study. Similarly, a more comprehensive amalgamation of the macroscopic insight we develop in this study with the microscopic observations made by Swapnil et al.⁴⁹ for grain-scale plate lengths can be attempted in the future studies.

3.3 Drag variations in the plowing of granular media

This case takes inspiration from the work of Gravish et al.⁶² who studied drag force fluctuations in the plowing of granular beds at different initial packing fractions. A Discrete Element Method (DEM) based study of the case was also performed by Kobayakawa et al.⁶³ Both of these studies observed increasing drag force fluctuations with an increasing initial packing fraction of the granular beds. Similar force fluctuations were observed by Kashizadeh and Hambleton⁶⁴ on reduced-order modeling of the plowing processes in sands and by Jin et al.⁶⁵ in the development of a new single-gravity (1-g) small scale testing methodology. As this phenomenon directly relates to the changing density of the media, we use the more detailed DPM model for this case.

A schematic representation of this case is given in Fig. 1C. Both of the reported studies were performed in 3D while our simulations are 2D plane-strain. Characterizing the effects of this difference in our studies is difficult, so we do not attempt an exact match. Nevertheless, we do expect the 2D simulations to capture the phenomena qualitatively and the drag forces to be similar in their magnitudes. We once again use the glass bead material properties used in the previous case. We use a grain density ρ_g of 2520 kg m⁻³, a critical packing fraction ϕ_c of 0.6, a steady-state critical internal friction μ_c of 0.4, and a scaling coefficient χ of 2.5. Similar to previous cases, the plates are modeled as elastic bodies with high elastic modulus to approximate rigid bodies. The media/plate interface friction (μ_f) was set to 0.35. We use a 2.4 m × 0.4 m granular bed, and a 4 × 10⁻³ m spatial resolution (Δx) for simulating these cases. The plowing plate dimensions are 0.03 × 0.08 m².

Fig. 7 shows the variation of horizontal drag forces from the Gravish et al.⁶² experiments alongside our continuum simulations. The continuum results are scaled proportionally to the out-of-plane width in Gravish et al.⁶² The mean drag force and force fluctuations from continuum modeling are plotted in Fig. 7D and E, showing smooth forces transition to larger, fluctuating forces as the initial packing fraction rises above ϕ_c. The same trends can be seen in experiments, cf.Fig. 7A and B. A visual free surface comparison between experiments and the simulations also shows the similarity of flows as ϕ_i varies, cf.Fig. 7C and F. In the over-compacted case, the continuum results show a stepped pattern forming on the surface but it is not as persistent as the wavy patterns observed in the experiments. This is due to the 2D nature of our simulations. The 2D nature restricts the material from flowing in an out-of-plane direction, which causes the material to flow over the existing waves and overrun the wavy patterns. Such patterns are otherwise not overrun in 3D experiments and are more visible in over-compacted cases.

Fig. 7 Case 3: force fluctuations during plowing (using DPM model): variation of drag forces at various initial packing fractions in Gravish et al.'s⁶² experiments: (A) time variation of drag forces, and (B) average drag forces. The average values of forces for three initial packing fractions in (A) are shown with corresponding colored arrows in (B). Corresponding continuum simulations: (D) time variation of drag forces (force plots for consecutive ϕ_i are shifted vertically by 0, 40, 80, and 120 N respectively for improved visualization), and (E) averaged drag forces. Visualisation of the free surface in under/over compacted granular media: top view from Gravish et al.'s⁶² experiments (C), and side view in continuum simulations (F). (G) Variation of material packing fraction from continuum simulations in initially under-compacted (ϕ_i = 0.57) and over-compacted (ϕ_i = 0.63) media cases considered in (F). The simulations are plane-strain and glass bead properties are used for simulating granular media (ρ_g = 2520 kg m⁻³, ϕ_c = 0.6, μ_c = 0.4, χ = 5.0, ϕ_min = 0.45). See Movie S3 (ESI†) for visualizing time evolution of material front during plowing of under- and over-compacted media from experiments and simulations. Figure (A)–(C) are modified from Gravish et al.⁶² [Photo credits: (C) N. Gravish, P. B. Umbanhowar and D. I. Goldman, Georgia Institute of Technology].

The observed shear patterns, free-surface profiles, and force fluctuations shown in Fig. 7A–C are in accord with the shear-softening (shear-strengthening) of over-compacted (under-compacted) granular media that is built-in to the DPM model. In a sample initially under-compacted (ϕ_i < ϕ_c), shear deformations cause compaction, resulting in higher densities along the sheared regions than in the bulk (see eqn (5)). This density increase results in higher shear strength along the shear zone (see eqn (4)). Thus, further loading in such systems induces shear to occur in the weaker material adjacent to the shear band, which effectively spreads the shearing in such systems. On the other hand, in the over-compacted case, shear deformations dilate the material, resulting in lower density in the sheared region than in the bulk, which results in lower shear strength there. Hence, continued loading causes shear to accumulate along the thin zone of initial failure causing the appearance of a strong shear band. This process continues until the total force requirement for shearing along the existing band exceeds that for creating a new shear band in the bulk (after which the same process repeats itself). Thus, in over-compacted media, a visually separable shear band formation pattern occurs (Fig. 7F). In initially over-compacted cases (ϕ_i > ϕ_c), increasing the initial packing fraction of the media results in an increased plate motion requirement between successive shear band formations (due to increased strength of the media in the bulk) and thus the force fluctuation magnitudes increase (and their spatial frequency decreases) with increasing packing fraction (observed also by Gravish et al.⁶²). We also plot the variation of the packing fraction in the media from continuum modeling in Fig. 7G. The figure provides a visualization of changing packing fraction ahead of the plate in accordance with the smooth versus banded mechanism explained above.

3.4 Development of a shear deformation zone in plate intrusions

This last case takes inspiration from the work of Aguilar and Goldman⁵⁹ and highlights the capability of the basic DPM model in capturing the development of the flow profile during vertical intrusion of single plates. A schematic representation of the case is given in Fig. 1D. Aguilar and Goldman⁵⁹ postulated that the formation of a rigid No Shear Zone (NSZ) ahead of a flat plate during vertical intrusion is an incremental growth process. In three dimensions, this growth takes the form of a rigid frustum shape which grows from a low height frustum with a fixed base (matching the shape of the intruder) to a fully developed cone/pyramid at the completion of the mechanism. For a circular base, the final shape is a cone. In two dimensions, like our case, this would translate to successive isosceles trapezoids (with larger parallel edges matching the intruder's intruding edge) leading to a wedge shape with the base as the intruder's leading edge.

Fig. 8 shows the variation of intrusion forces and successive velocity profiles in under-compacted granular intrusion experiments compared to our continuum simulations. We model the granular media, poppy seeds, with the DPM to incorporate the effect of density transitions. We use material properties for poppy seeds with a grain density ρ_g = 1100 kg m⁻³, a critical packing fraction ϕ_c = 0.60, a steady-state critical internal friction μ_c = 0.53, and a scaling coefficient χ = 5.0. The media/plate interface friction (μ_f) was set to 0.35. We use a 0.5 m × 0.2 m granular bed, and a 5 × 10⁻⁴ m spatial resolution (Δx) for simulating these cases. The intruding plate dimensions are 0.04 m × 0.02 m.

Fig. 8 Case 4: material front development during vertical intrusion (using DPM model): (A) experimental data from Aguilar and Goldman.⁵⁹ The experiments intruded a circular plate of 0.051 m diameter into poppy seeds (ρ_g = 1100 kg m⁻³) at various initial packing fractions (ϕ_i). (B) Schematic of Aguilar and Goldman's cone growth model.⁵⁹ (C) Continuum modelling force data from 2D plane strain simulations. We intrude 0.04 m wide flat plates into material with properties ρ_g = 1100 kg m⁻³, μ_c = 0.53, ϕ_c = 0.60, and χ = 5.0 as calibrated to poppy seeds. (D) Mechanism of solid zone development apparent from continuum modeling. (E) PIV analysis of constant speed plate intrusions: new set of experiments were conducted by intruding L = 40 mm wide rectangular flat plates next to a clear plexiglass wall at a constant low speed, v = 150 mm s⁻¹, in loosely packed poppy seeds. The initial three frames represent the cone growth phase and the last frame (t̃ ∼ 0.30) shows the stability of the cone at later stages of the intrusion. Time is non-dimensionalized (t̃ = t/t₀) by t₀ ≡ L/v. The velocity is normalized by intrusion speed, v, and strain rates by [small gamma, Greek, dot above]₀ ≡ 1/t₀. We plot these results in the intruder's frame of reference. (F) Continuum simulation results: vertical velocity and equivalent plastic strain rate fields during vertical intrusion in initially under-compacted media (ϕ_i = 0.55) for material properties discussed in (C). We also show the evolution of packing fraction (ϕ) from continuum modeling in (G). See Movie S4 (ESI†) for visualizing material flow over time for various initial packing fractions (ϕ_i) including the one considered in (F).

The force trends from simulation qualitatively match the trends from Aguilar's experiments (compare Fig. 8A–C), keeping in mind that Aguilar uses a 3D circular plate while our simulations are in 2D plane-strain. In a set of separate experiments using a rectangular plate intruder and PIV (Fig. 8E) we observed the flow zone development appears similar to continuum results (Fig. 8F). These trends also agree with general observations from a 3D DEM study done by Feng et al.⁶⁶ The Aguilar and Goldman⁵⁹ study presents a model whereby the rigid cone emerges from a growing rigid frustum of constant base angle that gets progressively longer until converging to the final cone shape (Fig. 8B). However, in our simulations, the zone actually starts as a ‘short’ wedge coincident with the plate bottom, having a larger apex angle (expected by slip-line theory to be approximately π/2 − tan⁻¹ μ_i for μ_i the friction at the initial density) due to the low initial packing fraction of the media. The rigid front then grows by ‘fanning out’ from the diagonal edges (see Fig. 8D), as the edges represent the zone experiencing maximum shear-compaction and hence the most strengthening. This growth can be observed as well from density variations, shown in Fig. 8G. The growing density of the region results in an increasing internal friction value in the zone below the intruder and results in the development of a progressively sharpening, quasi-rigid trapezoid-like shape under the plate. As the intruder moves deeper, the process converges to a final wedge shape (with a sharper half-angle equal to ‘π/2 − tan⁻¹ μ_c’). Thus the DPM model provides an apt description of the observed behavior in under compacted granular intrusion.

4 Approach limitations and their implications

While the continuum modeling and the numerical implementation used in this study are able to represent the considered cases to a sufficient degree, both the model and the method have their limitations. Clearly, as emphasized in the introduction, the phenomena incorporated in a constitutive model limits the behavior the model can capture, and the constitutive relations we use here intentionally exclude certain effects for the benefit of simplicity. Similarly, MPM has known accuracy limitations given by the choice of the grid resolution, material point density, shape functions, and the means of representing contact between domains. For instance, MPM inherently captures a locally volume averaged material response everywhere. If during the flow an element has a low number of interior material points, the accuracy of the integration is also diminished. But these issues can be overcome with an appropriate choice of shape functions and refined discretization. Specifically for these issues, use of more advanced methods such as the hybrid DEM-MPM approach (such as Yue et al.⁶⁷ and Chen et al.⁶⁸) or dynamic particle enrichment (such as Zhu et al.⁶⁹) could be used at the expense of computation time. The use of smoother and wider shape functions could also help decrease numerical fluctuations often observed in MPM.⁷⁰

5 Conclusions

In this work we demonstrated the efficacy of continuum modeling in four intrusion cases using two continuum descriptions of granular media—(1) depth-dependent force response in horizontal submerged intruder motion; (2) separation-dependent drag variation in parallel plate vertical intrusion; (3) initial density-dependent drag fluctuations in free-surface plowing; and (4) flow zone development in vertical plate intrusions in under compacted granular media (see Fig. 1). The study shows that relatively simple, friction-based plasticity models capture a large variety of granular intrusion phenomena. Moreover, the models provide a useful macroscopic understanding of granular intrusion processes, which are often primary interests in engineering applications, and remove the additional complexity of trying to determine large-scale physics from grain-scale observations. The simplicity of these continuum models also streamlines this understanding, both by exclusion – i.e. if such a model works, it implies that mechanisms or effects lying outside the model's formulation are not crucial to the outcome – and by admitting simple scaling analyses as we have utilized throughout. Such simplifications will certainly limit the accuracy of these models in a variety of cases, but an incremental approach of adding physical augmentations (such as micro-inertial effects, particle size effects, or evolving fabric variables) provides a systematic approach for exploring the underlying physics in diverse cases. For instance, we do not use μ(I) rheology in either of the models in this study; the fact that our modeling still captures the observed behaviors indicates micro-inertial effects are not a key mechanism in the observed behaviors. We also emphasize that the continuum approach is useful for determining further-reduced models of intrusion, due to the simplicity of its large-scale characterization of systems. One such study was done in Agarwal et al.,²⁰ which develops a generalized, rate-dependent, dynamic resistive force theory (DRFT) for rapid intrusion of generic intruders. In the future, the work could be extended to three dimensions to extensively compare the computational advantage of using such methods. Furthermore, the continuum treatment could help reconcile granular behavior with similar behaviors observed in other, more standard continua. For example, Minetti et al.⁷¹ reported that during swimming in water, slight separation between fingers increases propulsive thrust, similar to our observation for slightly separated granular intruders.^71,72 Comparing and analyzing continuum forms could provide insights into the rationale behind such similarities, as was done in other flow resistance studies.⁷³

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

SA and KK acknowledge support from Army Research Office (ARO) grants W911NF1510196, W911NF1810118, and W911NF1910431. AK and DG acknowledge support from ARO grant W911NF-18-1-0120. We also acknowledge support from the U.S. Army's Tank Automotive Research, Development, and Engineering Center (TARDEC) and thank Dr Christian Hubicki for his support in conducting vertical intrusion experiments.

Notes and references

1. A. Soliman, S. Reid and W. Johnson, Int. J. Mech. Sci., 1976, 18, 279–284.
2. M. Omidvar, M. Iskander and S. Bless, Int. J. Impact Eng., 2014, 66, 60–82.
3. D. Ortiz, N. Gravish and M. T. Tolley, IEEE Robot. Autom. Lett., 2019, 4, 2630–2636.
4. D. Van Der Meer, Annu. Rev. Fluid Mech., 2017, 49, 463.
5. B. Andreotti, Y. Forterre and O. Pouliquen, Granular media: between fluid and solid, Cambridge University Press, 2013.
6. D. M. Wood, Soil behaviour and critical state soil mechanics, Cambridge University Press, 1990.
7. R. P. Jensen, P. J. Bosscher, M. E. Plesha and T. B. Edil, Int. J. Numer. Anal. Methods Geomech., 1999, 23, 531–547.
8. R. P. Jensen, M. E. Plesha, T. B. Edil, P. J. Bosscher and N. B. Kahla, Int. J. Geomech., 2001, 1, 21–39.
9. N. V. Brilliantov, F. Spahn, J.-M. Hertzsch and T. Pöschel, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 53, 5382.
10. L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, D. Levine and S. J. Plimpton, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 051302.
11. K. Kamrin and G. Koval, Phys. Rev. Lett., 2012, 108, 178301.
12. S. Kim and K. Kamrin, Phys. Rev. Lett., 2020, 125(8), 088002.
13. D. Jajcevic, E. Siegmann, C. Radeke and J. G. Khinast, Chem. Eng. Sci., 2013, 98, 298–310.
14. W. Zhong, A. Yu, X. Liu, Z. Tong and H. Zhang, Powder Technol., 2016, 302, 108–152.
15. C. Kelly, N. Olsen and D. Negrut, Multibody Syst. Dyn., 2020, 5, 355–379.
16. N. Vriend, J. McElwaine, B. Sovilla, C. Keylock, M. Ash and P. Brennan, Geophys. Res. Lett., 2013, 40, 727–731.
17. K. E. Daniels and N. W. Hayman, J. Geophys. Res.: Solid Earth, 2008, 113, B11.
18. S. Agarwal, C. Senatore, T. Zhang, M. Kingsbury, K. Iagnemma, D. I. Goldman and K. Kamrin, J. Terramechanics, 2019, 85, 1–14.
19. C. Li, T. Zhang and D. I. Goldman, Science, 2013, 339, 1408–1412.
20. S. Agarwal, A. Karsai, D. I. Goldman and K. Kamrin, Sci. Adv., 2021, 7, eabe0631.
21. N. Mahabadi and J. Jang, Appl. Phys. Lett., 2017, 110, 041907.
22. N. V. Brilliantov and T. Pöschel, Kinetic theory of granular gases, Oxford University Press, 2010.
23. J. T. Jenkins and M. Richman, Phys. Fluids, 1985, 28, 3485–3494.
24. A. Schofield and P. Wroth, Critical state soil mechanics, McGraw-hill, 1968.
25. H.-S. Yu, Keynote lecture. Proc. 12th Int. Conf. IACMAG, Goa, India, 2008, pp. 361–378.
26. P. Jop, Y. Forterre and O. Pouliquen, Nature, 2006, 441, 727–730.
27. J. Wiacek, M. Molenda, J. Horabik, et al., Mechanical properties of granular agro-materials. Continuum and discrete approach, Polska Akademia Nauk, 2011.
28. D. L. Henann and K. Kamrin, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 6730–6735.
29. K. Kamrin, Front. Phys., 2019, 7, 116.
30. S. Dunatunga and K. Kamrin, J. Fluid Mech., 2015, 779, 483–513.
31. G. Daviet and F. Bertails-Descoubes, ACM Trans. Graph., 2016, 35, 1–13.
32. G. Daviet and F. Bertails-Descoubes, J. Non-Newtonian Fluid Mech., 2016, 234, 15–35.
33. C. Dávalos, J. Cante, J. Hernández and J. Oliver, Int. J. Solids Struct., 2015, 71, 99–125.
34. C. Peng, X. Guo, W. Wu and Y. Wang, Acta Geotech., 2016, 11, 1231–1247.
35. S. Dunatunga and K. Kamrin, J. Mech. Phys. Solids, 2017, 100, 45–60.
36. D. W. Taylor, Fundamentals of soil mechanics, LWW, 1948, vol. 66.
37. S. Agarwal, Continuum Modeling of Granular Flows and Intrusions, 2020, MATLAB Central File Exchange.
38. Osborne Reynolds, Philos. Mag., 1885, 20(127), 469–481.
39. D. Sulsky, Z. Chen and H. L. Schreyer, Comput. Methods Appl. Mech. Eng., 1994, 118, 179–196.
40. A. Jarray, H. Shi, B. J. Scheper, M. Habibi and S. Luding, Sci. Rep., 2019, 9, 1–12.
41. J. Aguilar, T. Zhang, F. Qian, M. Kingsbury, B. McInroe, N. Mazouchova, C. Li, R. Maladen, C. Gong and M. Travers, et al. , Rep. Prog. Phys., 2016, 79, 110001.
42. H. Li, Y. Lai, L. Wang, X. Yang, N. Jiang, L. Li, C. Wang and B. Yang, Cold Reg. Sci. Technol., 2019, 157, 171–186.
43. X. Zhang, D. Sheng, G. P. Kouretzis, K. Krabbenhoft and S. W. Sloan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 022204.
44. F. Guillard, Y. Forterre and O. Pouliquen, Phys. Fluids, 2014, 26, 043301.
45. J. Wu, G. Kouretzis, L. Suwal, Y. Ansari and S. W. Sloan, Can. Geotech. J., 2020, 57, 1472–1483.
46. Y. Ding, N. Gravish and D. I. Goldman, Phys. Rev. Lett., 2011, 106, 028001.
47. M. Klinkmüller, G. Schreurs, M. Rosenau and H. Kemnitz, Tectonophysics, 2016, 684, 23–38.
48. R. L. De La Cruz and G. Caballero-Robledo, J. Fluid Mech., 2016, 800, 248–263.
49. A. Merceron, A. Sauret and P. Jop, EPL, 2018, 121, 34005.
50. S. Pravin, B. Chang, E. Han, L. London, D. I. Goldman, H. M. Jaeger and S. T. Hsieh, 2020, arXiv preprint arXiv:2010.15172.
51. A. H. Clark and R. P. Behringer, EPL, 2013, 101, 64001.
52. L. K. Roth, E. Han and H. M. Jaeger, Phys. Rev. Lett., 2021, 126, 218001.
53. H. Katsuragi and D. J. Durian, Nat. Phys., 2007, 3, 420–423.
54. T. A. Brzinski III, P. Mayor and D. J. Durian, Phys. Rev. Lett., 2013, 111, 168002.
55. W. Kang, Y. Feng, C. Liu and R. Blumenfeld, Nat. Commun., 2018, 9, 1–9.
56. C. S. Bester and R. P. Behringer, Phys. Rev. E, 2017, 95, 032906.
57. V. V. Sokolovski, Statics of soil media, Butterworths Scientific Publications, 1960.
58. A. F. Bower, Chapter 6, Applied mechanics of Solids, CRC Press, 2009.
59. J. Aguilar and D. I. Goldman, Nat. Phys., 2016, 12, 278–283.
60. L. Prandtl, Nachr. Ges. Wiss. Gottingen, 1920, 74–85.
61. R. Hill, Q. J. Mech. Appl. Math., 1949, 2, 40–52.
62. N. Gravish, P. B. Umbanhowar and D. I. Goldman, Phys. Rev. Lett., 2010, 105, 128301.
63. M. Kobayakawa, S. Miyai, T. Tsuji and T. Tanaka, J. Terramechanics, 2020, 90, 3–10.
64. E. Kashizadeh, J. Hambleton and S. Stanier, Proc. 14th Int. Conf. IACMAG, Kyoto, Japan, 2015, pp. 159–164.
65. Z. Jin, Z. Shi and J. Hambleton, Spree Internal Report, 2020, 20-7/495S.
66. Y. Feng, R. Blumenfeld and C. Liu, Soft Matter, 2019, 15, 3008–3017.
67. Y. Yue, B. Smith, P. Y. Chen, M. Chantharayukhonthorn, K. Kamrin and E. Grinspun, ACM Trans. Graph., 2018, 37, 1–19.
68. P. Y. Chen, M. Chantharayukhonthorn, Y. Yue, E. Grinspun and K. Kamrin, J. Mech. Phys. Solids, 2021, 104404.
69. F. Zhu, J. Zhao, S. Li, Y. Tang and G. Wang, Computer Graphics Forum, 2017, pp. 381–392.
70. M. Steffen, R. M. Kirby and M. Berzins, Int. J. Numer. Meth. Eng., 2008, 76, 922–948.
71. A. E. Minetti, G. Machtsiras and J. C. Masters, J. Biomech., 2009, 42, 2188–2190.
72. N. Sidelnik and B. Young, Sports Eng., 2006, 9, 129–135.
73. H. Askari and K. Kamrin, Nat. Mater., 2016, 15, 1274–1279.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sm00130b

---