Force chains, friction, and flow: Behavior of

granular media across length scales

Thesis by

Ryan C. Hurley

In Partial Fulfillment of the Requirements

for the Degree of
Doctor of Philosophy

California Institute of Technology

Pasadena, California

2016
(Defended September 8, 2015)
 ii

c 2015
Ryan C. Hurley
All Rights Reserved
 iii

Acknowledgements

I would like to thank my advisor, José Andrade, for the many lessons I have learned
from him about being a successful faculty member. His excitement for science has
motivated me throughout my time at Caltech. I hope to put his many lessons to use
in the future. I would also like to thank the other members of my thesis committee –
Kaushik Bhattacharya, Guruswami Ravichandran, and Michael Lamb – for invaluable
conversations, feedback, and professional advice.
I would like to thank my fiancée, Margaret. Without her, I would be a different
person today and this thesis would likely not exist.
I would like to thank all past and present members of the Computational Geome-
chanics Group at Caltech with whom I have shared the ups and downs of research:
Keng-Wit Lim, Utkarsh Mital, Reid Kawamoto, Alex Jerves, Eloise Marteau, Daniel
Arthur, Ivan Vlahinic, and Jason Marshall. I would also like to thank Stephen Hall
of Lund University for allowing me to collaborate with him now and in the future. I
consider this a great honor.
Finally, I would like to acknowledge the various funding agencies and projects that
I have had the pleasure of providing research for: AFOSR (# FA9550-12-1-0091),
DTRA (HDTRA1-12-1-0041), JPL (“Technologies for Increased Landing Accuracy,
Payload Mass Fraction, and Landing Site Elevation”), and KISS (“xTerramechanics:
Characterization and Modeling of Spacecraft-Regolith Interactions”). These projects
and the great teams I have worked with have opened my eyes to new and exciting
frontiers in science and engineering.
 iv

Abstract

We study the behavior of granular materials at three length scales. At the small-
est length scale, the grain-scale, we study inter-particle forces and “force chains”.
Inter-particle forces are the natural building blocks of constitutive laws for granular
materials. Force chains are a key signature of the heterogeneity of granular systems.
Despite their fundamental importance for calibrating grain-scale numerical models
and elucidating constitutive laws, inter-particle forces have not been fully quantified
in natural granular materials. We present a numerical force inference technique for
determining inter-particle forces from experimental data and apply the technique to
two-dimensional and three-dimensional systems under quasi-static and dynamic load.
These experiments validate the technique and provide insight into the quasi-static
and dynamic behavior of granular materials.
At a larger length scale, the mesoscale, we study the emergent frictional behavior
of a collection of grains. Properties of granular materials at this intermediate scale
are crucial inputs for macro-scale continuum models. We derive friction laws for
granular materials at the mesoscale by applying averaging techniques to grain-scale
quantities. These laws portray the nature of steady-state frictional strength as a
competition between steady-state dilation and grain-scale dissipation rates. The laws
also directly link the rate of dilation to the non-steady-state frictional strength.
At the macro-scale, we investigate continuum modeling techniques capable of sim-
ulating the distinct solid-like, liquid-like, and gas-like behaviors exhibited by granular
materials in a single computational domain. We propose a Smoothed Particle Hydro-
dynamics (SPH) approach for granular materials with a viscoplastic constitutive law.
The constitutive law uses a rate-dependent and dilation-dependent friction law. We
 v
provide a theoretical basis for a dilation-dependent friction law using similar analysis
to that performed at the mesoscale. We provide several qualitative and quantitative
validations of the technique and discuss ongoing work aiming to couple the granular
flow with gas and fluid flows.
 vi

Contents

Acknowledgements iii

Abstract iv

1 Introduction 1
1.1 Granular materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Overview of thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Background on inter-particle forces . . . . . . . . . . . . . . . . . . . 5

2 Extracting inter-particle forces in opaque granular materials: Be-

yond photoelasticity1 9
2.1 Experimental validation and application of GEM . . . . . . . . . . . 11
2.1.1 Experimental setup and procedure . . . . . . . . . . . . . . . 11
2.1.2 Strain fields and force results . . . . . . . . . . . . . . . . . . 13
2.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.1 Equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.2 Average particle stress . . . . . . . . . . . . . . . . . . . . . . 19
2.2.3 Cohesionless and tangential forces . . . . . . . . . . . . . . . . 20
2.2.4 Structure of equations and connection to experiments . . . . . 22
2.3 Mathematical framework . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.1 Inverse problem formulation . . . . . . . . . . . . . . . . . . . 23
2.3.2 Measurement noise and alternative formulation . . . . . . . . 25
1
Adapted from R. Hurley, E. Marteau, G. Ravichandran, and J.E. Andrade. Extracting inter-
particle forces in opaque granular materials: Beyond photoelasticity. Journal of the Mechanics and
Physics of Solids, 63:154-166, 2014. DOI: 10.1016/j.jmps.2013.09.013
 vii
2.3.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4 Numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.1 Global force inference . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.2 Noise reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.3 Local force inference . . . . . . . . . . . . . . . . . . . . . . . 30
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 Dynamic inter-particle force inference in granular materials: Method

and application2 32
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 The method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2.1 Experimental measurements . . . . . . . . . . . . . . . . . . . 33
3.2.2 Numerical optimization . . . . . . . . . . . . . . . . . . . . . . 34
3.3 Example 1: Validation . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.1 Experiment setup . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.2 Experiment analysis . . . . . . . . . . . . . . . . . . . . . . . 44
3.3.3 Experiment results and FEM comparison . . . . . . . . . . . . 46
3.4 Example 2: Application . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.4.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . 49
3.4.2 Experiment analysis . . . . . . . . . . . . . . . . . . . . . . . 50
3.4.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . 50
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4 Friction in inertial granular flows: Competition between dilation and

grain-scale dissipation rates3 58
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2
Adapted from R.C. Hurley, K.W. Lim, G. Ravichandran, J.E. Andrade. Dynamic inter-particle
force inference in granular materials: Method and application. Experimental Mechanics, 2015. DOI:
10.1007/s11340-015-0063-8
3
Adapted from R.C. Hurley and J.E. Andrade. Friction in inertial granular flows: Competition
between dilation and grain-scale dissipation rates. Granular Matter, 17(3)287-295, 2015. DOI:
10.1007/s10035-015-0564-2
 viii
4.2 The friction law for simple shear . . . . . . . . . . . . . . . . . . . . . 61
4.3 Numerical simulations of simple shear . . . . . . . . . . . . . . . . . . 63
4.3.1 Description of code . . . . . . . . . . . . . . . . . . . . . . . . 64
4.3.2 Results: Validity of friction law . . . . . . . . . . . . . . . . . 67
4.3.3 Results: Grain-scale dissipation mechanisms . . . . . . . . . . 68
4.3.4 Influence of material properties . . . . . . . . . . . . . . . . . 69
4.3.5 Results: Dilation and dissipation rates . . . . . . . . . . . . . 71
4.3.6 Results: Rate-dependent friction . . . . . . . . . . . . . . . . 73
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

5 Modeling dilative viscoplastic granular flows using SPH4 77

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2 Balance and constitutive Laws . . . . . . . . . . . . . . . . . . . . . . 80
5.2.1 Balance law . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.2.2 Constitutive law . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.2.3 Theoretical link between friction and dilation rate . . . . . . . 82
5.3 SPH formulation and algorithm . . . . . . . . . . . . . . . . . . . . . 84
5.3.1 Kernel function and basic equations . . . . . . . . . . . . . . . 84
5.3.2 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.3.3 Equation of state . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.3.4 Equation of motion . . . . . . . . . . . . . . . . . . . . . . . . 87
5.3.5 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 88
5.3.6 Time integration . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.4 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.4.1 Example 1: Angle of repose . . . . . . . . . . . . . . . . . . . 92
5.4.2 Example 2: Infinite inclined plane flow . . . . . . . . . . . . . 96
5.4.3 Example 3: Granular column collapse . . . . . . . . . . . . . . 102
5.4.3.1 Scaling laws . . . . . . . . . . . . . . . . . . . . . . . 102
4
Adapted from R.C. Hurley and J.E. Andrade. Modeling dilative viscoplastic granular flows
using SPH. Under review.
 ix
5.4.3.2 Reduced bulk modulus, friction angle, and resolution
studies . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.4.4 Example 4: Column collapse down inclined planes . . . . . . . 108
5.5 Discussion and conclusion . . . . . . . . . . . . . . . . . . . . . . . . 115

6 Discussion and conclusion 117

6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
6.2 Ongoing and future work . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.2.1 Inter-particle force inference . . . . . . . . . . . . . . . . . . . 118
6.2.2 SPH modeling of coupled granular media and fluid flows . . . 119

Bibliography 124
 x

List of Figures

1.1 The three length scales studied in this thesis: the microscale, the mesoscale,
and the macroscale. The respective quantities of interest at each scale
include force chains, friction, and flow. . . . . . . . . . . . . . . . . . . 2

2.1 The GEM methodology for inter-particle force-inference presented in

this paper. Experimental imaging techniques provide rich data sets for
extracting intra-particle strain fields and material fabric. These ingredi-
ents are input into a mathematical framework which yields inter-particle
forces by solving an appropriate inverse problem. Variables and equa-
tions in panel 3 are described in section 2.2. . . . . . . . . . . . . . . . 10
2.2 Experimental setup: 1. CMOS camera ; 2. Canon lens ; 3. Digimeter ;
4. Loading device. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 (a) Sample under macroscopic loading. Top face (red) was prescribed a
vertical load of 215 N using a smooth, rigid wall. Bottom and lateral
faces (blue) were smooth, stationary walls. (b) Results of segmentation
process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4 Sample dimensions in millimeters and information on the fabric: particle
numbers, position of contact points and centroids. . . . . . . . . . . . . 15
2.5 Strain components xx , yy and xy obtained from DIC. . . . . . . . . . 16
2.6 Resulting inter-particle forces (in N) from application of GEM. . . . . 17
2.7 Illustration of particle-particle and particle-boundary contacts. . . . . . 18
2.8 Unit vectors for contact i. . . . . . . . . . . . . . . . . . . . . . . . . . 21
 xi
2.9 (a) Numerical odometric test setup. (b) Inter-particle forces computed
with DEM and (2.13). Length scale in meters, forces in Newtons. Line
thickness proportional to force magnitude. . . . . . . . . . . . . . . . . 28
2.10 (a) Solution to (2.13) with artificial noise, kδf k = 0.154. (b) Solution
to (2.14) with same noise, kδf k = 0.053 (right). . . . . . . . . . . . . . 29
2.11 Solution to (2.13) performed locally within granular material. . . . . . 30

3.1 Illustration of various terms used in the governing equations. . . . . . 35

3.2 Pareto front. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3 Experimental materials and setup. (a) Table setup with air duct con-
nected to fan, rigid blocks fastened to table, camera and light source ver-
tically over the area of interest. (b) Close-up view of table top and the
top surface of a polyurethane disk with speckle pattern. (c) The under-
side of a polyurethane disk, showing plenum chamber for compressed air.
(d) Experimental setup with 6 polyurethane disks and rigid impactor
propelled at initial conditions vx = −1.141m/s and vy = −0.66m/s.
Particle labels will be used later and are in no meaningful order. . . . 44
3.4 The xy component of strain from (a) DIC and (b) FEM at four times
during the impact event. All figures share a common scale. Strain values
above and below the extreme values of the scale occur in small areas. 46
3.5 (a) Forces found using Eq. (3.18). Length and width of lines are propor-
tional to force magnitude. (b) Forces found from the FEM model. (c)
Eight contact point locations, chosen to include the major force chain
and otherwise at random. (d) A comparison of force evolutions using
Eq. (3.18) and the FEM model. . . . . . . . . . . . . . . . . . . . . . 48
3.6 Experimental setup. (a) A window featuring a piece of Plexiglas and
a sheet of PVC plastic holding the 2D granular bed in place. A high
speed camera is used to view the granular bed as it is impacted by a
steel intruder. (b) A close-up of the steel intruder about to impact the
granular bed at -2.4m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . 51
 xii
3.7 Inter-particle forces in the granular bed at six times during an impact
event. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.8 Inter-particle forces in the granular bed at six times during an impact
event with an initial packing of the bed different from the experiment
shown in Fig. 3.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.9 The vertical force felt by the intruder for the experiments shown in Figs.
3.7 and 3.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.10 Maximum inter-particle force as a function of distance from the surface
of the intruder for three experiments with distinct initial packings. Exp.
1 refers to Fig. 3.7. Exp. 2 refers to Fig. 3.8. An inter-particle force
is considered to be a distance ndmin from the surface of the intruder if
the contact point that it belongs to falls within a range of ndmin and
(n + 1)dmin away from the surface of the intruder, where dmin is the
minimum particle diameter in the granular bed, 2.54cm. . . . . . . . . 55

4.1 Friction encodes contact-scale and grain-scale information in a single

parameter for continuum analysis. . . . . . . . . . . . . . . . . . . . . 59
4.2 The two timescales associated with the inertial number, I = Tc /Tγ̇ . . . 59
4.3 (a) A rendering of the simple shear flows featured in this paper. The top-
most and bottom-most particles are used as rough boundaries. Colors
indicate the magnitude of velocity, where vx is the imposed wall velocity.
(b) Coefficient of restitution e in a two-particle collision with normal
p
velocity vc ρg /kn for data set 1 (- -) and data set 2 (-). . . . . . . . . 65
4.4 A comparison of (a) effective friction and (b) solid fraction from our sim-
ulations and available data sets taken from the literature. Blue squares
are from contact dynamics simulations [1]. Black triangles are from 3D
annular shear cell experiments [2, 3]. Error bars indicate standard de-
viations in time of measured quantities. The dashed line shows the fit
µ = µ0 + (µ1 − µ0 )/(1 + I0 /I) from [4]. . . . . . . . . . . . . . . . . . 67
 xiii
4.5 A comparison of the effective friction coefficient calculated from pro-
posed friction relationship in Eq. (4.6) and the stress formula in Eq.
(4.7). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.6 (a) The two terms µn = ZφhΓn i/(I Γ̃) and µs = ZφhΓs i/(I Γ̃) in the
additive decomposition of effective friction given in Eq. (4.9) for the
primary data set (- -) labeled DS1 and the secondary data set (-) labeled
DS2. (b) The total effective friction for the two data sets as a function
of I. Error bars are omitted from inset plots for clarity. . . . . . . . . . 70
4.7 (a) The coordination number Z as a function of inertial number. The
dashed line is the fit from Eq. (4.10). (b) The solid fraction φ as a
function of inertial number. The dashed line is the fit from Eq. (4.11).
(c) The average grain-scale dissipation rates hΓn i/Γ̃ and hΓs i/Γ̃ as a
function of inertial number. The dashed lines are power-law fits, with
hΓn i/Γ̃ ∝ I 2.4 and hΓs i/Γ̃ ∝ I 1.87 (d) The quantity Zφ/I as a function of
inertial number. The dashed lines represent the two regimes of behavior
in which Zφ/I ∝ I −1 and Zφ/I ∝ I −2 . Error bars in (d) are negligibly
small. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.8 The inset of Fig. 4.7c, showing the average grain-scale dissipation rates
as a function of inertial number on a log-log scale. The dashed lines
illustrate the scaling of each dissipation rate discussed in the text, as
well as the scaling proportional to I at the transition to quasi-static flow. 74

5.1 (a) Illustration of the smoothing kernel W with compact support of ra-
dius 2h about particle a. (b) Interior particles interacting with boundary
particles, which are given artificial velocity. . . . . . . . . . . . . . . . 86
 xiv
5.2 (a) Initial conditions of the slump test through an orifice. In color, blue
represents internal SPH particles and red represents boundary particles.
Length l = 60 cm, h = 50 cm, and the size of the simulated material in
the z dimension is 6 cm. (b) Final collapsed side profile of the slump
test through an orifice for φ = 30◦ . The repose angle θ is measured in
the middle 75% of the domain between the left wall and the orifice. (c)
Results of the simulations with and without including β in the calcula-
tion of friction. The dashed line and circles represents the curve θ = φ.
Results without β are blue squares which provide a lower bound to the
θ = φ curve. Results with β are red inverted triangles which nearly
match the θ = φ curve. (d) Pressure in the collapsed granular pile for
φ = 30◦ showing a slight pressure drop in the layer of SPH particles
directly above the bottom surface of the container. . . . . . . . . . . . 94
5.3 (a) Initial conditions of a typical inclined plane flow experiment showing
internal SPH particles above y = 0 and boundary particles at or below
y = 0. (b) A typical velocity profile of an inclined plane flow. (c) Veloc-
itiy profiles in the x direction as a function of height y in the granular
material for simulations that do not include β in the calculation of fric-
tion. Symbols represent simulation results and dashed lines represent
the best fit of the Bagnold profile (Eq. (5.30)) to the simulation data.
(d) Same as (c) but for simulations that include β in the calculation of
friction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
 xv
5.4 (a) Bulk acceleration ā˙ of the granular material in the inclined plane
simulations, without the use of β in the calculation of friction, as a
function of time. The highest curve represents θ = 33◦ and each lower
curve represents a reduction in θ by 3◦ . All accelerations appear to
asymptote to zero except for θ = 33◦ . The inset shows raw data overlaid
with the same fits. (b) Same as (a) except for simulations that include
β in the calculation of friction. (c) Slip velocity in simulations not
including β in the calculation of friction, found by averaging the vx
velocity of the bottom-most layer of interior SPH particles. All slip
velocities appear to asymptote to finite values except for θ = 33◦ . (d)
Same as (c) but for simulations that include β in the calculation of friction.101
5.5 Typical initial conditions of the column collapse simulations for (a) grit
and (b) fine glass. The right-most vertical wall in each figure is the
containing wall that is deleted at t = 0. In the color figure, blue particles
represent internal particles and red particles represent boundary particles.104
5.6 (a) Scaling laws for collapsed height in simulations not using β in the
calculation of friction. Symbols are SPH simulations and dashed lines
are linear best-fit lines in logarithmic space for all data with x coordinate
above 1.7. (b) Scaling laws for collapsed runout distance in simulations
not using β in the calculation of friction. (c) Same as (a) but for sim-
ulations using β in the calculation of friction. (d) Same as (b) but for
simulations using β in the calculation of friction. . . . . . . . . . . . . 106
 xvi
5.7 Results of parametric studies of various model parameters on numer-
ical results. (a) Result of varying the bulk modulus over two orders
of magnitude does not indicate any significant change in results. This
also justifies our use of κ = 105 in most simulations in this chapter to
eliminate the shorter time step that may otherwise be needed to ensure
stability. (b) Result of varying bulk and basal friction coefficients within
the experimental error bars of [5] illustrates the expected changes, but
only minor changes. (d) Results of mesh-refinement study illustrate
the coarse mesh used throughout the chapter produces results that are
nearly the same as the finest mesh. . . . . . . . . . . . . . . . . . . . . 109
5.8 (a) Initial and (b) final configurations for simulations of column collapse
down inclined planes. The right-most near-vertical wall in (a) represents
the containing wall that is lifted in the y direction with vx = 2 m/s at
t = 0. In the color figure, blue particles represent interior particles and
red particles represent boundary particles. . . . . . . . . . . . . . . . . 110
5.9 Time-dependent collapse profiles down inclined planes for four inclina-
tion angles: θ = 0◦ , 10◦ , 16◦ , and 22◦ . Symbols represent SPH profiles
and lines represent experimental profiles from [6]. (a) and (b) illus-
trate results with periodic boundary conditions in the z direction and
no sidewalls. (c) and (d) illustrate results with sidewalls. (a) and (c)
show results for simulations not including β in the calculation of friction.
(b) and (d) show results for simulations including β in the calculation
of friction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
 xvii
5.10 Snapshots of simulations shown in Fig. 5.9d for (a) θ = 0◦ , (b) θ =
10◦ , (c) θ = 16◦ and (d) θ = 22◦ . These simulations employ β in the
calculation of friction. In the color figure, the color of each particle
corresponds to the magnitude of the strain rate tensor D. We clearly
see the influence of sidewalls (not rendered) at the first nonzero time
for each inclination angle. At these times, particles appear to undergo a
higher strain rate in the foreground and background than in the center
of the flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

6.1 (a) Macroscopic load curve for oedeometric compression of single-crystal

quartz grains. (b) Sample volumetric image of the assembly obtained
with XRCT during the load cycle. (c) Preliminary results of force in-
ference using Eq. (3.19) applied to a single load step. Lines connect
contacting grain centroids and are darkened and thickened linearly with
force magnitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.2 Qualitative example of gas-driven soil erosion modeled using the SPH
framework. The left column of figures illustrates the sand response only
in the first 75 ms of the gas impingement event. The right column of
images illustrates the gas response only. . . . . . . . . . . . . . . . . . 123
 1

Chapter 1

Introduction

1.1 Granular materials

Granular materials are materials composed of discrete macroscopic particles with size
greater than 1 micron [7]. Snow flakes, sand, soil, rice, pharmaceutical pills, and even
asteroids can be considered granular. These materials are ubiquitous in our universe
and their behavior plays a crucial role in natural and engineered processes such as ero-
sion, mining, pharmaceutical processing, material transport, and food production [8].
Unlike gases and liquids, the mechanical behavior of granular materials is governed
by nonlinear inter-particle interactions rather than thermal fluctuations [7]. Further-
more, the lack of scale separation in these materials for many problems of interest
prohibits development of a single governing mathematical framework (e.g., Naiver-
Stokes) [7, 8, 9]. Thus, despite their ubiquity, granular materials possess a complexity
that necessitates a unique multi-scale approach.
This thesis examines the behavior of granular materials at three length scales: the
grain-scale, the mesoscale, and the macro-scale. These three scales and the behav-
iors of interest that will be studied are illustrated graphically in Fig. 1.1. At the
grain-scale, we study inter-particle forces and “force chains”. These unique features
are natural building blocks of constitutive laws for granular materials and frequently
manifest themselves in macroscopically observable phenomena. At the mesoscale, we
study emergent frictional behavior of a collection of flowing grains. The effective
frictional strength of granular materials plays a governing role in mechanical and
 2
LENGTH SCALE
μm mm m

yy
yx

GRAIN-SCALE MESOSCALE MACRO-SCALE

Ex: Force chains Ex: Friction Ex: Flows

Figure 1.1: The three length scales studied in this thesis: the microscale, the
mesoscale, and the macroscale. The respective quantities of interest at each scale
include force chains, friction, and flow.

dynamical behavior, both in quasi-static and dynamic events. At the macro-scale,

we investigate continuum modeling techniques for capturing the distinct solid-like,
fluid-like, and gas-like behaviors exhibited by granular materials. Full-field numer-
ical methods capable of capturing these phases are essential for predicting and un-
derstanding the processes underlying the behavior of granular media in nature and
industry.
The three scales and topics discussed in this thesis are connected in a manner that
highlights the multi-scale nature of granular materials. For instance, the mesoscale
effective friction coefficient and the macro-scale stress tensor can be derived from inter-
particle forces and grain kinematics (e.g., see [1,10,11,12]). Force chains arising from
grain-scale interactions may also induce nonlocal behavior that is widely observed at
meso- and macro-scales in number of flow configuration (e.g., see [90]). An improved
understanding of force chain dynamics and nonlocal constitutive theories like those
discussed in [90] may therefore provide new insight into the multi-scale physics of
granular media. Furthermore, mesoscale properties such as effective friction play
an important role in solid-like, fluid-like, and gas-like behavior and transitions at
the macro-scale. An improved understanding of these properties and new full-field
 3
modeling tools may help elucidate macro-scale behavior in insightful ways.
Despite the stated connections between the scales studied in this thesis, we do
not focus heavily on investigating these connections. Rather, we primarily develop
technologies for studying each scale and leave extensive investigation of connections
for future work. The lack of rigorous theories connecting scales is an ongoing challenge
in the field of granular materials. We hope that the technologies advanced in the
present work aid researchers in overcoming this challenge in the future.
Throughout this thesis, we do not focus on a particular type of granular material.
Instead, many of the results are generally applicable to any conglomeration of par-
ticulates. This practice of seeking results generalizable to many particulate systems
when studying granular media has a history of elucidating the behavior of many sim-
ilar systems: foams, colloids, emulsions, and suspensions [9]. However, we often make
reference to sands and cohesionless soils in an attempt to understand these most basic
yet fundamental examples of granular media.

1.2 Overview of thesis

The layout of this thesis is as follows. Chapters 2, 3, and 4 are adapted from published
journal articles. Chapter 5 is adapted from a journal article currently undergoing
peer-review. Most of the content of these four chapters is reproduced directly (with
permission of the publishers) in order to preserve their stand-alone quality. The
titles for these chapters are also the same as the journal articles they represent.
However, the introduction sections of chapters 2 and 3 are modified heavily. Rather
than presenting the original article introductions in those chapters, we consolidate
background information in section 1.3 to eliminate redundancies. We then dedicate
the introduction sections of these chapters to briefly outlining the research they will
describe. The introduction sections of chapters 4 and 5 are mostly unmodified to
provide a more seamless reading experience. However, we provide a brief overview of
these chapters, as well as chapters 2 and 3, in the remainder of section 1.2.
Chapter 2, adapted from a paper appearing in the Journal of the Mechanics and
 4
Physics of Solids [13], describes a technique for inferring inter-particle forces in quasi-
static experiments on granular materials. The force inference technique involves nu-
merically solving an optimization problem using experimental measurements of grain
positions, grain contact points, and volume-averaged grain stresses as input. The
technique is applicable to opaque, arbitrarily shaped grains, and makes no assump-
tion of contact or constitutive law for the grain material. However, the technique
requires that grains are deformable and that deformations can be measured experi-
mentally. We present the numerical framework, discuss its solution space, and provide
two examples of application to real and synthetic experimental data.
Chapter 3, adapted from a paper appearing in Experimental Mechanics [14], ex-
tends the force inference technique introduced in chapter 2 to dynamic experiments.
The extended technique involves augmented governing equations to accommodate
dynamic effects and a modified multi-objective optimization problem. We validate
this dynamic technique by comparing results of force inference applied to an impact
experiment with an Abaqus/Explicit simulation of the same boundary value problem.
We also present a study of force chains in impacted granular beds that illustrates the
link between force chains and intruder dynamics. The study also highlights an appar-
ent limit to force chain length that occurs at identical impact velocities for different
initial bed packings.
Chapter 4, adapted from a paper appearing in Granular Matter [10], provides
a derivation of a new friction law for steady-state granular flows. This friction law
analytically relates friction to shear rate, coordination number, porosity, and average
grain-scale dissipation rates. The friction law illustrates that rate-dependent friction
in steady-state flows arises from a delicate competition between material dilation and
grain-scale dissipation rates. We perform numerical simulations to investigate this
competition and illustrate how different grain-scale dissipation mechanisms are re-
sponsible for the macroscopic friction coefficient in granular shear flows. The reader
is directed to an additional paper appearing in Procedia Engineering [15] that ex-
tends this friction law to the non-steady-state regime and briefly discusses numerical
simulations.
 5
Chapter 5, adapted from a journal manuscript currently undergoing peer-review,
discusses a continuum modeling technique for granular flows. The technique combines
a mesh-free Smoothed Particle Hydrodynamics (SPH) framework with a viscoplastic
constitutive law that contains a rate- and dilation-dependent friction law. The friction
law is derived in a similar manner to those found in chapter 4, but from the governing
equations of continuum mechanics. Several examples illustrate the ability of this
method to reproduce realistic angles of repose in slump tests, to match the Bagnold
velocity profile in infinite inclined plane flows, to produce accurate scaling laws for
runout and final height of granular column collapses on flat surfaces, and to generate
time-dependent profiles of granular columns collapsing down inclined planes that
agree with experimental results.
Chapter 6 offers concluding remarks and briefly highlights ongoing and future
work related to the topics covered in each chapter. In particular, we briefly discuss
ongoing force inference experiments on single-crystal quartz grains imaged using x-ray
computed tomography and x-ray diffraction. These experiments are a continuation of
the work in chapters 2 and 3 and represent an important step toward understanding
force transmission in real materials and calibrating more complex grain-scale models.
We also discuss an ongoing effort to extend the SPH framework discussed in chapter
5 to modeling granular materials coupled with gas and fluid flows.

1.3 Background on inter-particle forces

Many studies have investigated analytical or empirical links between grain-scale and
macro-scale properties of granular media. For instance, researchers have experimen-
tally linked the so-called “fabric” of granular materials (i.e. the directional distribu-
tion of grain contact points) to the mobilized stress ratio and secant modulus [16].
Other researchers have developed methods of coarse-graining inter-particle forces and
contact point locations for the macro-scale stress tensor [17]. Several researchers have
derived analytical relationships between the mobilized angle of friction of a granular
material and the radial distribution of inter-particle contacts and forces [1, 12, 18].
 6
The length and dynamics of force chains has been shown to play an important role in
macro-scale behaviors such as intruder penetration [43]. Force chains may also play an
important role in inducing long-range communication in granular media, a behavior
captured by recent nonlocal constitutive laws (e.g., see [90]). Several studies have ex-
amined the collective statistical properties and heterogeneity of inter-particle forces in
quasi-static and dynamic settings [11,19,20,21,22,23,24,24,25,26,27,28,29,30,31,32]).
These studies have proven useful in unlocking the tools of statistical physics for gran-
ular media, although significant progress remains to be made before a unifying frame-
work exists for granular materials as it does for gases. One important milestone in
developing this framework and understanding granular media has been the advent of
grain-scale numerical techniques.
Grain-scale numerical models provide a platform for simulating the behavior of
individual particles in a granular material. The Discrete Element Method (DEM) rep-
resents the seminal work on grain-scale numerical modeling [33]. DEM uses Newton’s
second law and linear or Hertzian contact laws to explicitly integrate the equations
of motion of individual particles. Contact dynamics has also been introduced as an
alternative implicit approach [34,35]. In the past two decades, numerous studies have
employed these modeling techniques to further investigate links between grain-scale
features such as forces and macro-scale behavior (e.g., see [27, 36, 37]). As these tech-
niques evolve to handle complex grain shapes (e.g., see [38, 39]) and are increasingly
exploited to study inter-particle forces and other grain-scale behavior, it is important
to validate their results using experimental data.
Experimental measurements of inter-particle forces in model granular materials
has historically been accomplished using photoelasticity [21, 40, 41]). While photoe-
lasticity has proven immensely useful, it has several restrictions. First, photoelasticity
is limited to quasi-two-dimensional experiments on birefringent grains. Second, in-
version techniques used to determine inter-particle forces from photoelasticity data
require grains to remain elastic. Finally, photoelasticity poorly resolves the dynamic
evolution of inter-particle forces at the grain-scale, instead providing forces averaged
over tens or hundreds of grains (e.g., see [42, 43]).
 7
