Advances in Modeling Dense Granular Media

This review focuses on how the modeling of dense granular media has advanced over the last 15 years. The jumping-off point of our review is the (cid:22) ( I ) rheology for dry granular flow, which opened the door to generic flow field modeling but was primarily geared toward problems involving small monodisperse grains of simple shapes. Our review focuses on advances in modeling more material types and behaviors including new approaches for modeling finite-grain-size effects or nonlocality, polydispersity and un-mixing, and nontrivial grain shapes. We also discuss growing application areas with tractable order-reduction strategies with a focus on intrusion and locomotion problems.


INTRODUCTION
Granular media, despite their ubiquity in industry and geophysical applications, remain one of the most complex materials to model at scale.They are the second-most-handled material by weight (de Gennes 1999) in global industry (second only to water), and yet their dynamic properties have remained an elusive scientific question going back to the days of Coulomb.The 125th anniversary issue of Science (AAAS 2005) listed a general theory for the motion of granular media as one of the 125 critical scientific knowledge gaps.
The difficulty in modeling granular material comes from several sources.First, there are many material properties at play in granular media since they comprise many distinct individual particles.Developing a model that accounts appropriately for the distribution of grain shape, size, roughness, elasticity, restitution, and other properties is a very daunting task.Granular media also tend to unmix during flow, so the property distributions can and do change over space and time.Even if they stay well mixed, the mechanical behavior changes with flow history, and the negligibility of thermal fluctuations render many equilibrium principles inapplicable.Moreover, the grain size is typically not well separated from the size scale of the experiment-one can often look at an entire apparatus and still make out the individual grains.This means that the role of the microscopic size scale may not vanish even in the continuum limit.A discrete particle approach can be used in place of a continuum picture.However, this is often not a tractable solution given the number of grains involved; a handful of sand is already around one billion grains.
Despite these complicating factors, in the early and mid-2000s, the advent of the inertial rheology constituted a major advance, thanks to the combined efforts of several laboratories in the south of France, as described in MiDi (2004).It was first demonstrated in three dimensions by Jop et al. (2006), and a thorough review of the model was given by Forterre & Pouliquen (2008).Also known as the µ(I) model, the inertial rheology describes granular flow in the case of small, hard, quasi-monodisperse particles flowing steadily under a finite pressure.The basic premise of the model is that there is a one-to-one relationship between the ratio of shear stress to pressure acting on a granular material element, µ ࣕ τ /P, and the normalized flow rate or inertial number, I ≡ γ d √ ρ s /P, where d is the mean grain diameter and ρ s is the density of a grain.In the dense (noncollisional) limit, this function is often found to obey a roughly linear form (da Cruz et al. 2005) for constants µ s and b: The goal of the current review is to discuss some of the crucial advances that have been made in the understanding of granular media since the inertial rheology.We focus on the dense (i.e., liquidlike to solid-like) response, as well as questions stemming from material limitations of the µ(I) model.First, what if the grains are not small?And what if the grains are not quasi-monodisperse?Advances on these two questions have come in the form of nonlocal continuum models and models of particle segregation, discussed in Sections 2 and 3, respectively.Then, in Section 4, we consider the question of particle shape and shape distribution: How can shapes be modeled even at the small scale and how do they produce the large-scale mechanical properties emergent from the media?Lastly, despite the complexity of granular media and the efforts over recent years to model an everincreasing number of flow phenomena, in Section 5, we conclude with a discussion of application areas in which predictive models have emerged that are remarkably simple.This section focuses on reduced-order models for problems of impact and intrusion.

ACCOUNTING FOR FINITE GRAIN SIZE
The µ(I) model is a mechanically local rheology.We use the term "local" to refer to a constitutive relation that, when a flow reaches steady state, poses a direct algebraic relation between the stress tensor σ and the velocity gradient v at the same point x.Of note, a local relation does not relate any higher gradients of stress or velocity.Despite the commonality of local models, the fact that grains have a finite mean size d that is not always well separated from the geometry's size means that granular flows often display behaviors indescribable with a local form.For example, even though the µ(I) rheology holds rather well in homogeneous shear flows, the same relationship fails if the same grains are made to flow inhomogeneously (Koval et al. 2009).The fact that the presence of a strain-rate gradient causes a departure from the µ(I) model implies that granular rheology must be influenced by an intrinsic length scale that competes against flow gradients, implicating the importance of d.Likewise, the rheology must have a grain-size effect, rendering it nonlocal.The issue of determining a complete nonlocal form has made the problem of predicting granular flow fields very challenging, even for simple spherical beads.

Nonlocal Phenomena in Granular Flows
While instances of nonlocal effects have been considered in the soil mechanics literature for many decades, primarily based on micropolar or Cosserat continua (Cosserat & Cosserat 1909, Mühlhaus & Vardoulakis 1987, Mohan et al. 2002), over the last 20 years or so the fluids and statistical mechanics literature has increasingly appreciated nonlocality in dense granular flowsthat is, problems going into the very large strain regime, approaching more of a fluid-like flow.One of the earlier observations of grain-size effects on developed flow fields was in the silo flow literature (Litwiniszyn 1963;Mullins 1972Mullins , 1974;;Tüzün & Nedderman 1979).It was observed that the downward velocity component of grains in a flowing silo appeared diffusive in character, with the horizontal distribution of downward velocity appearing to spread diffusively with height, i.e., v z (x, z) ∼ e −x 2 /4Bz .The fact that the spread factor B is proportional to grain size d points to a nonlocal flow mechanism governed by the grain scale.The kinematic model of granular drainage (Tüzün & Nedderman 1979) and its related mechanisms [the void model (Mullins 1972) and spot model (Bazant 2006)] were inspired by the observation of silo diffusion.This diffusive effect causes a creeping decay of shear rate deep into the bulk that is not correctly predicted by inertial rheology (see Figure 1a), which predicts sharp flow/no-flow transitions.
Similarly, later studies on steady 2D flows in other geometries [reviewed by, e.g., MiDi (2004)] described annular shear flows, shear flows with gravity, vertical chute flows, rotating drum flows, and heap flows as exhibiting diffusive-looking flow fields with similar characteristics (Komatsu et al. 2001) (e.g., Figure 1b).In all cases, the flow field resembles a solution to the µ(I) rheology qualitatively, but with creeping, decaying tails of flow extending deep into the zones that µ(I) rheology predicts to be static.The decay length scales with grain size, implicating a particlesize effect.Perhaps the most infamous 3D example of this effect is the split-bottom Couette cell (Fenistein & van Hecke 2003, Depken et al. 2006).Here, the surface velocity profile has a central wide shear band far from all walls with a striking resemblance to an error function (Figure 1c).Nonzero strain rate occurs everywhere, while the µ(I) rheology predicts a sharp central shear band surrounded by totally rigid zones.
The influence of nonlocality can be seen outside of steady flow fields.For example, it affects the flow arrest criterion of systems with thin features that compete with the d size scale.Perhaps the most famous example of this is hopper clogging, where the flow can arrest by forming an arch over the opening even if the opening size is nonzero (Beverloo et al. 1961, Zuriguel et al. 2014).This runs counter to the µ(I) rheology, which predicts nonzero flow unless the opening size is reduced to zero (Staron et al. 2012).Note the irony here: The previous paragraph discusses how nonlocality can cause more material to flow than µ(I) would predict, but now we point out that nonlocality can also cause material to stay static when µ(I) would predict flow.A somewhat more controllable experiment for this effect is flow stoppage of a thin layer of grains on a rough incline (Figure 2a) (Pouliquen 1999;Silbert et al. 2001Silbert et al. , 2003)).The µ(I) model, like other local granular flow models, always predicts a single repose angle for a material.While thicker layers above O(10d ) are observed to arrest at an incline angle relatively independent of thickness, H/d, in thinner layers the arrest angle increases continually as the layer thins.This H stop effect and the silo clogging phenomena indicate that as a geometric length scale gets thin relative to d, the granular material acts stronger.The same effect occurs in the arrest of thin vertical chutes and annular shear cells (Liu & Henann 2018).An analogous strengthening effect can be seen in rod intrusion; the force needed to sustain drag increases as the grain size to rod width ratio increases above ∼0.1 (Seguin et al. 2011, Guillard et al. 2014).
Another unusual nonlocal effect is referred to as secondary rheology.To wit, probes within a body of granular media can be moved with an arbitrarily small force as long as the body is flowing somewhere, even if the primary flow is far away from the probe.Conversely, if the body is static, a yield force exists that must be overcome in order to move the probe.For example, Reddy et al. (2011) inserted a rod into an annular shear cell and applied a force to it (Figure 2b).When the inner wall of the cell is stationary, the rod must exceed a finite yield force F c depending on the rod depth and submerged volume in order for it to move.However, if the inner wall is rotating, causing a wall-located shear flow, the rod will move if any force is applied to it, even if it is in the quiescent material far from the wall.The rod's speed obeys v rod (r) ∼ v wall F(r)sinh ((F − F c )/F c G(r/d)).The explicit dependence on grain size in the (decreasing) function G(r/d) implies that a size effect is involved.Secondary rheology has been observed in other cases where a faraway moving boundary removes the yield stress everywhere, allowing force probes to move even for small forces (Nichol et al. 2010, Wandersman & Van Hecke 2014).

Nonlocal Models for Granular Flow
Over the last 15 years or so, a variety of approaches have been proposed to extend the µ(I) model to account for the nonlocal phenomena described in the previous section.Most recently, proposed models have focused on capturing nonlocal creep.Fewer have also been demonstrated to capture thin-body strengthening and secondary rheology-many of these models are not of a mathematical form capable of representing these effects.Nonlocal models often cannot be solved with common codes, which means users must create their own general solvers, restrict attention to situations with simply calculable solutions, or attempt to draw inferences without fully constructed solutions.Hence the extent of testing of the various proposed models varies greatly.Granular kinetic theory (Haff 1983, Jenkins & Savage 1983, Lun et al. 1984, Garzó & Dufty 1999, Brilliantov & Pöschel 2010) could be considered a mechanically nonlocal form since it couples the granular viscosity to a collisional temperature T that obeys a diffusive energy balance equation, where ϵ ij are the strain components.Both the viscosity η and the conductivity K are proportional to ∼Pd/ √ T and the volumetric dissipation rate is proportional to (1/d)T 3/2 .The way kinetic theory causes nonlocal creep is a coupled process: Stress drives flow in high-stress regions, which increases T there; T then diffuses everywhere, increasing temperature even in zones of low stress; the increase in T in low-stress regions reduces the viscous resistance, thereby promoting flow there despite the low stresses.Hence, nonlocal creep emerges as a self-heating process, an idea common to a variety of the nonlocal models we discuss.The length scale for the self-heating is d due to its role as the length scale in K. Extensions of kinetic theory to the dense limit ( Jenkins & Berzi 2010, Berzi 2014) have emerged over the last decade by generalizing the length scale to a larger multiple of d.Recent empirical power-law correlations show that temperature sensitivity is key to describing nonlocal creeping (Kim & Kamrin 2020, 2023).
Diffusion-based self-heating models are promising since they can at least qualitatively reconcile nonlocal creeping and secondary rheology-through diffusion of fluidization and subsequent material weakening-as well as thin-body strengthening-through wall boundary conditions that absorb fluidization-in a compact mathematical form.Some models employ the self-heating idea holistically by coupling the rheology to a diffusing fluidization variable in an attempt to formulate a simpler and potentially more quantitative (albeit more empirical) model.The first proposed holistic model of this type is the partial fluidization model (Aranson & Tsimring 2001, 2002;Aranson et al. 2008), which utilizes an order parameter ρ (not density) and assumes Landau-like purely dissipative dynamics sourced by a stress-dependent term: 3.
The viscosity is then presumed to be ∼1/(1 − ρ), so the material solidifies as ρ tends to 1 and is fluid-like for ρ < 1.The term δ drives the solution toward a hysteretic frictional response with dynamic friction µ d as 2 ρ tends to 0. The Laplacian term, scaled by l ∼ d, then diffuses the order parameter in inhomogeneous flow cases, giving rise to nonlocality in a fashion analogous to that of kinetic theory.Partial fluidization theory has been shown to qualitatively produce nonlocal creep and thin-body strengthening, although the demonstrations have been limited to quasi-1D examples.
A more recent diffusive self-heating model is the nonlocal granular fluidity model (NGF).Introduced in 2012 (Kamrin & Koval 2012) and inspired by a related model for emulsions (Bocquet et al. 2009), NGF is designed to give quantitative outcomes and a broad representation of the three aforementioned families of nonlocal phenomena.These aims are achieved by proposing the existence of a nonnegative state field called the granular fluidity, g, which satisfies the following system: Notice in the first relation that g behaves as a pressure-weighted inverse viscosity; g = 0 implies no flow.Most parameters above are borrowed from µ(I) rheology (e.g., ρ s , b, µ s ).The only new material parameters are the timescale t 0 (which is irrelevant at steady state) and the dimensionless nonlocal amplitude A, which tends to be order-one in size [≈0.45 for spherical beads (Henann & Kamrin 2013) and ≈0.85 for 2D disks (Kamrin & Koval 2012)].In the absence of flow gradients (or negligible d), the two equations above, at steady state, produce a fluidity that replicates the µ(I) model (Equation 1).But when flow gradients are present, the diffusion term spreads fluidization based on the grain size.
The NGF partial differential equation (PDE) can be reduced to the following approximate form at steady state: where ⟨•⟩ is the Macaulay bracket, ξ is the apparent cooperativity length, and g loc is what the fluidity formula would be if the local µ(I) model were correct.This form elucidates some of NGF's key properties.First, if µ ≤ µ s holds everywhere, g loc vanishes and the resulting PDE, g = ξ 2 2 g, is linear in g.Hence, the rheology becomes rate independent in the slow, nonlocally dominated regime.This odd feature has been observed in various studies (MiDi 2004, Koval et al. 2009) and highlights a dichotomy between nonlocally dominated dynamics and local µ(I) rheology, which is markedly rate dependent.The blow-up dependence of ξ on µ has been observed in steady flows in annular shear cells (Kamrin & Koval 2012, Tang et al. 2018) and split-bottom Couette flows (Henann & Kamrin 2013).
NGF has been tested in many circumstances with the aid of multiple custom-written numerical solvers using finite-element analysis (Henann & Kamrin 2013, 2016), finite volume methods (Faroux et al. 2022), and meshless methods (Dunatunga & Kamrin 2022).Using glass beads and hard disks as common benchmark materials, researchers have used NGF to make quantitative predictions of steady flows in two dimensions, including annular cells, shear under gravity, and vertical chutes (Kamrin & Koval 2012); flow fields in three dimensions, including split-bottom cells (Henann & Kamrin 2013) and wall-bounded heaps (Liu & Henann 2017); the H stop curve (Kamrin & Henann 2015), silo stoppage (Dunatunga & Kamrin 2022), and other thin-body strengthening diagrams (Liu & Henann 2018); and secondary rheology phenomenology in Reddy et al.'s (2011) rod-probe geometry (Henann & Kamrin 2014) as well as other forced probe tests (Li & Henann 2020).While the underlying physics of NGF is not fully understood, a connection between g and temperature has been established (Zhang & Kamrin 2017, Kim & Kamrin 2020), suggesting that the NGF PDE acts as an empirically simplified heat equation for a conglomerated thermally related variable.
We mention in closing other recent nonlocal formulations based on approaches other than a diffusing state variable.The I-gradient model of Bouzid et al. (2013Bouzid et al. ( , 2015) ) supposes µ = µ(I) (1 − ν 2 I/I), a direct nonlocal extension of the µ(I) model.It assumes that friction increases when neighboring media flow slower and vice versa.It captures the emergence of nonlocal rate independence, quantitatively predicts nonlocal creep in 2D test cases, and has the benefit of not requiring an auxiliary state field.However, it cannot capture thin-body strengthening and it has been shown to be mathematically ill posed in the sense of Hadamard (Li & Henann 2019), which could complicate solution procedures.Another set of models are based on system-wide integral equations (strongly nonlocal models) rather than PDEs.These include a body-wide self-activation model where perturbations from flow everywhere influence shear rate at each point (Pouliquen & Forterre 2009), a model based on convolving local rheology over a kernel (Nott 2017), and an integrated eddy-viscosity formulation (Miller et al. 2013).As these models are based on integral equations, they are more complex to solve, but PDE semilocalized approximations could be considered.