To overcome the elasticity restriction of photoelasticity, some researchers have
laboriously reconstructed individual contact forces by matching experimentally ob-
served fringe patterns to those generated from explicit simulations (see [44]). How-
ever, this technique is cumbersome and unrealistic for packings of more than a few
grains. To overcome the two-dimensional restriction of photoelasticity, several re-
searchers have studied force transmission in liquid droplets and gel packings (e.g.,
see [45, 46, 47]). However, these studies assume a Hertz contact law a priori and
therefore cannot validate inter-particle forces found in grain-scale numerical models
that also employ a Hertz contact formulation.
In chapters 2 and 3 of this thesis, we discuss a new force inference technique that
overcomes the limitations of methods employed in the past. In particular, the new
force inference technique is applicable to two- and three-dimensional particle packings,
is not restricted to material linearity or elasticity, makes no assumption of contact
law, and provides individual inter-particle forces during both quasi-static and dynamic
events. The basic equations for this force inference technique were originally derived
for quasi-static particle packings in [48]. However, [48] employed an incorrect nu-
merical technique for obtaining inter-particle forces and did not provide experimental
examples. Chapter 2 (adapted from [13]) provides a new presentation of the governing
equations and a new numerical inverse problem relying on a constrained optimiza-
tion approach for obtaining inter-particle forces. This chapter also includes the first
experimental example of applying the force inference technique to an indeterminate
packing of rubber cylinders. Chapter 3 (adapted from [14]) provides a new derivation
of the governing equations for dynamic events and a modified multi-objective opti-
mization approach for obtaining inter-particle forces. The dynamic form is validated
experimentally and used to study inter-particle forces beneath an intruder impacting
a granular bed. Chapters 4 and 5 address different topics, as discussed above and in
their respective introduction sections.
In chapter 6, we provide a brief discussion of ongoing and future work on the force
inference technique. The ongoing work includes inferring forces in opaque three-
dimensional single-crystal quartz grains imaged using combined x-ray computed to-
 8
mography and x-ray diffraction. This work represents the first time inter-particle
forces have been inferred within the bulk of an opaque indeterminately-packed sand-
like granular material.
 9

Chapter 2

Extracting inter-particle forces in

opaque granular materials: Beyond
photoelasticity1

This chapter presents the first example of inter-particle force inference in real gran-
ular materials using an improved version of the methodology known as the Granu-
lar Element Method (GEM). GEM combines experimental imaging techniques with
equations governing particle behavior to allow force inference in cohesionless materials
with grains of arbitrary shape, texture, and opacity. This novel capability serves as a
useful tool for experimentally characterizing granular materials, and provides a new
means for investigating force networks. In addition to an experimental example, this
chapter presents a precise mathematical formulation of the inverse problem involving
the governing equations and illustrates solution strategies.
The Granular Element Method (GEM) presented in this chapter is an extended
form of the method originally proposed in [48]. The contribution of this chapter is
to present an experimental validation of GEM and a precise mathematical formula-
tion of the proposed inverse problem involving the governing equations for particle
mechanics. This paper also provides an additional inverse problem formulation that
may help practitioners reduce solution error when experimental noise is present.
The GEM methodology can be visualized in Figure 2.1. Experimental imaging
1
Adapted from R. Hurley, E. Marteau, G. Ravichandran, and J.E. Andrade. Extracting inter-
particle forces in opaque granular materials: Beyond photoelasticity. Journal of the Mechanics and
Physics of Solids, 63:154-166, 2014. DOI: 10.1016/j.jmps.2013.09.013
 10
techniques such as high-resolution photography, 3D X-ray diffraction [49, 50, 51], and
X-ray computed tomography (XRCT) [52, 53, 54] provide rich data sets from which
Digital Image Correction (DIC) [55], level-set methods [56], and other techniques can
extract intra-particle strain fields and material fabric (contact locations and normals).
Intra-particle strains fields and material fabric are used as input into a mathematical
framework which yields intra-particle stress fields using an appropriate constitutive
relation and numerically solves an inverse problem using the governing equations of
particle statics. The result of the inverse problem is inter-particle forces, a key com-
ponent in the development of theories regarding granular material behavior. Because
GEM can be coupled with any experimental techniques capable of extracting strain
fields and fabric, recent advances in imaging (e.g., [49, 50]) will soon allow inference
of inter-particle forces in natural materials with small grains, such as sands.

Figure 2.1: The GEM methodology for inter-particle force-inference presented in

this paper. Experimental imaging techniques provide rich data sets for extracting
intra-particle strain fields and material fabric. These ingredients are input into a
mathematical framework which yields inter-particle forces by solving an appropriate
inverse problem. Variables and equations in panel 3 are described in section 2.2.

The layout of this chapter is as follows. Section 2.1 presents an example of the
 11
methodology applied to a real experiment involving rubber particles, showcasing the
applicability of the method to real materials of any texture and its adaptability to a
number of experimental imaging and algorithmic data analysis techniques. Sections
2.2 and 2.3 detail the ingredients of the mathematical framework of GEM. Section
2.4 illustrates additional features of GEM with a numerical example. Finally, section
2.5 offers concluding remarks.
The notation and terminology used in this chapter are defined as follows: Rm×n
denotes the set of real matrices with m rows and n columns, “·” denotes an inner
product (e.g., a · b = ai bi ; c · d = cij djk ), “⊗” denotes a dyadic product (e.g., a ⊗ b =
p
ai bj ), k · k2 denotes the Euclidian norm of a vector (e.g., kak2 = a21 + ... + a2n ), and
“sym” is the symmetric operator defined as sym(·) = 1/2((·) + (·)T ). R(A) refers to
the range of A, the set of all possible linear combinations of the columns of A. N (A)
refers to the nullspace of A, the set of all vectors z such that Az = 0.

2.1 Experimental validation and application of GEM

This section presents an example of GEM applied to a real experiment. The current
experiment uses Digital Image Correlation (DIC) to measure full-field particle strains
from which stresses are deduced. Information about the location of contact points and
particle shapes is obtained from digital images using segmentation algorithms. This
data is then passed to the GEM algorithm, which reconstructs the force distribution.
The results of this experiment showcase the applicability of GEM to real materials
and its potential to be combined with advanced experimental techniques.

2.1.1 Experimental setup and procedure

The experimental setup is shown in Figure 2.2. A CMOS camera with a 3.0 megapixel
sensor (PL-B623, PixeLINK, Ottawa, Canada) and a Canon lens was used at a work-
ing distance of 20 cm to image the assembly. A specifically designed loading device
is used to apply axial compression on the specimen, as shown in Figure 2.2. The
specimen is placed between four faces: the bottom and lateral faces are stationary
 12
and the top face can move vertically. The applied force was carefully measured with
a 500 g load cell (LCFA-50G, Omega, Stamford, CT) and monitored with a digimeter
(MD-40, Newport, Irvine, CA).

Figure 2.2: Experimental setup: 1. CMOS camera ; 2. Canon lens ; 3. Digimeter ;

4. Loading device.

The specimen used in the experiment was composed of rubber cylindrical grains.
Grain diameters were 7 mm, 10 mm, or 14 mm, and the grains’ out-of-plane length
was 20 mm. The grains’ Young’s modulus was 5.5 MPa and Poisson’s ratio was
approximately 0.5. The specimen shown in Figure 2.3 was compressed with a 215 N
vertical load on the top wall, while the side walls were held rigid.
The software VIC-2D was used to perform Digital Image Correlation (DIC) on
the images in order to determine in plane full-field strain [57, 58]. Digital Image
Correlation (DIC) is an optical tool based on digital image processing and numerical
computing which provides full-field displacements and strains by comparing the gray
intensity changes of the object surface before and after deformation [55,59]. The DIC
 13
procedure consists of tracking the same pixels between reference and deformed images.
To perform this tracking, a correlation window, or subset, is chosen and deformed
until the pattern in the deformed image matches the pattern in the reference image
as closely as possible. To determine an adequate subset size, a compromise between
resolution and measurement error needs to be found. The measurement error is
evaluated by correlating two subsequent images of the specimen without applying
any deformation. The values of the resulting strain components (xx , yy and xy ) for
different subset sizes are then compared in order to identify suitable configurations.
The degree of similarity between the reference and deformed subsets is computed using
a correlation coefficient and the best fit is achieved when the correlation coefficient
reaches its maximum. The position of the deformed subset is determined and the
in-plane displacement is obtained by calculating the difference between the position
of each point in the reference subset and the position of the corresponding point in
the deformed subset. The strain field is then computed by numerical differentiation
of the displacement field.
Segmentation algorithms were used to extract the material’s fabric: particle con-
tact points, centroids and areas. In particular, a circular Hough transform [60,61] was
first performed to determine the number of particles Np and approximate positions
of centroids and maximum radii, followed by a snake, or active contour model algo-
rithm [62,63] to determine the true contours of the grain. The different parameters of
the Hough transform and active contour model algorithm are manually adjusted by
visual inspection of the resulting segmentation of the grains. Segmentation is finally
achieved by partitioning the digital images into sets of pixels, each set constituting
different particles. The result of this segmentation procedure is shown in Figure 2.3.

2.1.2 Strain fields and force results

The contact locations, particle areas, and average particle strain were used as input
into GEM’s mathematical framework to be described in Sections 2.2 and 2.2.
A simple algorithm was used to determine the position of the contact points.
 14

Figure 2.3: (a) Sample under macroscopic loading. Top face (red) was prescribed a
vertical load of 215 N using a smooth, rigid wall. Bottom and lateral faces (blue)
were smooth, stationary walls. (b) Results of segmentation process.

First, grain boundaries were delineated as shown in Figure 2.3. Next, the euclidean
distance between pixels on two adjacent boundaries was calculated and the two pixels
were considered to belong to a contact surface if this distance was less than or equal
to a certain pixel threshold Lpix . Finally, the average value of positions of all pixels
belonging to a contact surface was used to determine a single coordinate (x, y) of the
true contact point. The results of this algorithm are shown in Figure 2.4, where it
can be seen that the 34 particles share a total of 78 contact points.
Figure 2.5 shows the full-field strain distribution obtained from DIC. The DIC
procedure was applied grain by grain. For each grain, the subset size was manually
chosen in order to achieve a reliable correlation analysis. The full-field strain distri-
bution was obtained by calculating values of xx , yy and xy inside of each particle
 15

58

50

Figure 2.4: Sample dimensions in millimeters and information on the fabric: particle
numbers, position of contact points and centroids.

using numerical differentiation of the displacement field. The average strain tensor,
¯p , was obtained by taking the arithmetic average of each tensor computed within a
given particle. The average stress σ̄ p of a particle was deduced from the average strain
¯p using generalized Hooke’s law, σ̄ p = c : ¯p , where σ̄ p is the average elastic stress
tensor at a particle and c is the elastic stiffness tensor of the particle. The Young’s
modulus and Poisson’s ratio used in the elastic stiffness tensor were determined using
a separate experiment on a single grain.
Using average particle strain and and fabric information, the regularized inverse
problem in equation (2.14) was solved with a regularization parameter of λ = 0.0024.
The resulting inter-particle forces are shown in Figure 2.6. It should be emphasized
that the presented forces are forces distributed over the length of the cylinder. Fig-
ure 2.6 offers the first look at inter-particle forces inferred in a real, opaque material
using GEM. The forces form expected patterns of force chains throughout the mate-
rial. The particle-boundary forces along the top platen were measured by the load
 16

x
0

0.025

0.045

0.035
0.035
0.025

0.045
✏xx ✏yy ✏xy

Figure 2.5: Strain components xx , yy and xy obtained from DIC.

cell to be 215 N and the calculated resultant forces sum to 215 N.

The average macroscopic stress in the array was also computed using the expres-
sion [17]:
1X c N
hσi = sym(f α ⊗ dα ) (2.1)
Ω α=1

where Ω is the volume of the assembly, Nc is the number of inter-particle contacts in

the assembly, hσi is the macroscopic average stress, f α is the inter-particle contact
force, and dα is the branch vector between particles. The resulting macroscopic stress
was computed to be:  
−2.56 0.26
hσi =   (kPa) (2.2)
0.26 −4.42

where compressional stresses are negative. The hσi22 component of this result cor-
responds well with the applied stress of −215N/0.05m = −4.3 kPa, a motivational
result considering that equation (2.1) is approximate since branch vectors between
particles and walls were not included.
The result of this experiment underscores the power of GEM for inferring inter-
particle forces in real opaque materials. The experiment also illustrates the versatility
of GEM: it is adaptable to any experimental technique able to furnish the required
 17

Figure 2.6: Resulting inter-particle forces (in N) from application of GEM.

input. The following sections describe GEM’s mathematical framework in more detail
and describe precisely how the required ingredients are used to produce inter-particle
force estimates.

2.2 Governing equations

This section describes the equations governing particle behavior in cohesionless gran-
ular materials in static equilibrium. Underlying the formulation of all governing
equations is the assumption of point contact. This assumption is valid for many stiff
particles of interest in granular physics, including, for instance, sands, powders, and
many pharmaceuticals. In the case of more compliant particles, this approximation
may remain accurate for moderate load levels.
 18
2.2.1 Equilibrium

Consider the pth particle in a static granular material, interacting with other particles
through Ncp contact points labeled with index α (see Figure 2.7). Balance of forces
and moments yields the two equilibrium equations:

p
Nc
X
fα = 0 (2.3)
α=1
p
Nc
X
xα × f α = 0 (2.4)
α=1

where xα is a vector from a conveniently chosen origin to the contact point α, and
f α is a force vector acting at α.

fj
ry
bounda

⌦p
fi

⌦q

Figure 2.7: Illustration of particle-particle and particle-boundary contacts.

Equations (2.3) and (2.4) can be combined into a single matrix expression for an
entire assembly of particles: K eq f = 0. In two dimensions, the system takes the
 19
form:

i j
    
.. .. ..
. 0 ··· 0 ··· . .
     
   i  
p 0 K ieq 0 K jeq 0   f  0
     
 .. ... ..   ..   .. 
 . 0 0 .   .  = . (2.5)
     
   j  
q 0 −K ieq 0 0 0  f  0
     
.. .. .. .. ..
. 0 . 0 . . .
| {z } | {z } | {z }
K eq f beq
where
 
1 0  
  f1i
  i  
K ieq = 0 1 ; f =
  f2i
−xi2 xi1

and where p and q represent particles, i and j represent particle-particle and particle-
boundary contacts as shown in Figure 2.7.
The extension to three dimensions is straightforward and omitted for brevity. In
general, K eq will have dNp (d + 1)/2 rows and dNc columns where d is the dimension
(e.g., d = 2 for 2D), Np is the total number of particles in the assembly, and Nc is
the total number of contact points in the assembly.

2.2.2 Average particle stress

The average Cauchy stress for a particle in equilibrium under the action of discrete
boundary forces can be derived by considering the volume averaged stress equation
for a particle p: Z
p1
σ̄ = σ p dv (2.6)
Ωp Ωp

where Ωp indicates integration over the deformed volume (in 3D) or area (in 2D) of
the particle p. By considering balance of linear momentum, the divergence theorem,
and the symmetry of the Cauchy stress tensor, this expression takes the form (see [48]
 20
for more details):
p
Nc
1 X
σ̄ = p
sym(f α ⊗ xα ) (2.7)
Ωp α=1

Equation (2.7) can be written in matrix form for an entire assembly of particles
as K st f = bst . In two dimensions, the system takes the form:

i j
     
.. .. ..
. 0 ··· 0 ··· . .
     
   i  p 
p 0 K ist 0 K jst 0   f  bst 
     
 .. .. ..   ..   .. 
 . 0 . 0 .   . = .  (2.8)
     
   j  q 
q 0 −K ist 0 0 0  f  bst 
     
.. .. .. .. ..
. 0 . 0 . . .
| {z } | {z } | {z }
K st f bst
where
   
p
xi1 0 Ωp σ̄11
   
   p 
K ist =  0 xi2  ; bpst =  Ωp σ̄22 
   
i i p
x2 x1 2Ωp σ̄12

Extension to 3D is again straightforward and is omitted for brevity. Similar to

K eq , the matrix K st has dNp (d+1)/2 rows and dNc columns, where d is the dimension
(e.g., d = 2 for 2D), Np is the total number of particles in the assembly and Nc is the
total number of contact points in the assembly.

2.2.3 Cohesionless and tangential forces

Cohesion-less granular materials in equilibrium obey two additional laws: normal

forces are repulsive and tangential forces are governed by a Coulomb type friction
law. The latter constraint requires that |ft | ≤ µ|fn |, where ft is a tangential force
magnitude, fn is a corresponding normal force magnitude, and µ is the inter-particle
coefficient of friction. These two constraints can be written for any contact point
 21
acting on the pth particle as:

−ei · f i ≥ 0 (2.9)
1 i
−(ei + t ) · fi ≥ 0 (2.10)
µ
1
−(ei − ti ) · f i ≥ 0 (2.11)
µ

where ei and ti represent normal and tangential unit vectors at the contact point i
for a particular particle Ωp , as shown in Figure 2.8.

ti
fi
⌦p

ei
⌦q

Figure 2.8: Unit vectors for contact i.

Equations (2.9), (2.10), and (2.11) can be combined into a single matrix expression
for an entire assembly of particles: Bf ≥ 0. In 2D, the system takes the form:

i i j j
 
i −ei1 −ei2 0 0 0
 
..  .. 
. 0 0 . 0 0 
    
 
j 0 0 0 −ej1 −e2 j
 f1i 0
    
 i   i   .. 
i  −e1 − µ1 ti1 −ei2 − µ1 ti2 0 0 0   f2   . 
    
..  ..   ..   .. 
. 0 0 . 0 0   .  ≥ . (2.12)
    
 1 j    .
j 0 0 0 −ej1 − µ1 tj1 −e2 − µ t2  f1j   .. 
j
    
 i 1 i 
i  −e1 + µ t1 −ei2 + µ1 ti2 0 0 0  f2j 0
 
..  .. 
. 0 0 . 0 0 
 
j 0 0 0 −ej1 + µ1 tj1 j 1 j
−e2 + µ t2

where subscripts refer to vector components and superscripts refer to particular con-
 22
tact points. Extension to 3D is straightforward and is omitted for brevity. In general,
B ∈ R3Nc ×dNc .

2.2.4 Structure of equations and connection to experiments

The three matrix equations detailed in this section encapsulate equilibrium relations,
constitutive relations, and contact laws for each particle in a granular material in
terms of the unknown force components. When the sum of the number of linearly
independent rows in both K st and K eq exceeds the number of columns in one of these
matrices, there are more equations than unknown force components. Furthermore,
when K eq has more linearly independent columns than rows, K eq f = 0 has an infinite
number of solutions and forces can be inferred using the inverse problem formulation
detailed in Section 2.3. In many practical examples of interest, these conditions
will be met as indicated by observed coordination numbers in stable static granular
packings [64].
Experimentally measured quantities are required input to the three matrix equa-
tions (2.5), (2.8), (2.12). These quantities can be measured with any suitable pro-
cedure and the structure of the equations will remain unchanged. In particular, the
matrices K st , K eq , and B contain contact locations and normal and tangent vectors
to each contact plane, which can be extracted from experiments using high-resolution
photography or X-ray techniques in connection with DIC or level-set methods. Fric-
tion coefficients in B can be estimated using existing literature or experiments. Fi-
nally, the vector bst contains particle areas (or volumes) and average particle stresses.
Particle areas (or volumes) can be extracted using image (or data) processing tech-
niques. Particle stresses can be determined by first extracting strain fields using DIC
or other algorithms, and then by applying a suitable constitutive relation. Because
the mathematical framework incorporates only a single dataset corresponding to one
instant in time, it does not currently allow for history dependent constitutive relations
and therefore requires a linear or nonlinear elastic relation.
 23
2.3 Mathematical framework
This section presents the inverse problem formulation for inferring inter-particle forces
using the equations K eq f = 0, Bf ≥ 0, and K st f = bst .

2.3.1 Inverse problem formulation

When the conditions discussed in section 2.2 are satisfied, the unknown force compo-
nents can be inferred by satisfying equilibrium precisely, and minimizing an L2 -norm
cost function involving forces and average particle stresses:

f = arg minkK st f − bst k2 (2.13a)

f

subject to: K eq f = 0 (2.13b)

Bf ≥ 0 (2.13c)

The L2 -norm cost function (2.13a) is typical in inverse problems and represents
only one possible cost function for force-inference. In theory, the mathematical frame-
work may be extended to incorporate other cost functions when the structure of noise
in bst (or in other variables) is well known. In general, however, this noise cannot be
easily characterized for all possible experimental techniques adaptable to GEM, and
in fact is very difficult to characterize even for a single procedure such as DIC (e.g.,
see [65]). Therefore, use of the L2 -norm cost function is retained in the present chap-
ter as an example of a method with simple interpretation and implementation, and
simple conditions for existence and uniqueness, but not as the only possible method.
The inverse problem formulation in equation (2.13) has been reliable in all practical
implementations by the authors.
It is important to note that the constraint in equation (2.13c) is not needed in
most cases to obtain an accurate solution. In fact, this constraint is unnecessary and
does not influence the solution obtained using equation (2.13) except when signifi-
cant noise is present in experimental measurements. In the case when such noise is
present, the constraint (2.13c) plays the role of ensuring that the selected solution
 24
is physically admissible since noise may bias the result towards a force distribution
with unrealistically high tangential to normal force ratios or attractive contact forces.
The constraint (2.13c) can also be modified to eliminate the dependence on Coulomb
friction when a Coulomb friction law is unjustified or when the friction coefficient
cannot be estimated. When Coulomb friction is abandoned, the restriction will only
require that normal forces be repulsive.
Conditions for the existence and uniqueness of a solution to equation (2.13) are
easy to understand. It must be emphasized, however, that existence and uniqueness
of a solution to equation (2.13) does not imply that the solution corresponds to the
true inter-particle forces but rather to a solution that minimizes the error between
experimental observations and calculations made with the governing equations of the
problem in an L2 -norm sense.
Let S = {f |K eq f = 0, Bf ≥ 0}, the set of force vectors that satisfy the con-
straints (2.13b) and (2.13c). A solution to (2.13) exists when S is nonempty since
the cost function in (2.13a) is bounded below by 0, e.g., when K st is positive semi-
definite. Furthermore, S is nonempty when K eq has more columns than rows, and
when Bf ≥ 0 is solvable, which can be ensured in practice by choosing µ conserva-
tively. The solution is unique if there is no w 6= 0 such that w ∈ N (K st ) ∩ N (K eq )
(see proof and other conditions in Theorem 1 of [66]); such a w could be added to
any existing solution without changing the value of the cost function in (2.13a) or
violating the equality (2.13b).
From a physical perspective, such a w is unlikely to exist: nonzero forces in
N (K st ) satisfy K st f = 0 and must cause rigid body particle motion while forces
in N (K eq ) satisfy K eq f = 0 and must result in equilibrium. When K eq and K st
are rank deficient, such a w may exist and additional criteria must be satisfied to
ensure uniqueness; namely, the w ∈ N (K st ) ∩ N (K eq ) must violate the restrictions
of constraint (2.13c) to ensure that the solution to (2.13) remains unique. It is merely
stated here that the authors have never found such a w to exist in both numerical
simulations and experiments, ensuring the uniqueness of solution to (2.13).
 25
2.3.2 Measurement noise and alternative formulation

Experimental imaging techniques contain error, or measurement noise. In addition,

algorithms used to extract strain fields, contact locations, contact planes, and con-
stitutive law parameters introduce noise. This noise manifests itself in the matrices
K eq , K st , and B and in the vector bst . While high-fidelity imaging can typically
ensure negligibly small error in the point-wise quantities used to populate K eq , K st ,
and B, the vector bst requires accurate estimation of the particle constitutive model
and associated parameters, and involves a sum over all point-wise measured stresses
in a particle, potentially introducing significant error.
To account for the possibility of significant measurement error in bst a simple
alternative to the inverse problem (2.13) is proposed which incorporates knowledge
of boundary forces, quantities that are typically found by using load cells in experi-
ments. The alternative method is motivated by the experience that when the solution
to (2.13) is affected by measurement error in the constitutive law, the relative sizes
of forces remains relatively unchanged (i.e., all forces are generally over- or underes-
timated). The alternative inverse problem is given by:

f = arg minkK st f − bst k2 + λ2 kf k2 (2.14a)

f

subject to: K eq f = 0 (2.14b)

Bf ≥ 0 (2.14c)

where λ is a regularization parameter, to be discussed.

Problem (2.14) employs Tikhonov regularization, a common technique in the so-
lution of ill-posed or rank-deficient inverse problems (e.g., [67, 68]). Tikhonov regu-
larization can be interpreted as a method for incorporating prior knowledge of the
magnitude of f and as a technique for selecting a “smoother” solution. This can be
seen intuitively: as λ is increased, the solution will decrease until each value of f
approaches 0. Furthermore, the solution will typically decrease with some uniformity
in that each value of f will approach 0 at a rate proportional to its size. It is im-
 26
portant to note, however, that the solution to problem (2.14) will still satisfy particle
equilibrium.
The primary challenge in using the alternative form (2.14) is selecting λ to find
a tradeoff between noise reduction and loss of information [68]. When the structure
of noise in bst is well known, many methods exist to select λ, although all methods
have known issues and limitations [69]. The technique used in this chapter to select
λ is simple and has been proven effective by experience: gradually increase λ until a
solution f to (2.14) minimizes the difference between known and calculated boundary
forces. Example of implementing this procedure are presented in the next section and
the experimental example in section 2.1. In general, the mathematical framework
presented in this section may be extended to incorporate penalty functions other
than Tikhonov regularization, but it is beyond the scope of this chapter to address all
of these possibilities. It is important to note that solution of equation (2.14), as with
the solution of equation (2.13), will produce forces that satisfy particle equilibrium.

2.3.3 Implementation

Many solvers and optimization packages exist to solve the minimization problem
(2.13) or (2.14). In particular, implementations of SeDuMi [70] such as CVX [71] and
Yalmip [72] provide efficient ways to solve the problem in Matlab. In addition, any
numerical optimization package capable of solving quadratic programming problems
can be used to solve (2.13) and (2.14) when the cost functions in these problems
are expanded to quadratic form (e.g.: min{f T K Tst K st f − 2bT K st f }). Since the
matrices K eq and K st may not be full rank in practice, solvers capable of handling
rank deficient matrices may be preferred.

2.4 Numerical example

This section presents a numerical example of implementing GEM. This example is
used to demonstrate the salient features of the method in an environment where the
exact solution, or exact inter-particle forces, is known. A metric used throughout this
 27
example is relative error, kδf k, defined as:

kf Calculated − f Actual k
kδf k = (2.15)
kf Actual k

To generate a suitable data set for demonstrating the features of GEM, an odomet-
ric compression test was performed on a rectangular sample of 58 particles using the
Discrete Element Method (DEM) [33]. The particle radii were uniformly distributed
between 2.5 mm and 4 mm. The particle-particle and particle-wall coefficients of fric-
tion were 0.5 and zero, respectively. No gravity acted on the particles. The normal
and tangential stiffnesses were 67.5 KN/mm. The applied compressive stress was
30 kPa on the top face, as shown in Figure 2.9.

2.4.1 Global force inference

After the DEM simulation neared an equilibrium state, contact forces, normal and
tangential contact vectors, and average particle stresses determined by equation (2.7)
were used to form the matrices K st , K eq , B, and the vector bst . The inter-particle
friction coefficient for all contacts was assumed to be 0.5. The problem (2.13) was
solved and the exact inter-particle forces were recovered with relative error, kδf k ≈ 0.
Since forces from DEM exactly matched those found using (2.13), the forces are
plotted only once in Figure 2.9. Only forces with magnitude over 300 N are labeled.

2.4.2 Noise reduction

To demonstrate the application of the modified inverse problem (2.14), noise was
added to each element of bst as follows: each element was first increased by 15% to
simulate overestimating the Young’s modulus in a linearly elastic constitutive model,
and then each element was added to the product of its original value with a unique
number generated from a Gaussian distribution with mean 0 and standard deviation
0.1 to simulate noise generated from image resolution and strain field calculation
 28
algorithms. Each element of the new bst therefore took the form:

(bst )i = (bst )i + (0.15 + N (0, 0.1))(bst )i (2.16)

where N (0, 0.1) represents a number generated from a Gaussian distribution with
mean 0 and standard deviation 0.1.

(a) (b)

Figure 2.9: (a) Numerical odometric test setup. (b) Inter-particle forces computed
with DEM and (2.13). Length scale in meters, forces in Newtons. Line thickness
proportional to force magnitude.

Problem (2.13) was first solved, yielding results with kδf k = 0.154 as shown in
Figure 2.10. Despite the noticeable differences between the forces in Figure 2.10 and
the exact forces, the solution in Figure 2.10 does not violate Coulomb friction by
virtue of the inequality constraints used in the optimization problem (2.13).
Next, problem (2.14) was solved, yielding results with kδf k = 0.053 as shown in
Figure 2.10. The total force on the top boundary measured from the DEM simulation
was approximately 1170 N. The regularization parameter λ was chosen to minimize
the difference between this value and the sum of the particle-boundary reaction forces
 29

(a) (b)

Figure 2.10: (a) Solution to (2.13) with artificial noise, kδf k = 0.154. (b) Solution
to (2.14) with same noise, kδf k = 0.053 (right).

on the top boundary:

Nc b
X
Total Difference = |1170 + f2α | (2.17)
α

where Ncb is the number of contact points on the top boundary and f2α is the vertical
component of force acting at contact points along the top boundary. A value of
λ = 0.0714 was found to minimize this total difference by solving (2.14) with a
variety of values over an interval between 0.06 and 0.075.
From the reduced relative error and a qualitative comparison of Figures 2.9 and
2.10, it is clear that the formulation in equation (2.14) significantly improved the
accuracy of the solution to the inverse problem. The effect of overestimating Young’s
modulus was nearly eliminated and the relative error was reduced to levels more repre-
sentative of the Gaussian noise alone. Other numerical examples not presented here
demonstrate the same capability with other forms of artificial measurement noise.
This alternative formulation of the inverse problem represents an important exten-
sion of GEM for practitioners, providing a means for obtaining solutions with higher
accuracy. It is important to note, however, that this method is largely based on judg-
 30
ment: selection of λ is an interactive process, particularly when using the method of
minimizing a non-convex function like the total difference in (2.17).

2.4.3 Local force inference

The proposed formulation in equation (2.13) has the capability of extracting forces
locally within a material. To demonstrate this capability, force inference was per-
formed on six particles, as shown in Figure 2.11. The contact points, corresponding
normal and tangent vectors, and average particle stresses for these particles only were
used to form the quantities K st , K eq , bst , and B. The contact forces were inferred
using (2.13). Comparison with Figure 2.9 shows that the inferred local contact forces
exactly match those from the DEM simulation.

Figure 2.11: Solution to (2.13) performed locally within granular material.

The capability of the proposed method to infer contact forces locally represents an
important feature of GEM. This feature allows force inference in regions of interest
(e.g., shear bands) without the necessity of solving an inverse problem for the entire
assembly. Furthermore, this feature may allow GEM to be used in connection with
 31
larger experiments since only a small portion of the material’s fabric needs to be
imaged and processed.

2.5 Conclusion
The improved formulation of GEM presented in this chapter provides a powerful
methodology for investigating inter-particle forces in granular materials. The first ex-
perimental example using GEM has been showcased to demonstrate that the method
can and has been applied to real materials. The presentation of the mathematical
framework illustrates its simplicity and versatility, and offers many possibilities of
extending the framework to incorporate additional experimental measurements such
as boundary forces. With progress in experimental imaging and intra-particle strain
field extraction (e.g., [49, 50]), GEM will soon be able to extract inter-particle forces
in materials with smaller grains like sands, providing the first chance to validate
many theories regarding force networks in natural granular media. GEM will also
advance the boundaries of the micro-mechanical understanding of granular materials
by offering a characterization tool for a new class of opaque complex assemblies.
 32

Chapter 3

Dynamic inter-particle force

inference in granular materials:
Method and application1

3.1 Introduction
In this chapter, we extend GEM to the dynamic regime, providing a new tool for quan-
titative inter-particle force measurements in granular media. The method currently
is restricted only by the limitation of experimental techniques to supply measure-
ments of particle positions, contact points, and volume-averaged grain stresses in the
material of interest. However, the mathematical foundation and inverse problem for-
mulation that the method employs is independent of particle properties or contact
law.
In section 3.2 we discuss the necessary experimental measurements, image pro-
cessing techniques, and theory behind the numerical optimization problem used in
the method. Section 3.3 illustrates an experimental example in which the method is
used to infer forces between a group of two-dimensional rubber grains impacted by
a foreign intruder. A comparison of the experimental results with a finite-element
simulation is provided to validate the results of the method. Section 3.4 will feature
results from a preliminary application of the method to force inference during impact
1
Adapted from R.C. Hurley, K.W. Lim, G. Ravichandran, J.E. Andrade. Dynamic inter-particle
force inference in granular materials: Method and application. Experimental Mechanics, 2015. DOI:
10.1007/s11340-015-0063-8
 33
of an intruder on a model granular bed. Section 3.5 will discuss the current challenges
faced during application of the method and provide insight into the new physics the
tool can unravel. Finally, section 3.6 will offer concluding remarks.

3.2 The method

The method for inferring inter-particle forces involves two steps: (1) experimental
measurements, and (2) solving a numerical optimization problem. Item (1) is briefly
discussed to highlight the required measurements and some techniques capable of
providing them. Item (2) is discussed in more detail since the numerical optimization
problem is the cornerstone of the method.

3.2.1 Experimental measurements

The first step of the method involves performing a dynamic experiment of interest
and extracting the evolution of particle positions, inter-particle contact points, inter-
particle contact surface normals, and the evolution of volume-averaged grain stresses.
These measurements are required as input for the numerical optimization problem
and, once measured, provide enough information to accurately infer inter-particle
forces. No restrictions are inherently placed on the particle shape or dimensionality
of the experiment and no prior assumption is needed of the contact law between
grains.
Particle positions and contact points can be found in 2D experiments using high-
speed photography and particle segmentation algorithms. In 3D, XRCT provides
a suitable method for extracting particle positions and contact points. As stated
in the introduction, XRCT is still an area of active research and has not yet been
applied to large dynamic experiments [52,53,54]. The proposed force inference method
is therefore not yet suitable for application to 3D experiments. Nevertheless, we
will elaborate on the 3D formulation of the inverse problem for completeness. Once
particle positions and contact points have been obtained, normal and tangent vectors
can be found for each contact point using the local edge features.
 34
Volume-averaged grain stresses can be found from volume-averaging the full-field
grain stress computed point-wise. For linear elastic materials volume-averaged grain
stresses can also be computed directly from volume-averaged grain strains using σ̄ p =
C : ¯p , where C is the stiffness tensor. In either case, a constitutive law is required
to obtain volume-averaged grain stresses. Such constitutive laws - linear elastic,
hyperelastic, plastic, etc. - can be determined using experiments not discussed here.
Full-field grain strains are needed to obtain volume-averaged grain stresses when
experiments involve non-linear elastic grain materials or grains undergoing plastic
deformation. Full-field grain strains can be found using digital image correlation
(DIC) [55] in 2D or a combination of 3D X-ray diffraction (3DXRD) and XRCT in
3D [49, 50, 51, 73]. Digital volume correlation (DVC) can also be used with XRCT or
similar methods to obtain full-field strains in 3D. The methods for obtaining full-field
grain strains (DIC, DVC, XRCT, 3DXRD) typically rely on internal grain contrast
(e.g., see [55]). Like XRCT, 3DXRD is a field of active research and has not yet been
applied to dynamic experiments.

3.2.2 Numerical optimization

The numerical optimization problem used to infer inter-particle forces from experi-
mental measurements requires three sets of governing equations: momentum balance,
stress-force relations, and constraint equations.
The momentum balance equations for a deformable body are derived as follows.
The balance of linear momentum arises from noting that the total force on a de-
formable particle equals the time derivative of its linear momentum:

Z ! Z Z
d
ρv dv = ρb dv + t ds (3.1)
dt Vp Vp ∂Vp

where ρ is the particle density, b is the body force, t are surface tractions, v is velocity,
and Vp and ∂Vp represent deformed grain volumes and edges, respectively. Taking
 35

ti fk
fi Vp
ei
xi x̃
xcm
p
x
fj
x2
BOUNDARY OF
REGION OF INTEREST
x3 x1

Figure 3.1: Illustration of various terms used in the governing equations.

the derivative and ignoring body forces, one obtains

Nc
X
p Z
α
f = ρa dv (3.2)
α=1 Vp

where f α is a contact point as shown in Fig. 3.1 and Ncp is the total number of
contact points for particle p. We have made use of the fact that integrating tractions
on disjoint sections of a surface is equivalent to summing the total force on each of
those sections. We can further reduce Eq. (3.2) by making use of the definition of
center-of-mass acceleration to obtain

p
Nc
X
f α = mp acm
p (3.3)
α=1

The balance of angular momentum arises from noting that the total torque on a
deformable particle is equal to the time derivative of its angular momentum:

Z ! Z Z
d
x × ρv dv = x × ρb dv + x × t ds (3.4)
dt Vp Vp ∂Vp

where x is a position vector and × represents a cross product of two vectors. Taking
the derivative and ignoring body forces, one obtains

Nc
X
p Z
c c
x ×f = x × ρa dv (3.5)
c=1 Vp
 36
where we have made use of the same argument as above to turn the integral into
a sum while approximating the contact point xc as the center of the contact area.
Now define x = xcm cm cm cm
p + x̃ and a = ap + ã, where xp and ap are the position and

acceleration of the center-of-mass, respectively. Plugging these into Eq. (3.5) yields

Nc
X
p Z
c c
x ×f = mp (xcm
p × acm
p ) + ρ(x̃ × ã) dv (3.6)
c=1 Vp

Note that the second term on the right hand side of this equation is negligibly small
when accelerations within a grain remain small relative to its center-of-mass acceler-
ation.
Linear and angular momentum balance, Eqs. (3.3) and (3.6), can be combined
into a single matrix equation for a group of contacting particles, taking the form

i j
     
.. .. ..
. 0 ··· 0 ··· . .
     
   i  p 
p 0 K im 0 K jm 0   f  bm 
     
 .. ... ..   ..   .. 
 . 0 0 .   . = .  (3.7)
     
   j  q 
q 0 −K im 0 0 0  f  bm 
     
.. .. .. .. ..
. 0 . 0 . . .
| {z } | {z } | {z }
Km f bm

where p and q represent particles, i and j represent particle-particle and particle-

boundary contacts as shown in Fig. 3.1, and K im , f i , and bpm encompass Eqs. (3.3)
 37
and (3.6) and in 2D take the form
 
1 0  
  f1i
 
1 ; f =  
i
K im =  0 (3.8)
  f2i
i i
−x2 x1
 
mp acm
 1 
p  
bm =  mp a2cm  (3.9)
 
cm cm cm cm
mp (x1 a2 − x2 a1 )

where we have omitted the integral term from Eq. (3.6) for brevity. In general, 3Np
equations are contained in Eq. 5.1, where Np is the number of particles in the region
of interest. Extension of K im , f i , and bpm to 3D is straightforward.
The stress force relations are derived as follows. The volume-averaged stress in a
deformable particle can be written as
Z
1
σ̄ p = σ dv (3.10)
Vp Vp

Left-multiplying the σ within the integral by the identity matrix xi,k and switching
to index notation yields
Z
1
σ̄ij = [(xi σkj ),k − xi σkj,k ] dv (3.11)
Vp Vp

Using the divergence theorem on the first term inside the integral of this equation
and balance of linear momentum in the absence of body forces (σkj,k = ρaj ) on the
second term, we obtain

p
Nc
X Z
c c
x ⊗ f = Vp σ̄ Vp + x ⊗ ρa dv (3.12)
c=1 Vp

where ⊗ is a dyadic product (a ⊗ b = ai bj ). Finally, defining x = xcm

p + x̃ and
 38
a = acm
p + ã as above, we obtain

Nc
X
p Z
c c cm cm
x ⊗ f = Vp σ̄ + mp (x ⊗a )+ x̃ ⊗ ρã dv (3.13)
c=1 Vp

Once again, we note that the third term on the right hand side of this equation is
negligibly small when accelerations within a grain remain small relative to its center-
of-mass acceleration.
The stress-force relation in Eq. (3.13) can also be written into a matrix equation
for a group of contacting particles, taking the form

i j
   
 
... .. ..
0 ··· 0 ··· . .
     
   i   p
p 0 K is 0 K js 0   f  b s 
     
 .. .. ..   ..   .. 
 . 0 . 0 .   . =. (3.14)
     
   j  q
q 0 −K is 0 0 0  f   b s 
     
.. .. ... .. ..
. 0 . 0 . .
| {z } | {z } | {z }
Ks f bs
where f i is the same as in Eq. (3.8), and K is and bps encompass Eq. (3.13) and in
2D take the form

 
i
x 0
 1 
i  
K s =  0 x2  i (3.15)
 
i i
x2 x1
 
cm cm
Vp σ̄11 + mp (x1 a1 )
 
 
bps =  Vp σ̄22 + mp (xcm
2 a cm
2 )  (3.16)
 
cm cm cm cm
2Vp σ̄12 + mp (x1 a2 + x2 a1 )

where we have omitted the integral terms from Eq. (3.13) for brevity. In general,
3Np equations are also contained in Eq. 3.14. Extension of K is and bps to 3D is
 39
straightforward.
Constraint equations can also be written for a group of cohesionless grains obeying
a Coulomb-type friction law. These constraints in 2D require that (f i · ei ) ≥ 0 and
(µei ± ti ) · f i ≥ 0 where ei and ti are the normal and tangential to the contact surface
of contact point i as shown in Fig. 3.1 and µ is the inter-particle friction coefficient.
This friction coefficient can be measured from a separate experiment performed on the
particle material and may be either dynamic or static depending on the nature of the
experiment. In 3D, the tangential direction must be specified by two vectors and the
Coulomb friction constraint written above must be modified during implementation in
the numerical optimization problem. We will continue with the 2D form and discuss
the 3D extension briefly below.
The cohesionless force and Coulomb friction constraint equations can be expressed
as a matrix equation for a group of contacting particles, taking the form Bf ≥ 0
where

i i j j
 
i ei1 ei2 0 0 0
 
..  ... 
. 0 0 0 0 
 
 
j 0 0 0 ej1 ej2 
 
 
i  µei1 + ti1 µei2 + ti2 0 0 0 
 
..  .. 
B = . 0 0 . 0 0  (3.17)
 
 j 
j 0 0 0 µej1 + tj1 j
µe2 + t2 
 
 i 
i  µe1 − ti1 µei2 − ti2 0 0 0 
 
..  .. 
. 0 0 . 0 0 
 
j 0 0 0 µej1 − tj1 j j
µe2 − t2

and f = (f1i f2i . . . f1j f2j )> . In general, 3Ncp constraints are contained in Eq. (3.17),
where Ncp is the total number of contacts.
The mechanics of a group of contacting grains is fully governed by the three matrix
Eqs. (5.1), (3.14), and Bf ≥ 0. Every term in these equations except for the inter-
particle forces f must be measured during an experiment as described in section 3.2.1.
 40
Once measured, finding the set of forces f that best fits the measured data involves
solving a multi-objective optimization problem of the form:

f = arg min (||K s f − bs ||2 + λ||K m f − bm ||2 ) (3.18a)

f
subject to: Bf ≥ 0 (3.18b)

Here, the objective functions ||K s f − bs ||2 and ||K m f − bm ||2 represent the desire
to minimize the difference between inferred forces and experimental measurements of
stress and momentum, respectively. The value λ is a weight that conveys the relative
importance of the two objectives. No previous knowledge of boundary forces is needed
as in past implementations of this inverse problem [13].
In order to solve Eq. (3.18), a suitable value of λ must be selected. With no prior
knowledge of the type or magnitude of error in stress or momentum measurements,
the best strategy is to seek a solution which does not favor one set of measurements
(e.g., momentum or stress) over another. To do this, we solve Eq. (3.18) with many
values of λ and obtain a curve of ||K s f − bs ||2 versus ||K m f − bm ||2 , as shown in
Fig 3.2. All points on this curve for which neither objective can be decreased without
increasing the other are referred to as Pareto optimal points, and the resulting curve
is referred to as the optimal trade-off curve or the Pareto front [74].

OPTIMAL TRADE-OFF
||K s f bs ||2 CURVE / PARETO FRONT

KNEE POINT

||K m f bm ||2
1
1

Figure 3.2: Pareto front.

 41
In principle, any point on the optimal trade-off curve is a possible solution to
the multi-objective problem. However, in order to use as much information in the
two objective functions as possible, we wish to select a point where the slope of the
trade-off curve is steep. At a steep point, a small decrease in one objective function is
accompanied by a small increase in the other. This steep “knee” of the optimal trade-
off curve therefore represents the solution where we satisfy each objective without
making an unfavorably large sacrifice in the other.
Finally, we mention that when solving Eq. (3.18) for 2D experiments, special care
must be taken because of the nonzero thickness of the disks or cylinders involved. In
particular, the mass used in Eqs. (5.1) and (3.14) must be mass per unit thickness
(i.e., actual mass divided by thickness of the particles) and the volume in Eq. (3.14)
must be the area of the face of the particle. The result of solving Eq. (3.18) will be
force per unit thickness, from which actual forces can be obtained by multiplying by
the particle thickness.
In 3D, a contact surface between two particles must be specified by a normal vector
and two orthogonal tangent vectors, ti1 and ti2 , rather than one. Defining q1i and q2i
to be the magnitude of f projected onto each tangent vector, the multi-objective
optimization problem takes the form

f = arg min (||K s f − bs ||2 + λ||K m f − bm ||2 ) (3.19a)

subject to: ei · f i ≥ 0 (3.19b)

q1i = (ti1 · f ) (3.19c)

q2i = (ti2 · f ) (3.19d)

µ(ei · f i ) − ||q i || ≥ 0 (3.19e)

where ||q i || = ((q1i )2 + (q2i )2 )1/2 . Many of the constraints in Eq. (3.19) can be cast
in matrix form as in Eq. (3.18). However, the last constraint is a second-order cone
constraint which is not present in the 2D problem and cannot be written in matrix
form, instead requiring a norm [74]. For more details on this type of constraint,
 42
see [74].

3.3 Example 1: Validation

This section contains a validation of the proposed force inference technique. An
experiment is performed with a group of two-dimensional rubber grains impacted by
a foreign intruder. Force inference is performed with the proposed method and the
results are compared with those of a finite-element simulation of the experiment. The
section is presented in three steps: (1) experimental setup; (2) experiment analysis;
(3) experimental results and FEM comparison.

3.3.1 Experiment setup

The experimental setup is shown in Fig 3.3. The setup involves a small 0.4m x 0.6m
table containing an air chamber under its upper surface. The chamber is fed through
an air duct connected to a fan (Hydrofarm Active Air 6 inch in-line fan) and air is
permitted to escape the top surface of the table through 1.6mm diameter holes drilled
on a grid with 19mm spacing in perpendicular directions, as shown in Fig. 3.3a-b.
Six 44.45mm diameter by 6.35mm thick disks were cut out of 60A durometer
polyurethane (McMaster-Carr product number 8784K55) with measured Young’s
modulus E = 5.85 MPa and Poisson’s ratio ν ≈ 0.5 through compression testing
up to 13% strain. The friction coefficient between disks was measured by pulling a
block (of known weight) fastened to a sheet of polyurethane across another sheet of
the same material and measuring the horizontal force with a load cell. The mean
friction coefficient from this experiment is approximately 0.6. This coefficient is used
throughout the numerical analysis since neighboring disks appear to engage mostly
static friction throughout the test. However, future work should address the quan-
titative sensitivity of the force inference results to this parameter as well as the ap-
propriate selection between a static and dynamic friction coefficient. The underside
of each disk is endowed with a shallow 38.1mm diameter by 1.19mm thick concentric
cut to create a small plenum chamber, as shown in Fig. 3.3c. This allows the disks to
 43
“float” on a frictionless bed of air on the top surface of the table and interact with one
another through mechanical contact only. The top surface of each disk is endowed
with a painted speckle pattern to enhance DIC results, as shown in Fig. 3.3b.
Two stiff wooden blocks are fastened to the top surface of the table to form a corner
and the six disks are placed at rest in the corner in two rows of three, as shown in
Fig. 3.3d. A stiff wooden impactor is then manually propelled toward the disks with
a measured initial velocity of vx = −1.141m/s, vy = −0.66m/s and negligible initial
rotational velocity. Initial contact between the disks and the impactor occurs at an
angle of approximately 27o with respect to the center of the top right particle in Fig
3.3d. The entire impact event is captured with a Phantom v310 high-speed camera
at 5000 frames per second. The impact and subsequent rebound of the impactor
and all disks are monitored for 0.028 seconds to obtain the necessary input to the
multi-objective optimization problem in Eq. (3.18).
A finite-element (FEM) model with the same initial conditions discussed above
was created and analyzed in Abaqus/Explicit [75] to compare inferred inter-particle
forces with the results of the multi-objective optimization problem. FEM was chosen
to model the grains rather than DEM because FEM uses only material properties and
the governing equations of continuum mechanics, making no assumption of contact
law. Since our goal is to validate a method that makes no assumption of contact law in
an experiment where we do not necessarily know the contact law, it is logical to use a
validation method with the same features. DEM instead relies on a prescribed contact
law (e.g., Hooke’s or Hertz’s law) and would likely require extensive calibration with
experimental results to produce quantitatively accurate dynamics. The grain material
is modeled using a 2D plane stress formulation with grain properties discussed above.
The walls are modeled as rigid boundaries. The impactor is modeled as a stiff 2D
linear elastic body. The Young’s modulus of the impactor does not affect the final
results so long as it is an order of magnitude higher than the grains’ Young’s modulus.
 44
44.45mm

19mm

19mm
CAMERA LIGHT
SOURCE
(b)

RIGID
BLOCKS

AIR DUCT
(a) (c)
RIGID RIGID
BLOCK IMPACTOR

2 5 1

y
4 6 3
RIGID
x BLOCK

(d)

Figure 3.3: Experimental materials and setup. (a) Table setup with air duct con-
nected to fan, rigid blocks fastened to table, camera and light source vertically over the
area of interest. (b) Close-up view of table top and the top surface of a polyurethane
disk with speckle pattern. (c) The underside of a polyurethane disk, showing plenum
chamber for compressed air. (d) Experimental setup with 6 polyurethane disks and
rigid impactor propelled at initial conditions vx = −1.141m/s and vy = −0.66m/s.
Particle labels will be used later and are in no meaningful order.

3.3.2 Experiment analysis

The DIC software VIC-2D [57] is used to calculate full-field intra-particle strains
and velocities in each of the six grains in the experiment. Particle diameters are
approximately 220 pixels in the 800x600 images captured at each time step. Subset
size is chosen to be 15 pixels and subset step size is selected as 2 pixels, giving a
ratio of particle diameter to subset size of approximately 15. Other DIC parameters
such as prediction margin, confidence interval, and matchability are selected with
default values of 0.02, 0.1, and 0.1, respectively. However, it is nearly impossible
 45
to determine how the selection of these parameters influences the calculated strain
values [76]. This is an area of ongoing research in the DIC community. Nevertheless,
our parameter selection appears to give smooth, visually reasonable results for the
displacement field within each grain. The influence of uncertainty in strain values on
the results of a similar optimization problem was discussed briefly in [13], where it
was shown that relative errors in strain values cause similar relative errors in force
values. Future work on further assessing how changing some parameter values, such
as subset size, influence the results presented here would be useful.
The xy component of the strain field as measured using DIC and FEM is shown
in Fig. 3.4 at four times during an impact event. The comparison illustrates an excel-
lent agreement between DIC measurements and FEM results. Intra-particle stresses
are obtained by using the full-field strains in a linear elastic plane stress formulation
at each pixel. It is important to note that we could alternatively use grain boundary
deformations to obtain volume-averaged grain strains (applying the divergence theo-
rem to the definition of volume-averaged strain) before computing volume-averaged
grain stresses with the linear elastic constitutive law. Material properties are taken
as E = 5.85 MPa and ν = 0.5 as mentioned above. Stresses are then averaged over
each pixel within a particular grain, yielding the average stress needed in Eq. (3.14).
Particle edges and centroids are calculated at each time step (i.e., each movie
frame) using a circular Hough transform in Matlab’s image processing toolbox [77].
In particular, the Matlab function imfindcircles is used, which makes use of the
techniques developed in [78] and [79].
The resulting centroids and radii are then used to find contacting particles: if the
distance between two centroids is less than two times the radii of the corresponding
particles, the particles are taken to be contacting with appropriate normal and tangent
vectors. While the circular geometry of our grains makes this step of the analysis
simple, the method is not restricted to circular particles.
Accelerations of the center-of-mass of each particle can be calculated in two ways:
by summing and numerically differentiating intra-grain velocity fields calculated using
VIC2D, or by twice numerically differentiating centroids obtained from the circular
 46
✏xy ✏xy
EXPERIMENT FEM
Printed using Abaqus/CAE on: Thu Nov 21 22:22:30 Pacific Standard Time 2013

t = 2ms

Printed using Abaqus/CAE on: Thu Nov 21 22:22:39 Pacific Standard Time 2013

t = 5ms

Printed using Abaqus/CAE on: Thu Nov 21 22:22:50 Pacific Standard Time 2013

t = 8ms

Printed using Abaqus/CAE on: Thu Nov 21 22:23:00 Pacific Standard Time 2013

t = 11ms

(a) (b)

✏ij  0.04 ✏ij = 0 ✏ij 0.04

Figure 3.4: The xy component of strain from (a) DIC and (b) FEM at four times
during the impact event. All figures share a common scale. Strain values above and
below the extreme values of the scale occur in small areas.

Hough transform. We use the latter approach to obtain the results presented in this
chapter, although both methods produce very similar center-of-mass accelerations.
While using the latter approach, we also employ a low-pass filter on the centroid and
acceleration time series. The filter is a finite-impulse response filter with a Hamming
window, provided by Matlab’s signal processing toolbox. Filter order and normalized
cutoff frequency are selected to maintain a signal to noise ratio greater than 10 for
each acceleration time series, as done in other works [42].

3.3.3 Experiment results and FEM comparison

The optimization problem, Eq. (3.18), is solved at each time step using CVX [71]
in MATLAB. The result of solving Eq. (3.18) is shown in Fig. 3.5 next to forces
 47
calculated using the FEM model. It is important to use caution when comparing the
forces calculated using the proposed method and the FEM model. The FEM model
does not represent the exact forces in the experiment because the material model
and initial and boundary conditions are approximate. However, the FEM model
does provide an example of the force evolution at each contact point in an idealized
experiment and should therefore capture major features of the force evolution (i.e.,
large peaks) and the appropriate order of magnitude. Integral terms in Eqs. (3.8)
and (3.15) are neglected as they were assumed to be small for the observed particle
accelerations. The results validate this decision.
The evolution of forces obtained using Eq. (3.18) shows excellent agreement with
the FEM results in Fig. 3.5. Figs. 3.5a and 3.5b illustrate that force magnitudes
and directions obtained using Eq. (3.18) qualitatively match those obtained from the
FEM model. Fig 3.5d illustrates that most major features of the force evolution at
each contact are captured using Eq. (3.18). The order of magnitude of all forces
is also in excellent agreement, often differing by less than 20% from those obtained
using FEM. Eq. (3.18) is clearly capable of extracting quantitative force information
representative of the true physical state of the system of particles.
We also compared the calculated impulse imparted to each particle using the
forces from Eq. (3.18) with those obtained from the high-speed images. The im-
pulses obtained from Eq. (3.18) are calculated by numerically integrating all forces
on a particular particle over the period of time in which they are in contact with
neighboring bodies. The results are presented in Table 3.1 in terms of ∆v obtained
from Eq. (3.18) and those measured from high-speed images. Again, the agreement
between the results of Eq. (3.18) and actual measurements is excellent in most cases.
There are some slight deviations in the calculated and measured changes in velocity,
likely due to experimental measurement error and the simplifying assumptions of the
constitutive law.
 48

(c)
1
4 3 2

FORCE MAGNITUDE FORCE MAGNITUDE

6 5
EXPERIMENT FEM 8 7
t = 2ms

CONTACT 1 CONTACT 2
40 40 40 40
F (N) FEM
FEM
FEM F (N) FEM
FEM
FEM
Eq.(16) Eq.(16)
Eq.(16)
GEM Eq.(16)
GEM
30 30
30 30

Force (N)

Force (N)
20 20

Force (N)
Force (N)
t = 5ms 20 20
10 10

Time (s) Time (s)

0 0
10 0 0.005 0.01 0.015 10 0.02 0 0.005 0.01 0.015 0.02
Time (s) Time (s)
CONTACT 3 CONTACT 4
40 40 40 40
F (N) FEM
FEM
FEM
0 F (N)
FEM
FEM
FEM
0 Eq.(16) Eq.(16)
0 0.005 0.01 0 0.02 0.005
0.015 Eq.(16)
GEM 0.01 0.015 Eq.(16)
GEM0.02
30 Time (s) 30 Time (s)
30 30

Force (N)
Force (N)

t = 8ms 20 20
Force (N)

Force (N)
20 20
10 10

Time (s) Time (s)