ACCOUNTING FOR GRAIN PROPERTY DISTRIBUTION
Up to this point in the review, our discussion has focused on how size affects the response of a granular material in a relatively monodisperse system (monodisperse in size, density, shape, roughness, elasticity, etc.).Most natural particulate systems comprise particles differing in these properties.When subjected to shear, jostling, or other excitation, the constituents of such systems tend to unmix by particle properties, resulting a system nonuniform not only in characteristic particle properties but also in local response to applied stresses and other excitations (Barker & Mehta 1993, Knight et al. 1993).Such apparent sorting or segregation gives rise to a range of striking patterns in nature, from fingering in pyroclastic flows to sorting boulders to the front in debris flows (Figure 3).Typically, as in the previous section, local and global constitutive behavior is not linearly related to particle size or other properties even for small changes in particle size distribution.Thus, to capture the response to applied stresses or other excitation, we need to simultaneously capture the spatial variability of particle properties and the feedback between segregation and constitutive behaviors over time and space.
Here, we focus on segregation associated with sheared granular flows, as these are the conditions under which the most significant advances have been made in the last decade or so.In addition, we focus on what we might call the mechanisms, or drivers, responsible for dynamically segregating particles, sending them preferentially in different directions according to their properties.Other factors such as advection associated with complex boundary conditions play a role in the final segregation pattern obtained in a relatively dense system; associated problems saw significant advances in the 1990s and early 2000s (Ottino & Khakhar 2000).A commonly reported observation is that when mixtures of particles differing only in size are sheared, shaken, or otherwise disturbed, the larger particles segregate to the top-the so-called Brazil nut problem (Gray 2018).The situation is rarely this simple (Hill 2016, Hill & Fan 2016).Larger particles may segregate upward, downward, neither, or both, depending on boundary conditions and other particle properties (Félix & Thomas 2004).We use this relatively simple context to examine some of the recent developments in segregation modeling.

Segregation and Kinetic Theory
Kinetic theory applied to granular materials can capture segregation spontaneously.The evolution of local concentrations is derived based on first principles based on interparticle interactions and external forces such as gravity.Model predictions result from deriving transport coefficients through moment analysis of the Boltzmann equation.Simplifying assumptions that are typical include the form of the velocity distribution functions, two-particle interaction models, and simplifications around collision integrals.Over the last 30 years, most kinetic theory formulations for granular flows involve assumptions of Maxwellian distribution of particle velocities and assumptions of binary collisions and collisional chaos.In other words, particles only contact one other particle at a time with velocities independent of one another and of any previous interparticle interactions.Complete pair distribution functions at collisions are approximated by the product of single particle velocity distribution functions for each disk and a factor akin to a radial distribution function.In theory, these models can be completely accurate if the assumptions hold.
Before considering these problems at higher density in detail, we consider the strengths in the context of segregating particles.We can think of the segregating dynamics in this framework as strictly the dynamics associated with differing collisional interactions and external forces (e.g., gravity) for the particle constituents.In the context of continuum theories, we consider averaged macroscopic dynamics associated with the continuum expressions derived from kinetic theory.For an example, we consider segregation associated with variations in both particle size and density of constituents A and B, based on work by Yoon & Jenkins (2006), in a shear flow in the direction of gravity.For the expression to be somewhat easily represented, size and mass ratios are relatively small, and dissipation is small but nonnegligible.To make the segregation drivers more apparent, we express a segregation flux of particles of type B relative to the mixture The four terms on the right-hand side of Equation 6 describe the effect of gravity, partial pressure, granular temperature gradient, and diffusion, respectively, as emphasized in Equation 7. K AB is a radial distribution function of contacting particle pairs, which includes the local solid volume fraction of the mixture; T is a local granular temperature measuring average kinetic energy associated with velocity fluctuations relative to the barycentric velocity; and π i is a (local) partial pressure component associated with species i, functions of T, and the granular temperatures of each species i, T i .We note that the third term (often called a thermal diffusion term owing to the dependence on granular temperature gradient) does not result in segregation because of an underlying assumption by Jenkins and Yoon of "equipartition of energy" between the species.Other theoretical derivations (e.g., Garzó 2009) predict an importance of this term in moderately dense mixtures.The diffusion term depends on the number densities of the two species relative to one another (n A and n B ).Many of the prefactors in Equation 7 are complicated functions of relative particle sizes and densities.However, it is notable that the relative segregation flux direction depends directly and solely on the relative mass of the two particles.All else being equal, the more massive particles segregate toward decreasing temperature.Under certain boundary conditions, predictions of this type successfully capture segregation dynamics produced in experiments and particle-scale simulations (Conway et al. 2006, Fan & Hill 2011b).But at higher solid fractions, the assumptions above typically overpredict the intensity of segregation, and at the highest solid fractions, segregation of different-sized particles could oppose that predicted from the above kinetic theory (Conway et al. 2006).Discrete element method (DEM) simulations demonstrate that theory and results start to disagree when multiple particle contacts grow to a scale that dominates the system (Fan & Hill 2011b).In the last decade, kinetic-theory-based models for shear flow of mixtures that start with assumptions similar to those outlined above have been extended to include representations of correlations of forces of a longer range than the particle scale and whose average pressure and temperature fields more accurately reflect boundary conditions (Larcher & Jenkins 2019).Yet, to date, the mixture dynamics of the interspecies interactions such as those described by Equations 6 and 7 remain based on the assumptions provided by the theory with no changes in the response to segregation drivers such as temperature gradients.