10 0 0
0 0.005 0.01 0.015 10 0.02 0 0.005 0.01 0.015 0.02
Time (s) Time (s)
CONTACT 5 CONTACT 6
40 40 40 40
FEM
FEM
FEM FEM
FEM
FEM
0 F (N) 0
Eq.(16) F (N) Eq.(16)
0 0.005 0.01 0.015Eq.(16)
GEM0 0.02 0.005 0.01 0.015Eq.(16)
GEM 0.02
30 30
t = 11ms Time (s) Time (s)
30 30
Force (N)

Force (N)
Force (N)

20 20
Force (N)

20 20
10 10
Time (s)
Time (s)
10 0 0
0 0.005 0.01 0.01510 0.02 0 0.005 0.01 0.015 0.02
Time (s) Time (s)
CONTACT 7 CONTACT 8
40 40 40 40
FEM
FEM
FEM FEM
t = 14ms 0 F (N) 0
Eq.(16) F (N) FEM
FEM
Eq.(16)
0 0.005 0.01 0.015 Eq.(16)
0 0.02 0.005
GEM 0.01 0.015 Eq.(16)
GEM0.02
30 Time (s) 30 Time (s)
30 30
Force (N)

Force (N)

20
Force (N)

20
Force (N)

20 20
10 10
Time (s) Time (s)

10 0
0 0.005 0.01 0.01510 0.02
0
0 0.005 0.01 0.015 0.02
Time (s) Time (s)
(a) (b) (d)
0 0
0 0.005 0.01 0.015 0 0.02 0.005 0.01 0.015 0.02
Time (s) Time (s)
Figure 3.5: (a) Forces found using Eq. (3.18). Length and width of lines are
proportional to force magnitude. (b) Forces found from the FEM model. (c) Eight
contact point locations, chosen to include the major force chain and otherwise at
random. (d) A comparison of force evolutions using Eq. (3.18) and the FEM model.
 49

Table 3.1: Comparison of expected change in velocity obtained by dividing impulse

obtained from Eq. (3.18) by particle mass and by direct measurement from high-speed
video. All results are in m/s and particle numbers are those in Fig. 3.3d.
Particle ∆vx ∆vx Particle ∆vy ∆vy
Eq. (3.18) Measured Eq. (3.18) Measured
1 0.76 0.54 1 0.28 0.25
2 0.16 0.04 2 0.12 0.12
3 0.38 0.40 3 0.07 0.04
4 0.26 0.00 4 -0.01 0.00
5 0.51 0.47 5 0.18 0.09
6 0.23 0.21 6 0.11 0.00

3.4 Example 2: Application

This section contains a preliminary application of the proposed force inference method
to the study of impact of an intruder on a model granular bed. This example is
intended to demonstrate how the method may be used in practice in experiments
more complex than the one illustrated in section 3.3. The section is presented in two
steps: (1) experimental setup and (2) experimental results.

3.4.1 Experimental setup

The experimental setup is shown in Fig. 3.6. A 1.22m tall by 0.81m wide granular
bed is prepared between a Plexiglas window and a sheet of PVC plastic. A small
1.11cm gap is provided between the Plexiglas and PVC to allow the 0.95cm thick
grains to move without significant interference. The granular bed is composed of
several hundred rubber grains made from polyurethane (E = 12MPa, ν ≈ 0.5, friction
coefficient 0.6) dropped randomly into the gap. Two grain diameters were used to
prevent crystallization: 2.54cm and 3.175cm.
Each dynamic experiment begins by dropping a 20cm steel intruder into the gran-
ular bed. Under the weight of gravity, the intruder accelerates to an impact velocity
of approximately 2.4m/s in each of the experiments presented here. The entire im-
pact event is captured with a high-speed Phantom v310 camera at 4500 frames per
second.
 50
3.4.2 Experiment analysis

As in the previous example, Vic2D is used to extract intra-particle strain fields and
MATLAB’s image processing toolbox is used to determine particle positions at each
time step. The smallest particle diameters are approximately 72 pixels in the 800x600
images captured at each time step. Subset size is chosen to be 13 pixels and the
subset step size is chosen as 2 pixels, giving a ratio of particle diameter to subset size
of approximately 5.5 in the smallest particles. This ratio is lower than in the first
example, implying that there may be larger uncertainty in the strains calculated in
this case. However, without a systematic way of determining the influence of subset
size on calculated strains, there is no way to quantify this additional uncertainty [76].
Nevertheless, parameter selection in this example also appears to give smooth, visually
reasonable results for the displacement field in each grain. Other DIC parameters are
chosen to be the same as in the previous example.
It is important to note that by using images of the prepared granular bed as
reference images in Vic2D, we are implicitly ignoring prestress in the grains. The
result of force inference therefore does not account for forces existing in the granular
bed before impact. However, in dynamic experiments like those presented here, such
pre-existing forces are likely negligible.

3.4.3 Experimental results

The inter-particle forces are once again obtained by solving Eq. (3.18) at each time
step. An example of inter-particle forces inferred as a function of time are shown
for one experiment in Fig. 3.7. A force-chain develops beneath the intruder at
approximately 5ms into the experiment before it collapses at approximately 28ms.
The emergence of a persistent force chain immediately beneath the intruder in these
experiments agrees with past experiments that use photoelastic disks [42, 43]. In
contrast to past experiments, the present work resolves individual force values at
grain-grain and grain-intruder contact points.
Figure 3.8 shows inter-particle forces in the granular bed during an experiment
 51

STEEL INTRUDER
WINDOW

STEEL
INTRUDER

CAMERA RUBBER GRAINS

(a) (b)

Figure 3.6: Experimental setup. (a) A window featuring a piece of Plexiglas and a
sheet of PVC plastic holding the 2D granular bed in place. A high speed camera is
used to view the granular bed as it is impacted by a steel intruder. (b) A close-up of
the steel intruder about to impact the granular bed at -2.4m/s.

with a different initial packing than that shown in Fig. 3.7. The response of the
granular bed is clearly different in the two cases. The force response shown in Fig.
3.8 demonstrates the emergence of force chains, as in Fig. 3.7, but through different
particles and branching in different directions. Such sensitivity to initial packing
conditions is characteristic of granular media. This behavior is the reason that past
research on force transmission and force chains has focused so extensively on the
statistics of forces [31,80]. Similar statistics are expected to arise regardless of distinct
local features of the response and may elucidate critical intrinsic behaviors of the
material.
Figure 3.9 shows the vertical force on the intruder as a function of time in the
experiments shown in Figs. 3.7 and 3.8. Both force responses demonstrate several
peaks during the time in which a persistent force chain exists beneath the intruder.
We attribute these peaks to repeated excitations of the existing force chain as the
intruder’s energy is sent into the granular bed along this path [42]. The excitations
 52

1.77ms 5.33ms

11.1ms 17.7ms

23.3ms 28.8ms
60N
30N

Figure 3.7: Inter-particle forces in the granular bed at six times during an impact
event.

may be caused by repeatedly collapsing and restrengthening of the primary force

chain as particles rearrange beneath the intruder. The primary force chain and its
collapse is labeled for the curve pertaining to Fig. 3.7. In both experiments, the
force response at early times is similar, while the force response at later times differs
significantly.
A guiding question in the ongoing research into granular bed impact response is:
Does the response of a granular bed to low-velocity impact admit any repeatable fea-
tures indicative of intrinsic collective behavior of the grains? For instance, past work
has indicated that the distribution of inter-particle forces evolves similarly during an
impact event regardless of initial packing conditions [32]. In Fig. 3.10 we demon-
strate another feature of the granular bed response that appears to be independent
 53

1.77ms 5.33ms

11.1ms 17.7ms

23.3ms 28.8ms
60N
30N

Figure 3.8: Inter-particle forces in the granular bed at six times during an impact
event with an initial packing of the bed different from the experiment shown in Fig.
3.7.
 54

250
Force chain emerges

Vertical Force (N)

Force chain collapses
200

150

100

50 Fig. 7 Exp.
Fig. 8 Exp.
0
0 0.005 0.01 0.015 0.02 0.025 0.03
Time (s)

Figure 3.9: The vertical force felt by the intruder for the experiments shown in Figs.
3.7 and 3.8.

of initial packing conditions. The three curves in Fig. 3.10 represent the maximum
force between two grains as a function of distance to the surface of the intruder in
three experiments with three distinct initial packing conditions. Experiment 1 refers
to the response shown in Fig. 3.7 and experiment 2 refers to the response shown in
Fig. 3.8. All experiments have an intruder impact velocity of approximately 2.4m/s.
Despite the emergence of force chains along different paths in each experiment, the
maximum force as a function of distance from the intruder surface is very similar
in each experiment. This suggests that, regardless of initial packing configuration,
force chains are only able to transmit large loads a certain distance into a granular
medium before the forces split into different paths and decay, a theory proposed in
past work [21, 42]. This phenomenon, and other phenomena involving inter-particle
forces are the subject of ongoing research and will be discussed in more depth in
future work. GEM offers a valuable tool for investigating these questions regarding
intrinsic material response.

3.5 Discussion
The force inference method proposed in this chapter has a number of advantages and
disadvantages. Compared to photoelasticity [42, 43], the proposed method does not
 55
250
Exp. 1

Force (N) Force (N)

Exp. 2
200 Exp. 3

Maximum Inter-Particle
150

100

50
0 2 4 6 8 10 12
Distance
Distance (d) (d )
from Intruder min

Figure 3.10: Maximum inter-particle force as a function of distance from the surface
of the intruder for three experiments with distinct initial packings. Exp. 1 refers to
Fig. 3.7. Exp. 2 refers to Fig. 3.8. An inter-particle force is considered to be a
distance ndmin from the surface of the intruder if the contact point that it belongs to
falls within a range of ndmin and (n + 1)dmin away from the surface of the intruder,
where dmin is the minimum particle diameter in the granular bed, 2.54cm.

place any inherent restriction on the grain material used in experiments. The material
can be opaque and, for materials with linear constitutive laws, without any internal
contrast. The grains can undergo plastic deformation so long as a constitutive law
exists to obtain volume-averaged stresses from strain measurements. However, the
grains’ deformation must be of the appropriate order of magnitude to be captured by
a technique for calculating strain. While the governing equations in photoelasticity
can be modified for force inference in arbitrarily shaped grains, the proposed method
requires essentially no modification, making it straightforward to apply when each
grain has a different shape. A possible disadvantage of the proposed method is the
requirement of volume-averaged strain calculation. For grains with non-linear consti-
tutive laws, this requires calculating full-field grain strains, which remains a challenge
in 3D dynamic experiments. However, in static experiments, strain fields can now be
measured in realistic 3D quartz grains by using 3DXRD, offering the opportunity to
apply the proposed method to such data sets in the near future [81].
While photoelasticity has been used to infer an approximate force averaged over
tens of grain diameters in dynamic experiments [42], the proposed method has been
 56
used to provide individual inter-particle forces. Thus, whereas photoelasticty has the
advantage of being able to infer approximate forces across larger domains with more
grains, the proposed method provides finer force resolution in local areas of interest.
This resolution has previously only been demonstrated in numerical simulations using,
for instance, the discrete element method [33] or FEM simulations [82]. The interest
in resolving individual inter-particle forces is driven by a desire to understand the
grain scale origin of the macroscopic response in many environments. For instance,
in impact into a granular bed, researchers have long been interested in understanding
frictional drag forces between the intruder and grains [32], the importance of inter-
particle friction in energy dissipation within the granular bed [83], and the dynamics
of force chains [42]. While numerical simulations can provide some of this information
in idealized systems, experimental methods are the only way to study such physics
in real systems. To date, the proposed method is the only one that has been used to
infer quantitative inter-particle force values in such dynamic experiments.
Both photoelasticity and the proposed method have strengths and weaknesses.
We are not advocating the use of one over the other in all scenarios. Rather, the
proposed method may provide increased resolution and the ability to study different
materials in certain experiments, while photoelasticity may prove to be more practical
in other experiments.

3.6 Conclusion
In this chapter, we have proposed a new method for quantitatively inferring the
evolution of inter-particle forces during a dynamic experiment on granular materials.
The proposed method has been validated using a 2D experimental example in which
inferred forces were compared to an FEM simulation. The results illustrate that
the method can accurately infer the evolution of inter-particle forces at individual
contact points in granular assemblies undergoing dynamic deformation. Preliminary
results from experiments in which the method is applied to study impact in a granular
bed demonstrate a potential application of the method: for unraveling the physics
 57
of model granular materials in a variety of environments. Future applications of the
method may include the study of granular flow, impact, and intruder dynamics with
a variety of complex shaped grains and grain materials. The ongoing development of
3D imaging techniques may soon allow dynamic force inference in real powders and
sands.
 58

Chapter 4

Friction in inertial granular flows:

Competition between dilation and
grain-scale dissipation rates1

4.1 Introduction
Granular flows are ubiquitous in nature and technology [84]. Geologic events such
as landslides and earthquakes occur because granular materials are able to transition
from a solid state to a flowing state. Industrial processes such as hopper flows and
powder transport involve the flow of food and pharmaceutical particles. Defense
applications of brittle ceramics rely on the flow of a pulverized bulk material for
energy dissipation. All of these applications demonstrate the need to understand
granular flows at a fundamental level. Effective friction describes the shear resistance
of a flowing granular medium and, in continuum simulations of these events, encodes
information about grain-scale and contact-scale processes in a single parameter (see
Fig. 4.1).
Granular flows can be classified as quasi-static, inertial (also referred to as dense),
or rapid based on a dimensionless shear rate known as the inertial number [3,11]. The
quasi-static behavior of granular media is typically modeled using critical state soil
mechanics [85]. Rapid granular flows are similar in some respects to gases and have
1
Adapted from R.C. Hurley and J.E. Andrade. Friction in inertial granular flows: Competition
between dilation and grain-scale dissipation rates. Granular Matter, 17(3)287-295, 2015. DOI:
10.1007/s10035-015-0564-2
 59
LENGTH SCALE
μm mm m

contact scale grain scale rve scale field scale

Figure 4.1: Friction encodes contact-scale and grain-scale information in a single

parameter for continuum analysis.

P ˙

Tc
I=
T˙

d 1
Tc = p T˙ =
P/⇢g ˙

Figure 4.2: The two timescales associated with the inertial number, I = Tc /Tγ̇ .

therefore been extensively modeled using kinetic theories (e.g., [86, 87]). The inter-
mediate regime of inertial granular flows has, unlike the quasi-static and rapid cases,
eluded a unified modeling approach. Nevertheless, researchers have made important
progress in understanding inertial granular flows in recent years.
p
The inertial flow regime corresponds to flows with an inertial number, I = γ̇d/ P/ρg ,
between approximately 0.001 and 1. Here, γ̇ is the shear rate (|γ̇| in 3D), d is the
grain diameter, P is the confining pressure, and ρg is the grain density. The inertial
p
number is the ratio of the particle relaxation time d/ P/ρg to the macroscopic shear
time γ̇ −1 , as illustrated in Fig. 4.2 [11, 88]. This interpretation will be revisited when
we develop a new friction relationship in section 4.2.
Researchers studying the inertial flow regime have developed empirical relation-
ships between I and the steady state effective friction coefficient µ, the solid fraction
φ, and the coordination number Z [4, 11]. For instance, da Cruz and co-workers pro-
 60
posed a linear relationship between µ and I given by µ = a + bI for 2D simple shear
flows, where a and b are empirical constants [11]. Jop and co-workers have proposed
the nonlinear relationship µ = µ1 +(µ2 −µ1 )/(I0 /I +1) for 3D flows where µ1 , µ2 , and
I0 are empirical constants [4]. Jop and co-workers have also developed a constitutive
law for predicting the stress distribution and flow profile for well developed granular
flows, using the empirical friction law described above [4]. Other local and non-local
continuum models have been developed for inertial flows. Each model takes advantage
of one of the empirical friction laws described above [89, 90, 91]. These models have
provided promising tools for predicting the behavior of granular media in a variety of
flow configurations. However, investigative studies of inertial granular flows continue
on a more basic level in an effort to understand the processes underlying frictional
rate-dependence and the microstructure that develops during an inertial flow [2].
Da Cruz and co-workers studied the evolution of forces and anisotropy in iner-
tial granular flows, showing that the anisotropy of the contact network can be ex-
plicitly related to friction [11]. Azema and Radjai similarly showed that a classical
stress-force-fabric relation holds for inertial flows, demonstrating another link between
friction and contact network anisotropy [1, 92]. Hatano and Kuwano [93] provided
another interpretation of friction, using an energy balance equation to derive a steady
state friction law very similar to that of rate-and-state theory. Jenkins has provided
interesting links between friction and various attributes of inclined plane flows by
extending hydrodynamic equations valid for rapid flows to the inertial regime [94].
Sun and co-workers have studied the energy characteristics of inertial granular flows
and revealed a number of correlations between the friction coefficient and energy ra-
tios [95]. All of these studies have provided valuable insight into the nature of effective
friction in inertial granular flows.
In this chapter, we intend to contribute an additional interpretation of effective
friction in granular flows by explicitly relating it to the inertial number, the coordi-
nation number, the solid fraction, and grain-scale dissipation rates. In section 4.2, we
develop the friction relationship for steady state simple shear flows by performing an
energy balance and a simple statistical analysis. We discuss the resulting picture of
 61
friction as a competition between dilation and grain-scale dissipation rates. In section
4.3, we discuss numerical simulations of simple shear flows and present results showing
the accuracy of the proposed friction relation. Simulation results are used to illus-
trate how the effective friction coefficient can be decomposed into contributions from
grain-scale dissipation mechanisms. An analysis of the scaling of each of each term in
the friction relationship elucidates the mechanisms controlling rate-strengthening. In
section 4.4, we briefly compare our friction law with others proposed in past research.
Finally, section 4.5 offers concluding remarks.

4.2 The friction law for simple shear

This section provides a derivation of our steady state friction relationship for simple
shear flows. Fig. 4.3a illustrates such a well-developed flow in which the velocity
profile in the direction of flow is quasi-linear.
Similar to past analyses of simple shear flows [93, 96], our starting point is the
energy balance relationship

d
(T + U ) = Dij σji − Γ (4.1)
dt

where T and U are the kinetic and potential energy densities, respectively, Dij =
∂ui /∂xj is the velocity gradient tensor, and Γ is the dissipation rate per unit volume.
In steady state simple shear flows, only one component of the velocity gradient tensor
is nonzero (Dxy = γ̇ in our case), and a time average of Eq. (5.10) yields

γ̇σyx = Γ (4.2)

where the (·) indicates a time-average. Defining the effective friction coefficient as
µ = σ yx /σ yy and assuming dissipation occurs only at grain contact points, Eq. (4.2)
can be rewritten as
Nc hΓc i
γ̇µσyy = (4.3)
V
 62
where Nc is the number of grain contacts in the system, V is the volume of the
system, and hΓc i is the average dissipation rate at grain contacts in the system. We
P c
have made use of the fact that N c=1 Γc = Nc hΓc i.

Numerical simulations to be described in the next section demonstrate that there

is less than a 1% correlation in the fluctuations in the terms in Eq. (4.3). We will
therefore assume averages of products can be written as products of averages and we
will drop all time-averaging symbols. Variables in all following equations should be
assumed to be time-averaged unless otherwise noted. Equation (4.3) becomes

Nc hΓc i
µ= (4.4)
σyy γ̇V

We can further simplify Eq. (4.4) by noting: (1) the number of contact points
is related to the coordination number by Nc = ZNp /2 where Np is the number of
particles in the flow; (2) the number of particles can be related to the solid volume
of grain material Vs by defining d such that Np (4/3)πd3 /8 = Vs ; (3) Vs /V = φ where
φ is the solid fraction. The definition of d to satisfy (2) is consistent with d being
the grain diameter in the case of monodisperse spheres, an average grain diameter in
the case of polydisperse spheres, and a characteristic grain size in flows of complex
shaped grains.
Combining the simplifications described above, Eq. (4.4) becomes

√ !
Zφ 3 ρg
µ= 3/2
hΓc i (4.5)
I πd2 σyy

3/2 √
where the quantity πd2 σyy /(3 ρg ) is a pressure dependent term with units of energy
dissipation rate. We therefore call this quantity Γ̃ and rewrite Eq. (4.5) as

Zφ hΓc i
µ= (4.6)
I Γ̃

Equation (4.6) is the most general form of our friction relationship. This expres-
sion makes no assumption of contact law or grain properties and only imposes the
 63
restrictions that the flow is in steady state and energy dissipation occurs at contact
points. This assumption does not prohibit the incorporation of material plasticity or
fracture so long as such processes are assumed to arise because of contact between
grains. In section 4.3, Eq. (4.6) will be applied to a specific contact law to discuss
the results of numerical simulations.
Before discussing numerical simulations, we can provide a physical interpretation
of Eq. (4.6). The coordination number Z and solid fraction φ convey the connec-
tivity and compactness of the granular material. Both of these quantities, as well
as Zφ/I taken together, can be assumed to decrease during shearing dilation, a pro-
cess by which the packing expands at higher shear rates. In contrast, the grain-scale
dissipation rates hΓc i may be expected to increase with shear rate due to higher
inter-particle forces and collision velocities. The term Γ̃ remains constant when the
confining pressure is held fixed.
The friction law in Eq. (4.6) therefore conveys a competition between dilation
and microscopic dissipation rates. At low shear rates, Zφ/I is large and hΓc i/Γ̃ is
small. At high shear rates, Zφ/I is small and hΓc i/Γ̃ is large. The result of this
competition dictates whether the material is rate-strengthening or rate-weakening
and highlights the role of dilation in effective friction. This interpretation of effective
friction is closely related to the interpretation of the inertial number given in Fig.
4.2. The same time scales at work in the inertial number, those of confinement and
macroscopic shear, are at work in determining the effective friction coefficient. The
confinement time scale dictates the value of Zφ/I while the macroscopic shear time
scale dictates the frequency and intensity of particle collisions.

4.3 Numerical simulations of simple shear

In this section, we discuss numerical simulations of simple shear flows. The simula-
tions explore the behavior of the variables in Eq. (4.6) and elucidate the competing
roles of dilation and microscopic dissipation rates during simple shear at a variety of
imposed shear rates. The numerical code used for the simulations is first discussed,
 64
followed by a discussion of results.

4.3.1 Description of code

We use a discrete element code [33] to study the various components of Eq. (4.6) in the
inertial flow regime. Our simulations use a modified version of the granular module
from the molecular dynamics code LAMMPS [97, 98]. Grains are modeled as spheres
and interact with a Hertzian contact model. The normal force, F n = F m v
n + F n , has a
√
mechanical portion, F m ∗
n = R kn δ
3/2
nij , and a viscous portion F vn = δR∗ meff γn v n ,
p
where R∗ = Ri Rj /(Ri + Rj ), kn is a spring constant, δ is the particle overlap,
nij is a vector from the centroid of particle j to the centroid of particle i, meff =
mi mj /(mi + mj ), γn is a damping coefficient, and v n is the normal component of
the relative velocity vector. The tangential force, F t = min(F m
t , µp |F n |), has a
√
mechanical portion, F m ∗
t = − δR kt ∆s , and a Coulomb slider enforcing |F t | <
m

µp |F n | where kt is a spring constant, ∆s is the accumulated tangential displacement of

the grains, and µp is the inter-particle friction coefficient. The accumulated tangential
displacement of the grains is frozen when grains are sliding. The constants kn and kt
have units of force per area, consistent with the model presented in [99] and discussed
in [97]. Thus, kn and kt are material properties that are independent of grain size and
can be explicitly linked to grain properties if desired [100]. The constant γn prescribes
a velocity-dependent coefficient of restitution of a binary collision, consistent with
experiments [101, 102].
Simple shear is achieved by compressing approximately 10,000 bidisperse spheres
between rough boundaries made of grains and moving the boundaries at a specified
velocity in opposite directions, as shown in Fig. 4.3a. Grain radii are (1±0.2)d˜ where
d˜ is specified. The height of the flow h is chosen such that h/d˜ ≥ 20. The rough
boundaries are moved in the y direction to maintain a constant confining pressure
throughout each simulation. Periodic boundary conditions are used in the x and z
directions.
Confining pressure and grain stiffness is chosen such that (kn /σyy )2/3 ≥ 104 , mak-
 65

yy

(a) vx
|v| = vx

x vx |v| = 0
z
1

0.9
e DS1
DS2
0.8

(b)
0.7
0 0.01 q 0.02 0.03
vc ⇢g /kn

Figure 4.3: (a) A rendering of the simple shear flows featured in this paper. The
top-most and bottom-most particles are used as rough boundaries. Colors indicate
the magnitude of velocity, where vx is the imposed wall velocity.p (b) Coefficient of
restitution e in a two-particle collision with normal velocity vc ρg /kn for data set 1
(- -) and data set 2 (-).
 66
ing the grains “rigid” as described in [11]. The parameter kt is chosen to be 1/2 of kn .
Our primary data (DS1) set features 26 simulations across the inertial flow regime
in which µp = 0.3 and γn is set to prescribe the coefficient of restitution e shown as
the dashed line in Fig. 4.3b. The inter-particle friction coefficient of µp = 0.3 is cho-
sen to provide a balance between lower values found in recent experiments [103] and
slightly higher values used in recent simulations [1, 11]. We found that changing this
inter-particle friction coefficient has minimal qualitative influence on the results, and
mainly acts to shift the µ(I) curve up or down as reported in [11]. A secondary data
set (DS2) features 18 simulations throughout the inertial flow regime with µp = 0.3
and γn set to prescribe e as the solid line in Fig. 4.3b. This data set is only referred to
in order to illustrate how grain viscoelasticity influences the material response. Unless
otherwise specified, data should be assumed to belong to the primary data set. We
leave an in-depth study of the effects of varying µp , γn , particle size distribution, and
contact laws for future work.
Stress is measured using the equation

1 X
Nc
σij = lc f c (4.7)
V c=1 i j

where V is volume, c are contact point labels, Nc is the number of contacts in the
material, lic is a branch vector pointing from the centroid of particle j to the centroid
of particle i, and fjc is the force vector from particle j to i [104]. The effective friction
coefficient is computed using µ = σ yx /σ yy , where averages are carried out over several
thousand stress calculations once a steady state velocity profile has been reached.
Figure 4.4 compares µ and φ found in our simulations (using DS1) with available
data from contact dynamics simulations [1] and 3D annular shear cell experiments [2,
3]. Our simulations show an excellent collapse with the other data sets both in terms
of effective frictional response and solid fraction. The effective friction coefficient
increases from its quasi-static value throughout the inertial flow regime, approaching
a plateau at the transition to the rapid flow regime. The solid fraction decreases
approximately linearly throughout the inertial flow regime from a maximum quasi-
 67
(a) (b)
0.7 0.6
DEM Data DEM Data
Azema and Radjai (2014) Azema and Radjai (2014)
Savage and Sayed (1984) 0.56
0.6
0.52
0.7 0.6

0.5
µ

φ
0.56
0.6 0.48
0.52
0.5

φ
0.4 0.44 0.48
0.4
0.44

0.3 −3
10 10
−2 −1
10
0
10
0.4 0.4 −3
10 10
−2 −1
10
0
10
I I
0.3
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
I I

Figure 4.4: A comparison of (a) effective friction and (b) solid fraction from our
simulations and available data sets taken from the literature. Blue squares are from
contact dynamics simulations [1]. Black triangles are from 3D annular shear cell
experiments [2,3]. Error bars indicate standard deviations in time of measured quan-
tities. The dashed line shows the fit µ = µ0 + (µ1 − µ0 )/(1 + I0 /I) from [4].

static value of 0.59.

In Fig. 4.4, as in all figures in this paper, plotted quantities are obtained as follows.
First, a quantity of interest (e.g., coordination number) is computed at periodic times
(approximately 5 × 104 times) once steady state flow has been achieved. The total
strain over which quantities are extracted is taken such that averages over larger
strains do not change the results. Next, the average of these quantities is used to
obtain the plotted data points. Finally, the sample standard deviation is used to
obtain the error bars. Error bars are typically omitted from inset plots for clarity.

4.3.2 Results: Validity of friction law

Figure 4.5 displays the effective friction coefficients for our primary data set computed
using the friction relationship in Eq. (4.6) and the stress formula in Eq. (4.7).
The figure demonstrates that the proposed friction law excellently approximates the
effective friction coefficient throughout the inertial flow regime. We have confirmed
that a similarly accurate fit exists for other grain properties.
 68
0.7
Stress Formula, Eq. (7)
Friction Law, Eq. (6)

0.6

0.7

0.5

µ
0.6

0.5

µ
0.4 0.4

0.3 −3 −2 −1 0
10 10 10 10
I
0.3
0 0.2 0.4 0.6 0.8 1
I

Figure 4.5: A comparison of the effective friction coefficient calculated from proposed
friction relationship in Eq. (4.6) and the stress formula in Eq. (4.7).

4.3.3 Results: Grain-scale dissipation mechanisms

The contact law discussed in section 4.3.1 implies that Eq. (4.6) can be written as

Zφ hΓn i + hΓs i
µ= (4.8)
I Γ̃

where hΓn i = hF n · v n i is the average viscoelastic dissipation rate, averaged over

all contacts, and hΓs i = h|F t ||v t |i is the average dissipation rate from grain sliding,
averaged only over sliding contacts. An additive decomposition of Eq. (4.8) into
µ = µn + µs yields

Zφ hΓn i Zφ hΓs i
µn = and µs = (4.9)
I Γ̃ I Γ̃

The effective friction coefficient can thus be written in a form that clearly decouples
contributions from the two grain-scale dissipation mechanisms, viscoelasticity and
grain sliding.
Figure 4.6a illustrates how the two terms in Eq. (4.9) evolve throughout the
inertial flow regime. At low shear rates, dissipation from grain sliding is the primary
 69
contributor to effective friction. At higher shear rates, the contribution from grain
sliding remains constant or declines as the contribution from viscoelastic dissipation
becomes increasingly prominent.

4.3.4 Influence of material properties

In order to highlight how grain properties influence the relative contributions of mi-
croscopic dissipation mechanisms, results from the secondary data set are also shown
in Fig. 4.6 (the solid lines). We recall that the primary data set (DS1) features the
same grain properties as the secondary data set (DS2) except for a lower coefficient
of restitution, as shown in Fig. 4.3b.
Figure 4.6a compares how the two terms in Eq. (4.9) evolve as a function of
shear rate for each data set. Compared to DS1, DS2 features a larger effective fric-
tion contribution from grain sliding and a smaller contribution from viscoelasticity
throughout most of the inertial flow regime. This occurs because grain viscoelastic-
ity dissipates less energy for a given particle collision in DS2, leaving more kinetic
energy in the system to be dissipated by grain sliding. Despite the difference in the
grain-scale contributions to effective friction, both data sets feature similar values for
µ until I ≈ 0.4, as shown in Fig. 4.6b.
In past work [1,11], some researchers have ignored the influence of the coefficient of
restitution e in shear flows because of the similarity of the effective friction coefficient
when measured in systems using different values of e. Simulations are therefore often
carried out using e = 0 [1]. The finding in Fig. 4.6 illustrates that although effective
friction may be similar, the grain-scale dissipation mechanism responsible for friction
is different in systems with different values of e. It may be interesting to explore the
range of grain-scale behaviors that emerge from varying µp and e and the reason that
they have such a minor influence on effective friction below I ≈ 0.4. We leave such
an investigation as future work.
Increasing or decreasing the inter-particle friction coefficient has the effect of shift-
ing the effective friction curves in Fig. 4.6b up or down, respectively, but does not
 70

(a)
0.4

0.3 µn
µs
DS1
0.2 DS2
7i

0.4

0.3

0.1 0.2

7i
0.1

0
-3 -2 -1 0
10 10 10 10

0
I

0 0.2 0.4 0.6 0.8 1

(b) I
0.7
DS1
DS2
0.6

0.5
7

0.7

0.6

0.5
0.4
7

0.4

0.3
-3 -2 -1 0
10 10 10 10
I
0.3
0 0.2 0.4 0.6 0.8 1
I

Figure 4.6: (a) The two terms µn = ZφhΓn i/(I Γ̃) and µs = ZφhΓs i/(I Γ̃) in the
additive decomposition of effective friction given in Eq. (4.9) for the primary data set
(- -) labeled DS1 and the secondary data set (-) labeled DS2. (b) The total effective
friction for the two data sets as a function of I. Error bars are omitted from inset
plots for clarity.
 71
significantly affect their shape or the magnitude of viscoelastic dissipation. Thus,
the inter-particle friction coefficient primarily sets the baseline magnitude of effective
friction while the coefficient of restitution controls the grain-scale contributions.

4.3.5 Results: Dilation and dissipation rates

The evolution of each term in Eq. (4.6) is shown in Fig. 4.7. Shearing dilation is
captured in the evolution of Z, φ, and Zφ/I in Figs. 4.7a, 4.7b, and 4.7d, respectively.
At low shear rates, Z and φ maintain maximum quasi-static values that depend upon
properties such as the inter-particle friction coefficient and particle shape. For our
primary data set, these quasi-static values are approximately 4 and 0.59 for Z and φ,
respectively. As shear rates increase throughout the inertial flow regime, both Z and
φ decrease as the material dilates. At all shear rates investigated, Z is well described
by
Z2
Z ≈ Z1 + (4.10)
b+I
and φ is well described by
φ ≈ φmax − mI (4.11)

where Z1 , Z2 , b, φmax , and m are constants. These approximations hold for both data
sets set discussed here and for simulations with different values of µp and γn , and
different particle size distributions which we do not discuss here.
When combined, Eqs. (4.10) and (4.11) suggest the two scaling regimes of Zφ/I
shown in Fig. 4.7d. The scaling Zφ/I ∝ I −1 arises because Z and φ maintain quasi-
static values at the low end of the inertial regime. The scaling Zφ/I ∝ I −2 arises
because Z and φ decrease in agreement with Eqs. (4.10) and (4.11) at higher shear
rates. This decreasing contribution of Zφ/I in Eq. (4.6) reflects a decrease in both
number of contact points and total solid fraction as shearing dilation increases. From
an energy perspective, this decay conveys the decrease in internal surface area over
which the material can dissipate energy.
While Zφ/I decreases with shear rate due to shearing dilation, average grain-scale
 72

(a) (b)
4 4 0.6 0.6
3

2
0.56

φ
0.56

Z
1
0.52

0.48 −3 −2 −1 0
−3 −2 −1 0
10 10 10 10 10 10 10 10
2 0.52
Z

φ
I

1 0.48

0 0.44
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
I I
(c) (d)
2
101
104
0
10
10-1
10
-2
103
10-3
h!i i=!
~

10-4 1
I
10-5
102
h!i i=!
~

Z?=I

-6
10
1 10-7
10-8 -3
10 10-2 10-1 100 101 2
I I

˜ 100
h n i/
h s i/
˜
0 10-1 -3
0 0.2 0.4 0.6 0.8 1 10 10-2 10-1 100
I I

Figure 4.7: (a) The coordination number Z as a function of inertial number. The
dashed line is the fit from Eq. (4.10). (b) The solid fraction φ as a function of inertial
number. The dashed line is the fit from Eq. (4.11). (c) The average grain-scale
dissipation rates hΓn i/Γ̃ and hΓs i/Γ̃ as a function of inertial number. The dashed
lines are power-law fits, with hΓn i/Γ̃ ∝ I 2.4 and hΓs i/Γ̃ ∝ I 1.87 (d) The quantity
Zφ/I as a function of inertial number. The dashed lines represent the two regimes
of behavior in which Zφ/I ∝ I −1 and Zφ/I ∝ I −2 . Error bars in (d) are negligibly
small.
 73
dissipation rates increase as shown in Fig. 4.7c. Both viscoelastic and grain sliding
dissipation rates approximately follow a power-law dependence on I throughout the
inertial flow regime, with hΓn i/Γ̃ ∝ I 2.4 and hΓs i/Γ̃ ∝ I 1.87 , as shown in the inset
of Fig. 4.7c. We generally expect hΓn i/Γ̃ to scale at least as fast as I 2 for our
chosen contact law because collision velocities scale approximately with I and viscous
normal contact forces scale with collision velocity. Surprisingly, we also find hΓn i/Γ̃
to scale at least as fast as I 2 when we impose sub-linear dependence of viscous normal
forces on grain collision velocities. This likely occurs because the correlation between
collision velocities and viscous normal contact forces implies that hΓn i ∝ hF vn · v n i >
hF vn i hv n i ∝ I 2 .
We have not found a similar argument for the scaling of hΓs i/Γ̃ with shear rate,
but we have always observed this term to scale slower than hΓn i/Γ̃.

4.3.6 Results: Rate-dependent friction

When combined, the competing processes of shearing dilation and grain scale dissi-
pation rates give rise to a rate-strengthening effective friction coefficient in our data
sets, as shown in Fig. 4.6b. Rate-strengthening seems to occur because Zφ/I never
decays faster than I −2 while hΓn i/Γ̃ increases at least as fast as I 2 in the inertial
regime.
Data sets using other values for µp and γn , as well as other simple contact laws
(linear springs or nonlinear dependence of viscous normal force on collision velocity)
have been investigated and yield similar results: Zφ/I never decays faster than I −2
and hΓn i/Γ̃ always increases at least as fast as I 2 in the inertial regime. Velocity-
strengthening therefore appears to be a generic system response for bidisperse spheres
interacting with many viscoelastic contact models in the inertial regime. We have
not observed a transition to rate-weakening friction at low shear rates as observed in
some recent experiments [93]. Given the variety of grain properties and particle size
distributions we have studied (but not discussed here), we suspect such a crossover,
if it exists, to be caused by processes not captured by the current viscoelastic contact
 74

101
100
10-1 I 1.87
10-2 I 2.4
10-3

h!i i=!
~
10-4 I
10-5 I
-6 ˜
10 h n i/
-7 ˜
10 h s i/

10-8 -3
10 10-2 10-1 100
I

Figure 4.8: The inset of Fig. 4.7c, showing the average grain-scale dissipation rates
as a function of inertial number on a log-log scale. The dashed lines illustrate the
scaling of each dissipation rate discussed in the text, as well as the scaling proportional
to I at the transition to quasi-static flow.

model, such as flash heating [93].

In Fig. 4.7c, dissipation rates approach 0 as I → 0. This seems to suggest effective
friction also approaches 0 as I → 0; however, I appears in the denominator of Eq.
(4.6). The scaling of dissipation rates with I is therefore the quantity controlling
the approach to a quasi-static value of effective friction, not the absolute value of
dissipation rates. To clarify this, we enlarge the inset of Fig. 4.7c in Fig. 4.8. At the
transition from inertial to quasi-static flow, close to I = 10−3 , both dissipation rates
trend toward scaling as hΓi i/Γ̃ ∝ I. In the quasi-static regime below I = 10−3 , both
dissipation rates reach this scaling. When this occurs, the effective friction coefficient
reaches a quasi-static value and the strength is rate-independent. This transition
is significant for understanding the onset of rate-dependent behavior and deserves
further investigation in future work.
 75
4.4 Discussion
The primary finding of sections 4.2 and 4.3 are: (1) effective friction in inertial gran-
ular shear flows can be interpreted as a competition between shearing dilation and
grain-scale dissipation rates; and (2) rate-strengthening effective friction occurs in the
inertial flow regime for many viscoelastic contact models and grain properties because
grain-scale dissipation rates win the competition with shearing dilation. Although the
first finding may seem intuitive, Eq. (4.6) provides a quantitative means of studying
it. Equation (4.9) also provides a method of tracing macroscopic frictional energy
loss down to the grain-scale.
The interpretation of steady state effective friction given in this chapter comple-
ments interpretations offered by past work [1, 11]. In these past works, researchers
explicitly linked friction to the anisotropy of the contact network. An increase in
normal force anisotropy was implicated in causing rate-strengthening. This inter-
pretation offers an understanding of how the structural organization of grains may
influence macroscopic frictional response. The interpretation of friction in the current
chapter, however, makes no mention of structural organization and rather relies upon
energy dissipation to explain changes in the macroscopic frictional response. That
each interpretation can independently explain the change in effective friction as a
function of shear rate is interesting and deserves further investigation.
Several other interesting behaviors were observed during the simulations carried
out for this work. Many of these behaviors warrant future investigation. First,
a similarity in the effective friction coefficient was observed for simulations using
different values for the coefficient of restitution e. This similarity persisted even
while differences were observed in the grain-scale contributions to friction (see Fig.
4.6). This behavior has been noted by previous authors [1, 11] and has not yet been
explained. A number of interesting changes in behavior were also observed to occur
near I ≈ 0.1. Authors in [1] noted a number of topological transitions occurring
at this shear rate. In this chapter, we observe the transition from Zφ/I ∝ I −1 to
Zφ/I ∝ I −2 close to I ≈ 0.1. We also observe a dramatic increase in the fraction
 76
of sliding contacts at I ≈ 0.1, regardless of the inter-particle friction coefficient µp
used in simulations. Finally, we have observed that Eq. (6) holds for other contact
laws and holds locally in layers of roughly constant inertial number in systems and
flow configurations where the inertial number varies spatially. It would be useful
to investigate such systems in more detail in future work to determine whether the
quantitative results of this chapter persist.

4.5 Conclusion
We have presented a relationship between steady state effective friction, the inertial
number, coordination number, solid fraction, and grain-scale dissipation rates in a
granular shear flow. This relationship elucidates the rate- and porosity-dependent
nature of effective friction in granular flows. Numerical simulations of simple shear
flows have been used to illustrate how effective friction is furnished by grain-scale
dissipation mechanisms. Rate-strengthening was seen to occur because terms encom-
passing shearing dilation decay more slowly with shear rate than terms encompassing
grain-scale dissipation rates increase. We discussed how our findings compare with
other interpretations of effective friction and mentioned several observations that
warrant future investigation.
 77

Chapter 5

Modeling dilative viscoplastic

granular flows using SPH1

5.1 Introduction
Many geologic and industrial processes involve large deformation and flow of granular
media such as soils, sands, powders, and pulverized brittle materials. Examples of
such processes include landslides, debris flows, asteroid impact and cratering, bulk
food transport, pharmaceutical processing, and ballistic impact of ceramics. Numeri-
cal models that can accurately simulate these processes are particularly important as
predictive tools. However, the lack of well-established rheological relations for flowing
granular media and the difficulty of representing arbitrarily large deformations with
traditional numerical methods makes predictive modeling of these processes difficult.
The rheological behavior of flowing granular media is complex and has been a
subject of study for many decades. Dry granular materials can exhibit solid-like,
liquid-like, and gas-like behavior depending upon their excitation [8]. The ability
of granular media to accommodate additional liquids and gases in their pore space
further complicates their flow behavior [105]. Through decades of experimental and
numerical investigation, the constitutive behavior of granular materials in several
flow configurations has been elucidated (e.g., see [3, 4, 86, 88]). Particularly relevant
to the present chapter, the rheology of dry granular flows has been examined in flows
1
Adapted from R.C. Hurley and J.E. Andrade. Modeling dilative viscoplastic granular flows
using SPH. Under review.
 78
down inclined planes ( [4,106]), simple shear flows [11], slow flows [85,107], and rapid
flows [86]. Authors in [4] proposed a rate-dependent constitutive law for steady-
state granular flows down inclined planes that has since proven to be an accurate
framework for modeling large-scale processes such as impact cratering and granular
column collapses (e.g., see [89, 108]). The rate-dependent nature of friction in this
constitutive law has been examined analytically and shown to have grain-scale origins
(e.g., see [10]). This rheological picture provides one of the most unifying frameworks
for granular flows to-date. Below, we will adopt part of this constitutive law to model
granular media in both their solid-like and liquid-like states.
Numerous numerical methods have been used for modeling granular materials.
Depth-averaged finite difference methods have been a popular means of predicting
geological flows for many decades (e.g., see [109]). The finite element method (FEM)
has also been a popular tool for modeling both the quasi-static and dynamical flow
of granular media using classical constitutive laws for soils and newer viscoplastic
constitutive laws similar to the one developed by [4] (e.g., see [110] and [111]). More
recently, the material point method (MPM) has been used to model the dynamic
behavior of granular media using a Matsuoka-Nakai constitutive framework (e.g.,
see [112]). Smoothed Particle Hydrodynamics (SPH) has also been used to model
the dynamic behavior of granular media using a Drucker-Prager plasticity law (e.g.,
see [113,114]). While each of these methods has strengths, mesh-free methods (MPM
and SPH) are particularly attractive for very large and rapid material deformations
since they do not require mesh-refinement. Furthermore, Lagrangian mesh-free meth-
ods that do not require a background grid (SPH) eliminate the convective effects of
classical Eulerian formulations. Compared to Eulerian methods, SPH also offers pre-
cise interface definition and arbitrary resolution of flow details.
In this chapter, we present an SPH approach to modeling the flow of dry granular
media. The novelty of our approach comes from combining an SPH formulation
with a rate-dependent constitutive law similar to that proposed by [4] and used by
[107]. We augment the rate-dependent constitutive law used in [4] and [107] with the
rate of dilation and justify this approach with an analysis of the continuum energy
 79
balance equation. We then present four qualitative and quantitative examples of the
method. In particular, we model angle of repose tests, steady-state rate-dependent
flows down inclined planes, granular column collapses on flat surfaces, and granular
column collapses down inclined planes. These examples provide a partial validation of
the method against experimental results and illustrate that, with some future work,
it may be used as a predictive tool for many processes involving granular media.
Compared to past mesh-free simulations of granular media (e.g., [112, 113, 114]), our
method employs a new constitutive law well-suited for rate-dependent behavior, and
is implemented three dimensionally.
The layout of this chapter is as follows. In section 5.2, we describe the balance and
constitutive laws we use to model granular materials. We present our rate-dependent
constitutive law in this section. We also show in this section that analyzing the
continuum energy balance equation in the context of our constitutive law suggests
a relationship between the rate of dilation and the effective friction coefficient. Sec-
tion 5.3 is dedicated to describing the SPH framework used in this chapter and the
algorithm we employ in our simulations. Four numerical examples of our method
are described in section 5.4. We first discuss simulations of slumping granular media
that demonstrate the ability of our method to produce angles of repose consistent
with the friction angle used in our model. We then discuss simulations of inclined
plane flows that illustrate the rate-dependent feature of our numerical method and
the ability to match Bagnold velocity profiles. Next, we simulate granular column
collapses on flat surfaces to illustrate the ability of our method to produce scaling
laws of slump and runout consistent with those found experimentally. Finally, we
compare time-dependent profiles of granular column collapses down inclined planes
to demonstrate that our method accurately models the dynamic structure of granular
flows. Section 5.5 offers a discussion of future work and concluding remarks.
 80
5.2 Balance and constitutive Laws

5.2.1 Balance law

The governing balance law we will solve is the equation of momentum balance given
by

∇ · σ + ρb = ρa (5.1)

where σ is the Cauchy stress tensor, ρ is the mass density, b is a body force, and a
is the material acceleration.

5.2.2 Constitutive law

Defining the spatial velocity gradient as L = ∇v, we introduce the strain rate tensor
as

D = L + LT (5.2)

where the T indicates a transpose. This definition differs by a factor of two from the
classical definition in order to be consistent with [4], whose constitutive law we now
adopt. We propose a stress tensor similar to that used for yield-stress (e.g., Bingham)
fluids

(µp + c)D
σ = −pI + (5.3)
||D||
where p is pressure, I is the identity tensor, µ is the friction coefficient, c is cohesion,
and ||D|| = ( 21 D : D)1/2 is the second invariant of the strain rate tensor. This form
of the stress tensor is similar to that proposed by [4], the difference being the addition
of a cohesion term.
The second term in Eq. (5.3) can be identified as the shear stress tensor τ . When
D is nonzero, the material flows with shear stress magnitude ||τ || = µp + c where
 81
||τ || = ( 21 τ : τ )1/2 . Thus, the yield criterion is given by

||τ || = µp + c (5.4)

This form is again similar to the Drucker-Prager-type yield criterion in [4], but with
the addition of cohesion.
The pressure in the granular material p is derived from an equation of state relating
pressure to density. The general form of our equation of state is reminiscent of Tait’s
equation of state for nearly-incompressible fluids and is given by
   
γ

κ ρ
ρ0
− 1 ρ ≥ ρ0
p= (5.5)

0 ρ < ρ0

where κ and γ are parameters and ρ0 is the critical “jamming” density of the granu-
lar material, or the loosest packing density at which the material supports a nonzero
stress state. We choose γ = 3/2 to produce a pressure-density relationship consistent
with that found at the jamming transition for granular solids in [115]. The param-
eter κ is then constrained by the desired bulk modulus of the material through the
relationship
  
3/2

3κ ρ
dp 2 ρ0
ρ ≥ ρ0
K=ρ = (5.6)
dρ  
0 ρ < ρ0

We use a rate-dependent and dilation-dependent friction coefficient given by

µh − µl
µ = µl + ∗
+β (5.7)
D /||D|| + 1
where µl is the friction coefficient at low shear rates, µh is the asymptotic friction
coefficient at high rates, D∗ is a hardening or softening strain rate scale, and β is
related to the rate of dilation by

1
2
tr(D)
β= (5.8)
||D||
 82
where tr(D) is the trace of D. Note that the friction coefficient given by Eq. (5.7)
can be either rate-strengthening or rate-weakening depending upon the chosen values
of µl and µh . In granular materials, the classical choice is for a rate-independent
or rate-strengthening friction coefficient. The theoretical link between the friction
coefficient and the rate of dilation β is discussed next.

5.2.3 Theoretical link between friction and dilation rate

Several classical stress-dilatancy relationships link the frictional strength of a granular

material to expansion or contraction (e.g., see [85]). One example is that of [116] who
proposed a stress-dilatancy relationship of the form

tan φ0m = µ + tan ψ (5.9)

where φ0m is the externally mobilized friction angle, µ is a frictional constant represent-
ing local internal material strength, and ψ is the dilation angle. Similar relationships
have been used extensively in modeling quasi-static soil behavior (e.g., see [107]).
Recently, a friction law linking dynamic friction and the rate of dilation has also
been derived for simple shear of a granular material ( [15]). As with the the deriva-
tion in [116], this friction law emerges from analysis of the governing energy balance
equation. Below, we give a similar derivation but for three-dimensional flow of a
material with a stress tensor given by Eq. (5.3). This derivation is meant to establish
a theoretical justification for including the dilation rate β in Eq. 5.7.
Consider a granular material undergoing deformation at a rate D. In the absence
of heat sources, local energy balance requires that

1
ρė = σ : D − ∇ · q (5.10)
2
where e is the specific internal energy, q is the heat flux vector, and the factor of 1/2
accounts for the difference between our definition of the strain rate tensor given in
Eq. (5.2) and the classical definition. The strain-rate tensor can be decomposed into
 83
volumetric and deviatoric components such that

D = Dv + Ds (5.11)

where D v = tr(D)I/3 and D s = D − D v . Note that tr(D s ) = 0 by definition and

thus D v : D s = 0. Substituting Eqs. (5.11) and (5.3) into (5.10) and rearranging,
while ignoring cohesion, therefore yields

1
− p tr(D) + µp||D|| = ρė + ∇ · q (5.12)
2
Solving this equation for µ yields

1
ρė ∇·q 2
tr(D)
µ= + + (5.13)
p||D|| p||D|| ||D||
At steady-state, granular materials are observed to deform at constant volume (
[4,85]) and constant average internal energy. Therefore, sufficiently long time averages
of the first and third terms in Eq. (5.13) are zero at steady-state. The second term
in Eq. (5.13) is then the rate-dependent steady-state friction coefficient represented
by the first two terms in Eq. (5.7). The first term in Eq. (5.13) is a transient term
associated with changes in internal energy of the material that accompany non-steady
deformations (e.g., see [15]). In this work, we neglect representing this term in our
friction law Eq. (5.7), making the assumption the time required for local internal
energy equilibrium is sufficiently short that it does not play a significant role in the
problems we wish to solve. Relaxing this assumption in future work will be discussed
in section 5.5. The third term in Eq. (5.13) is identified as the rate of dilation of the
material and is represented directly in Eq. (5.7).
In the examples presented in this chapter, the value of β is limited to 0 ≤ β ≤ 0.5.
This choice is meant to accomplish three things. First, these limits eliminate the
influence of contraction of SPH particles (β < 0) near boundaries that can reduce
friction and result in excessive slip. Second, we wish to capture the role of dilation in
a similar manner to past models (e.g., see [107]), which use it primarily as a means
 84
of enhancing the ultimate strength of a yielding granular material. Third, we wish to
eliminate the effects of free-floating particles (i.e., those with p = 0) on the material
dilation. One way of accomplishing this is with the current limits, which restrict β
to values we perceive as more physically realistic and important as highlighted by
past work ( [107]). We will highlight the consequences of limiting β in this manner in
section 5.3.4 and examples 1 and 2 in sections 5.4.1 and 5.4.2. We will discuss future
extensions of our treatment of dilation in sections 5.3.4 and 5.5.

5.3 SPH formulation and algorithm

5.3.1 Kernel function and basic equations

We use SPH to model granular materials. Lagrangian particles of fixed mass repre-
sent a fixed mass of granular material and move according to the governing balance
and constitutive laws. All field quantities are expressed at these particle locations us-
ing summation interpolants and a kernel function W with smoothing length h. The
theory and history of SPH, a description of various kernel functions and summation
interpolants, and a discussion of special considerations such as free surfaces and in-
stabilities can be found in various monographs on the topic (e.g., see [117, 118, 119]).
Here, we merely state that classical SPH has a number of attractive properties, in-
cluding zero intrinsic dissipation, exact conservation of mass, momentum, angular
momentum, energy, and entropy. We also note that resolution follows the chosen
SPH particle mass.
We use the classical cubic spline kernel in all summation interpolants ( [120])




 1 − 32 q 2 + 34 q 3 0≤q<1


1
W (r, h) = 3 1
(2 − q)3 1≤q<2 (5.14)
h π
 4



0 q≥2

where q = r/h and r = |xa − xb | is the distance between two particles labeled a and
b. Equation (5.14) is normalized for three dimensions since all simulations in this
 85
chapter are three-dimensional. One-dimensional and two-dimensional normalizations
can be found elsewhere ( [117]). Alternative higher-order kernels may be implemented
in the future to reduce instabilities (e.g., see [121]), but were not required for the
examples presented here. The smoothing length h is kept constant in this work since
the modeled granular media is fairly incompressible. The length h is chosen to be
1.2 times the spacing of SPH particles when the granular material is in its loosest
packing state. The fairly incompressible nature of the material prevents significant
particle concentrations that might cause instabilities discussed in other works (e.g.,
see [122]).
The basic SPH interpolants used throughout this work include the summation
interpolant for a field variable f at a particle location a given by (e.g., see [117])

X mb
fa = fb W (r, h) (5.15)
b
ρb

and the interplant for the first derivative of a field variable, such as the pressure,
given by (e.g., see [123])

X  
pb pa
∇pa = − mb ρa + ∇W (r, h) (5.16)
b
ρ2b ρ2a

where m is mass and ∇W (r, h) = (xa − xb )F , F being a scalar symmetric function

in a and b.

5.3.2 Density

We use the basic SPH equation for density with a Shepard filter applied at each time
step. For particle a, we first compute ( [117])

X
ρ̃a = mb W (r, h) (5.17)
b∈N

where the set N includes neighboring particles in the support domain of a as shown
 86
W Interior particles
va
da a

a 2h Boundary db
b vb b
Boundary particles
(a) (b)

Figure 5.1: (a) Illustration of the smoothing kernel W with compact support of radius
2h about particle a. (b) Interior particles interacting with boundary particles, which
are given artificial velocity.

in Figure 5.1a. We then apply the Shepard filter to obtain

P
mb W (r, h)
ρa = Pb∈N mb (5.18)
b∈N ρ̃b W (r, h)

This filter corrects particle deficiencies near free surfaces and boundaries and is equiv-
alent to modifying the kernel to satisfy the partition of unity

X mb
W (r, h) = 1 (5.19)
b∈N
ρb

We note that we have also employed the more common and efficient free surface
density formulation proposed by [124]. While we obtain similar results in many cases
using this formulation, we use Eq. (5.18) instead in order to retain high accuracy in
density near the free surfaces. Future work may find this to be unnecessary.

5.3.3 Equation of state

After the density ρa at each particle has been computed using Eq. (5.18), pressure is
calculated by direct application of Eq. (5.5).
 87
5.3.4 Equation of motion

To update SPH particle locations, we discretize the momentum balance relationship

Eq. (5.1) to give

X  σa σb

aa = ρa + ∇W (r, h) + b (5.20)
b∈N
ρ2a ρ2b

where σ is given by Eq. (5.3). In computing Eq. (5.3) before the calculation of Eq.
(5.20), the strain rate tensor D must be calculated using Eq. (5.2) and the velocity
gradient computed by

X mb
La = (v b − va ) ⊗ ∇W (r, h) (5.21)
b∈N
ρb

The value of p used in σ is given by applying the equation of state to the current
density computed using Eq. (5.18). When rate-dependent friction is employed, the
value of µ used in σ is computed with Eq. (5.7) applied to the particle velocities at
the last time step.
When the rate of dilation β is included in calculating the value of µ from Eq.
(5.7), β should be calculated by computing the third term in Eq. (5.13) for each
particle. Unfortunately, doing this from SPH continuum fields in a simulation with
arbitrary resolution (particle size and time step) may introduce artifacts due to tem-
poral fluctuations in D that are time-step-dependent and spatial fluctuations in D
that are particle-size-dependent. A more appropriate approach that recognizes β as
arising from processes occurring “within” each SPH particle (i.e. inherently below the
resolution of the continuum method) would be to evolve the rate of dilation β of each
SPH particle by it’s own evolution equation. This could be done with a source-decay
differential equation that increases β due to a source (i.e. rapid change in strain rate
of the particle) and allows β to decay to zero as a function of strain when strain rate
is roughly constant. For example, such an evolution equation could take the form

dβ β
=− ∗ +S (5.22)
ds s
 88
where s is total strain, ∗s is a scale of decay, and S is a source term. This method
for evolving β may be the focus of future work and may embrace a hierarchical
multiscale approach like that proposed in [125]. However, since data to calibrate this
method are lacking in the literature and since the present chapter is more focused on
qualitative demonstration of the SPH framework and constitutive laws, we use the
former approach of computing β directly from the SPH D field. We note that this
choice introduces a few numerical artifacts into the solution. Most importantly, β
is not seen to be identically zero when averaged during steady-state flows like those
encountered in example 2 in section 5.4.2, primarily because of fluctuations in D due
to SPH particle disorder, the choice of ignoring data outside the limits 0 ≤ β ≤ 0.5,
and the finite non-local smoothing effect employed in SPH. We will discuss this topic
and needed future work more in sections 5.4.2 and 5.5.