Phenomenological Segregation Models and Mixture Theory
Over the same time period, researchers developed complementary macroscale models for segregation in denser flows where assumptions of kinetic theory break down.These have increased the breadth of boundary conditions that could be represented mathematically.Several original details of segregating dynamics are relatively phenomenological and, unlike the models associated with kinetic theory, are essentially superposed with a representation of the average dynamics.For simplicity, we contextualize this framework in simple shear, where the system is uniform in all directions except one, and define both the direction of nonuniformity and the direction of segregation in the −y direction.Early examples took the form of a gravity-driven term and a diffusion (mixing) term for segregation flux: For example, Khakhar et al. (1997) proposed a buoyancy mechanism for gravity-driven segregation that models segregation of particles of similar size and different material density akin to the pressure on a submerged object in a surrounding fluid of different density: Superscripts ld and d refer to less dense and dense species, respectively, and subscript m refers to the material density of the particles (as opposed to a local mass density of the mixture).C KMO and D KMO are drag and diffusion coefficients, respectively, which are both determined phenomenologically.
For different-sized particles, Gray and colleagues (Gray & Thornton 2005, Gray & Chugunov 2006) captured effects of kinetic sieving proposed by Savage & Lun (1988) in a continuum mixture-theory framework.Statistically, in a mixture of different-sized particles, all particles are pulled in the direction of gravity, and small particles are more likely to find an opening they can fit below a large particle rather than the other way around.Hence, they are sieved downward, and the large particles are pushed upward.Gray and colleagues captured this effect using partial pressures; the large particles bare a higher fraction of the pressure (associated with gravity) than the smaller ones and are pushed upward under the pressure gradient.Segregation fluxes of large and small particles (superscripts l and s, respectively) are expressed according to ϕ l where b i (i = s or l), C GTA , and D GTA are coefficients of segregation, drag, and diffusion, respectively, which are all initially set on a case-by-case basis.Slightly more specifically, b i scales with the magnitude of the deviation from the equipartition of pressure; b l > 0 and b s < 0. We note the approximate mathematical equivalence with Equation 9, particularly for the cases where the phenomenologically determined coefficients are all constant.However, additional work is needed to clarify the underlying physics.
In the last decade advances to these frameworks have been made, thanks in part to advances in computer modeling and physical experiments.Some researchers have found fluid-like analogs useful for gaining insights into segregation.Specifically, not unlike internal viscous and turbulent stresses responding to applied stresses in fluids, applied stresses associated with gravity in granular mixtures may evoke internal stresses associated with interparticle contacts and stresses associated with granular temperatures.Each of these may be partitioned independently between species.Hill and colleagues (Fan & Hill 2011a, Hill & Tan 2014) showed that including this in the mass and momentum balance equations of the mixtures gives rise to the following form for segregation fluxes: i THF = i g,THF + i T ′ ,THF + i D,THF , which is similar to Equation 8, with an additional temperature gradient term.The form of each term on the right-hand side is as follows: 10.
The partitioning of the stresses is represented by ψ k,i = σ k,i yy /σ k yy and ψ c,i = σ c,i yy /σ c yy ; σ k yy is the normal stress associated with the fluctuations around the barycentric velocities (and scales with local granular temperature); σ c yy is the stress associated with the contacts only; ϕ i m is the local mass fraction associated with species i; and ρ i = ϕ i ρ, where ρ is the local mass density of the mixture.In mixtures of spherical particles, typically the kinetic stress gradient term dominates over the gravity term because the relative particle density scales with the partitioning of the contact stress.
Other recent work has improved other details of the physics such as diffusion and drag behaviors (Umbanhowar et al. 2019).Inspired by other fluid-like analogs, Sarkar & Khakhar (2008) and Tripathi & Khakhar (2013) derived representative drag, diffusion, and segregation coefficients from first principles for different density particle mixtures.Another example involves Saffman-like segregation forces on an intruding particle (van der Vaart et al. 2018), not unlike the squeezeexpulsion accompanying the kinetic sieving first modeled by Savage (1998).Other advancements apply to stress partitioning and particle asphericity: Hernández-Delfin et al. (2022) showed that for ellipsoids, both kinetic stress and gravity terms can dominate, depending on the aspect ratio of the species in the mixture.