5.3.5 Boundary conditions

In all simulations shown in this chapter, SPH particles are used to represent solid
boundaries. Figure 5.1b shows a configuration of these “boundary” particles rela-
tive to the “interior” particles comprising the material of interest. The density of
boundary particles is evolved using Eq. (5.18) with a summation over all neighboring
particles, including those on the boundary and interior. Interior particles similarly
use all neighboring particles in their density calculation. Interior particles also use all
neighboring particles in their velocity gradient and stress tensor calculation. A basal
friction coefficient µb may be used in place of the first two terms in Eq. (5.7) when
an interior particle interacts with a boundary particle.
Boundary particles are given an artificial velocity when interacting with interior
particles in order to simulate a smooth velocity gradient across the boundary. As
discussed in [121], when interacting with interior particle a, a boundary particle b
will have an artificial velocity of

v b = (1 − β)v a + βv w (5.23)
 89
where v w is the imposed velocity of the boundary and

 
db
β = min βmax , 1.0 + (5.24)
da
We choose a value of βmax = 1.5, consistent with past work ( [113, 121]). The lengths
db and da are shown in Fig. 5.1b.
As discussed in [113] and [126], the artificial velocity of boundary particles is
insufficient to provide a no-slip condition; artificial stress values must also be assigned
to boundary particles. In the present work, we do not assign artificial stress values
to boundary particles and instead directly compute the stress values of boundary
particles as we do with interior particles. The effect of this is to reduce the shear stress
gradient for interior particles close to boundaries, thus allowing a certain about of slip
across the boundary. We note that this choice does not reduce the density, pressure,
or velocity gradient of interior particles near the boundary since these quantities are
computed by directly including boundary particles in the appropriate summations.
Boundary slip is a key feature of granular systems that distinguishes them from
idealized fluid systems. Future work may compare our boundary condition with that
obtained by imposing a no-slip condition in various scenarios. An alternative to both
the present approach and the no-slip approach used in [113] would be to assign an
artificial stress value to boundary particles, the value of which is between the value
computed by evolving this quantity directly and the average stress value that would
be assigned artificially from interior particles. In the present work, we simply quantify
the amount of boundary slip in an example of inclined plane flow.
There are several disadvantages of representing boundaries using real SPH parti-
cles. One disadvantage is the inability to represent curved surfaces when the simula-
tion resolution is low. Another disadvantage is the potential for interior SPH particles
to penetrate a boundary composed of SPH particles. We have found that modeling
boundaries using real SPH particles with the same parameters as interior particles
sometimes permits penetration of SPH particles into the boundary. One solution to
this problem was proposed by [127], who modified particle velocities by an average of
 90
velocities nearby. To avoid modifying bulk dynamics, we instead choose to alter the
parameter ρ0 used in calculating the pressure of boundary particles only. We set ρ0
to approximately 3% below the initial packing density of boundary SPH particles in
order to create a higher pressure p of stationary boundary particles and thus a higher
pressure gradient for interior particles that interact with the boundary. This modifi-
cation is not performed on interior particles, whose value of ρ0 is maintained at the
initial packing density. The modification of the boundary ρ0 value slightly decreases
the pressure of the interior particles closest to the boundary. Bulk dynamics are not
noticeably affected; the modification only serves to increase the pressure gradient of
interior particles closest to the boundary and therefore prevent penetration of interior
particles through boundary particles. We illustrate this effect in section 5.4.1.

5.3.6 Time integration

We use a Velocity Verlet time integration scheme equivalent to the time-synchronized

Leapfrog method. This symplectic integrator has been shown to conserve momenta
and energy, making it particularly attractive with SPH which possesses conservation
properties ( [117]). This approach involves the following updates for particle a’s
position and velocity

xn+1
a = xna + v na ∆t + 21 ana ∆t2 (5.25)
n+1
an
a +aa
v n+1
a = v na + 2
∆t (5.26)

As an intermediate step, the half-step velocity can be written as

ana
v n+1/2
a = v na + ∆t (5.27)
2
The acceleration an+1
a is then computed using Eq. (5.27) and xn+1 . The precise
update steps are described in Algorithm 1 below. We note that all field quantities
ρ, p, D, and σ are computed at each time step rather than their values being updated
 91
rom previous time steps. In this manner, the procedure is simpler than traditional
implementations of plasticity, which require a return mapping algorithm. The absence
of a mesh eliminates the need for a mesh-updating, re-meshing, or mesh-mapping step
characteristic of finite element or material point method implementations.
The time-step ∆t must be chosen to satisfy a Courant condition, a limit imposed
by maximal force, and a viscous diffusion condition (see [121, 123]). These conditions
amount to

 
h h h2 ρa
∆t ≤ min 0.25 , 0.25 , 0.125 (5.28)
ca |f a | µ a p a + ca
p
where ca = Ka /ρa is the bulk sound speed and the last condition emerges from
analyzing the viscous condition discussed in [123] and [121]. In practice, some trial and
error may be necessary to ensure the time-step satisfies these conditions, particularly
the second condition since maximum particle accelerations may not be known a priori.

Algorithm 1 Velocity Verlet time step procedures for SPH simulation.

1. Begin with fixed particle mass ma , particle positions xna , velocities v na , and
accelerations ana from last time step, and constants β, ρ0 , γ, µl , µh , D∗ , c, µb , h.

2. Compute updated particle positions xn+1 a and partially updated particle veloc-
n+1/2 n
ities v a using aa and Eqs. (5.25) and Eqs. (5.27).

3. Perform intermediate steps to obtain v n+1

a .

(a) Compute density ρn+1

a using Eq. (5.18) and xn+1
a .
(b) Compute pressure pn+1
a using Eq. (5.5) and xn+1
a .
n+1/2
(c) Compute strain rate tensor D n+1
a using xn+1
a , va and Eqs. (5.2) and
(5.21).
n+1/2
(d) Compute friction coefficient µn+1
a using D n+1
a , xn+1
a , va and Eqs. (5.7)-
(5.8).
n+1/2
(e) Compute stress tensor σ n+1
a using xn+1
a , va and Eq. (5.3).
n+1/2
(f) Compute particle acceleration an+1
a using xn+1
a , va and Eq. (5.20).
(g) Compute updated velocity v n+1
a using Eq. (5.26).

4. Set xn+1
a , v n+1
a , an+1
a to xna , v na , ana and n + 1 to n and return to Step 2 for next
time step.
 92
5.4 Examples
This section provides several demonstrations of the SPH algorithm applied to mod-
eling granular media. Qualitative and partial quantitative comparisons with existing
experimental data are provided. The major findings of each example are: (1) for a
given friction coefficient µ, the method produces an angle of repose consistent with
classical theory; (2) the rate-dependent feature of the friction law produces flow pro-
files in inclined plane settings that are consistent with Bagnold’s theory; (3) the
method produces scaling laws for column collapse runout and slump that agree well
with experimental results; and (4) the method produces time-dependent column col-
lapse profiles down inclined planes that agree remarkably well with those observed
experimentally. In addition to these major findings, we also illustrate with these
examples the role of β in our friction law (Eq. 5.7), the influence of our boundary
condition on slip velocities and interior pressures, the effect of scaling the bulk modu-
lus on simulation results, and the influence of numerical resolution (i.e. SPH particle
mass).

5.4.1 Example 1: Angle of repose

This example illustrates how the friction coefficient µ that appears in Eq. (5.3)
influences the simulated angle of repose of the material. The friction coefficient µ
appearing in Eq. (5.3) is readily related to the internal friction angle φ of the material
by comparing Eq. (5.4) to the classical Mohr-Coulomb failure criterion to obtain

φ = tan−1 (µ) (5.29)

[128] noted that the angle of repose of a dry granular material has a roughly
constant value close to the angle of internal friction of the material in the loosest
packing state. This definition suggests that the angle of repose should be a lower
bound for the internal friction angle of the material. To test our model’s ability to
approximate this effect, we simulate a slump test through an orifice using various
 93
internal friction coefficients and compare the resulting angles of repose θ with the
internal friction angle computed using Eq. (5.29). As the first example of this method,
we also use the results to illustrate the pressure-drop in the bottom-most layer of the
slumped granular piles (mentioned in section 5.3.5).
The initial simulation conditions are shown in Fig. 5.2a. The simulated bulk of
granular material measures 60 cm in the x dimension, 50 cm in the y dimension, and
6 cm in the z dimension (into the page). Periodic boundary conditions are used in the
z direction to produce the effect of an infinitely wide slope. The granular material is
permitted to leave the domain through a 5 cm wide orifice in the bottom boundary.
The final simulation conditions are shown in Fig. 5.2b. Once the granular media
has flowed through the orifice and come to rest (kinetic energy density of less than
1 × 10−5 J/m3 ), the resulting angle of repose is measured by a linear least-squares fit
to the surface of the material in the middle 75% of the domain between the left wall
and the orifice. The simulation parameters used for this example are shown in Table
5.1. The SPH particle spacing is chosen to be ∆x = 0.01m and the SPH particle
mass is chosen to produce the loosest packing density of the granular material ρ0 at
this particle spacing. Five simulations are run, each with a constant internal friction
coefficient between 0.364 (φ = 20◦ ) and 0.839 (φ = 40◦ ). The parameter κ is set
to 105 to reduce the wave speed and permit a large time step. It will be shown in
example 3 that this choice, rather than a more realistic choice (e.g., 107 or higher),
does not significantly influence the material response in simulations where the wave
speed does not control the quantities of interest or the quasi-static deformation field
is not of interest. In simulations in which wave speed plays a significant role or the
quasi-static deformation must be captured accurately, κ should be chosen to enforce
a more realistic bulk modulus. The dilation β is limited to 0 ≤ β ≤ 0.5 as discussed
in section 5.2.3.
Simulation results are shown in Fig. 5.2c. Results are shown both for simulations
with and without β included in the calculation of friction (see Eq. (5.7)). The simu-
lations that do not use β in the calculation of friction (β = 0) yield an angle of repose
slightly lower than the internal friction angle of the material. This is consistent with
 94

h
✓

y
x
Internal Boundary 5 cm
(a) (b)

40 p = 0Pa p = 4kPa
Repose Angle (Degrees)

35

30

25
✓=
20 0   0.5
=0
15 (d)
20 25 30 35 40
Friction Angle (Degrees)
(c)

Figure 5.2: (a) Initial conditions of the slump test through an orifice. In color,
blue represents internal SPH particles and red represents boundary particles. Length
l = 60 cm, h = 50 cm, and the size of the simulated material in the z dimension is 6
cm. (b) Final collapsed side profile of the slump test through an orifice for φ = 30◦ .
The repose angle θ is measured in the middle 75% of the domain between the left
wall and the orifice. (c) Results of the simulations with and without including β in
the calculation of friction. The dashed line and circles represents the curve θ = φ.
Results without β are blue squares which provide a lower bound to the θ = φ curve.
Results with β are red inverted triangles which nearly match the θ = φ curve. (d)
Pressure in the collapsed granular pile for φ = 30◦ showing a slight pressure drop in
the layer of SPH particles directly above the bottom surface of the container.
 95

Table 5.1: Model parameters used in the angle of repose simulations. Note that the
friction law is taken as rate-independent, so µh is set to µl and the value of D∗ is
irrelevant. Basal friction coefficient µb is also set to µl . The initial SPH particle
spacing of ∆x = 0.01 m corresponds to the loosest packing state of the material (i.e.
in which ρ = ρ0 .)

ρ0 (kg/m3 ) µl µh D∗ c (Pa) µb κ (Pa) γ ∆x (m)

1550 0.364-0.839 µl - 0 µl 105 1.5 0.01

the classical theory of the angle of repose discussed by [128]. Simulations including
β in the calculation of friction (0 ≤ β ≤ 0.5) yield an angle of repose that nearly
matches the internal friction angle of the material, again consistent with the classical
theory. Both simulation results deviate from the classical theory at large friction
angles. This deviation may suggest a limitation of using a viscoplastic model for high
friction angle materials but has not yet been fully investigated.
These results suggest that the constitutive law proposed in this work, both with
and without the inclusion of dilation, is able to reproduce a key feature of granular
systems: a steady-state angle of repose θ that is a lower bound for the internal
friction angle φ of the material. We can see from Fig. 5.2c that including β in
the calculation of friction in Eq. (5.7) has the effect of slightly increasing the yield
strength of the granular material. This feature is consistent with the classical theory
that β enhances the ultimate strength of a granular material leading up to steady-state
(e.g., see [85,107]). However, including β in our model also appears to enhance steady-
state strength slightly. This feature of the constitutive law and chosen parameters is
discussed in more detail in section 5.4.2.
Fig. 5.2d illustrates the pressure distribution in the collapsed granular pile (for
φ = 30◦ ). This distribution is qualitatively similar to our simulations with other
material parameters. Interior SPH particles in the layer directly above the bottom
surface of the container exhibit a pressure slightly below that of the second layer. This
occurs because of the reduced value of ρ0 taken for boundary particles, as mentioned
in section 5.3.5. Using an identical value of ρ0 for boundary and internal particles
 96
eliminates this effect and restores the linear pressure increase as a function of depth
up to and including the bottom layer of interior SPH particles. However, using
identical values of ρ0 for boundary and internal particles can lead to internal particles
penetrating the boundary in some simulations. We have verified that aside from
this artifact in the bottom-most layer, the pressure distribution in a pile under self-
weight agrees with that expected from the overburden load. Therefore, we continue to
employ this reduced value of ρ0 for boundary particles and discuss possible extensions
in section 5.5.

5.4.2 Example 2: Infinite inclined plane flow

Inclined plane flow experiments are a classical means of studying the rate-dependent
frictional strength of granular media (e.g., see [4]). Experiments in this flow configu-
ration have been used to derive rate-dependent frictional parameters, study Bagnold-
type velocity profiles, investigate boundary-slip in granular flows, and examine non-
local and finite-size effects of granular media (see [4, 88, 129, 130]). We use inclined
plane flow simulations to qualitatively demonstrate the rate-dependent behavior of
our friction law and to illustrate the slip generated by our boundary conditions. We
also compare our results with a Bagnold-type velocity profile for completeness. Fi-
nally, we discuss the effect of including β and limiting it to 0 ≤ β ≤ 0.5 on the
frictional strength.
[105] presented a scaling argument for relating shear stress to strain rate in inertial
granular flow that has since been used to derive a velocity profile for inclined plane
flows. [130] and [88] derived Bagnold velocity profiles for inclined plane flows and
found a close comparison when fitting to numerical and experimental results. The
general form of the Bagnold profile is (see [88])

p (h3/2 − (h − y)3/2 )
vx (z) = A(θ) gd (5.30)
d3/2
where A(θ) is a θ dependent constant, θ is the slope inclination angle, g is the mag-
nitude of gravity, d is the grain diameter, h is the height of the flow, and y is the
 97

Table 5.2: Model parameters used in the infinite inclined-plane flow simulations. Note
that the basal friction coefficient µb is taken as equal to the value of µ and thus is
also rate-dependent.

ρ0 (kg/m3 ) µl µh D∗ (s−1 ) c (Pa) µb κ (Pa) γ ∆x (cm)

1550 0.268 0.577 20 0 µ 105 1.5 0.01

height above the inclined boundary. Close to the inclined plane surface, results do
not fit well due to a clearly observed nonzero boundary slip velocity (see [130] and
discussion in [88]). In addition, the Bagnold profile does not fit the results at the
free-surface, where [88] notes that there is a nonzero shear rate exhibited numerically
and experimentally but not modeled by Bagnold’s theory.
Figure 5.3a illustrates the setup of our inclined-plane flow simulations. The sim-
ulated bulk of granular material measures 10 cm in the x dimension, 15 cm in the y
dimension and 5 cm in the z dimension (into the page). Periodic boundary conditions
are used in the x and z directions to simulate the effect of an infinitely wide and long
slope. Table 5.2 lists the simulation parameters. The material is modeled as rate-
dependent with a friction coefficient varying from µl = 0.268 (φ = 15◦ ) to µh = 0.577
(φ = 30◦ ) with parameter D∗ = 20 s−1 . The bottom boundary is modeled using five
layers of fixed SPH particles with the same properties as those used to represent the
granular material, including frictional properties µl , µh , D∗ . Often a constant basal
friction coefficient is used in simulations of granular flows (e.g., see [109]). We find
that such a choice results in an unbounded boundary slip velocity when the inclination
angle exceeds the basal friction angle. Therefore, a constant basal friction coefficient
does not capture the classical nature of rate-dependent inclined plane flows in which
the material can find a steady state at a variety of inclination angles. The granular
material is permitted to settle under self weight before tilt is instantaneously applied
by changing the direction of gravity. The velocity profile is monitored by spatially
averaging SPH particle velocity in bins of various heights (in 0.01m increments from
the bottom boundary).
Simulation results are shown in Figs. 5.3c-d and 5.4a-d. Figs. 5.3c-d illustrate
 98

y y
h

x x

Internal Boundary
vx = 0 vx = max
(a) (b)
8 6
✓ = 27 =0 ✓ = 27 0  0.5
7
✓ = 24 5 ✓ = 24
6 ✓ = 21 ✓ = 21
✓ = 18 4 ✓ = 18
5
vx (m/s)

vx (m/s)

4 ✓ = 15 3 ✓ = 15
vxBag vxBag
3
2
2
1
1
0 0
0 5 10 15 0 5 10 15
y (cm) y (cm)
(c) (d)

Figure 5.3: (a) Initial conditions of a typical inclined plane flow experiment showing
internal SPH particles above y = 0 and boundary particles at or below y = 0. (b)
A typical velocity profile of an inclined plane flow. (c) Velocitiy profiles in the x
direction as a function of height y in the granular material for simulations that do
not include β in the calculation of friction. Symbols represent simulation results and
dashed lines represent the best fit of the Bagnold profile (Eq. (5.30)) to the simulation
data. (d) Same as (c) but for simulations that include β in the calculation of friction.
 99
the steady-state velocity profiles (symbols) for inclination angles angles ranging from
θ = 15◦ to θ = 27◦ in increments of 3◦ for simulations both with and without β
included in the calculation of friction. We clearly see a nonzero steady-state velocity
profile for angles θ = 18◦ through θ = 27◦ , as expected from the friction law. For
simulations with θ = 12◦ (not shown) and θ = 15◦ , the steady-state velocity profile
is zero because the imposed shear stress on the material equals but does not exceed
the yield strength (incipient flow). For simulations using θ = 30◦ (also not shown),
the velocity profile appears to approach a nonzero steady-state for simulations both
with and without the inclusion of β in the calculation of friction. However, relaxation
is very slow and we do not simulate the response to steady-state because of the
immense computational cost. For simulations using θ = 33◦ , the velocity profile does
not approach a steady-state profile for simulations either with or without β included
in the calculation of friction because the imposed shear stress on the material exceeds
maximum material strength, leading to a constant downslope acceleration.
Also shown on Figs. 5.3c-d as dashed lines are the Bagnold velocity profiles,
found using a least-squares fit of Eq. (5.30) to the data. These profiles demonstrate
an excellent match with the simulation data (symbols) within the bulk of the flow.
The mismatch at the boundary and free-surface is consistent with the discussion
in [88], highlighting the nonzero slip at the boundary and the nonzero strain rate at
the free surface.
Fig. 5.4a-b illustrates the average bulk acceleration ā˙ of the granular material as
a function of time for all simulations. This average bulk acceleration is computed by
averaging the velocities of all interior SPH particles in each simulation and taking
the time derivative using backward finite difference with a time step of 0.01s. We
note that due to finite numerical precision, all simulations eventually contain minor
particle disorder and therefore fluctuations in the ā˙ versus time curve. We thus fit
each the data with a curve of the form

b−a
ā˙ = a + (5.31)
c/t + 1
 100
where a, b, and c are best-fit parameters, and t is time in seconds. These curves are fit
to the data from the point of maximum acceleration, which occurs at roughly 0.02s
after gravity is rotated in each simulation. These curves are the ones shown in Figs.
5.4a-b, with the raw data and fits shown superimposed in the inset plots. It is clear
to see that for inclination angles from θ = 12◦ to θ = 27◦ , a nonzero steady-state
velocity profile is reached for simulations with and without β. This is confirmed by
the best-fit parameter b in Eq. (5.31) being zero or very close to zero for each of these
data sets. For inclination angle θ = 30◦ , the bulk acceleration approaches zero but
does not reach zero in the time frame of our simulation, as mentioned above. For
inclination angle θ = 33◦ , the bulk acceleration asymptotes to a nonzero value for
simulations both with and without β being included in the calculation of friction.
Figs. 5.4c-d show the slip velocities for all simulations, computed by averaging
the velocity of the bottom-most layer of SPH particles parallel to the stationary
boundary. The raw data is shown for each inclination angle as well as a dashed-line
found by least-squares fit of the data to Eq. (5.30). Consistent with other numer-
ical simulations and experiments, we observe a nonzero slip velocity that increases
monotonically with inclination angle (see insets for slip velocity versus θ). The slip
velocities asymptote toward a constant value for inclination angles θ = 12◦ through
θ = 30◦ for simulations both including and excluding β. Slip velocities undergo un-
bounded increase for inclination angles above θ = 30◦ in simulations with and without
β.
The slip velocities are a general feature both of the boundary conditions imposed in
our simulations and of the constitutive law. Future work could involve implementing
no-slip or reduced-slip boundary conditions by imposing an artificial stress tensor
on boundary SPH particles, with a value falling between the stress tensor computed
directly from Eq. (5.3) and that of the neighboring interior particles. This reduced-
slip boundary condition could be calibrated to experimental data on boundary slip
in real granular materials. We have verified that applying an artificial stress tensor
to boundary SPH particles does indeed produce a no-slip condition, as discussed
by [113, 121]. More details on our boundary condition are found in section 5.3.5.
 101

20 12
20 12
=0 0  0.5 10
15
10 8
15

7_ (m/s2 )

7_ (m/s2 )
10 6
8

a
4
5
7_ (m/s2 )

7_ (m/s2 )
2
10 0
6 0
0 10 20 30 40 0 10 20 30 40
Time (s) Time (s)
a

a
4
5
2

0 0
0 10 20 30 40 0 10 20 30 40
Time (s) Time (s)
(a) (b)
10 3 10
3
=0 2.5 0  0.5 2.5
Slip Velocity (m/s)

Slip Velocity (m/s)

8 2
8 2
Slip Velocity (m/s)

Slip Velocity (m/s)

1.5
1.5
1
1
6 0.5 6 0.5
0
10 15 20 25 30 35 0
3 (Degrees) 10 15 20 25 30 35
4 4 3 (Degrees)

2 2

0 0
0 10 20 30 40 0 10 20 30 40
Time (s) Time (s)
(c) (d)

Figure 5.4: (a) Bulk acceleration ā˙ of the granular material in the inclined plane
simulations, without the use of β in the calculation of friction, as a function of time.
The highest curve represents θ = 33◦ and each lower curve represents a reduction in
θ by 3◦ . All accelerations appear to asymptote to zero except for θ = 33◦ . The inset
shows raw data overlaid with the same fits. (b) Same as (a) except for simulations
that include β in the calculation of friction. (c) Slip velocity in simulations not
including β in the calculation of friction, found by averaging the vx velocity of the
bottom-most layer of interior SPH particles. All slip velocities appear to asymptote
to finite values except for θ = 33◦ . (d) Same as (c) but for simulations that include
β in the calculation of friction.
 102
5.4.3 Example 3: Granular column collapse

Granular column collapse is an important test of the flow and yield behavior of gran-
ular materials. This test is significant for a number of processes such as landslides,
granular avalanches, hill-slope stability, and granular dam-break scenarios. Many
researchers have studied this problem numerically and experimentally in order to un-
derstand the dynamics of the yielding and flow processes, the timescales associated
with various regimes of flow, the local rheology within the flow, and the final scaling
laws that describe the slumped profile of the medium (see [5, 108, 112, 131, 132, 133,
134, 135, 136]). Scaling laws describing the slumped runout length and final height
of a collapsed granular medium are typically studied using an axisymmetric column
or quasi-two-dimensional rectangular column as an initial state. These studies are
designed to predict the basic physics of a sudden collapse event such as a landslide.
This example examines the ability of our numerical model to produce scaling laws
for quasi-two-dimensional column collapses consistent with those found in existing
experimental literature. We also use this flow configuration to illustrate the influence
of resolution (SPH particle mass and spacing), bulk modulus, and friction parame-
ters on the final results of a column collapse simulation. The results of this latter
study establish the resolution and reduced bulk modulus used in other examples as
acceptable for accurately modeling processes such as granular avalanches.

5.4.3.1 Scaling laws

In order to compare our scaling laws with experimental data, we simulate column
collapse using two materials with properties reported by [5]. In particular, we model
“grit” (a fine-grained sand) and fine glass beads. These two materials have distinct
bulk friction angles and loose-packing densities. Furthermore, both experience basal
friction angles φb that are different than their bulk friction angles φ, as reported by [5].
Table 5.3 tabulates all properties of both materials, as well as the parameters used
in modeling each for this particular analysis. Both materials are modeled as rate-
independent and cohesionless, using the nominal parameters reported by [5]. The
 103
basal friction coefficient is µb = tan−1 (φb ).
Figs. 5.5a and 5.5b illustrate the initial setup of a typical SPH simulation of grit
and fine glass, respectively. The most significant difference between the two initial
configurations is the different initial length of the columns L0 . Simulations of grit
use L0 = 9 cm while simulations of fine glass use L0 = 2.5 cm. These values roughly
reflect those used for fitting scaling laws in [5] (i.e. see Figs. 10 and 11 of that
reference). A finer resolution is used for the fine glass simulations (∆x = 0.005 cm)
than for the grit simulations (∆x = 0.01 cm) in order to accommodate this difference.
Various initial heights H0 (19 for grit simulations, 13 for fine glass simulations) are
used in order to examine the final runout L and final height H as a function of initial
column aspect ratio a = H0 /L0 . In all simulations, the simulated bulk of grit or fine
glass measures 10 cm in the z dimension. Periodic boundaries are employed in the z
dimension to produce the effect of an infinitely wide slope. We note that this differs
from experiments presented by [5], who used sidewalls spaced 20 cm apart. We ignore
the effect of sidewalls in this analysis of scaling laws since we are primarily interested
in a scaling exponent. Sidewalls would have the effect of increasing the apparent
bulk friction φ slightly, as discussed in [137]. This would influence the runout and
final height of individual runout simulations but would not significantly influence the
derived scaling laws from a number of simulations using the same conditions. We let
the simulated grit or fine glass settle under self weight before instantaneously deleting
the right-most vertical wall and allowing the material to flow until it comes to rest.
After the material comes to rest, we measure the final height H by finding the y
coordinate of the highest interior SPH particle in the domain. We measure the final
runout length L by finding the largest x coordinate of an interior SPH particle in
the domain. Using materials with distinct bulk and basal friction coefficients and
loose-packing densities suggests that our method is likely able to accurately capture
the behavior of various materials flowing over various substrates.
Figs. 5.6a-d illustrate the our column collapse results for simulations with (Figs.
5.6c and 5.6d) and without (Figs. 5.6a and 5.6b) including the rate of dilation β in
the calculation of the friction coefficient. As reported in [5], data for each material
 104

L0 Grit L0 Fine glass

H0 H0

Internal Boundary Internal Boundary

(a) (b)

Figure 5.5: Typical initial conditions of the column collapse simulations for (a) grit
and (b) fine glass. The right-most vertical wall in each figure is the containing wall
that is deleted at t = 0. In the color figure, blue particles represent internal particles
and red particles represent boundary particles.

Table 5.3: Model parameters used in the granular column collapse simulations. Both
materials are modeled as rate-independent. Bulk and basal friction angles are taken
as the nominal values reported in [5].

Material ρ0 (kg/m3 ) µl µh D∗ s−1 c (Pa) µb κ (Pa) γ ∆x (cm)

Grit 1404 0.74 0.74 - 0 0.335 105 1.5 0.01
Glass 1450 0.456 0.456 - 0 0.263 105 1.5 0.005
 105

Table 5.4: Comparison of scaling laws derived from simulation data and those re-
ported in [5]. Note that [5] does not report values for λ1 or λ2 .

Source λ1 αa λ2 α2
SPH Grit (β = 0) 1.24 0.55 1.99 0.72
SPH Glass (β = 0) 2.07 0.61 3.49 0.72
SPH Grit (0 ≤ β ≤ 0.5) 0.97 0.61 0.79 1.01
SPH Glass (0 ≤ β ≤ 0.5) 1.38 0.68 1.36 0.96
[5] Grit - 0.58 - 0.88
[5] Glass - 0.6 - 0.92
[5] All - ∼0.6 - 0.9±0.1

with a > 1.7 collapse well onto a linear line in logarithmic space indicating scaling
laws of the form

H0 (L − L0 )
∼ λ1 aα1 and ∼ λ2 aα2 (5.32)
H L0
where λ1 , λ2 , α1 and α2 are scalars. The symbols in Fig. 5.6 represent the results of
our SPH simulations while the dashed line indicates a least-squares linear fit to the
data in logarithmic space. The reported slopes represent the scaling exponents α1
and α2 . We find similar scaling exponents α1 for final heights for grit and glass using
β = 0 (Fig. 5.6a) and 0 ≤ β ≤ 0.5 (Fig. 5.6c). As shown in table 5.4, these scaling
exponents agree well with the experimental results of [5] in wide slots (most similar to
infinitely wide slope modeled here; see Fig. 10 in that reference), indicating that our
method, both with and without β, accurately captures scaling of the granular column
collapse scenario. We find different values for final runout scaling when using β = 0
(Fig. 5.6b) and 0 ≤ β ≤ 0.5 (Fig. 5.6d). The final runout scaling for 0 ≤ β ≤ 0.5
agrees more closely with results of [5] for their gate experiments (most similar to the
infinitely wide slope modeled here; see Fig. 11 in that reference), as shown in table
5.4. Individual simulations of laboratory-scale experiments using β = 0 tend to over-
predict runout length, as discussed more in section 5.4.4. These results collectively
indicate that our present simulations that include β in the calculation of friction can
accurately reproduce experimental data. Section 5.4.4 further supports this claim.
 106

20 30
Grit (SPH) Grit (SPH) 1
12 Glass (SPH) 1 15 Glass (SPH) 0.72
8 Best-fit 0.61 Best-fit

(L ! L0 )=L0
8
0.72
H0 =H

4 4 1
0.55
1
2
2
1010
0 =0 =0
101
0.5
0.25 1 0 1.7 3 5 8
10 15 30 0.25 1010 2 3 5 8 15 30
a = H0 =L0 a = H0 =L0
(a) (b)
20 30
Grit (SPH) Grit (SPH)
12 Glass (SPH) 15 Glass (SPH) 1
8 Best-fit Best-fit 0.96
1
(L ! L0 )=L0

8
0.68
H0 =H

4 1.01
4
0.61 1
2
2 1
10 01
0 0  0.5 0  0.5
101 0.5
0.25 1 0 1.7 3 5 8
10 15 30 0.25 10
10 2 3 5 8 15 30
a = H0 =L0 a = H0 =L0
(c) (d)

Figure 5.6: (a) Scaling laws for collapsed height in simulations not using β in the
calculation of friction. Symbols are SPH simulations and dashed lines are linear best-
fit lines in logarithmic space for all data with x coordinate above 1.7. (b) Scaling laws
for collapsed runout distance in simulations not using β in the calculation of friction.
(c) Same as (a) but for simulations using β in the calculation of friction. (d) Same
as (b) but for simulations using β in the calculation of friction.
 107
5.4.3.2 Reduced bulk modulus, friction angle, and resolution studies

We briefly explore the influence of resolution (SPH particle spacing and mass), re-
duced bulk modulus, and employed friction angle on the numerical results. For this
study, we simulate the experimental setup of [5] with grit in a wide slot. The sim-
ulation is identical to those discussed in section 5.4.3.1 and listed in the first row of
table 5.3 except the initial length L0 is set to 20 cm, the initial height H0 is set to 13
cm, and the z dimension is periodic and measures 18 cm. The simulation parameters
for three parametric studies are outlined in table 5.5. Since the goal of this particular
parametric study is not to investigate the influence of β on the results, we present
some comparisons entirely with β = 0 and others entirely with 0 ≤ β ≤ 0.5.
The first study explores the effect of varying κ over two orders of magnitude
in the EOS presented in Eq. (5.5). The results of this study are shown in Fig.
5.7a which portrays simulation results with symbols and the experimentally observed
collapse profile with a dashed black line. This figure illustrates that the bulk modulus
has only a minor influence on the final results for this type of simulation. This is
further supported by the results of example 4, which illustrate an excellent agreement
between simulated time-dependent dynamic data using a reduced bulk modulus and
experimental results. The second parametric study explores the effect of varying the
bulk and basal friction coefficients within the error bounds provided in [5] for grit:
φ = 36.5◦ ±4.5◦ and φb = 18.5◦ ±1.5◦ . As shown in Fig. 5.7b, there is a clear difference
in the response of the material, although the difference is minor. The difference in
response is also as expected: a lower friction angle increases runout and reduces the
slumped height of the material. This topic will be discussed in more detail in section
5.4.4. Finally, the third parametric study explores the effect of varying the resolution
of the simulation. In SPH, the resolution simply follows SPH particle spacing and
mass. We vary initial particle spacing over almost an order of magnitude, using
4,680 particles in the coarsest simulation with ∆x = 0.01 m, 37,440 particles in the
intermediate resolution simulation with ∆x = 0.005 m, and 299,250 particles in the
finest simulation with ∆x = 0.0025 m. Fig 5.7c illustrates “mesh-convergence” in
 108

Table 5.5: Model parameters used in parametric studies of granular column collapse.
All models are rate-independent.

Test ρ0 (kg/m3 ) µl µh D∗ s−1 µb c (Pa) κ (Pa) γ ∆x (cm)

Modulus 1404 0.74 0.74 - 0.335 0 105 1.5 0.01
Modulus 1404 0.74 0.74 - 0.335 0 106 1.5 0.01
Modulus 1404 0.74 0.74 - 0.335 0 107 1.5 0.01
Friction 1404 0.74 0.74 - 0.335 0 105 1.5 0.01
Friction 1404 0.87 0.87 - 0.364 0 105 1.5 0.01
Friction 1404 0.625 0.625 - 0.305 0 105 1.5 0.01
Resolution 1404 0.74 0.74 - 0.335 0 105 1.5 0.01
Resolution 1404 0.74 0.74 - 0.335 0 105 1.5 0.005
Resolution 1404 0.74 0.74 - 0.335 0 105 1.5 0.0025

the sense that progressive refinements in resolution appear to approach the results of
the finest resolution (∆x = 0.0025 m). Fig. 5.7c also demonstrates that the coarsest
resolution of ∆x = 0.01 m used in most of the examples in this chapter yield nearly
identical results to the finest resolution. Further studies into mesh refinement along
with the results presented in Fig. 5.7c can elucidate the appropriate mesh to select
for field scale simulations of events such as dry granular avalanches. Future studies
on resolution are discussed in section 5.5.

5.4.4 Example 4: Column collapse down inclined planes

The ability of a numerical method to reproduce time-dependent two or three-dimensional

granular collapse profiles down inclined planes is important for several reasons. First,
granular column collapses down inclined planes reach significantly higher velocities
and strain rates than those on flat surfaces ( [6, 109]), making comparison with such
experiments a more challenging test of the flow rheology. Second, reproducing profiles
of individual flows is a more challenging test of the flow rheology than reproducing
scaling laws since many numerical methods such as DEM can accurately reproduce
scaling laws despite strongly overestimating the runout of individual experiments
( [134, 135]). Finally, capturing the full two or three-dimensional profile of a collapse
is crucial since many complex features that govern the flow, such as convection cells,
only occur below specific length scales ( [111]). This last point further motivates the
 109
0.15 0.15 0.15
k = 105 Pa = 36.5 , b = 18.5 x = 0.01 m
k = 106 Pa = 41 , b = 20 x = 0.005 m
k = 107 Pa = 32 , b = 17 x = 0.0025 m
0.1 Experiment 0.1 Experiment 0.1 Experiment
h (m)

h (m)

h (m)
0.05 0.05 0.05

0  0.5 =0 =0
0 0 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
x (m) x (m) x (m)
(a) (b) (c)

Figure 5.7: Results of parametric studies of various model parameters on numerical

results. (a) Result of varying the bulk modulus over two orders of magnitude does
not indicate any significant change in results. This also justifies our use of κ = 105 in
most simulations in this chapter to eliminate the shorter time step that may otherwise
be needed to ensure stability. (b) Result of varying bulk and basal friction coefficients
within the experimental error bars of [5] illustrates the expected changes, but only
minor changes. (d) Results of mesh-refinement study illustrate the coarse mesh used
throughout the chapter produces results that are nearly the same as the finest mesh.

need for the three-dimensional full-field modeling techniques like the one presented
in this chapter. One-dimensional depth-averaged finite difference models commonly
employed in modeling geological flows are unable to resolve these complex features
( [111]).
In this example, we compare time-dependent two-dimensional profiles of simulated
column collapses down inclined planes with those reported at [6]. This comparison
illustrates the capability of our method to reproduce the complex three-dimensional
time-dependent dynamical processes occurring during the collapse event. This com-
parison also demonstrates the importance of including the rate of dilation β in the
calculation of the friction coefficient for such processes.
We study four sets of numerical simulations of quasi-two-dimensional column col-
lapses down inclined planes with parameters listed in table 5.6. Each row of table 5.6
represents the parameters used in one set of simulations that includes the rate of dila-
tion β in the calculation of friction and one set of simulations that does not include β.
The distinction between the rows is the use of sidewalls. Two simulations employ pe-
riodic boundary conditions in the z direction to produce the effect of an infinitely wide
slope, while the other two explicitly model the sidewalls present in experiments. Figs.
 110
Initial Configuration Final Configuration
L0
H0 y L
H y
x x
✓ ✓

Internal Boundary Internal Boundary

(a) (b)

Figure 5.8: (a) Initial and (b) final configurations for simulations of column collapse
down inclined planes. The right-most near-vertical wall in (a) represents the con-
taining wall that is lifted in the y direction with vx = 2 m/s at t = 0. In the color
figure, blue particles represent interior particles and red particles represent boundary
particles.

5.8a and 5.8b illustrate the typical initial and final configurations, respectively, of all
simulations. Sidewalls are not shown in these figures. When used, sidewalls sand-
wich the interior SPH particles in the foreground and background between boundaries
with the same properties as the basal boundary. In all cases, the initial length of the
simulated granular material in the x dimension is 20 cm, the initial height in the y
dimension is 14 cm, and the initial width in the z dimension (into the page) is 10 cm.
Simulations progress by allowing the simulated material to settle under self weight
before lifting the containing wall in the y direction with velocity 2 m/s starting at
t = 0s. We use rate-independent friction in all simulations, leaving full paramet-
ric fitting of the rate-dependent rheology to existing experimental data for a future
study. We note that this contradicts the work of [111] who used a rate-dependent
model when simulating these collapses and [4] who derived rate-dependent parameters
from laboratory experiments with similar flow rates. Nevertheless, we employ rate-
independent friction and use the four sets of numerical simulations to illustrate how
including β in the calculation of friction can increase the ultimate strength enough to
accurately match experimental results. We note that the simulated granular material
with parameters described in [6] and presented in table 5.6 is a collection of glass
beads of diameter d between 600 and 800 µm.
 111

Table 5.6: Model parameters used in column collapse simulations for deriving final
runout and height scaling laws.

Sidewalls ρ0 (kg/m3 ) µl µh D∗ (s−1 ) µb c (Pa) κ (Pa) γ ∆x (cm)

No 1550 0.435 0.435 - 0.435 0 105 1.5 0.01
Yes 1550 0.435 0.435 - 0.435 0 105 1.5 0.01

Figs. 5.9a and 5.9b illustrate the results of the simulations without sidewalls using
β = 0 and 0 ≤ β ≤ 0.5, respectively, for four inclination angles: θ = 0◦ , 10◦ , 16◦ , and
22◦ . Symbols in these figures represent the height of the highest SPH particle at
three times after the containing wall is lifted. The lines represent the experimentally
observed collapse profiles at these times reported in [6]. A comparison between Figs.
5.9a and 5.9b illustrates that simulations using β in the calculation of friction exhibit
more accurate slumping in the bulk of the simulated material near the left-most
boundary than simulations using β = 0 in Fig. 5.9a. Furthermore, simulations
using β demonstrate reduced runout compared to those without β. These results are
consistent with the simulations of [111], who found excessive slumping and runout
even when using a rate-strengthening model without β. We have found in other
simulations, not shown here, that increasing µ as suggested in [137] and [111] to
model the effect of sidewalls does not correct the excessive slumping. The role of β in
our simulations therefore plays the unique role of increasing the ultimate strength of
the material to enhance yield strength and prevent slumping. It also appears to play
a strengthening role near the flow front, where material is constantly turned over and
subjected to rapid changes in strain rate.
Figs. 5.9c and 5.9d illustrate the results of the simulations with sidewalls using
β = 0 and 0 ≤ β ≤ 0.5, respectively, for the same four inclination angles as those in
Figs. 5.9a and 5.9b. These simulations are performed in an effort to exactly repro-
duce the experimental conditions of [6]. In previous simulations where comparison to
individual collapse profiles was not performed, simulating an infinitely wide slope was
sufficient for qualitative and partial quantitative validation. The sidewall boundaries
have the same properties as the bottom boundary. Fig. 5.9c illustrates that the pres-
 112
ence of sidewalls, even without including β in the calculation of friction, improves the
slumped profile and runout. However, the simulated profile near the left boundary
is still significantly below that found in experiments. Furthermore, the numerical
simulations continue to strongly overestimate the runout. Fig. 5.9d illustrates that
including both sidewalls and β in the calculation of friction produces results that
agree remarkably well at nearly all times for both the slump and runout. The ex-
ception to this remarkable fit occurs only for the highest inclination angle θ = 22◦ ,
where the final runout is overestimated. This overestimation likely occurs because a
rate-strengthening friction law is truly needed at the higher strain rates present at
this inclination angle. The excellent agreement at nearly all times suggests that the
presented numerical method, when modeling the same initial and boundary condi-
tions as those found experimentally, and including β in the calculation of friction,
can quantitatively reproduce collapse profiles and dynamics of granular flows down
inclined planes.
For brevity, we leave a full investigation of the three-dimensional continuum fields
and flow structures of our simulated incline plane flows for future work. Such an
investigation will elucidate the physics behind the dynamical processes observed in
simulations and experiments. In Fig. 5.10, we merely present snapshots of the 3D
inclined plane flow with sidewalls at the same times as those shown in Fig. 5.9d.
These snapshots are for the simulation employing β in the calculation of friction.
SPH particles are rendered as spheres and colored by the magnitude of their strain
rate tensor D. We clearly see at all inclination angles that the flow progresses over
a nearly stationary wedge of material. The influence of the sidewalls is evidenced
by higher shear rates in the foreground and background where the sidewalls (not
rendered) are located. It may be instructive in future work to investigate the flow/no-
flow transition, the basal stresses, and the full three-dimensional structure of other
continuum fields in these simulations.
 113

16 16
=0 ✓=0 0   0.5 ✓=0
12 12

h (cm)
t = 0.18s

h (cm)
No sidewalls No sidewalls t = 0.18s
8 t = 0.50s 8 t = 0.50s
t = 1.06s t = 1.06s
4 4
0 0
16-20 0 20 40 60 80 100 120 16-20 0 20 40 60 80 100 120
=0 x (cm) ✓ = 10 0   0.5x (cm) ✓ = 10
12 12

h (cm)
h (cm)
No sidewalls t = 0.18s No sidewalls t = 0.18s
8 t = 0.54s 8 t = 0.54s
4 t = 1.32s 4 t = 1.32s
0 0
16-20 0 20 40 60 80 100 120 16-20 0 20 40 60 80 100 120
=0 x (cm) ✓ = 16 0   0.5x (cm) ✓ = 16
12 12

h (cm)
h (cm)

No sidewalls t = 0.36s No sidewalls t = 0.36s

8 t = 1.12s 8 t = 1.12s
4 t = 1.62s 4 t = 1.62s
0 0
16-20 0 20 40 60 80 100 120 16-20 0 20 40 60 80 100 120
=0 x (cm) ✓ = 22
12 0   0.5x (cm) ✓ = 22
12

h (cm)
h (cm)

No sidewalls t = 0.32s No sidewalls t = 0.32s

8 t = 1.12s 8 t = 1.12s
t = 2.30s 4 t = 2.30s
4
0 0
-20 0 20 40 60 80 100 120 -20 0 20 40 60 80 100 120
x (cm) x (cm)
(a) (b)
16 16
=0 ✓=0 0   0.5 ✓=0
12 t = 0.18s 12 t = 0.18s
h (cm)

h (cm)

Sidewalls Sidewalls
8 t = 0.50s 8 t = 0.50s
t = 1.06s t = 1.06s
4 4
0 0
16-20 0 20 40 60 80 100 120 -20
16 0 20 40 60 80 100 120
=0 x (cm) ✓ = 10 0   0.5 x (cm) ✓ = 10
12 12
h (cm)

h (cm)

Sidewalls t = 0.18s Sidewalls t = 0.18s

8 t = 0.54s 8 t = 0.54s
4 t = 1.32s 4 t = 1.32s
0 0
16-20 0 20 40 60 80 100 120 16-20 0 20 40 60 80 100 120
=0 x (cm) ✓ = 16 0   0.5 x (cm) ✓ = 16
12 12
h (cm)

h (cm)

Sidewalls t = 0.36s Sidewalls t = 0.36s

8 t = 1.12s 8 t = 1.12s
4 t = 1.62s 4 t = 1.62s
0 0
16-20 0 20 40 60 80 100 120 16-20 0 20 40 60 80 100 120
=0 x (cm) ✓ = 22 0   0.5 x (cm) ✓ = 22
12 12
h (cm)
h (cm)

Sidewalls t = 0.32s Sidewalls t = 0.32s

8 t = 1.12s 8 t = 1.12s
t = 2.30s 4 t = 2.30s
4
0 0
-20 0 20 40 60 80 100 120 -20 0 20 40 60 80 100 120
x (cm) x (cm)
(c) (d)

Figure 5.9: Time-dependent collapse profiles down inclined planes for four inclination
angles: θ = 0◦ , 10◦ , 16◦ , and 22◦ . Symbols represent SPH profiles and lines represent
experimental profiles from [6]. (a) and (b) illustrate results with periodic boundary
conditions in the z direction and no sidewalls. (c) and (d) illustrate results with
sidewalls. (a) and (c) show results for simulations not including β in the calculation
of friction. (b) and (d) show results for simulations including β in the calculation of
friction.
 114

✓=0 ✓ = 10
t = 0s
t = 0s

t = 0.18s
t = 0.18s

t = 0.50s
t = 0.54s

t = 1.06s
t = 1.32s
1
||D||( s ) 1
||D||( s )
0 25
(a) 0 45
(b)

✓ = 16 ✓ = 22

t = 0s
t = 0s

t = 0.36s
t = 0.32s

t = 1.12s
t = 1.12s

t = 1.62s
t = 2.30s
1
||D||( s )
1
||D||( s )
0 50
(c)
0 (d) 50

Figure 5.10: Snapshots of simulations shown in Fig. 5.9d for (a) θ = 0◦ , (b) θ = 10◦ ,
(c) θ = 16◦ and (d) θ = 22◦ . These simulations employ β in the calculation of
friction. In the color figure, the color of each particle corresponds to the magnitude
of the strain rate tensor D. We clearly see the influence of sidewalls (not rendered)
at the first nonzero time for each inclination angle. At these times, particles appear
to undergo a higher strain rate in the foreground and background than in the center
of the flow.
 115
5.5 Discussion and conclusion
We have presented an SPH method for simulating the flow of granular materials. A
rate-dependent and dilation-dependent friction law has been employed and justified
from an analysis of the continuum energy balance equation. Numerical examples
have demonstrated that the SPH approach with the proposed constitutive law can
accurately reproduce an angle of repose consistent with the employed friction angle,
rate-dependent flows down inclined planes, scaling laws for granular column collapses
on flat surfaces, and the dynamic structure of column collapses down inclined planes.
We have seen through numerical examples that future work is required before
the numerical method can reliably be applied to a broad range of granular flows.
In particular, four major aspects of the numerical method or constitutive law need
additional study. First, the slip condition created by the imposed boundary condition
needs additional calibration to experimental data in order to be reliable for many
granular flows. Although it appears that the method is able to reproduce the runout
of granular column collapses down inclined planes (e.g., see section 5.4.4), the inclined-
plane flows in section 5.4.2 exhibit significant boundary slip that may not be realistic
in all scenarios. Second, the method for calculating the rate of dilation β used in
the friction coefficient needs to be refined and possibly given its own evolution law
as discussed in section 5.3.4. This will involve calibrating the evolution of β as a
function of numerous rheological parameters to experimental data as was done for the
frictional parameters µl , µh , and Ds∗ in [4]. Thirdly, we have not provided examples
of simulating granular materials with cohesion (i.e. c 6= 0). We have done this
purposefully since there exists less experimental validation data for flows of cohesive
materials. Simulating cohesive flows poses an additional challenge in SPH of the
so-called tensile instability. Methods exist for remedying the tensile instability (e.g.,
see [119, 121, 138]) and preliminary simulations of our method with cohesion do not
exhibit a tensile instability, likely due to the real viscosity we employ. Nevertheless,
future work should address any needed additions to our SPH framework that will
permit simulation of highly cohesive flows in order to expand the range of problems
 116
the method can address. Fourthly and finally, the friction law used in this chapter
may need additional work before it can be applied to a range of very rapid granular
flows. We have employed a friction law originally intended for steady-state flows in
the present chapter. Eq. (5.13) and [15] have shown that additional terms related to
changes in internal energy may be needed in flows exhibiting rapidly changing strain
rates. Furthermore, using the inertial number rather than Ds∗ in Eq. (5.7) as in [4]
may be more accurate for many flows.
In addition to the future work described above, it would be instructive to use
the method presented in this chapter to study the full three-dimensional structure
of continuum fields in column collapses down inclined planes such as those presented
in [6] and [136]. Future validation must also establish the the ability of the method
to accurately reproduce stress fields at the base of granular avalanches, debris flows,
and under sand piles. This test would elucidate the important role of basal stresses
in erosion and deposition of sediment during geophysical flows [139]. Finally, we are
currently implementing a two-phase approach for coupling granular flows with gas
and fluid flows. This extension will permit full-field simulation of geophysical debris
flows and other natural and industrial processes.
 117

Chapter 6

Discussion and conclusion

6.1 Summary
In chapters 2 through 5, we have provided four studies of granular materials across
length scales. We have focused on inter-particle forces at the micro-scale, illustrating
how these quantities influence macro-scale behavior such as intruder dynamics. The
inter-particle force inference technique detailed in chapters 2 and 3 provides a method
of studying force chains and validating grain-scale numerical modeling techniques.
The appeal of these methods is their applicability to arbitrarily shaped grains, non-
linear materials and non-Hertzian contacts, and dynamic events. Ongoing and future
work with these techniques is detailed in section 6.2.1.
In chapter 4, we have studied friction in granular flows at the mesoscale. We de-
rived a new friction law relating the friction coefficient to the shear rate, coordination
number, porosity, and grain-scale dissipation rates. This friction law allows us to
examine, through grain-scale simulations with arbitrary contact laws (or potentially
through experiments in the future), the grain-scale dissipation pathways responsi-
ble for causing macroscopic friction. This relationship also illustrates friction as a
delicate competition between material dilation and grain-scale dissipation. When
these two processes do not balance one another, rate-dependence is observed. Rate-
strengthening friction is common in granular materials and has been studied using
this relationship and numerical simulations. We direct the reader to another paper
in which we extend this friction law to non-steady-state flows [15].
 118
We have focused on continuum modeling at the macro-scale in chapter 5. Here, we
developed a new technique for capturing the solid-like, liquid-like, and gas-like states
of granular media in a single computational domain. The technique combines SPH
and a viscoplastic constitutive law containing a strain-rate and dilation-dependent
yield criterion. This yield criterion is justified theoretically using an analysis of the
continuum energy balance equation analogous to that provided in [15]. Numerous
examples demonstrate that the modeling technique accurately captures theoretically
predicted and experimentally observed behavior of granular media in quasi-static
and dynamic settings. Ongoing and future work will extend this technique to allow
coupling of the granular material with gases and fluids, as described in section 6.2.2.

6.2 Ongoing and future work

6.2.1 Inter-particle force inference

We are currently applying the inter-particle force inference technique developed in

chapters 2 and 3 to opaque three-dimensional single-crystal quartz grains imaged
using combined X-ray tomography and x-ray diffraction. This work is in collaboration
with Stephen Hall at Lund University and follows work described in [49].
The data set currently being studied comes from an oedeometric compression test
on 77 nearly-spherical single-crystal quartz grains contained in a 1.5 mm diameter by
1.5 mm tall quartz oedeometer. The grains are subjected to one-dimensional load-
unload-reload up to 70 N on the load curve and up to failure on the reload curve.
We note that failure occurs by crushing of a single grain in a strong force chain.
Fig. 6.1a illustrates the macroscopic load as measured by the load cell during this
experiment. X-ray computed tomography (XRCT) and X-ray diffraction (XRD) are
used periodically during the compression cycle to extract a volumetric image of the
assembly and volume-averaged grain strains (e.g., see [49, 81]). Fig. 6.1b shows a
sample volumetric image of the assembly obtained with XRCT during the load cycle.
Various in-house algorithms are being applied to these volumetric images to extract
 119
the evolution of particle position and contact point location and direction throughout
the loading cycle. Preliminary results of force inference using Eq. (3.19) at a single
load step are illustrated in Fig. 6.1c. Force chains are clearly seen to develop in the
bulk of the material, splitting and recombining at various locations. While we will
defer presenting most quantitative results until a future publication, we will mention
that the results agree remarkably well with expected statistical trends. In particular,
inferred forces exhibit an exponential probability distribution above the mean force
at each load step. Forces also divide into a strong (40% of total contacts) and weak
(60% of total contacts) network in close agreement with the predictions of [27].
The goal of the work described above is to provide the first three-dimensional
example of our force inference technique applied to opaque grains with sand-like
stiffness. In addition, the ongoing work will determine whether forces in such systems
exhibit features and statistics consistent with past theory and simulations. Force
distributions and force chains will also be studied leading up to the crushing event
that occurs during the reload cycle. This crushing event is a characteristic means of
dissipating energy during failure for granular media. A deeper understand of force
arrangement leading up to crushing could improve plasticity models and predictions of
energy dissipation for granular materials. However, the grain crushing event cannot be
accurately modeled with available grain-scale numerical techniques. Thus, the force
inference technique can provide information currently only experimentally available,
and can also serve as a calibration tool when grain-scale models for particle crushing
are developed.

6.2.2 SPH modeling of coupled granular media and fluid

flows

We are currently extending the SPH technique described in chapter 5 to handle gas
and fluid flows coupled to granular flows. This work was originally motivated by
a need to model gas impinging on a porous surface of sand or soil. This problem
has important applications for planetary exploration, where the current technique of
 120

(a) (b)
80

60
Load (N)

1.5mm
40

20

0
0 0.05 0.1 0.15 0.2
Displacement (mm)
1.5mm

(c)

Figure 6.1: (a) Macroscopic load curve for oedeometric compression of single-crystal
quartz grains. (b) Sample volumetric image of the assembly obtained with XRCT
during the load cycle. (c) Preliminary results of force inference using Eq. (3.19) ap-
plied to a single load step. Lines connect contacting grain centroids and are darkened
and thickened linearly with force magnitude.
 121
landing on celestial bodies involves firing rockets downward to arrest motion during
the final stages of vehicular descent. Predicting erosion during this descent would
tighten the bounds on mission risk and provide a quantitative means of comparing
landing sites. While extensive work has addressed part or all of this problem (e.g.,
see [140, 141, 142, 143]), the problem involves numerous coupled phase interactions
that are not easily examined separately. For example, viscous erosion, fluidized soil
erosion, and soil diffusion are all processes involving the coupled interaction of a gas
and granular material. Thus, a full-field modeling technique that captures the coupled
behavior of these two phases is required for a truly predictive simulation of erosion.
Furthermore, such a technique would enhance our understanding of debris flows, gas
or fluid driven terrestrial erosion, and industrial processes involving porous media.
The proposed approach in chapter 5 involved combining SPH with a viscoplastic
constitutive law. SPH is well-suited for the problem described above because it pro-
vides precise interface definition and allows one to calculate ballistic soil trajectories.
Eulerian techniques do not provide this precise interface definition or the ballistic soil
trajectories needed to predict the potentially hazardous motion of eroded soil.
Currently, we have implemented a two-phase saturated porous media formulation
in SPH. We have provided porosity-based switches for drag between phases that per-
mit one to model the Darcy-like flow regime characteristic of low porosity states and
the dusty-gas flow regime characteristic of high porosity states (see [144, 145]). Com-
bining these capabilities in a single computational domain will allow us to perform
predictive simulations of gas-driven soil erosion in the future.
A qualitative example of predicting gas-driven soil erosion is provided in Fig. 6.2.
This figure depicts a jet of gas impinging on a sandy surface. The domain is 0.5 m
wide in the x dimension and 0.1 m thick in the y dimension (into the page). Periodic
boundary conditions are used in the x and y dimensions. The sand has density
ρ0 = 1600 kg/m3 (using the conventions of chapter 5), and the gas is modeled as
ideal with a viscosity of 10 µPa·s. The frictional strength of the granular material
is a constant µ = 0.8, permeability is 2×10−10 m2 , and gravity is set to 3.711 m/s2 .
The soil layer is approximately 0.25 m thick in the z dimension, sitting above an
 122
impermeable solid substrate (i.e. fixed SPH particles). The gas is driven at 100 m/s
downward and the switch between Darcy drag and dusty-gas drag is at a porosity of
0.38.
The images in the left column of Fig. 6.2 illustrate the sand response only in the
first 75 ms following the gas impingement event. The images in the right column
illustrate the gas response only. At early times, rapid diffusion of the gas into the
granular medium results in a high-pressure region beneath the impingement zone.
This high-pressure region is diminished as gas diffuses deeper into the sand. During
this diminishing process, a crater begins to form on the surface of the granular material
due to the high-pressure bulb of gas forcing the granular medium upward. At later
times, the sand crater is seen further developing, with sand mass accelerating outward
and upward as it is entrained in the background dusty-gas flow.
Future work at the continuum scale will also address modeling other porous ma-
terials using the SPH framework described above and in chapter 5. Adding the
appropriate balance and constitutive laws will allow us to address a wide range of
problems involving large deformation and flow of materials. Some examples include
post-fracture flow of brittle ceramics and comminution of brittle rocks.
 123