ACCOUNTING FOR NONTRIVIAL GRAIN SHAPES
Advances in grain-scale measurements and models have brought the concept of a virtual laboratory for granular rheology increasingly in reach, the goal of which is the production of high-fidelity micro-and macroscopic data that can be correlated to derive either next-generation constitutive models that account for microscopic physical variables as part of their internal variables or datadriven models that bypass the constitutive model and evaluate the evolution of stress in space and time.This section explores these ideas in the context of particle morphology as an important micromechanical feature in the description of granular solid-like and fluid-like behavior.It is clear that grain-scale modeling and observations should serve to complement continuum models rather than replace them.Physics-based predictive grain-scale models allow us to do what experiments can do and more since they give us access to classic macroscopic measures such as stress and strain, while also providing micromechanical information such as particle kinematics and kinetics that have thus far been elusive in most experimental settings.For example, grain-scale models can routinely provide force chains deep inside a 3D wet or dry system, which is only possible with advanced experimental techniques (e.g., X-ray computed tomography, diffraction, neutrons) that are not readily available and have numerous limitations.In addition, these models enable the reproducibility of initial conditions and controllability of boundary conditions that are not possible in physical experiments.Hence, virtual experiments accounting for particle shape offer an opportunity to redefine granular continuum descriptions.
The continuum description of deformable materials is a topic that has seen significant progress over the past century (Maugin 2013), which has primarily been driven by constitutive models developed in the framework of continuum thermomechanics (Gurtin et al. 2010).This modeling approach requires fitting phenomenological functions and parameters whose link to the underlying physics is often obscure or unavailable.In the context of granular mechanics, one possible alternative to continuum modeling lies in DEM simulations (Cundall & Strack 1979), a framework developed in the 1970s for direct micromechanical simulations.Initial DEM formulations only considered circular or spherical particles.Consequently, models based on discrete elements were strongly limited to descriptive rather than predictive investigations of granular materials (Christoffersen et al. 1981).This limitation naturally stems from the fundamental influence of particle morphology on the salient macroscopic properties of granular media, including stiffness, strength (Cho et Overview of (a) standard spherical, (b) sphere-clumping, (c) polyhedral, (d) potential particle, and (e) level-set discrete element methods.f(x) and ϕ(x) are example level sets.
Several improvements of the DEM were later proposed to provide more realistic representations of grains, focusing on the pressing issue of particle shape.One of the first efforts employed ellipsoidal particles (Ting et al. 1993), while a more versatile approach relied on the clumping of spheres (Garcia et al. 2009b) (Figure 5b).The latter method is perhaps the most popular and pragmatic but requires a very large number of spheres to approximate the local curvature of actual grain shapes.Another popular approach considers polyhedral elements (Cundall 1988, Alonso-Marroquin & Herrmann 2002) (Figure 5c), potentially providing a more accurate approximation of the local morphological features.However, these methods require many vertices to approximate the local curvature, rendering the contact interaction computationally demanding and requiring cumbersome algorithms, particularly for 3D simulations.
More recent variants of DEM have tackled the issue of particle shape from different fronts.An example is the combined finite-discrete element method (Munjiza 2004), where each particle is treated as a continuum and discretized into a finite element mesh to capture the (generally finite) deformation field.This method provides capabilities of particular relevance for modeling deformable grains but implies a prohibitive computational cost for large numbers of particles.Focusing again on approximately rigid particles, other efforts have focused on dedicated mathematical representations of contact between arbitrary particle shapes.One method of this nature that deviates from the soft-contact formulation of DEM is the convex complementarity approach (Tasora & Anitescu 2010), where nonoverlapping contacts are enforced geometrically through differential variational inequalities.The method does not rely on contact stiffness parameters, allowing for large simulation time steps, but each time step requires the solution of a computationally demanding and rather complex optimization problem that may lead to force indeterminacy in underconstrained granular systems (typical in three dimensions).Conversely, the potential particle approach (Houlsby 2009) (Figure 5d) remains faithful to the overlapping particle approach of DEM but describes general particle shapes through convex functions whose level surfaces characterize the contact interaction.This method assumes convex shapes and the knowledge of explicit shape-dependent morphological functions.More recently, the granular element method (Andrade et al. 2012) was proposed to capture general particle morphologies accurately, including nonconvex shape interactions (Lim et al. 2014), making use of nonuniform rational basis splines at the cost of relatively high computational requirements.
Despite this progress, a notable gap was observed between the capabilities of experimental and computational studies on granular media.This gap motivated the development of imaging techniques to characterize particle shapes extracted directly from X-ray computed tomography (XRCT), achieving a seamless transition from experimental images to virtual grain avatars (Vlahinić et al. 2014).This inspired the development of the level set discrete element method (LS-DEM) (Kawamoto et al. 2016) (Figure 5e), a framework that overcomes many limitations of the aforementioned DEM techniques, providing a simple contact interaction procedure at a computational cost that generally lies between standard DEM and other shape-enhanced variants.LS-DEM defines a level set function for each particle, which convects with the particle kinematics.The level set enables the evaluation of distance to the boundary of the particle (the zero-th level set), which is used to readily evaluate any interpenetration between two particles using classic leader-follower approaches in contact mechanics.LS-DEM has been used to conduct high-fidelity virtual experiments of real granular materials extracted from XRCT and has accurately captured, for instance, shear band formations (Kawamoto et al. 2018).Its capabilities have been extended beyond grains to, for example, structured fabrics that exhibit a transition from fluidlike to solid-like behavior (Wang et al. 2021).In a recent study (Lai et al. 2022), a unified framework has been presented that encapsulates LS-DEM and other DEM variants as particular cases.
Still, large-scale simulations of discrete models remain intractable with current computational resources.How can we then exploit the use of DEM for practical engineering problems?One possibility is to use modern DEM models to guide the development of new constitutive models or the calibration of existing ones.In this context, constitutive models for granular media featuring a direct link to grain-scale mechanisms can be found in the literature (Nicot & Darve 2005, Einav 2007, Misra & Poorsolhjouy 2015).Unfortunately, the crucial effect of shape is persistently absent in micromechanics-based constitutive models, likely due to the difficulty of arriving at functional forms that embed such complex characteristics.The constitutive model presented in a recent study by Buscarnera & Einav (2021) is an exception; it considers the effect of coevolving grain shape and size distribution, providing some hints for future steps in this direction.Likely, a key step is to provide a robust mapping between shape parameters and continuum models through a reliable quantification of shape distribution.
Another possibility is to employ shape-enhanced DEM in a multiscale modeling framework.The multiscale or upscaling approach dates back to the so-called FE 2 approach, which uses FE simulations at the macroscale that call lower-level (micro) FE simulations at the material point level.The idea was simple and powerful and it established the field of multiscale computational mechanics.In granular mechanics, variants of the FE 2 approach have been proposed within the socalled FEM (finite-element method)×DEM approach, where a discrete element model is called as a surrogate for constitutive models.This approach can be as accurate as the underlying DEM but it has a few shortcomings: (a) It is unable to account for data stemming from physical experiments, (b) it requires new simulations for every material point for stress-strain data, and (c) it may lose predictive accuracy in cases where the flow gradients are large, making the underlying micromechanical DEM simulation nonrepresentative at a given Gauss point in the macroscopic computation.To ameliorate these shortcomings, researchers have recently proposed data-driven approaches where data stemming from physical and virtual experiments can be amalgamated and used to replace constitutive responses on demand.While this field is at its infancy, it offers the attractiveness of purely using data to model (granular) mechanics problems.Whether one chooses to upscale the information directly using data-driven approaches or by constructing micromechanics-informed constitutive models, it is clear that the accuracy of the microscale data is of paramount importance.Morphological data, including shape, grain size distributions, and their coevolution, provide critical parameters to reckon with for the next-generation of constitutive models.It is encouraging to see that physical and virtual experiments are paving the road to accurate and tractable models, as shown in the next section.

LOCOMOTION STUDIES PROVIDE AND MAKE USE OF NEW ADVANCES IN MODELING
In this section we discuss how studies of locomotion of animals and robots on and within granular media have led to new understanding of granular intrusion processes.A key principle is that reduced fluid-like models quantitatively capture the performance of such systems; many organisms at the macroscale locomote via similar principles to microscopic swimmers.
For decades many of the basic science studies of bulk granular flow have drawn inspiration from common grain-handling devices; these include chute flows, silos, and fluidized beds.One characteristic of such problems is that forcing tends to come from the motion of, or fluxes through, boundaries.Less attention has been paid to problems that involve local forcing, such as the intrusion of an object along a path (e.g., the locomotion of an animal or a robot in dry desert sand).Problems involving the impact of objects into granular media under gravity are the exception, which were first studied in relation to gunnery and later as models of meteor impact.Scientific studies on this topic began in the seventeenth century and continued intermittently until becoming of much interest in the 1990s and 2000s.Diverse phenomenological models were proposed and debated, and DEM simulations provided guidance and insight into the grain-scale mechanics responsible for impact in uniformly shaped objects (e.g., spheres, disks, cylinders, rods).Essentially, a phenomenological form for the impact force was obtained that contained terms quadratic in intruder speed and increasing in depth (Katsuragi & Durian 2007).Later modifications based on volume fraction, object shape, and so on were proposed (Goldman & Umbanhowar 2008, Umbanhowar & Goldman 2010).While impact is amenable to experimental simplicity, other studies have examined more controlled scenarios (Albert et al. 1999(Albert et al. , 2001)), largely horizontal and vertical intrusion of simple geometrical shapes (cones, disks, plates, etc.).Typical features of such intrusions are the dependence of the drag force on cross-sectional area, the increasing force with increasing depth, and the weak dependence of force on drag speed.Localized intrusions are also highly relevant in biological and robotic locomotion, and recently these have motivated a broader look at intrusion phenomena.Initially studied phenomenologically, a stronger theoretical foundation has emerged in recent years thanks to advances in continuum modeling.
We now discuss broad categories of granular dynamics relevant to locomotion on and within granular media; more details can be found in the cited literature and in a recent review (Astley et al. 2020).Diverse organisms live in dry desert sand and move on the surface as well as below (Mosauer 1932, Gans 1973, Ding et al. 2013a, Hosoi & Goldman 2015, Irschick & Jayne 1999).Extraterrestrial surfaces like that of the moon and Mars are composed of so-called regolith, which must be navigated by rovers.However, study of organisms and robots in natural media is a challenge.Starting in the mid-2000s, Goldman and coworkers embarked on a program to systematically and comparatively study locomotion in dry granular media in various forms and conditions (e.g., Figure 6).Using robots as robophysical models, Goldman's group developed  hypotheses for form and function in organisms that use limbless, legged, and mixed (body and leg) approaches for propulsion.The first such studies were conducted on limbless locomotors: the sandfish lizard and the shovel-nosed snake, two organisms that have adapted to desert conditions via specialized morphology (e.g., nose shape, low-friction scales).Using high-speed X-ray imaging and air fluidized beds to control the initial bed conditions, they discovered that these organisms propelled themselves using head-to-tail waves of body curvature (lateral undulation) (Maladen et al. 2009, Sharpe et al. 2015).Organisms controlled for wave shapes and parameters and, remarkably, resembled animals at much smaller scales (e.g., nematodes).At roughly the same time, Goldman's group began study of legged systems moving on the surface of granular media (Li et al. 2012).Here it became clear that adaptations led to locomotor performance seen in seemingly nonspecialized animals like hatchling sea turtles (Mazouchova et al. 2010), which use flexible wrist flippers to propel themselves forward without significantly yielding nearby material.This principle of minimal slip carried over into the locomotion of the highly specialized sidewinder rattlesnake (Marvi et al. 2014), one of three species of snake that moves effectively 90°to the direction one would expect a snake to slither in.As the slope angle increases from level to 20°, the sidewinders use more of their body surface to engage the granular material.Limbless locomotors that do not deploy this strategy tend to fail to climb granular slopes of even modest incline angles.Robophysical models (Aguilar et al. 2016, Aydin et al. 2019) have also been used to elucidate important features associated with different locomotor modes and regimes.Using such tools to construct models of the sand-swimmers, Maladen et al. (2011b) discovered that a multimotor linear limbless robot propagating an approximately sinusoidal wave of body curvature from head to tail performed optimal granular swimming in a local frictional fluid in terms of displacement per undulatory wave cycle at a wave amplitude to wavelength ratio comparable to that of the sandfish lizard.In addition, modifications to the robot in the form of heads of varying angular shape produced varying lift forces.Wedge-shaped heads with angles comparable to the actual animal incurred comparable forward drag but generated negative lift, allowing the robot to maintain approximately fixed depth beneath the surface (Maladen et al. 2011b).In the solid-like regime of locomotion, a robophysical model of the sea turtle with flexible flippers displayed robust performance over a range of limb intrusion parameters; the low slip intrusion scheme was analogous to that used by high-performing hatchling sea turtles (Mazouchova et al. 2013).When the robot used rigid limbs instead, it fluidized the surrounding material upon retraction of the limbs such that the robot sank into the material and incurred large belly drag.It was observed that hatchling sea turtles similarly avoided reinteracting with the disturbed ground.Robophysical experiments also revealed the importance of sidewinding in granular media (Astley et al. 2015), particularly on inclines.The multisegmented limbless robophysical models (Marvi et al. 2014) were able to locomote over a wide range of body contact fractions on level granular media, but as the slope angle increased, the fraction of body required for regular slope ascent increased while the range of contacts that led to successful climbing without slipping or pitching downhill narrowed.
When confronted with the above complexity and diversity of locomotor modes, one way to proceed is via DEM modeling with approximately spherical particles.Indeed, DEM has quantitatively captured the swimming performance of a lizard (Maladen et al. 2011a) and the running dynamics of a lightweight robot (Zhang et al. 2013) (e.g., Figure 6).However, these experiments were conducted with large-diameter particles (3-mm glass), which create unrealistic conditions.Thus, insight into locomotor tasks in more realistic granular materials by DEM would require orders of magnitude more grains than could possibly be simulated numerically.Seeing the benefits of a potentially reduced-order model, a scheme called resistive force theory (RFT) was hypothesized for granular locomotors (Zhang & Goldman 2014).RFT was initially developed in the mid-twentieth century by Gray & Hancock (1955) and Lighthill (1976) to model the locomotion of tiny organisms moving in viscous fluids, with the key feature that inertia is not important relative to dissipative forces.This leads to a feature as described by Purcell (1977): Time-symmetric self-deformation patterns produce no net displacement.Another major assumption in early formulations of RFT (e.g., Gray & Hancock 1955) is that local forces on the locomotor can be represented as elemental resistive forces produced by segments as if they translated alone through open fluids.These are typically solved assuming Stokes drag, and for long, slender elements, the ratio of the forces perpendicular to the surface to forces parallel (the so-called drag anisotropy) is ∼1.5 in free fluids but can range from 1.4 to 4.1 depending on conditions (larger values occur near solid boundaries) (Sznitman et al. 2010).The final assumption is that the net forward force on the locomotor is given by the summation of thrust minus drag from each elemental force.This scheme produces reasonable agreement with observations for locomotor performance of undulating flagella on nematode worms when wavelengths are long enough so that interaction between neighboring elements can be neglected.Given the Stokeslet form, these interactions can be long ranged, so RFT in fluids is often supplanted by so-called Stokeslet and slender body theories ( Johnson 1980, Cortez et al. 2005), which produce better agreement with experimental data at the cost of additional computational complexity; readers are referred to Rodenborn et al. (2013) for a comprehensive study of the different techniques.
To apply RFT to swimming in granular media, one needs a granular equivalent of Stokes laws to infer drag on elemental body units as a function of velocity, orientation, and angle of attack of each element.These can be found experimentally.Drag anisotropy was observed to be not a simple ratio but rather a function that was consistently greater than the drag ratio of low-Reynolds number fluids on slender objects (Maladen et al. 2009).Remarkably, when wave dynamics from the sand-swimmer were used as input, granular RFT found that the swimmer uses the optimal amplitude-to-wavelength ratio for forward motion and mechanical energy use.RFT was also successful in capturing the swimming of a robot in granular media (Maladen et al. 2011a), as well as the subsurface swimming of an elongate snake (Sharpe et al. 2015).The theory rationalized the wave shape on the body and agreed well with observations once the tangential forces were scaled to the appropriate grain-body friction.
Next, RFT was applied in legged locomotion.Elemental stresses on vertical descending and ascending plates were measured and used to model the locomotion of a legged robot with varying limb shapes (Li et al. 2013).Again, the key feature was that RFT assumed a linear superposition of elemental forces, but here it did not assume a force and torque balance on the robot so that inertia of the locomotor is important but material inertia is not.Remarkably, this works quantitatively as well over a variety of grain sizes and shapes (spherical glass to kidney-shaped poppy seeds).RFT began to break down with more complex-shaped and polydisperse grains, which are a frontier for future modeling.Since then, RFT has been applied to diverse situations of granular locomotion, including sidewinding (Rieser et al. 2021, Chong et al. 2022a), quadruped and hexapedal locomotion (Chong et al. 2021(Chong et al. , 2022a,b) ,b) and wheels (Agarwal et al. 2021, Shrivastava et al. 2020), and it has been used to study the neuromechanical control of sandfish (Ding et al. 2013b, Sharpe et al. 2013).
The original drag relations in RFT for low-Reynolds number swimmers were calculable from solutions of the Stokes equations.In a similar vein, one can infer elemental drag relations to be used in granular RFT given an appropriate continuum model for granular flow.As long as the intruding objects are much larger than the grain size, a local, tension-free frictional rheology has been shown in simulated drag tests to produce results that are very close to experimentally obtained data in cohesionless grains (Askari & Kamrin 2016).This constitutive model invokes only two material parameters: a constant friction coefficient µ s during dense flow, and a connected density ρ c , below which the material is deemed to be separated and stress-free.The intruding body may also have a surface friction µ w that differs from µ s .Hence, by dimensional analysis alone, on the continuum model, we can infer that the drag force on a buried plate element of size L × L oriented with normal n and moving slowly in direction v in a material at ambient pressure P must have the form F = L 2 P ψ( v, n; µ s , µ w ) for some dimensionless vector-valued function ψ.One reason RFT works so well in a granular media is the fact that this basic drag force relation gives perfect scaling with plate area-if the plate area quadruples (L doubles), the force quadruples.RFT would predict exactly this-doubling L makes a plate that looks like four copies of the original plate and thus RFT's independent superposition concept would say that the drag force should quadruple.This is actually an advantage over the original viscous fluid RFT, which does not have this perfect areal scaling and is thus known to be only approximate even in the simplest cases.
Granular RFT has proven to be powerful and is most studied in the regime of slightly polydisperse, dry, flat, nonhysteretic flows.This is a very small subset of the important regimes that locomotors encounter.Some progress has been made in modifying RFT to scenarios of rapid intrusion or interaction.In these situations, changes in the material state can then be incorporated in the RFT.For example, in accelerating intrusions, Aguilar & Goldman (2016) demonstrated that the RFT intrusion surface could be formed via the solidifying cone of media beneath a flat plate.In work studying wheel locomotion, Agarwal et al. (2021) demonstrated that the surface height behind a rapidly spinning wheel could be used to scale the stresses to obtain predictions and that material inertia could also be accounted for by quadratic terms in the RFT; these effects were dubbed dynamic RFT.
Despite the success in the wheel study, surface effects present the greatest challenge to RFT largely because of the memory of the material state.For example, slithering organisms and robots on the surface of sufficiently frictional sand encounter material that remains solidified in new structures [as observed by Schiebel et al. (2020)], which could presumably be accounted for by modified RFT coefficients.A path to modeling success was hinted at in studies of a different situation: wheel/legged locomotion on flat ground in which there was a spinning wheel, which allowed the theory to work since it broke the solid-like structures and allowed RFT to function.The situation becomes more challenging when locomotors ascend slopes.In their study of sidewinding, Marvi et al. (2014) discovered that the normal force (due to downhill plate drag) could be used via a scaled stress relationship (by angle of slope relative to avalanching angle for that volume fraction); here, however, the animals' careful control of the environment (minimal slipping) contributed to the success of the theory since material was not transported downhill.While the simple force scaling rationalized the forces generated by variable-width contacts, an improved RFT model would measure nondownhill forces on slopes, as well as account for interactions during downward-flowing avalanches, and could help rationalize the inability of most snake species to ascend sandy slopes, thereby providing clues to the evolutionary pressure that generated this style of locomotion.

SUMMARY POINTS
1. Predictive models for granular rheology that incorporate particle-size effects have now been introduced that account for a length scale directly within the rheology, and thereby capture more subtle behaviors that elude the µ(I) rheology and other mechanically local models.In particular, models based on a diffusing fluidization field can describe different manifestations of the nonlocal effect.
2. New models have been posed for predicting granular unmixing and, in particular, modeling grain-size-and mass-based segregation fluxes.These models construct segregation fluxes under the additive influence of gravity and gradients of different fields, such as concentration and partitions of the stress.
3. New particle-scale methods have been developed to represent nonspherical grain shapes and enable virtual experiments that, combined with multiscale and machine learningbased algorithms, reveal insights into how shape influences flow response.
4. Intrusion force modeling is an application area where relatively simple, reduced-order approaches have proven quite useful for generic grains.One such model, granular resistive force theory (RFT), has led to many new advances in modeling biological, robotic, and vehicular locomotion.Its key assumption is that the intrusion force on each small element of a complex body obeys a simple, localized formula.

FUTURE ISSUES
1. Is there a tractable way to merge the modeling ingredients discussed here into a universal model that can simultaneously predict flow and unmixing in granular media?This question becomes even more complex when including effects not discussed in this review, such as the transient role of evolving structural state variables.
2. If nonlocality is mediated by diffusing fluidization, what precise physical mechanism underlies this grain-size-dependent diffusion?The fluidization appears directly linked to grain fluctuations, so this question could relate to determining a precise heat equation for granular temperature.Other nonlocal mechanisms have been proposed but the degree of testing is not consistent-a comprehensive set of benchmark tests against which to compare proposed models should be determined.
3. What governs the apparent transition in segregation dynamics from collisional energetic flows to dense overdamped flows?A model framework that can capture these physics could be useful for a broader range of segregation problems.Perhaps the same mechanisms that give rise to the complicated rheology of dense flows-sensitive to the particle shape and size and possibly interstitial fluids (often non-Newtonian in nature)-also play a first-order role in dense segregation phenomena.
4. Why does RFT work so well in granular media-better in many cases than it does in viscous fluids, for which the theory was originally derived?The precision of granular RFT is derivable from plasticity concepts in simple cases, but it is not clear why it works so well for general geometries.
5. Can a predictive formula be developed that relates grain shape-and size-distribution properties to continuum rheology, or is the most realistic solution a model-free approach based purely on data?New methods would be needed to rapidly calibrate the dependence of rheology on the myriad material parameters.
Figure 1 Nonlocal effects causing spatially decaying creep.(a) Flat-bottom silo flows.(Top) Simulation of the µ(I) continuum model shows sharp cutoffs between zones undergoing shearing (color scale) and static zones (in black).(Bottom) Corresponding discrete element method simulation showing the zones of shearing diffusing out with no truly static zones anywhere.Panel adapted from Kamrin (2010) (MIT license).(b) Heap flows: Long-exposure photos of a heap flow with displayed exposure times showing that the boundary between what appears to be static media and flowing media descends as exposure time increases, indicative of a decaying slow creep, and that there is no truly static zone.Panel adapted with permission from Komatsu et al., Physical Review Letters, 86(9):1757, 2001; copyright 2001 American Physical Society.(c) The split-bottom Couette cell: This fully 3D flow geometry produces a wide shear zone on the free surface whose flow profile, when plotted as normalized rotation rate ω = v θ /V wall versus normalized position λ, collapses almost exactly onto an error function for many different choices of filling height, channel width W, and split radius R split .Panel adapted with permission from Fenistein et al., Physical Review Letters, 92(9):094301, 2004; copyright 2004 American Physical Society.
Figure 2 (a) Nonlocal effects causing thin-body strengthening.Experimental data showing that when layers are thick (large H/d, where d is the mean grain diameter) there is a clear individual angle of arrest; however, for thin layers, the thinner the layer, the steeper the angle θ at which flow arrests.Panel adapted from Pouliquen (1999) (MIT license).(b) Nonlocal effects causing secondary rheology.The Reddy et al. (2011) experimental setup, consisting of an annular shear cell and a force probe in the form of a submerged rod.When > 0, the rod can glide through the material for any F > 0, but otherwise it experiences a yield force.The dashed line marks the apparent boundary between the wall-induced flow zone and the quiescent zone.Panel adapted with permission from Reddy et al., Physical Review Letters, 106(10):108301, 2011; copyright 2011 American Physical Society.
Figure 3 Examples of segregation.(a) Fingering in pyroclastic flows.Panel adapted from Pouliquen & Vallance (1999) (MIT license).The two arrows indicate two fingers, and the helicopter provides a sense of scale.(b) Coarsening of the snout of a debris flow.Panel adapted from Bardou (2002).(c-e) Examples of somewhat complex patterns in rotating drums from (c,d) Hill et al. (1999) and (e) Hernández-Delfin et al. (2022).(c) Segregation pattern of steel (dark) and glass (light) spheres arising from a combination of chaotic advection in the square mixer combined with the more universally present segregation and diffusion.Panel adapted with permission from Hill et al. (1999); copyright National Academy of Sciences, USA.(d) Segregation pattern of small (blue) and large (white) glass spheres arising from a positive feedback mechanism driven by fill level of drum, independent of shape (e.g., Hill et al. 2004).Panel adapted with permission from Hill et al. (1999); copyright National Academy of Sciences, USA.(e) Illustrations of the particle-shape dependence of segregation behaviors.Figure adapted with permission from Hernández-Delfin et al. (2022), Physical Review Letters, 106(5):054614, 2022; copyright 2022 American Physical Society.Abbreviation: AR, aspect ratio.

Figure 4
Figure 4 Effect of particle shape on (a) the stress-strain response and (b) the plastic behavior of granular assemblies.Panel a adapted with permission from Athanassiadis et al. (2014).Panel b adapted from Murphy et al. (2019) (CC BY 4.0).
Figure 6 (a) Animals and robophysical models interacting with dry granular media.(i) The sandfish lizard (left), which swims beneath the surface of the sand, and its multisegmented model (right).(ii) The hatchling loggerhead sea turtle (left) and a robophysical model (right) using a passive flexible wrist.(iii) The sidewinder rattlesnake (left) and multisegmented model (right) (Marvi et al. 2014).(b) Approaches to model locomotion in and on granular media.(i) DEM simulation of the sandfish lizard showing the local surrounding frictional fluid.(ii) Characteristic RFT calculation scheme.(iii) The speed and energy use predicted during sandfish swimming; the gray region in the top plot and the PDF in the bottom plot indicate the region of the amplitude to wavelength ratio A/λ used by animals.(iv) Applying RFT in an inertial material regime by calculating coefficients on transient structures.Geometric and added-mass model of granular cone evolution during vertical flat foot intrusion (a plate of area A foot ).In subpanel iv C indicates a region (yellow) of virtual mass surrounding the solidified conical core (brown).The orange and blue arrows pictured on the growing cone sequence represent the stresses on the flat and rounded portions of the cone, each of which can be modeled accurately using RFT-based intrusion formulas from Aguilar & Goldman (2016).Abbreviations: CoT, cost of transport; DEM, discrete element method; PDF, probability density function; RFT, resistive force theory; STL, snout-tail length.Panel a, subpanel ii adapted from Mazouchova et al. (2013) (MIT license).Panel a, subpanel iii adapted with permission from Marvi H, Gong C, Gravish N, Astley H, Travers M, et al. 2014.Sidewinding with minimal slip: snake and robot ascent of sandy slopes, Science 346(6206):224-29; reprinted with permission from AAAS.Panel b, subpanel ii adapted with permission from Zhang & Goldman (2014).Panel b, subpanel iii adapted from Sharpe et al. (2013) with permission from the Journal of Experimental Biology.Panel b, subpanel iv adapted from Aguilar & Goldman (2016).