p (Pa) p (Pa)
Initial setup Walls
z 1000 2000 3000 1000 2000 3000

x 0 4e+03 0 4e+03
Initial Surface

Soil response Gas response

p (Pa) p (Pa)
t=4ms
1000 2000 3000 1000 2000 3000

0 4e+03 0 4e+03

Soil response Gas response

p (Pa) p (Pa)
t=6ms
1000 2000 3000 1000 2000 3000

0 4e+03 0 4e+03

Soil response Gas response

p (Pa) p (Pa)
t=75ms 1000 2000 3000 1000 2000 3000

0 4e+03 0 4e+03

Soil response Gas response

Figure 6.2: Qualitative example of gas-driven soil erosion modeled using the SPH
framework. The left column of figures illustrates the sand response only in the first
75 ms of the gas impingement event. The right column of images illustrates the gas
response only.
 124

Bibliography

[1] E. Azéma and F. Radjaı̈, “Internal structure of inertial granular flows,” Physical
Review Letters, vol. 112, no. 7, p. 078001, 2014.

[2] Y. Forterre and O. Pouliquen, “Flows of dense granular media,” Annual Review
of Fluid Mechanics, vol. 40, pp. 1–24, 2008.

[3] S. Savage and M. Sayed, “Stresses developed by dry cohesionless granular ma-
terials sheared in an annular shear cell,” Journal of Fluid Mechanics, vol. 142,
pp. 391–430, 1984.

[4] P. Jop, F. Yoel, and O. Pouliquen, “A constitutive law for dense granular flows,”
Nature, vol. 441, pp. 727–730, 2006.

[5] N. Balmforth and R. Kerswell, “Granular collapse in two dimensions,” Journal

of Fluid Mechanics, vol. 538, pp. 399–428, 2005.

[6] A. Mangeney, O. Roche, O. Hungr, N. Mangold, G. Faccanoni, and A. Lu-

cas, “Erosion and mobility in granular collapse over sloping beds,” Journal of
Geophysical Research: Earth Surface (2003–2012), vol. 115, no. F3, 2010.

[7] P. G. de Gennes, “Granular matter: a tentative view,” Reviews of Modern

Physics, vol. 71, pp. S374–S382, Mar 1999.

[8] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, “Granular solids, liquids, and

gases,” Reviews of Modern Physics, vol. 68, no. 4, p. 1259, 1996.
 125
[9] I. S. Aranson and L. S. Tsimring, “Patterns and collective behavior in granular
media: Theoretical concepts,” Reviews of Modern Physics, vol. 78, no. 2, p. 641,
2006.

[10] R. C. Hurley and J. E. Andrade, “Friction in inertial granular flows: competition

between dilation and grain-scale dissipation rates,” Granular Matter, vol. 17,
no. 3, pp. 287–295, 2015.

[11] F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, “Rheophysics

of dense granular materials: Discrete simulation of plane shear flows,” Physical
Review E, vol. 72, no. 2, p. 021309, 2005.

[12] L. Rothenburg and R. J. Bathurst, “Analytical study of induced anistropy in

idealized granular materials,” Géotechnique, vol. 39, no. 4, pp. 601–614, 1989.

[13] R. Hurley, E. Marteau, G. Ravichandran, and J. E. Andrade, “Extracting inter-

particle forces in opaque granular materials: Beyond photoelasticity,” Journal
of the Mechanics and Physics of Solids, vol. 63, pp. 154–166, 2014.

[14] R. Hurley, K. Lim, G. Ravichandran, and J. Andrade, “Dynamic inter-particle

force inference in granular materials: Method and application,” Experimental
Mechanics, 2015.

[15] R. C. Hurley and J. E. Andrade, “Strength of granular materials in transient

and steady state rapid shear,” Procedia Engineering, vol. 103, pp. 237–245,
2015.

[16] M. Oda, “Initial fabrics and their relations to mechanical properties of granular
material,” Soils and Foundations, vol. 12, no. 1, pp. 17–36, 1972.

[17] J. Christoffersen, M. M. Mehrabadi, and S. Nemat-Nasser, “A micromechani-

cal description of granular material behavior,” Journal of Applied Mechanics,
vol. 48, no. 2, pp. 339–344, 1981.
 126
[18] R. J. Bathurst and L. Rothenburg, “Observations on stress-force-fabric rela-
tionships in idealized granular materials,” Mechanics of Materials, vol. 9, no. 1,
pp. 65–80, 1990.

[19] M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, and P. Claudin, “Jamming, force

chains, and fragile matter,” Physical Review Letters, vol. 81, pp. 1841–1844,
Aug 1998.

[20] J. F. Peters, M. Muthuswamy, J. Wibowo, and A. Tordesillas, “Characterization

of force chains in granular material,” Physical Review E., vol. 72, p. 041307, Oct
2005.

[21] T. S. Majmudar and R. P. Behringer, “Contact force measurements and stress-

induced anisotropy in granular materials,” Nature, vol. 435, pp. 1079–1082, 06
2005.

[22] S. N. Coppersmith, C. h. Liu, S. Majumdar, O. Narayan, and T. A. Wit-

ten, “Model for force fluctuations in bead packs,” Physical Review E., vol. 53,
pp. 4673–4685, 1996.

[23] N. Guo and J. Zhao, “The signature of shear-induced anisotropy in granular

media,” Computers and Geotechnics, vol. 47, no. 0, pp. 1 – 15, 2013.

[24] C.-h. Liu, S. R. Nagel, D. A. Schecter, S. N. Coppersmith, S. Majumdar,

O. Narayan, and T. A. Witten, “Force fluctuations in bead packs,” Science,
vol. 269, no. 5223, pp. 513–515, 1995.

[25] S. Ostojic, E. Somfai, and B. Nienhuis, “Scale invariance and universality of

force networks in static granular matter,” Nature, vol. 439, pp. 828–830, 02
2006.

[26] F. Radjai, M. Jean, J. J. Moreau, and S. Roux, “Force distributions in dense

two-dimensional granular systems,” Physical Review Letters, vol. 77, pp. 274–
277, Jul 1996.
 127
[27] F. Radjai, D. E. Wolf, M. Jean, and J. J. Moreau, “Bimodal character of stress
transmission in granular packings,” Physical Review Letters, vol. 80, no. 1,
pp. 61–64, 1998.

[28] M. Satake, “Fabric tensor in granular materials,” in Deformation and Failure

of Granular Materials (P. Vermeer and H. J. Luger, eds.), pp. 63–68, 1982.

[29] P. Claudin, J.-P. Bouchaud, M. E. Cates, and J. P. Wittmer, “Models of stress

fluctuations in granular media,” Physical Review E., vol. 57, pp. 4441–4457,
Apr 1998.

[30] J.-P. Bouchaud, P. Claudin, D. Levine, and M. Otto, “Force chain splitting in
granular materials: A mechanism for large-scale pseudo-elastic behaviour,” The
European Physical Journal E, vol. 4, no. 4, pp. 451–457, 2001.

[31] D. M. Mueth, H. M. Jaeger, and S. R. Nagel, “Force distribution in a granular

medium,” Physical Review E, vol. 57, no. 3, p. 3164, 1998.

[32] M. P. Ciamarra, A. H. Lara, A. T. Lee, D. I. Goldman, I. Vishik, and H. L.

Swinney, “Dynamics of drag and force distributions for projectile impact in a
granular medium,” Physical Review Letters, vol. 92, no. 19, p. 194301, 2004.

[33] P. A. Cundall and O. D. L. Strack, “A discrete numerical model for granular

assemblies,” Géotechnique, vol. 29, no. 1, pp. 47–65, 1979.

[34] M. Jean, “The non-smooth contact dynamics method,” Computer Methods in

Applied Mechanics and Engineering, vol. 177, no. 3, pp. 235–257, 1999.

[35] J. J. Moreau, “Unilateral contact and dry friction in finite freedom dynamics,”
in Nonsmooth Mechanics and Applications, pp. 1–82, Springer, 1988.

[36] A. Taboada, K.-J. Chang, F. Radjaı̈, and F. Bouchette, “Rheology, force trans-
mission, and shear instabilities in frictional granular media from biaxial numeri-
cal tests using the contact dynamics method,” Journal of Geophysical Research:
Solid Earth (1978–2012), vol. 110, no. B9, 2005.
 128
[37] E. Azéma, F. Radjai, R. Peyroux, and G. Saussine, “Force transmission in a
packing of pentagonal particles,” Physical Review E, vol. 76, no. 1, p. 011301,
2007.

[38] X. Garcia, J.-P. Latham, J. Xiang, and J. Harrison, “A clustered overlap-

ping sphere algorithm to represent real particles in discrete element modelling,”
Geotechnique, vol. 59, no. 9, pp. 779–784, 2009.

[39] A. A. Peña, P. G. Lind, and H. J. Herrmann, “Modeling slow deformation of

polygonal particles using dem,” Particuology, vol. 6, no. 6, pp. 506–514, 2008.

[40] A. Drescher and G. de Josselin de Jong, “Photoelastic verification of a me-

chanical model for the flow of a granular material,” Journal of Mechanics and
Physics of Solids, vol. 20, pp. 337–340, 1972.

[41] D. Howell, R. P. Behringer, and C. Veje, “Stress fluctuations in a 2d granular

couette experiment: A continuous transition,” Physical Review Letters, vol. 82,
pp. 5241–5244, Jun 1999.

[42] A. H. Clark, L. Kondic, and R. P. Behringer, “Particle scale dynamics in gran-

ular impact,” Physical Review Letters, vol. 109, p. 238302, Dec 2012.

[43] A. H. Clark, A. J. Petersen, and R. P. Behringer, “Collisional model for granular

impact dynamics,” Physical Review E, vol. 89, no. 1, p. 012201, 2014.

[44] A. Shukla, M. H. Sadd, R. Singh, Q. Tai, and S. Vishwanathan, “Role of particle

shape and contact profile on the dynamic response of particulate materials,”
Optics and Lasers in Engineering, vol. 19, no. 1, pp. 99–119, 1993.

[45] J. Zhou, S. Long, Q. Wang, and A. D. Dinsmore, “Measurement of forces inside

a three-dimensional pile of frictionless droplets,” Science, vol. 312, no. 5780,
pp. 1631–1633, 2006.
 129
[46] M. Saadatfar, A. P. Sheppard, T. J. Senden, and A. J. Kabla, “Mapping forces
in a 3d elastic assembly of grains,” Journal of the Mechanics and Physics of
Solids, vol. 60, no. 1, pp. 55 – 66, 2012.

[47] N. Brodu, J. A. Dijksman, and R. P. Behringer, “Spanning the scales of granular

materials through microscopic force imaging,” Nature Communications, vol. 6,
2015.

[48] J. Andrade and C. Avila, “Granular element method (gem): linking inter-
particle forces with macroscopic loading,” Granular Matter, vol. 14, pp. 51–61,
2012.

[49] S. Hall, J. Wright, T. Pirling, E. Andò, D. Hughes, and G. Viggiani, “Can in-
tergranular force transmission be identified in sand?,” Granular Matter, vol. 13,
pp. 251–254, 2011.

[50] R. V. Martins, L. Margulies, S. Schmidt, H. F. Poulsen, and T. Leffers, “Si-

multaneous measurement of the strain tensor of 10 individual grains embedded
in an al tensile sample,” Materials Science and Engineering: A, vol. 387–389,
no. 0, pp. 84 – 88, 2004.

[51] H. F. Poulsen, X. Fu, E. Knudsen, E. M. Lauridsen, K. Margulies, and

S. Schmidt, “3dxrd - mapping grains and their dynamics in 3 dimensions,”
Materials Science Forum, vol. Recrystallization and Grain Growth, pp. 1363–
1372, 2004.

[52] K. A. Alshibli and A. H. Reed, Applications of X-ray Microtomography to Ge-

omaterials. Wiley-ISTE, 1st ed., 2010.

[53] J. Desrues, G. Viggiani, and P. Bésuelle, eds., Advances in X-ray Tomography

for Geomaterials. Wiley-ISTE, 2006.

[54] J. Y. Wang, L. Park and Y. Fu, “Representation of real particles for dem simu-
lation using x-ray tomography,” Construction and Building Materials, vol. 21,
no. 2, pp. 338–346, 2007.
 130
[55] M. A. Sutton, J. J. Orteu, and H. Schreier, Image Correlation for Shape, Mo-
tion and Deformation Measurements: Basic concepts, theory and applications.
Springer, 2009.

[56] I. Vlahinic, E. Andò, G. Viggiani, and J. E. Andrade, “Towards a more accurate

characterization of granular media: extracting quantitative descriptors from
grain-scale images,” Granular Matter, 2012. Under review.

[57] Correlated Solutions, Vic-2D, Reference Manual.

[58] Correlated Solutions, Vic-2D, Testing Guide.

[59] B. Pan, K. Qian, H. Xie, and A. Asundi, “Two-dimensional digital image corre-
lation for in-plane displacement and strain measurement: a review,” Measure-
ment Science and Technology, vol. 20, no. 6, p. 062001, 2009.

[60] D. H. Ballard, “Generalizing the Hough transform to detect arbitrary shapes,”

Pattern Recognition, vol. 13, pp. 111–122, Jan. 1981.

[61] T. Peng, “Detect circles with various radii in grayscale image via hough trans-
form.” http://www.mathworks.com/matlabcentral/fileexchange, 2010.

[62] M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: Active contour models,”

International Journal of Computer Vision, vol. 1, pp. 321–331, Jan. 1988.

[63] D.-J. Kroon, “Snake: Active contour.” http://www.mathworks.com/matlabcentral/fileexchang

2011.

[64] J. C. Santamarina, K. A. Klein, and M. A. Fam, Soils and Waves. J. Wiley &
Sons, 2001.

[65] M. Bornert, F. Brémand, P. Doumalin, J.-C. Dupré, M. Fazzini, M. Grédiac,

F. Hild, S. Mistou, J. Molimard, J.-J. Orteu, L. Robert, Y. Surrel, P. Vacher,
and B. Wattrisse, “Assessment of digital image correlation measurement errors:
Methodology and results,” Experimental Mechanics, vol. 49, pp. 353–370, 2009.
 131
[66] P. Lötstedt, “Perturbation bounds for the linear least squares problem subject
to linear inequality constraints,” BIT Numerical Mathematics, vol. 23, pp. 500–
519, 1983.

[67] C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems. Society for Indus-
trial and Applied Mathematics, 1995.

[68] V. C. R., Computational Methods for Inverse Problems. Society for Industrial
and Applied Mathematics, 2002.

[69] F. Bauer and M. A. Lukas, “Comparing parameter choice methods for regu-
larization of ill-posed problems,” Mathematics and Computers in Simulation,
vol. 81, no. 9, pp. 1795 – 1841, 2011.

[70] J. F. Sturm, “Using sedumi 1.02, a MATLAB toolbox for optimization over
symmetric cones,” Optimization Methods and Software, vol. 11–12, pp. 625–
653, 1999.

[71] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex program-
ming, version 1.21,” Apr. 2011.

[72] J. Lofberg, “Yalmip : a toolbox for modeling and optimization in matlab,” in

Computer Aided Control Systems Design, 2004 IEEE International Symposium
on, pp. 284 –289, 2004.

[73] J. Oddershede, S. Schmidt, H. F. Poulsen, H. O. Sørensen, J. Wright, and

W. Reimers, “Determining grain resolved stresses in polycrystalline materials
using three-dimensional X-ray diffraction,” Journal of Applied Crystallography,
vol. 43, pp. 539–549, Jun 2010.

[74] S. Boyd and L. Vandenberghe, Convex Optimization. New York, NY, USA:
Cambridge University Press, 2004.

[75] Abaqus, Version 6.11. D S Simulia, Dassault Systèmes.

 132
[76] M. Bornert, F. Brémand, P. Doumalin, J.-C. Dupré, M. Fazzini, M. Grédiac,
F. Hild, S. Mistou, J. Molimard, J.-J. Orteu, et al., “Assessment of digital
image correlation measurement errors: methodology and results,” Experimental
Mechanics, vol. 49, no. 3, pp. 353–370, 2009.

[77] MATLAB, version 8.1.0 (R2013a). Natick, Massachusetts: The MathWorks

Inc., 2013.

[78] E. R. Davies, Machine vision: theory, algorithms, practicalities. Elsevier, 2004.

[79] T. J. Atherton and D. J. Kerbyson, “Size invariant circle detection,” Image and
Vision Computing, vol. 17, no. 11, pp. 795–803, 1999.

[80] S. Coppersmith, C.-h. Liu, S. Majumdar, O. Narayan, and T. Witten, “Model

for force fluctuations in bead packs,” Physical Review E, vol. 53, no. 5, p. 4673,
1996.

[81] S. A. Hall and J. Wright, “Characterisation of 3d force transmission in real

granular media.” Presented at the International Conference on Experimental
Mechanics, July 2014.

[82] S. Dwivedi, R. Teeter, C. Felice, and Y. Gupta, “Two dimensional mesoscale

simulations of projectile instability during penetration in dry sand,” Journal of
Applied Physics, vol. 104, no. 8, pp. 083502–083502, 2008.

[83] L. Tsimring and D. Volfson, “Modeling of impact cratering in granular media,”

in Powders and Grains 2005 (R. Garca-Rojo, H. J. Herrmann, and S. McNa-
mara, eds.), vol. 2, pp. 1215–1223, A. A. Balkema, Rotterdam, 2005.

[84] H. M. Jaeger and S. R. Nagel, “Physics of the granular state,” Science, vol. 255,
no. 5051, pp. 1523–1531, 1992.

[85] D. M. Wood, Soil behaviour and critical state soil mechanics. Cambridge uni-
versity press, 1990.
 133
[86] C. S. Campbell, “Rapid granular flows,” Annual Review of Fluid Mechanics,
vol. 22, no. 1, pp. 57–90, 1990.

[87] I. Goldhirsch, “Rapid granular flows,” Annual Review of Fluid Mechanics,

vol. 35, no. 1, pp. 267–293, 2003.

[88] G. MiDi, “On dense granular flows,” The European Physical Journal E, vol. 14,
no. 4, pp. 341–365, 2004.

[89] M. Jutzi and E. Asphaug, “Forming the lunar farside highlands by accretion of
a companion moon,” Nature, vol. 476, no. 7358, pp. 69–72, 2011.

[90] K. Kamrin and G. Koval, “Nonlocal constitutive relation for steady granular
flow,” Physical Review Letters, vol. 108, no. 17, p. 178301, 2012.

[91] M. Tankeo, P. Richard, and C. Édouard, “Analytical solution of the µ(i)-

rheology for fully developed granular flows in simple configurations,” Granular
Matter, vol. 15, no. 6, pp. 881–891, 2013.

[92] L. Rothenburg and R. J. Bathurst, “Analytical study of induced anistropy in

idealized granular materials,” Géotechnique, vol. 4, no. 1, pp. 601–614, 1989.

[93] T. Hatano and O. Kuwano, “Origin of the velocity-strengthening nature of

granular friction,” Pure and Applied Geophysics, vol. 170, no. 1-2, pp. 3–11,
2013.

[94] J. T. Jenkins, “Dense inclined flows of inelastic spheres,” Granular Matter,

vol. 10, no. 1, pp. 47–52, 2007.

[95] Q. Sun, F. Jin, and G. G. Zhou, “Energy characteristics of simple shear granular
flows,” Granular Matter, vol. 15, no. 1, pp. 119–128, 2013.

[96] M. Babic, H. H. Shen, and H. T. Shen, “The stress tensor in granular shear
flows of uniform, deformable disks at high solids concentrations,” Journal of
Fluid Mechanics, vol. 219, no. 10, pp. 81–118, 1990.
 134
[97] S. Plimpton, “Lammps.”

[98] S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,”

Journal of Computational Physics, vol. 117, no. 1, pp. 1–19, 1995.

[99] H. Zhang and H. Makse, “Jamming transition in emulsions and granular mate-
rials,” Physical Review E, vol. 72, no. 1, p. 011301, 2005.

[100] A. D. Renzo and F. P. D. Maio, “Comparison of contact-force models for the

simulation of collisions in dem-based granular flow codes,” Chemical Engineer-
ing Science, vol. 59, no. 3, pp. 525–541, 2004.

[101] F. G. Bridges, A. Hatzes, and D. N. C. Lin, “Structure, stability and evolution

of saturn’s rings,” Nature, vol. 309, pp. 333–335, 1984.

[102] N. V. Brilliantov, F. Spahn, J.-M. Hertzsch, and T. Pöschel, “Model for colli-
sions in granular gases,” Physical Review E, vol. 53, no. 5, p. 5382, 1996.

[103] K. Senetakis, M. R. Coop, and M. C. Todisco, “The inter-particle coefficient

of friction at the contacts of leighton buzzard sand quartz minerals,” Soils and
Foundations, vol. 53, no. 5, pp. 746–755, 2013.

[104] K. Bagi, “Stress and strain in granular assemblies,” Mechanics of Materials,

vol. 22, no. 3, pp. 165–177, 1996.

[105] R. A. Bagnold, “Experiments on a gravity-free dispersion of large solid spheres

in a newtonian fluid under shear,” in Proceedings of the Royal Society of London
A: Mathematical, Physical and Engineering Sciences, vol. 225, pp. 49–63, The
Royal Society, 1954.

[106] O. Pouliquen and Y. Forterre, “Friction law for dense granular flows: applica-
tion to the motion of a mass down a rough inclined plane,” Journal of Fluid
Mechanics, vol. 453, pp. 133–151, 2002.

[107] J. E. Andrade, Q. Chen, P. H. Le, C. F. Avila, and T. M. Evans, “On the

rheology of dilative granular media: Bridging solid-and fluid-like behavior,”
 135
Journal of the Mechanics and Physics of Solids, vol. 60, no. 6, pp. 1122–1136,
2012.

[108] L. Lacaze and R. R. Kerswell, “Axisymmetric granular collapse: a transient 3d

flow test of viscoplasticity,” Physical Review Letters, vol. 102, no. 10, p. 108305,
2009.

[109] A. Lucas, A. Mangeney, and J. P. Ampuero, “Frictional velocity-weakening in

landslides on earth and on other planetary bodies,” Nature Communications,
vol. 5, 2014.

[110] J. E. Andrade and R. I. Borja, “Capturing strain localization in dense sands with
random density,” International Journal for Numerical Methods in Engineering,
vol. 67, no. 11, pp. 1531–1564, 2006.

[111] I. R. Ionescu, A. Mangeney, F. Bouchut, and O. Roche, “Viscoplastic modeling

of granular column collapse with pressure-dependent rheology,” Journal of Non-
Newtonian Fluid Mechanics, vol. 219, pp. 1–18, 2015.

[112] C. M. Mast, P. Arduino, P. Mackenzie-Helnwein, and G. R. Miller, “Simulating

granular column collapse using the material point method,” Acta Geotechnica,
vol. 10, no. 1, pp. 101–116, 2014.

[113] H. H. Bui, R. Fukagawa, K. Sako, and S. Ohno, “Lagrangian meshfree parti-

cles method (sph) for large deformation and failure flows of geomaterial using
elastic-plastic soil constitutive model,” International Journal for Numerical and
Analytical Methods in Geomechanics, vol. 32, no. 12, p. 1537, 2008.

[114] W. Chen and T. Qiu, “Numerical simulations for large deformation of granu-
lar materials using smoothed particle hydrodynamics method,” International
Journal of Geomechanics, vol. 12, no. 2, pp. 127–135, 2011.

[115] H. P. Zhang and H. A. Makse, “Jamming transition in emulsions and granular

materials,” Physical Review E, vol. 72, p. 011301, Jul 2005.
 136
[116] D. W. Taylor, “Fundamentals of soil mechanics.,” Soil Science, vol. 66, no. 2,
p. 161, 1948.

[117] J. J. Monaghan, “Smoothed particle hydrodynamics,” Reports on progress in

physics, vol. 68, no. 8, p. 1703, 2005.

[118] G.-R. Liu and M. B. Liu, Smoothed particle hydrodynamics: a meshfree particle
method. World Scientific, 2003.

[119] M. Liu and G. Liu, “Smoothed particle hydrodynamics (sph): an overview

and recent developments,” Archives of Computational Methods in Engineering,
vol. 17, no. 1, pp. 25–76, 2010.

[120] J. J. Monaghan and J. C. Lattanzio, “A refined particle method for astrophys-

ical problems,” Astronomy and Astrophysics, vol. 149, pp. 135–143, 1985.

[121] J. P. Morris, P. J. Fox, and Y. Zhu, “Modeling low reynolds number incom-
pressible flows using sph,” Journal of Computational Physics, vol. 136, no. 1,
pp. 214–226, 1997.

[122] D. J. Price and J. Monaghan, “Smoothed particle magnetohydrodynamics–ii.

variational principles and variable smoothing-length terms,” Monthly Notices
of the Royal Astronomical Society, vol. 348, no. 1, pp. 139–152, 2004.

[123] J. J. Monaghan, “Smoothed particle hydrodynamics,” Annual Review of As-

tronomy and Astrophysics, vol. 30, pp. 543–574, 1992.

[124] J. J. Monaghan, “Simulating free surface flows with sph,” Journal of Compu-
tational Physics, vol. 110, no. 2, pp. 399–406, 1994.

[125] J. E. Andrade and X. Tu, “Multiscale framework for behavior prediction in

granular media,” Mechanics of Materials, vol. 41, no. 6, pp. 652–669, 2009.

[126] P. Randles and L. Libersky, “Smoothed particle hydrodynamics: some recent

improvements and applications,” Computer Methods in Applied Mechanics and
Engineering, vol. 139, no. 1, pp. 375–408, 1996.
 137
[127] J. Monaghan, “On the problem of penetration in particle methods,” Journal of
Computational physics, vol. 82, no. 1, pp. 1–15, 1989.

[128] K. Terzaghi, Theoretical soil mechanics. Wiley, 1943.

[129] O. Pouliquen and Y. Forterre, “A non-local rheology for dense granular flows,”
Philosophical Transactions of the Royal Society of London A: Mathematical,
Physical and Engineering Sciences, vol. 367, no. 1909, pp. 5091–5107, 2009.

[130] L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, D. Levine, and S. J. Plimpton,

“Granular flow down an inclined plane: Bagnold scaling and rheology,” Physical
Review E, vol. 64, no. 5, p. 051302, 2001.

[131] G. Lube, H. E. Huppert, R. S. J. Sparks, and M. A. Hallworth, “Axisymmetric

collapses of granular columns,” Journal of Fluid Mechanics, vol. 508, pp. 175–
199, 2004.

[132] G. Lube, H. E. Huppert, R. S. J. Sparks, and A. Freundt, “Collapses of two-

dimensional granular columns,” Physical Review E, vol. 72, no. 4, p. 041301,
2005.

[133] R. Kerswell, “Dam break with coulomb friction: A model for granular slump-
ing?,” Physics of Fluids (1994-present), vol. 17, no. 5, p. 057101, 2005.

[134] L. Staron and E. Hinch, “Study of the collapse of granular columns using two-
dimensional discrete-grain simulation,” Journal of Fluid Mechanics, vol. 545,
pp. 1–27, 2005.

[135] P.-Y. Lagrée, L. Staron, and S. Popinet, “The granular column collapse as a
continuum: validity of a two-dimensional navier–stokes model with a µ (i)-
rheology,” Journal of Fluid Mechanics, vol. 686, pp. 378–408, 2011.

[136] M. Farin, A. Mangeney, and O. Roche, “Fundamental changes of granular

flow dynamics, deposition, and erosion processes at high slope angles: Insights
 138
from laboratory experiments,” Journal of Geophysical Research: Earth Surface,
vol. 119, no. 3, pp. 504–532, 2014.

[137] P. Jop, Y. Forterre, and O. Pouliquen, “Crucial role of sidewalls in granular sur-
face flows: consequences for the rheology,” Journal of Fluid Mechanics, vol. 541,
pp. 167–192, 2005.

[138] J. J. Monaghan, “Sph without a tensile instability,” Journal of Computational

Physics, vol. 159, no. 2, pp. 290–311, 2000.

[139] R. M. Iverson, “The physics of debris flows,” Reviews of Geophysics, vol. 35,
no. 3, pp. 245–296, 1997.

[140] H.-Y. Ko and R. Scott, “Transient rocket-engine gas flow in soil.,” AIAA Jour-
nal, vol. 6, no. 2, pp. 258–264, 1968.

[141] L. Roberts, “The interaction of a rocket exhaust with the lunar surface.,” 1966.

[142] P. T. Metzger, J. E. Lane, and C. D. Immer, “Modification of roberts’ theory

for rocket exhaust plumes eroding lunar soil,” Earth and Space, pp. 1–8, 2008.

[143] P. T. Metzger, J. Smith, and J. E. Lane, “Phenomenology of soil erosion due

to rocket exhaust on the moon and the mauna kea lunar test site,” Journal of
Geophysical Research: Planets (1991–2012), vol. 116, no. E6, 2011.

[144] O. C. Zienkiewicz, A. Chan, M. Pastor, B. Schrefler, and T. Shiomi, Computa-

tional geomechanics. Wiley Chichester, 1999.

[145] G. Laibe and D. J. Price, “Dusty gas with smoothed particle hydrodynamics–i.
algorithm and test suite,” Monthly Notices of the Royal Astronomical Society,
vol. 420, no. 3, pp. 2345–2364, 2012.