# Role of gravity or confining pressure and contact stiffness in granular
rheology

Abhinendra Singh, Vanessa Magnanimo, Kuniyasu Saitoh and Stefan Luding  Multi
Scale Mechanics (MSM), MESA + , CTW, University of Twente, Post Office Box
217, 7500 AE Enschede, The Netherlands<EMAIL_ADDRESS>
###### Abstract

The steady shear rheology of granular materials is investigated in slow quasi-
static states and inertial flows. The effect of the gravity field and contact
stiffness, which are conventionally trivialized is the focus of this study.
Series of Discrete Element Method simulations are performed on a weakly
frictional granular assembly in a split-bottom geometry considering various
gravity fields and contact stiffnesses. While traditionally the inertial
number, i.e., the ratio of stress to strain-rate timescales describes the flow
rheology, we find that a second dimensionless number, the ratio of softness
and stress timescales, must also be included to characterize the bulk flow
behavior. For slow, quasi-static flows, the density increases while the
macroscopic friction decreases with respective increase in particle softness
and gravity. This trend is added to the $\mu(I)$ rheology and can be traced
back to the anisotropy in the contact network, displaying a linear correlation
between macroscopic friction and deviatoric fabric in the steady state.
Interestingly, the linear relation holds when the external rotation rate is
increased for a given gravity field and contact stiffness.

###### pacs:

81.05.Rm, 45.70.Mg, 47.57.Gc

##  1 Introduction

Gravity is a critical factor in many natural (granular) phenomena like
avalanches, landslides, sand-piles and even in some industrial applications
[1, 2]. Avalanches and debris flows play an important role in the
transportation of mass existing at the surface of earth. Gravity-driven flows
have also been observed on other planetary bodies of our Solar System and are
of particular importance in understanding the geology of other planets and
asteroids as well as for the human exploration of the Moon and Mars in the
coming decades [3]. Currently, surface features found on Mars [4], Venus [5],
and the Moon [6] are hypothesized to be the results of avalanches of granular
material.

One of the important aspects of granular shear flows is the dependence of
stress on external driving. Various experimental and numerical studies have
shown that for slow–dense, quasi-static flows, the ratio of shear to
compressive stress (effective friction coefficient) is independent of the
imposed driving rate [7, 8, 9]. However, very little is known regarding the
same in the presence of very weak gravity fields or low confining stress.
Shear tests performed on parabolic flights have shown an increase in the
friction coefficient at low confinement [10, 11, 12, 13, 14, 15]. Brucks et
al[16] also obtained similar trend using centrifuge experiments at gravity
levels larger than Earth’s gravity. Despite these studies, the effect of
external compression (gravity) on granular flows is still poorly understood,
which leads to issues like exploration vehicles getting stuck in the Martian
soil.

Soft materials like hydrogel and elastomer, which can support large
deformation are of increasing importance in engineering and biological
applications such as tissue scaffolding, biosepration and micro-and-nano–
printing [17]. While inertial number has been relatively successful in
understanding the dynamics of rigid particles [18], elasticity becomes
relevant for soft particles [8, 19]. The deformability of the soft particles
has been shown to affect the force network close to the jamming transition
[20]. Recent study by Vaart et al[21] has shown different rheological behavior
for hard and soft particle suspensions. Despite the increasing importance, the
models for soft deformable particles have been largely ignored.

We claim here that these two factors, i.e., gravity and softness are two
aspects of the same phenomenon. We aim to test this claim by answering the
following questions: (1) How does gravity and softness of particles affect the
bulk flow behavior? (2) Is there a unique law that can describe the flow
behavior on Earth, Mars and the Moon for both soft and rigid particles?

In this paper, we address the above questions with a focus on dense,
frictional, quasi-static granular flow. Using Discrete Element Method (DEM),
we simulate cohesionless frictional granular material in a split-bottom ring
shear cell. An important aspect of this setup is that the shear rate is given
solely by external rotation rate and the geometry. At the same time, in this
geometry the local strain rate does not depend strongly on the external
compression [22], unlike the inclined plane and rotation drum where gravity
has a strong effect on the local strain rate [16, 23, 24]. To study the effect
of gravity and particle softness, we independently vary both parameters by two
orders of magnitude. A change in particle softness provides an adjustment on
the microscopic scale, while gravity is a macroscopic (field) modification. We
find that they have similar effect at the mesoscopic (local) scale. The bulk
behavior can be described well using a dimensionless parameter, defined as the
ratio between the time scales due to gravitational compression and contact
stiffness. Further, by increasing the external rotation rate, we study the
dependence of macroscopic friction and contact network anisotropy (deviatoric
fabric) on the inertial number. The dependence of macroscopic friction and
deviatoric fabric on pressure is added to $\mu(I)$ rheology. Additionally, we
find some non-local effects in our results due to the presence of gradients in
both stress and strain rate, which are quantified by following an approach
similar to Koval et al[25].

The outline of this paper is as follows: In section 2 we explain our numerical
setup and methodology. We present our results for quasi-static state in
section 3. In section 4, we provide results on the rheology and combine it
with the results from quasi-static state to present new rheological laws. We
then close the paper with a discussion and conclusion of our results in
section 5 along with a possible outlook for future research.

##  2 Discrete Element Method (DEM)

We present our numerical simulation scheme and setup in sections 2.1 and 2.2
respectively. Section 2.3 briefly presents our averaging methodology and
definitions of the tensorial quantities. We summarize various time scales
associated with the system in section 2.4.

###  2.1 Model

Our computational techniques are based on the soft-sphere DEM simulations as
developed by Cundall and Strack [26], Walton [27] and Luding [28, 29, 30]. The
normal force between the particles in contact is modeled by a Hookean spring
as $f_{n}=-k_{n}\delta_{n}-\eta_{n}v_{n}$, where $k_{n}$, $\delta_{n}$,
$\eta_{n}$, and $v_{n}$ are the normal stiffness, particle overlap, normal
viscous damping coefficient, and relative velocity in the normal direction,
respectively. Similarly, the tangential force is modeled as
$f_{t}=-k_{t}\delta_{t}-\eta_{t}v_{n}$, where $k_{t}=2k_{n}/7$, $\delta_{t}$,
$\eta_{t}=\eta_{n}/4$, and $v_{t}$ are the tangential stiffness, relative
displacement in tangential direction, tangential viscous damping coefficient,
and relative velocity in tangential direction, respectively. We also introduce
Coulomb’s friction between the particles, where the tangential force $f_{t}$
is switched to the sliding force $f_{s}=-\mu_{p}|f_{n}|$, $\mu_{p}$ being the
particle friction coefficient when $f_{t}$ exceeds the critical value, i.e.
$|f_{t}|>\mu_{p}|f_{n}|$ (with $\mu_{p}=0.01$) [30].

To study the effect of particle softness on macroscopic behaviors, we explore
a range of normal contact stiffness $k_{n}$, from $10\ \mathrm{Nm^{-1}}\leq
k_{n}\leq 10^{4}\ \mathrm{Nm^{-1}}$. While changing $k_{n}$, the time step for
numerical integration $\delta t$ is adjusted to be $\frac{1}{50}$ times the
contact duration to ensure accurate dynamic integration [30]. When $k_{n}$ is
changed, $k_{t}$ and $\eta_{n}$ are changed as well, to ensure that the
coefficient of restitution remains unchanged.

###  2.2 Split-bottom ring shear cell

Figure 1: (Color online) A sketch of our numerical setup consisting of a fixed
inner part (light blue shade) and a rotating outer part (white). The white
part of the base and the outer cylinder rotate with the same angular velocity,
$\Omega$, around the symmetry axis (dot-dashed line). The inner, split, and
outer radii are given by $R_{i}=8\langle a\rangle$, $R_{s}=40\langle
a\rangle$, and $R_{o}=52\langle a\rangle$, respectively, where each radius is
measured from the symmetry axis. The gravity, $g$, points downwards as shown
by an arrow.

Figure 1 is a sketch of our numerical setup [31, 32, 33, 34]. In this figure,
the inner, split, and outer radii are given by $R_{i}$, $R_{s}$, and $R_{o}$,
respectively, where the concentric cylinders rotate relative to each other at
a rate $\Omega$ around the symmetry axis (the dot-dashed line). The ring
shaped split at the bottom separates the moving and static parts of the
system, where a part of the bottom and the outer cylinder rotate at the same
rate. The system is filled with $N\approx 3.7\times 10^{4}$ polydispersed
spherical particles with density $\rho=2000\ \mathrm{kgm^{-3}}=2\
\mathrm{gcm^{-3}}$ up to height $H$. The average size of particles is $\langle
a\rangle=1.1$ mm, and the width of the homogeneous size-distribution (with
$a_{\rm min}/a_{\rm max}=1/2$) is $1-{\cal A}=1-\langle a\rangle^{2}/\langle
a^{2}\rangle=0.18922$. The cylindrical walls and the bottom are roughened due
to some (about $3\%$ of the total number) attached/glued particles [33, 35,
36]. The initial state of the system is prepared frictionless. Friction is
switched on at the onset of shear.

When there is relative motion at the split, a shear band initiates at the
bottom from the split position $R_{s}$ and propagates upwards and inwards,
remaining far away from the cylindrical walls in most cases [32, 34]. The
qualitative behavior is governed by the ratio $H/R_{s}$ and three regimes can
be observed as reported in [32, 34]. We keep $H/R_{s}<0.5$, such that the
shear band reaches the free surface and stays away from inner wall.

To understand the effect of gravity on the bulk behavior, gravity is varied in
the range $0.5\ \mathrm{ms^{-2}}\leq g\leq 50\ \mathrm{ms^{-2}}$. The details
regarding rotation rate of the system are reported in table 1. The total
simulation time is chosen such that the cylinder completes half a rotation in
that time.

###  2.3 Local averaging

One of the goals of current research in the granular community is to derive
macroscopic continuum theory based on the given micro-mechanical properties
[37, 38, 39]. We assume translational invariance in the azimuthal
$\theta-$direction. The averaging is thus performed over toroidal volumes over
many snapshots in time, leading to fields $Q(r,z)$ as function of the radial
and vertical positions. Here, the averaging is performed with spacings $\Delta
r$ and $\Delta z$ around two particles diameter in radial and vertical
directions (averaging procedure for two and three dimensions is discussed in
detail in Refs. [35, 38] and hence not discussed further here).

####  2.3.1 Macroscopic (tensorial) quantities

Here, we present general definitions of the averaged microscopic and
macroscopic quantities – including strain rate, stress and fabric (structure)
tensors.

The strain rate is calculated by averaging velocity gradients of particles
over cells and is given by

$\dot{\epsilon}_{\alpha\beta}(\mathbf{r})=\frac{1}{2}\displaystyle\sum_{p\in
V}\left(\bigtriangledown_{\beta}{\rm{v}}^{p}_{\alpha}+\bigtriangledown_{\alpha}{\rm{v}}^{p}_{\beta}\right)~{},$
(1)

where Greek letters represent components radial distance $r$, azimuthal angle
$\theta$, and height $z$, while $\bigtriangledown$ represents the gradient in
cylindrical coordinate [35].

The stress tensor is given by

${{\sigma}_{\alpha\beta}}(\mathbf{r})=\displaystyle\frac{1}{V}\left[\displaystyle\sum_{p\in
V}{{m^{p}}({\mathrm{v}^{p}_{\alpha}})({\mathrm{v}^{p}_{\beta}})}-\displaystyle\sum_{c\in
V}{{l^{c}_{\alpha}}{f^{c}_{\beta}}}\right]~{},$ (2)

with particles $p$, contacts $\mathrm{c}$, mass $\mathrm{m}^{p}$, velocity
$\mathrm{v}^{p}$, force $f^{c}$ and branch vector $l^{c}$, while Greek letters
represent components $r$, $\theta$, and $z$ [35]. The first term is the sum of
kinetic energy fluctuations, and the second term involves the dyadic product
of the contact-force with the contact-branch vector.

The quantity which describes the local configuration of a granular assembly is
the fabric tensor [40, 39] and is given by

${F}_{\alpha\beta}(\mathbf{r})=\displaystyle\frac{1}{V}\displaystyle\sum_{p\in
V}V^{p}\displaystyle\sum_{c\in p}{n^{c}_{\alpha}}{n^{c}_{\beta}}~{},$ (3)

where $V^{p}$ represents the particle volume which lies inside the averaging
volume $V$, $n^{c}$ is the normal unit branch-vector pointing from the center
of particle $p$ to contact $c$.

####  2.3.2 Isotropic and deviatoric parts

Any given tensor $\bm{Q}$ can be uniquely decomposed into isotropic and
deviatoric parts as

$\bm{Q}=Q_{V}\bm{I}+\bm{Q_{D}}$ (4)

with $Q_{V}=\frac{1}{3}\Tr{\bm{Q}}$ and the traceless deviator $\bm{Q_{D}}$.
The latter contains information about the eigensystem of $\bm{Q}$, that is
identical to the eigensystem of $\bm{Q_{D}}$ itself.

Let us assume $Q_{1}$, $Q_{2}$, and $Q_{3}$ are the eigenvalues of
$\bm{Q_{D}}$ sorted such that $Q_{1}\geq Q_{2}\geq Q_{3}$. Based on the
normalization, we use following definition to quantify the anisotropy of any
tensor $\bm{Q_{D}}$ using a single scalar quantity:

$Q_{\rm
dev}=\sqrt{\frac{{((Q_{1}-Q_{2})^{2}+(Q_{2}-Q_{3})^{2}+(Q_{3}-Q_{1})^{2})}}{6}}~{},$
(5)

the deviators $\dot{\epsilon}_{\mathrm{dev}}$, $\sigma_{\mathrm{dev}}$, and
$F_{\mathrm{dev}}$ refer to strain rate $\dot{\epsilon}_{\alpha\beta}$, stress
$\sigma_{\alpha\beta}$ and fabric $F_{\alpha\beta}$ respectively. Other
definitions of $Q_{\rm dev}$ are reported in [41].

$\dot{\gamma}=\dot{\epsilon}_{\mathrm{dev}}$ quantifies the local strain rate
magnitude, which is very close to the form defined in cylindrical coordinates
[42] as tested in [35]. The pressure $p$ is the isotropic hydrostatic stress,
while $\tau=\sigma_{\rm dev}$ quantifies objectively the shear stress. The
deviatoric stress ratio, $\mu={\tau}/{p}$, quantifies the “stress
anisotropy”or the macroscopic friction coefficient. The volumetric fabric
$3F_{\rm v}$ represents the contact number density, while the deviatoric
fabric $F_{\rm dev}$ quantifies the anisotropy of the contact network (as
studied in detail in [43]).

Table 1: Table showing the values of $g$ (units of $\rm ms^{-2}$) and particle stiffness $k_{n}$ (units of $\rm Nm^{-1}$) used in our simulations, and various time scales associated with the system, as discussed in the main text (in units of $\rm s$). The values of $T_{\dot{\gamma}}$ and $T_{p}$ are the average values reported at $z=2\langle d\rangle,H/2,H-2\langle d\rangle$ in the center of the shear band. $g$ | $\frac{\Omega}{2\pi}$ | $k_{n}$ | $T_{c}\times 10^{-3}$ | $T_{\eta}\times 10^{-3}$ | $T_{g}\times 10^{-2}$ | $T_{\dot{\gamma}}$ | $T_{p}\times 10^{-3}$
---|---|---|---|---|---|---|---
$0.5$ | 0.005 | 100 | $2$ | $5.6$ | $6.6$ | 25, 20, 10 | ($1.7$, $2.5$, $5$)
$1$ | 0.01 | 100 | $2$ | $5.6$ | $4.7$ | 10.9, 7.8, 2.7 | ($32$, $15.3$, $12.5$)
$2$ | 0.01 | 100 | $2$ | $5.6$ | $3.3$ | 10.7, 7.5, 2.7 | ($9$, $11$, $22$)
$5$ | 0.01 | 100 | $2$ | $5.6$ | $2.1$ | 10.3, 7.4,2.6 | ($5.9$, $7$, $14.6$)
5 | 0.01 | 500 | $0.1$ | $2.5$ | $2.1$ | 10.6, 7.5, 2.1 | ($5.1$, $7$, $14.1$)
$20$ | 0.01 | 100 | $2$ | $5.6$ | $1.05$ | 9.7, 7.0, 2.6 | ($2.9$, $3.8$, $8$)
20 | 0.01 | 400 | $10$ | $11.2$ | $1.05$ | 10, 7.1, 2.7 | ($2.9$, $3.6$, $7.4$)
$50$ | 0.01 | 100 | $2$ | $5.6$ | $0.66$ | 8.7, 6.6, 2.5 | ($1.8$, $2.2$, $4$)
50 | 0.01 | 1000 | $0.66$ | $18$ | $0.66$ | 10.1, 7.1, 2.6 | ($1.9$, $2.5$, $7$)
$10$ | 0.01 | 100 | $2$ | $5.6$ | $1.5$ | 9.9, 7.0, 2.6 | ($4$, $5.6$, $24$)
10 | 0.01 | 1000 | $0.66$ | $18$ | $1.5$ | 9.1, 8.1, 2.6 | ($4$, $5.6$, $27$)
10 | 0.01 | 10000 | $0.2$ | $0.56$ | $1.5$ | 10.7, 7.3, 2.8 | ($4$, $5.4$, $31$)
10 | 0.1 | 100 | $2$ | $5.6$ | $1.5$ | 1.1, 0.7, 0.23 | ($4$, $6$, $9$)
10 | 0.5 | 100 | $2$ | $5.6$ | $1.5$ | 0.2, 0.15, 0.05 | ($4$, $5$, $10$)
10 | 1.0 | 100 | $2$ | $5.6$ | $1.5$ | 0.1, 0.07, 0.02 | ($4$, $5$, $20$)
10 | 2.0 | 100 | $2$ | $5.6$ | $1.5$ | 0.02 0.03, 0.008 | ($4$, $6$, $18$)
Table 2:  Table showing the values of $g$ (units of $\rm ms^{-2}$) and particle stiffness $k_{n}$ (units of $\rm Nm^{-1}$) used in our simulations, and various dimensionless numbers as introduced in the text, as discussed in the main text (in units of $\rm s$). The average values of $I$ are reported at $z=2\langle d\rangle,H/2,H-2\langle d\rangle$ in the center of the shear band. $g$ | $\frac{\Omega}{2\pi}$ | $k_{n}$ | $I$ | $\kappa^{2}_{g}$
---|---|---|---|---
$0.5$ | 0.005 | 100 | ($7$, $12$, $50$) $\times 10^{-3}$ | $2\times 10^{-5}$
$1$ | 0.01 | 100 | ($1.3$, $2$, $4$) $\times 10^{-3}$ | $5\times 10^{-4}$
$2$ | 0.01 | 100 | ($0.75$, $2$, $1.4$) $\times 10^{-3}$ | $3.4\times 10^{-3}$
$5$ | 0.01 | 100 | ($2$, $0.7$, $0.9$) $\times 10^{-3}$ | $1\times 10^{-3}$
5 | 0.01 | 500 | ($2.5$, $0.8$, $1.4$) $\times 10^{-3}$ | $2\times 10^{-4}$
$20$ | 0.01 | 100 | ($1$, $0.5$, $0.9$) $\times 10^{-3}$ | $1\times 10^{-3}$
20 | 0.01 | 400 | ($2$, $4$, $7$) $\times 10^{-4}$ | $1\times 10^{-4}$
$50$ | 0.01 | 100 | ($8$, $3.4$, $7.2$) $\times 10^{-5}$ | $2.5\times 10^{-3}$
50 | 0.01 | 1000 | ($5$, $3$, $6$) $\times 10^{-5}$ | $2.5\times 10^{-4}$
$10$ | 0.01 | 100 | ($1.1$, $0.6$, $1.5$) $\times 10^{-3}$ | $5\times 10^{-4}$
10 | 0.01 | 1000 | ($1.4$, $0.5$, $1.6$) $\times 10^{-3}$ | $5\times 10^{-5}$
10 | 0.01 | 10000 | ($1.2$, $0.6$, $0.9$) $\times 10^{-3}$ | $5\times 10^{-6}$
10 | 0.1 | 100 | ($1.5$, $0.6$, $0.9$) $\times 10^{-2}$ | $5\times 10^{-4}$
10 | 0.5 | 100 | ($7.5$, $2.6$, $6$) $\times 10^{-2}$ | $5\times 10^{-4}$
10 | 1.0 | 100 | ($1.5$, $0.5$, $10.2$) $\times 10^{-1}$ | $5\times 10^{-4}$
10 | 2.0 | 100 | 0.3, 0.16, 1.51 | $5\times 10^{-4}$

###  2.4 Time scales

We characterize the dynamics of the system looking at different time scales.
At first, we define two microscopic time scales related to the contact
duration and the viscous damping coefficient between two particles in contact,
respectively, as

$T_{k}=\sqrt{\frac{\langle m\rangle}{k_{n}}}~{},\hskip
14.22636ptT_{\eta}=\frac{\langle m\rangle}{\eta_{n}}~{},$ (6)

where $\langle m\rangle$ is the mass of a particle with mean diameter. $T_{k}$
and $T_{\eta}$ can be related to contact time
$T_{c}=\frac{\pi}{\sqrt{\frac{1}{T^{2}_{k}}-\frac{1}{T^{2}_{\eta}}}}$. Next,
two time scales associated with external forces, i.e. the gravity and external
rotation rate, can be introduced as

$T_{g}=\sqrt{\frac{\langle d\rangle}{g}}~{},\hskip
14.22636ptT_{\Omega}=\frac{2\pi}{\Omega}~{},$ (7)

respectively, where $T_{g}$ is the time taken by a particle with zero initial
velocity to fall a distance $\langle d\rangle/2$.

The time scales, (6) and (7), are functions of material constants and applied
external forces, hence, are constants throughout the system. In this sense,
the time scales, $T_{c}$, $T_{\eta}$, $T_{g}$, and $T_{\Omega}$, are _global_.
On the other hand, two local time scales related to the local shear rate
$\dot{\gamma}$ and pressure $p$, as introduced in [18, 44] are:

$T_{\dot{\gamma}}=\frac{1}{\dot{\gamma}}~{},\hskip 14.22636ptT_{p}=\langle
d\rangle\sqrt{\frac{\rho}{p}}~{}.$ (8)

As shown in the following sections, the spatial distributions of pressure and
shear rate are inhomogeneous due to gravity and shear band localization.
Therefore, in contrast to the global time scales, $T_{\dot{\gamma}}$ and
$T_{p}$ are _local_ field variables that depend on spatial position.

The time scales can be combined to formulate dimensionless numbers that give
indications of dominance of one of the time scales. For example, the inertial
number, as introduced by [44, 45, 46],

$I\equiv T_{p}/T_{\dot{\gamma}}=\dot{\gamma}\langle
d\rangle/\sqrt{p/\rho}~{},$ (9)

provides an estimate of the local rapidity of the flow. For $I\ll 1$, the flow
is _quasi-static_ , where particles interact via enduring contacts and
inertial effects are negligible. For $I\sim 1$, the flow is in the _dense
inertial regime_ , and for $I\gg 1$, the flow is in the rapid, collisional gas
like state. A variant of inertial number $I_{k}=\dot{\gamma}\langle
d\rangle/\sqrt{p_{kin}/\rho}$ can be defined by using only the kinetic
pressure instead of the total pressure. Other expressions such as Savage or
Coulomb number ${\dot{\gamma}}^{2}{\langle d\rangle}^{2}\rho/{p}$ is simply
the square of the inertial number $I$ [47].

A dimensionless parameter _global compressibility_ can be introduced as the
ratio between $T_{g}$ and $T_{k}$:

$\kappa_{g}\equiv\frac{T_{k}}{T_{g}}=\sqrt{\frac{\langle m\rangle
g}{k_{n}\langle d\rangle}}~{},$ (10)

providing a global measure of compressibility of the bulk material. A high
$\kappa_{g}$ signifies that the bulk material is compressible, which comes
from either very high confinement by strong gravity or from low contact
stiffness at particle level. On the other hand, when $\kappa_{g}$ is small,
the average overlap is very small compared to the particle diameter, which
means that the bulk material is closer to being the rigid limit. A similar
dimensionless parameter can be introduced as:

$\kappa_{p}\equiv\frac{T_{k}}{T_{p}}=\sqrt{\frac{p\langle
d\rangle}{k_{n}}}~{},$ (11)

which estimates the compressibility of the material at the local scale.

Table 1 shows typical values of timescales in our simulations, and table 2
reports various dimensionless numbers. We observe that for flows with a
rotation rate $\Omega/2\pi=0.01\ \rm s^{-1}$ and the gravity $g\geq 1\
\rm{ms^{-2}}$, the inertial number $I$ is well below $1$ for all values of
$k_{n}$. For lower values of gravity $g=0.5\rm ms^{-2}$, the rotation rate is
chosen to be $\Omega/2\pi=0.005\ s^{-1}$, such that $I$ stays in the same
range. From this table we infer that the systems for wide a range of $g$ and
$k_{n}$ can be safely assumed to be in the rate-independent quasi-static
state. However, we observe that with increase in rate of rotation
$\Omega/2\pi$, $I$ begins to increase and the flow behavior enters into the
rate dependent inertial regime.

##  3 Quasistatic rheology

In this section, we present our results on the analysis of macroscopic
rheology in the quasi-static state. At first, in section 3.1 we study the
steady state and critical state and the amount of time system requires to
reach the critical state. We explore the effect of gravity and stiffness on
the macroscopic friction coefficient in section 3.2. We also show the results
of local volume fraction in section 3.3, and connect the rheology to the
microscopic structure tensor in section 3.4. We will extend our analysis to
dense inertial flows in section 4.

###  3.1 Steady state and critical state

Shearing leads to dilation and a build-up of shear stress and anisotropy in
fabric. Additionally, anisotropy in the fabric of the granular medium needs a
finite amount of strain to build up. Since, we are interested in the steady
state flow properties, we need to ask: what is the minimum time or
equivalently strain required to reach the steady state flow regime?

We define $\varphi=\frac{180}{2\pi}\Omega\Delta t$ (where $\Delta t$ is the
simulation time) as a shear parameter, which is the same as the amount, in
degrees, through which the system rotates. In general, the time required for
stabilization can depend on the considered variable [25]. In the following, we
consider the relaxation of various quantities to the steady state (both
locally and globally averaged). As we expect this relaxation to be slowest for
small $\Omega$, we perform these tests for $\Omega/2\pi=0.01\ \rm s^{-1}$ that
is the smallest $\Omega$ we explore in our simulations.

At first, we analyze the global quantities (averaged over the whole system)
like the averaged kinetic energy and average number of contacts. We find that
they relax very fast ($\varphi\sim 5$) and are not reliable since they only
represent the fastest relaxation in the system.

Next, we analyze the relaxation of local quantities. In order to estimate the
relaxation time for local quantities, we analyze local velocity profiles. We
find that such a relaxation requires ($\varphi\sim 30$ or more) which is in
agreement with Ries et al[22]. Consequently we consider that the condition to
approach the steady state is $\varphi\geq 30$. Experimental results from
Wortel et al[48] also show that the system needs to be rotated by an angle of
approximately 30 degrees in order to reach a reasonable steady state. When the
steady state is reached, the local averaging is performed over almost 600
snapshots distributed over $\varphi\geq 30$.

The consistency of the local averaged quantities also depends on the
accumulated local shear strain during the procurement of data. We concentrate
our interest in the region where the system can be considered in a critical
state, which occurs at large enough shear strain $\gamma=\dot{\gamma}\Delta
t\geq 1$. In other words, the system can be assumed to be in the critical
state. The critical state is a unique state reached after long shear, where
material deforms with applied strain without any change in normal stress,
shear stress and volume fraction, and the system forgets its preparation
history [49].

Figure 2 shows the local shear stress $\tau(r,h)$ plotted against the local
pressure $p(r,h)$. We observe that for a given pressure, $\tau$ is higher for
larger $\dot{\gamma}$, however for $\dot{\gamma}>\dot{\gamma}_{c}$ (with
$\dot{\gamma}_{c}\approx 0.08\ {\rm s^{-1}}$), $\tau$ becomes almost
independent of the local strain rate. This means that $\tau/p$ is almost
constant for all data points with strain rate $\dot{\gamma}>\dot{\gamma}_{c}$.
In other words, if the dimensionless shear length
$l_{\gamma}=t_{av}\dot{\gamma}$ [35] exceeds unity (where $t_{av}$ is the time
over which averaging of the data is performed), layers can be assumed to be
sheared by almost a particle diameter. For $l_{\gamma}\geq 1$, the shear
deformation can be considered to be fully established, and the system is
assumed to be in the critical state [35].

In the same setup, Ries et al[22] and Szabó et al[50] also found that after
long enough shear, the material inside the shear band reaches the critical
state and can be characterized by estimating local accumulated strain
$\gamma\geq 1$. Our previous works [35, 51] also showed that for rotation rate
$\Omega/2\pi=0.01\ \rm s^{-1}$, $\dot{\gamma}_{c}\approx 0.1\ {\rm s^{-1}}$ is
the shear rate above which the shear band is well established. Since we are
interested in the flow behavior of the material, as default in the rest of the
paper we focus only on the data well inside the shear band with local strain
rate,

$\dot{\gamma}(r,h)>\dot{\gamma}_{c}(\Omega)\equiv\frac{10\Omega}{2\pi}~{},$
(12)

unless specified otherwise.

Figure 2: (Color online) The local shear stress, $\tau(r,h)$, plotted against
the local pressure, $p(r,h)$, for different values of the local shear rate,
$\dot{\gamma}(r,h)$ as given in the legend (in units of $\mathrm{{s}^{-1}}$),
with gravity, $g=10\ \rm{ms^{-2}}$ and $\Omega/2\pi=0.01\ \rm s^{-1}$.

###  3.2 Friction coefficient

We will now turn our attention towards the effect of gravity and softness on
the macroscopic friction coefficient $\tau/p$ in the quasi-static state. In
previous studies, the friction coefficient has been assumed to be independent
of both the particle stiffness and gravity. However, particles used in these
were extremely rigid. Few studies were performed systematically investigating
the dependence of the flow behavior on gravity [10, 11, 12, 13, 14, 15].

Figures 3 and 3 show shear stress-pressure curves for different values of
normal stiffness $k_{n}$ and external gravity $g$, respectively. In these
plots, for the sake of clarity, both shear stress and pressure are plotted
only at the center position $R_{c}$ of the shear band ($R_{c}$ being the the
position at which the local shear rate is maximum). For a better comparison,
both shear stress and pressure are normalized with the maximum pressure
$p_{\rm max}$ reached in the simulation with particular $k_{n}$ and $g$ (so
that both abscissa and ordinates are of the same order). We observe that both
softness of the particles (interpreted as opposite of contact stiffness) and
external gravity drastically affect the shear stress. Moreover, they act in
the same direction as for a given pressure $\tau$ decreases with increase in
either particle softness or the external gravity. We also see that with
increasing softness or gravity, the relation between $\tau$ and $p$ becomes
non-linear.

Figure 3: (Color online) (Left) The local shear stress plotted against the
local pressure for different values of the normal stiffness as given in the
legend in units of $\mathrm{Nm^{-1}}$. Here, the gravity is fixed to $g=10\
\mathrm{ms^{-2}}$. (Right) The local shear stress plotted against the local
pressure for different values of the gravity as given in the legend in units
of $\mathrm{ms^{-2}}$. Here, the normal stiffness is fixed to $k_{n}=100\
\mathrm{Nm^{-1}}$. Both $\tau(r,h)$ and $p(r,h)$ are scaled by the maximum
pressure $p_{\mathrm{max}}(r,h)$, respectively.

####  3.2.1 Linear approximation

To understand the dependence of the macroscopic friction coefficient in a
quasi-static state on the softness and gravity, we estimate it as the slope of
a linear fitting function for the shear stress against pressure in the same
fashion as the Mohr-Coulomb failure criterion, i.e.

$\tau(r,h)\simeq\mu_{0}^{\mathrm{global}}p(r,h)~{},$ (13)

where $\mu_{0}^{\mathrm{global}}$ is the _global friction coefficient_ which
depends neither on the shear rate nor on the pressure.

Figure 4 displays the global friction coefficient $\mu_{0}^{\mathrm{global}}$
plotted against gravity $g$ for different values of the normal stiffness
$k_{n}$, as given in the legend. We observe that $\mu_{0}^{\mathrm{global}}$
decreases with increasing gravity, while it increases with increasing $k_{n}$.
Figure 4 shows the global friction coefficient with different values of the
normal stiffness $k_{n}$ and gravity $g$, where all results of
$\mu_{0}^{\mathrm{global}}$ are collapsed if we plot them against
$\kappa=\kappa^{2}_{g}$ ($\kappa_{g}$ is given by (10)).

In figure 4, the solid line is given by

$\mu_{0}^{\mathrm{global}}=\mu_{r}^{\mathrm{global}}-\left(\frac{\kappa}{\kappa_{0}}\right)^{\alpha}~{},$
(14)

where $\mu_{r}^{\mathrm{global}}=0.17\pm 0.01$ is the global friction
coefficient in the rigid particle limit and the exponent and characteristic
global compressibility are given by $\alpha\simeq 0.5$ and $\kappa_{0}\simeq
2.01$, respectively.

Previous microgravity parabolic flight and centrifuge experiments [10, 12, 13,
15, 16, 14] showed a similar decreasing trend of the macroscopic friction
coefficient on gravity. Few authors [12, 15] attributed this dependence to the
fact that at low gravity, the body forces become weak and the electrostatic
cohesive forces begin to dominate. Klein et al[12] also argued that a load-
dependent interparticle friction coefficient might be responsible for this
behavior. However, no cohesive force or load-dependent friction was
implemented in any of the DEM simulation data presented here. Hence, we claim
that there should be an additional mechanism responsible for this interesting
behavior. In order to gain a better understanding in the following, we have a
closer look by studying the system locally.

Figure 4: (Color online) The global friction coefficient,
$\mu_{0}^{\mathrm{global}}$, plotted against (left) gravity $g$, and (right)
the global compressibility, $\kappa=mg/(k_{n}\langle d\rangle)$, on a log-
linear scale for different values of the normal stiffness as shown in the
legend. The solid line represent (14).

####  3.2.2 Local compressibility

Since our system is heterogeneous in both stress and strain rate fields, a
local description of the system is highly desirable. In the shear stress-
pressure ($\tau-p$) curves for different softness and gravity (figure 3), the
dependence of shear stress on pressure slightly “bends” with increasing
softness and gravity, i.e., the friction coefficient has a dependence on the
pressure and the shear stress becomes a nonlinear function of pressure as:

$\tau(r,h)=\mu_{0}^{\mathrm{local}}(p,k_{n})p(r,h)~{},$ (15)

where $\mu_{0}^{\mathrm{local}}(p,k_{n})$ is a _local friction coefficient_
which depends on pressure and contact stiffness.

Figure 5 shows the local friction coefficient with different values of the
normal stiffness and gravity, where all results of
$\mu_{0}^{\mathrm{local}}(p,k_{n})$ collapse if we introduce the _local
compressibility_ ,

$p^{\ast}\equiv\kappa^{2}_{p}=\frac{p\langle d\rangle}{k_{n}}~{},$ (16)

defined as square of the ratio between two time scales, $T_{k}$ and $T_{p}$.
$p^{\ast}$ can also be interpreted as non-dimensional average overlap (scaled
with mean particle diameter). In this figure, we scanned through a wide range
of $p^{*}$ by systematically varying $g$ and $k_{n}$. We observe that for
$p^{*}<5\times 10^{-4}$, $\mu_{0}^{\mathrm{local}}(p^{\ast})$ is almost
constant, while for higher values, $\mu_{0}^{\mathrm{local}}(p^{\ast})$
decreases with $p^{*}$ up to $p^{\ast}\approx 0.1$.

This dependence can be written in the form,

$\mu_{0}^{\mathrm{local}}(p^{\ast})=\mu_{\mathrm{r}}^{\mathrm{local}}-\left(\frac{p^{*}}{p^{*}_{\sigma}}\right)^{\beta_{1}},$
(17)

where $\mu_{\mathrm{r}}^{\mathrm{local}}=0.172\pm 0.01$ is the value of
macroscopic friction in the rigid limit, which is in fair agreement with
contact dynamics shear simulations for the same particle friction coefficient
[52]. The exponent is found to be $\beta_{1}\approx 0.5\pm 0.04$ and the
characteristic local compressibility is ${p^{*}_{\sigma}}=10.1\pm 0.2$. As one
extreme of $p^{*}$, for $p^{*}=0.1$ the average overlap is almost $10\%$
relative to the mean particle diameter, i.e., the _soft particle limit_. The
upper bound of $\mu_{0}^{\mathrm{local}}(p^{\ast})$ is the low compression
case, i.e., _the rigid limit_ , where the average overlap is much smaller
(1/10000) compared to the particle diameter and
$\mu_{0}^{\mathrm{local}}(p^{\ast})$ is almost double as large as for
$p^{*}\approx 0.1$.

Figure 5: (Color online) The local friction coefficient,
$\mu_{0}^{\mathrm{local}}(p^{\ast})$, plotted against the local
compressibility, $p^{\ast}$, on a log-linear scale. Different symbols
represent different values of $\kappa$ as given in the legend of figure 6,
while the solid line represents (17).

###  3.3 Local volume fraction

In figure 6, the local volume fraction $\nu(r,h)$ is plotted against the local
compressibility, $p^{*}$. Because of slow quasi-static flows, no strong
dilation is observed, i.e., no strong dependence of $\nu$ on local shear rate
is observed. We observe that the packing is rather loose for lower $p^{*}$ and
tends to a critical value $\nu_{c}=0.642\pm 0.002$, as observed in [53]. The
data can be fitted well by the functional form

$p^{*}=a^{*}\left(\nu-\nu_{\rm{c}}\right),$ (18)

with $a^{*}=0.48\pm 0.02$ ($a^{*}$ can be further expressed in terms of
volumetric fabric as reported in [40, 54]). Interestingly, no significant
difference in volume fraction $\nu$ is observed for $p^{*}<10^{-3}$, while for
$p^{*}>10^{-3}$ within the fluctuations, $\nu$ increases almost linearly with
$p^{*}$ (the curvature is due to the logarithmic $p^{*}$ axis). The relation
between $\nu$ and $p^{*}$ is well established in the case of static packings
[40, 54, 55]. Here we show that the same relation holds for a slow granular
flow, involving considerable finite but small strain rates.

Figure 6: (Color online) The local volume fraction, $\nu(r,h)$, in the system
plotted against the local compressibility, $p^{*}$, on a log-linear scale.
Different symbols represent different values of $\kappa$ as given in the
legend. The solid line represents (18).

###  3.4 Local structure

Shearing of a granular assembly always leads to the buildup of a contact
anisotropy in the system [56, 57, 58]. To study contact anisotropy in our
system, we analyze the deviatoric fabric as defined in (3). We use the second
invariant or $F_{\rm dev}$ (5) of the deviatoric fabric tensor to quantify
anisotropy of the contact network in the system.

####  3.4.1 Anisotropy

Figure 7 displays the local deviatoric fabric, $F_{\mathrm{dev}}(r,h)$,
plotted against the local compressibility $p^{\ast}$, where
$F_{\mathrm{dev}}(r,h)$ for different values of the particle stiffness and
gravity is found to collapse on a unique curve (solid line). This dependence
can be written in a similar fashion to (17),

$F_{\rm dev}(p^{\ast})=F^{r}_{\rm
dev}-\left(\frac{p^{\ast}}{p^{\ast}_{F}}\right)^{\beta_{2}},$ (19)

where $F^{r}_{\rm dev}$ is the anisotropy of contact network in the rigid
limit, the exponent is found to be $\beta_{2}\approx 0.5\pm 0.03$, and
$p^{\ast}_{F}\approx 26.3\pm 0.6$. The decrease in $F_{\rm dev}$ with
increasing $p^{\ast}$ can be explained in terms of the increasing volume
fraction $\nu(r,h)$ with increase in $p^{\ast}$. In the case of a denser
packing, particles have less free space to re-arrange (and thus align along
preferential direction) and build up anisotropy in response to the local
shear, compared to a relatively loose packing (at low $p^{\ast}$).

Figure 7: (Color online) The local deviatoric fabric, $F_{\rm dev}(r,h)$,
plotted against the local compressibility, $p^{*}$, on a log-linear scale.
Different symbols represent different values of $\kappa$ as given in the
legend of figure 6. The solid line represents the corresponding fit to (19).

In figures 5 and 7, we observe that the the local shear resistance and the
local anisotropy of the contact network in the quasi-static state show a
similar trend. In figure 8, we plot $\mu_{0}^{\mathrm{local}}(p^{\ast})$
against $F_{\rm dev}(p^{\ast})$ for different values of $\kappa$, where a
linear correlation can be inferred as,

$\mu_{0}^{\mathrm{local}}(p^{\ast})=\mu_{\rm{iso}}+bF_{\rm dev}(p^{\ast})~{},$
(20)

where $\mu_{\rm{iso}}=0.01\pm 0.01(\approx 0)$ is the friction coefficient in
the (extrapolated) limit of the isotropic contact network $(F_{\rm dev}=0)$
and $b=1.38\pm 0.02$ is a constant of proportionality. This clearly shows that
the shear resistance seems to accompany the anisotropy in the contact network
in the critical state. It is worthwhile to mention that no information about
such a relation in the transient regime can be inferred from our analysis,
while it is strongly supported by our data in the critical state. It also
links the increase in the friction coefficient with decreasing gravity (as
observed in figure 4) to the change in the contact network configuration of
the material at lower $g$.

Figure 8: (Color online) The local friction coefficient,
$\mu_{0}^{\mathrm{local}}(p^{\ast})$, plotted against the local deviatoric
fabric, $F_{\rm dev}(r,h)$, for different values of $\kappa$. Different
symbols represent different values of $\kappa$ as given in the legend of
figure 6. The solid line represents the corresponding fit to (20).

####  3.4.2 Shape factor

In this section, we compare the shape factor ($\frac{Q_{2}}{Q_{1}}$) for
stress and fabric tensors, where $Q_{2}$, and $Q_{1}$ are the eigenvalues of
the deviatoric tensors as defined in section 2.3.2. The shape factor
quantifies the out of shear plane neutral eigenvalue, compared to the maximum
eigenvalue. Note that for strain rate tensor
$\dot{\epsilon}_{2}/\dot{\epsilon}_{1}\equiv 0$ due to geometry and symmetry,
i.e., we have plain-strain shear.

In figure 9, we plot the shape factor for the stress tensor. Different symbols
represent different values of compressibility $\kappa$ as given in the legend
of figure 6, while black circles show the data in the center of shear bands
for these simulations. We observe that the shape factor is highly fluctuating
for the two extremes of $p^{*}$, while in the range $10^{-4}<p^{*}<10^{-2}$,
it is significantly below zero. This implies that the stress in the shear
plane is higher as compared to the axial stress along the neutral direction
orthogonal to the shear plane. The sign means that this axial stress is
reduced perpendicular $-1$ and enhanced $+1/2$ within the shear plane [39].
The dashed line is the data from [39], where authors studied the flow behavior
on an inclined plane. We observe that the sign for both the cases is negative,
while the magnitude is different, which might be due to the difference in
interparticle friction. In figure 9, the shape factor of the fabric tensor
fluctuates (strongly) around zero. It is important to mention that the
fluctuations in the data are from a single simulation.

These two observations suggest that the fabric and stress tensors behave
differently even though they are proportional in magnitude (norm) as shown in
figure 8. The fabric tensor is in a planar state like the strain rate tensor,
whereas the induced stress state is triaxial, as expected for a solid-like
material [39]. $F_{2}/F_{1}$ tends to positive values for larger $p^{*}$,
further establishing the difference between structure and stress tensors.
However, in order to have a clear picture for the fabric tensor, the strong
and weak subnetworks should be studied separately, since only the strong
subnetwork carries almost all of the fabric anisotropy [43, 59].

As discussed in section 3.1 the cutoff shear rate $\dot{\gamma}_{c}$ can
depend on the simulation time or the averaging time. In this section, we
focused on the data only inside the shear band, which are in the critical
state and have forgotten their initial configuration due to large strain.
However, the velocity gradients in the setup are smooth, which implies that
part of the system outside the shear band is also flowing, but only very
slowly. If the simulation runs longer, the cutoff can be lowered, eventually
if the simulation would run infinitely long, no cutoff on the local strain
rate is needed. If we lower the cutoff on local strain rate (see next
section), by setting $\dot{\gamma}_{c}(\Omega)\equiv\frac{\Omega}{2\pi}$, we
observe much wider variation of data in our results, especially for the local
friction coefficient, the deviatoric fabric and the volume fraction. However,
the shape factors are not strongly influenced by a change in
$\dot{\gamma}_{c}$, by nearly an order of magnitude.

Figure 9: (Color online) Shape factor for (left) stress, and (right) fabric
plotted against dimensionless pressure $p^{*}$. Different symbols represent
different values of $\kappa$ as given in the legend of figure 6. Black circles
represent the data in the center of the shear band. Solid line represents
zero, while dashed line is the prediction from [39].

##  4 Complete rheology

The previous section showed that in the quasi-static state the friction
coefficient and deviatoric fabric are strongly correlated in the critical
state. Motivated by this, we check if this correlation also works in the rate-
dependent inertial regime. To test the correlation, high inertial number data
are generated by varying the external rotation rate $\Omega$ for a fixed
gravity and contact stiffness. In the following, we will explore the evolution
of the local macroscopic friction coefficient $\mu$ and deviatoric fabric with
$I$ and combine it with the dependence of both $\mu$ and deviatoric fabric on
$p^{*}$, to propose the complete rheological law. It is important to mention
that in this section the cutoff on local strain rate is lowered to
$\dot{\gamma}_{c}(\Omega)\equiv\frac{\Omega}{2\pi}$, so as to capture the
maximum data and present a unique rheology law outside and inside the shear
band.

###  4.1 Friction law

The local friction coefficient $\mu$ is plotted against inertial number $I$ in
figure 10. Different symbols show data from different rates of rotation as
given in the inset; the black solid circles represent the data in the center
of the shear band. The solid black line shows the friction law, as proposed in
[45]:

$\mu(I,p^{*})=\mu_{0}^{\mathrm{local}}(p^{*})+\frac{\mu_{2}^{\mathrm{local}}(p^{*})-\mu_{0}^{\mathrm{local}}(p^{*})}{1+I^{\sigma}_{0}/I}~{},$
(21)

with $\mu_{0}^{\mathrm{local}}(p^{*})$, as given in (17). We observe that data
from simulation using a single gravity and contact stiffness does not give a
wide variation in $\mu$ and $\mu_{0}^{\mathrm{local}}=0.14$,
$\mu_{2}^{\mathrm{local}}=0.5$ and $I^{\sigma}_{0}=0.1$ fits well the data. A
Taylor expansion (in the range $I<I^{\sigma}_{0}$ ) to above equation is
$\mu(I)\cong\mu_{0}^{\mathrm{local}}+({\mu_{2}^{\mathrm{local}}-\mu_{0}^{\mathrm{local}}})\frac{I}{I^{\sigma}_{0}}$,
which is similar to the frictional law proposed in [18, 44]. Two different
trends emerge, i.e., the shear band center data can be very well fitted by
(21) and for $I\geq 0.01$ data collapse on a unique curve. On the other hand,
for lower values of $I$, deviations from this relation are observed, depending
on the external rotation rate. The friction coefficient in slow flows (steady
state) becomes smaller than $\mu_{0}^{\mathrm{local}}$, i.e. in our system the
granular material is able to flow below $\mu_{0}^{\mathrm{local}}$. The
deviation of our data from the main law (21) is consistent with observations
in [25, 60], where this deviation was explained based on the heterogeneity in
the stress field (due to strain rate). In our system, we have gradients in
stress arising due to gradients in both strain rate and pressure (due to
gravity).

In order to quantify the deviation from (21), we fit the data with:

$\mu(I<I^{*},p^{*})=\mu_{0}^{\mathrm{local}}(p^{*})\left[1-\alpha\ln{(I/I^{*})}\right]~{},$
(22)

where $\alpha$ is a constant and $I^{*}$ is the characteristic inertial number
when $\mu\cong\mu_{0}^{\mathrm{local}}$. This relation is similar as proposed
in [25] for two-dimensional ring shear cell data. As the above relation was
derived for a two-dimensional setup with constant pressure, we fit it to our
data at three different heights (i.e. pressure levels), close to top, at mid-
height and close to bottom. In figure 10, different colored dashed lines
represent this fit at the mid-height of the system. We observe that the
prediction is in close agreement with the data, even though our system has
different dimension and boundary conditions. A shows the data and
corresponding fits for different heights (pressures). We find that both
$\alpha$ and $I^{*}$ do not show a clear dependence on pressure.

Figure 10: (Color online) The friction coefficient plotted against the
inertial number $I$ for results from simulations with different rates of
rotation. Different symbols represent different rates of rotation as given in
the legend, lines with the same color represent solution of (22) with
$\mu_{0}^{\mathrm{local}}=0.14$, $\mu_{2}^{\mathrm{local}}=0.3$ and
$I^{\sigma}_{0}=0.026$. Dot-dashed line shows the Taylor expansion of (22).
Black circles represent the data in the center of the shear band.

###  4.2 Fabric anisotropy

In order to look for the connection between anisotropic fabric and macroscopic
friction coefficient in the inertial regime, here we explore the dependence of
$F_{\rm dev}$ on $I$. In figure 11, we plot the local $F_{\rm dev}$ as
obtained by simulations with different rates of rotation against $I.$ We
observe that like $\mu$, $F_{\rm dev}$ varies strongly against $I$ and its
dependence on $I$ can be represented as:

$F_{\rm dev}(I,p^{*})=F_{\rm dev}^{0}(p^{*})+\frac{F_{\rm
dev}^{(2)}(p^{*})-F_{\rm dev}^{0}(p^{*})}{1+I^{F}_{0}/I}~{},$ (23)

with $F_{\rm dev}^{0}$ being the fabric anisotropy in the quasi-static state,
$F_{\rm dev}^{(2)}$ is the limiting fabric anisotropy, and $I^{F}_{0}$ is the
typical inertial number, which is an order of magnitude different from
$I^{\sigma}_{0}$. Green, red and black lines show the fit to above relation at
pressure levels 100, 200 and 400 $\rm{Nm^{-2}}$, respectively, with points
inside the shear band highlighted (black circles). Fit parameters to these
results are presented in table 3. It is noticeable that unlike $\mu$, $I$
alone is not able to describe $F_{\rm dev}$, with the effect of pressure being
prominent in case of slow flows i.e., low $I$. In contrast, for fast flows,
the anisotropy seems to become independent of pressure.

The increase in the contact anisotropy with inertial number is in accordance
with the recent studies [39, 58]. It is important to mention that for even
higher rates of rotation of the system, i.e. inertial number $I>0.1$, $F_{\rm
dev}$ shows a different behavior as predicted by (23) and a decreasing trend
is observed (as reported in [61]), which is beyond the scope of this work.
This might be due to the fact that for $I>0.2$ the packing becomes very loose
($\nu\leq 0.55$). Also for such high rates of rotation, the centrifugal force
on grains due to rotation becomes comparable to the gravitational force. As a
result, the top surface is not flat anymore, instead the surface develops a
dip in the middle, as observed previously [61, 62, 63]. In this range, also
the kinetic and contact contributions of the local macroscopic friction $\mu$
become comparable.

Figure 11: (Color online) The local fabric anisotropy $F_{\rm dev}$ plotted against the inertial number $I$ for results from simulations with different rates of rotation. Different symbols represent different rates of rotation as given in the legend. Black circles represent the data in the center of the shear band. Green, red, and black lines are fit to (23) for pressure levels $p=100,200,$ and $400$ $\rm{Nm^{-2}}$ respectively, with fit parameters given in Table 3. The arrow shows increasing pressure.  $p$ | $F_{\rm dev}^{0}$ | $F_{\rm dev}^{(2)}$ | $I^{F}_{0}$
---|---|---|---
100 | 0.1$\pm$ 0.0005 | 0.17$\pm$ 0.0005 | 0.012
200 | 0.095$\pm$ 0.0008 | 0.17$\pm$ 0.0001 | 0.011
400 | 0.085$\pm$ 0.0001 | 0.17$\pm$ 0.0004 | 0.009
Table 3: Table showing the fit parameters $F_{\rm dev}^{0}$, $F_{\rm
dev}^{(2)}$, and $I^{F}_{0}$ in (23) for different values of pressure $p$ (in
units of ${\rm Nm^{-2}}$).

Starting from both variations of macroscopic friction and fabric anisotropy as
a function of inertial number $I$, it is tempting to ask the question if the
correlation in (20) holds for the inertial regime as well. The result is
displayed in figure 12. Solid line shows (20), which fits well the shear band
center data being highlighted by black circles. Here again, we find data
outside the shear band showing a different trend. Red dashed line separates
the data into quasi-static and dense inertial regimes. It is interesting to
note that (20) (as shown by black solid line) very well captures the trend at
the onset of dense inertial regime. However, for faster flows, where rate
dependence becomes prominent, a different trend is observed which can also be
fitted by a slightly different linear relation (dot-dashed line). This shows
that the macroscopic friction and fabric anisotropy are correlated even in the
dense inertial regime.

Figure 12: (Color online) $\mu$ plotted against $F_{\rm dev}$ for data in
previous section and different rates of rotation. Different symbols represent
different rates of rotation as given in the legend. The solid line represent
(20), while the dashed line (with slope $\backsimeq 3.5$) is fit by the eye.
Black circles represent the data in the center of the shear band.

##  5 Discussion and Conclusion

To summarize, this work is an exploration of the flow behavior for 3D granular
shear flow using DEM simulations. We particularly focused on the effect of
external compression (gravity) and the particle softness on the flow behavior,
considering both stress and structure.

##### Quasi-static flows

Our study shows that the shear strength (macroscopic friction coefficient
$\mu$) of the material decreases with increase in either gravity or particle
softness for the quasi-static flows. We find that the data for different
gravity and particle softness can be expressed as a unique power law, when
analyzed in terms of a control parameter called local compressibility
$p^{\ast}$ (as defined in (16)). $p^{\ast}$ can be interpreted as the non-
dimensional local average overlap (scaled with mean particle diameter) and can
be used to quantify the softness of the bulk material relative to the local
compression level. Low values of $p^{\ast}$ signify rigid particles, while a
high $p^{\ast}$ implies soft, deformable particles. Both softness and gravity
are also found to affect the local microstructure, i.e., the anisotropy of the
contact network, which is quantified by the deviatoric fabric ($F_{\rm dev}$).
We show that the deviatoric fabric can also be expressed as a power law of
$p^{\ast}$ with an exponent similar to that of the shear strength. This points
out that the local anisotropy of the contact network (deviatoric fabric) and
the shear strength of the material are highly correlated in the slow, quasi-
static flows and the shear strength follows the anisotropy of the contact
network.

These results are highly significant for planetary studies regarding the shear
strength of the granular material in extraterrestrial bodies such as Moon or
Mars. Significant amount of experimental work using parabolic flight have
shown the increase in shear strength of the material with decreasing gravity
[10, 11, 12, 13, 14, 15] without proper explanation. We propose that the
anisotropy, i.e., the rearrangement in the contact network, is the key
parameter for this dependence if no new forces are involved. With decreasing
gravity, the packing becomes loose (due to decrease in body force acting on
the particles), which in turn provides more free space to the particles to
rearrange (and thus align) in response to the local shear. The fact that the
particle softness and gravity have been shown to have similar effects on the
local flow behavior, makes this work equally relevant for soft particles, that
find their applications in many engineering and biological systems [17].

Since it is extremely difficult and expensive to perform in situ experiments
on the Moon (or the parabolic flight), the ‘compensation’ effect we find with
the ratio of gravity $g$ and particle stiffness $k_{n}$ allows to mimic
variable gravity by tweaking/tuning the particle stiffness.

##### Inertial flows

Further, to study the rheology of the system for gravity $g=10\ \rm{ms^{-2}}$,
the rotation rate of the system is increased. For faster flows the system
enters into a rate dependent inertial regime, consistent with previous studies
[18, 44, 46]. We find that both the macroscopic friction $\mu$ and contact
anisotropy show a similar increasing behavior with inertial number $I$. This
shows that macroscopic friction and contact anisotropy are correlated also in
the dense inertial regime. The increase of $\mu$ with $I$ as observed in the
inertial regime is accompanied by the evolution of the microstructure (contact
anisotropy) with increasing inertial number. This picture is consistent with
the recent study of Azéma et al [58].

##### Open issues

We find that the frictional laws obtained from homogeneous shear flows [44]
can be applied locally in the inertial regime, while they fail to predict the
behavior of the material in the slower, quasi-static regime. The local
rheology laws can be applied to our data in the center of the shear bands,
where the strain-rate and stress gradients are zero, hence the material can be
assumed to be homogeneously sheared. However, away from the center of the
shear bands, in the quasi-static regime, we observe a nearly identical range
of $\mu$ values corresponding to a completely different range of $I$. We find
that in our system the material is able to flow even below
$\mu_{0}^{\mathrm{local}}$, but only very slowly. We have shown that some
observations can be explained by using an approach similar to Koval et al[25].
These deviations might be well captured using the non-local models by Kamrin
and coworkers [60, 64, 65]; this work is in progress. Another related issue
which remains untouched is the effect of particle softness and external
compression (gravity here) on the non-locality. A study of effect of gravity
on primary and secondary velocity fields, as done recently in [66, 67] also
deserve a further study.

##### Conclusion

The macroscopic friction (shear strength) of the material is found to be
affected not only by the local shear rate, but also by external compression
(gravity) and softness of the particles. While traditionally the inertial
number, the ratio of stress to strain-rate time-scales, is dominating the flow
rheology, we find that a second dimensionless number, the ratio of softness
and stress time scales, must be involved to characterize the bulk flow
behavior. For very slow shear rate the former can be ignored, while the latter
two affect the shear strength by decreasing it with increase in either gravity
(and thus local pressure) or particle softness. For faster flows, the
macroscopic friction is found to increase in general with increasing shear
rate. However, the tails of shear bands feature an anomalously small
macroscopic friction - as observed previously [25, 60, 65]. For the dependence
of macroscopic friction on above three quantities, the change in local
microstructure (contact anisotropy) is found to be a key parameter, that was
not often considered yet.

Looking towards the future, we are now in a position to address various
important issues, such as unexpectedly high shear strength of the material at
low (confining) stress or reduced gravity and a direct relation between the
contact anisotropy and the shear strength of the material. These issues are
vital for a better explanation of the macroscopic behavior of the granular
systems from a microscopic observation. The current study dealt with a dense
system with small interparticle friction ($\mu_{\rm{p}}=0.01$), where the
effect of softness on the macroscopic behavior is more direct than for large
$\mu_{\rm{p}}$. However, an issue which remains unanswered and will be an
extension of this study is whether the same effect can also be observed for
relatively loose system (with higher interparticle friction). Further, the
question whether the correlation between contact anisotropy and shear strength
is just a consequence of relatively low interparticle friction or if it will
also hold for a more realistic material (with higher interparticle friction)
remains to be answered.

##  6 Acknowledgements

We would like to thank C R K Windows-Yule for his careful revision. We thank
W. Losert, D. van der Meer, M. Sperl, D. Lohse, B. Tighe, P. Jop, H. Hayakawa,
A. Ikeda, O. I. Imole, and L. Silbert for stimulating discussions. Financial
support through the “Jamming and Rheology” project of the Stichting voor
Fundamenteel Onderzoek der Materie (FOM), which is financially supported by
the “Nederlandse Organisatie voor Wetenschappelijk Onderzoek” (NWO), is
acknowledged.

##  Appendix A Pressure dependence of local macroscopic friction

In this section, we explore the pressure dependence of our rheological laws as
presented in section 4. Figure 13 shows the fits for three different pressure
levels (height) and we find for all of them the fitting is in good agreement
with the data. Figure 14 shows fit parameters for different pressure levels,
where we observe that none of them show a clear strong dependence on pressure.

Figure 13: (Color online) $\mu$ plotted against $I$ for different local
pressures in the system (a) $p=100$, (b) $p=200$, and (c) $p=400$
$\mathrm{Nm^{-2}}$.

Figure 14: (Color online) (a) $\alpha$, and (b) $I^{*}$ plotted against the
external rotation rate for different local pressures in the system.

##  Appendix B Pressure dependence of correlation

In this section, we test the correlation between $\mu$ and $F_{\rm dev}$ that
was presented in section 4. Figure 15 shows the data for three different
pressure levels, we observe that the correlation holds very well for all
rotation rates, except for the fastest rotation rate, which seems to fall off
from the prediction of (20).

Figure 15: (Color online) $\mu$ plotted against $F_{\rm dev}$ for different
local pressures in the system (a) $p=100$, (b) $p=200$, and (c) $p=400$
$\mathrm{Nm^{-2}}$. The solid line represent the corresponding fit to (20),
while the dashed line is the best fit to the data.

## References

## References

  * [1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer.  Granular solids, liquids, and gases.  Reviews of Modern Physics, 68(4):1259–1273, 1996.
  * [2] J. Duran.  Sands, Powders, and Grains-An Introduction to the Physics of Granular Materials.  Springer, New York, 2000.
  * [3] K Krohn, R Jaumann, K Otto, T Hoogenboom, R Wagner, DL Buczkowski, B Garry, DA Williams, RA Yingst, J Scully, et al.  Mass movement on vesta at steep scarps and crater rims.  Icarus, 2014.
  * [4] Troy Shinbrot, N-H Duong, L Kwan, and MM Alvarez.  Dry granular flows can generate surface features resembling those seen in martian gullies.  Proceedings of the National Academy of Sciences of the United States of America, 101(23):8542–8546, 2004.
  * [5] Michael C Malin.  Mass movements on venus: Preliminary results from magellan cycle 1 observations.  Journal of Geophysical Research: Planets, 97(E10):16337–16352, 1992\.
  * [6] Keith A Howard.  Avalanche mode of motion: Implications from lunar examples.  Science, 180(4090):1052–1055, 1973.
  * [7] Gabriel I. Tardos, Sean McNamara, and Ilkay Talu.  Slow and intermediate flow of a frictional bulk powder in the Couette geometry.  Powder Technology, 131(1):23–39, 2003.
  * [8] Charles S. Campbell.  Granular shear flows at the elastic limit.  Journal of Fluid Mechanics, 465:261–291, 2002.
  * [9] Bruno Andreotti, Yoël Forterre, and Olivier Pouliquen.  Granular media: Between Fluid and Solid.  Cambridge University Press, 2013.
  * [10] Emir J Macari-Pasqualino, Stein Sture, and Kenneth Runesson.  Analysis of low effective stress characteristics of granular materials in reduced gravity.  In Geotechnical Engineering Congress 1991, pages 1222–1233. ASCE, 1991.
  * [11] Nicholas C Costes, VC Janoo, and S Sture.  Microgravity experiments on granular materials.  Material Research Society Symposium Proceedings, 87:203–212, 1987\.
  * [12] S. P. Klein and B. R. White.  Dynamic shear of granular material under variable gravity conditions.  AIAA Journal, 28:1701–1702, 1990.
  * [13] Khalid A Alshibli, Stein Sture, and Nicholas C Costes.  Constitutive and stability behavior of soils in microgravity environment.  In Space Technology and Applications International Forum-2000, volume 504, pages 246–252. AIP Publishing, 2000.
  * [14] Khalid A Alshibli, Nicholas C Costes, and Ronald F Porter.  Mechanics of granular materials (mgm).  In SPIE’s 1996 International Symposium on Optical Science, Engineering, and Instrumentation, pages 303–310. International Society for Optics and Photonics, 1996.
  * [15] Stein Sture, Nicholas C Costes, Susan N Batiste, Mark R Lankton, Khalid A AlShibli, Boris Jeremic, Roy A Swanson, and Melissa Frank.  Mechanics of granular materials at low effective stresses.  Journal of Aerospace Engineering, 11(3):67–72, 1998.
  * [16] Antje Brucks, Tim Arndt, Julio M. Ottino, and Richard M. Lueptow.  Behavior of flowing granular materials under variable $g$.  Phys. Rev. E, 75:032301, Mar 2007.
  * [17] Allan S Hoffman.  Hydrogels for biomedical applications.  Advanced drug delivery reviews, 54:3–12, 2002.
  * [18] G. D. R. MiDi.  On dense granular flows.  The European Physical Journal E: Soft Matter and Biological Physics, 14:341–365, 2004.
  * [19] Michio Otsuki, Hisao Hayakawa, and Stefan Luding.  Behavior of pressure and viscosity at high densities for two-dimensional hard and soft granular materials.  Progress of Theoretical Physics Supplement, 184:110–133, 2010.
  * [20] Hernán A. Makse, David L. Johnson, and Lawrence M. Schwartz.  Packing of Compressible Granular Materials.  Phys. Rev. Lett., 84, 2000.
  * [21] Kasper van der Vaart, Yasser Rahmani, Rojman Zargar, Zhibing Hu, Daniel Bonn, and Peter Schall.  Rheology of concentrated soft and hard-sphere suspensions.  Journal of Rheology (1978-present), 57:1195–1209, 2013.
  * [22] A. Ries, D. E. Wolf, and T. Unger.  Shear zones in granular media: Three-dimensional contact dynamics simulations.  Phys. Rev. E, 76:051301, 2007.
  * [23] Leonardo E Silbert, Gary S Grest, Robert Brewster, and Alex J Levine.  Rheology and contact lifetimes in dense granular flows.  Physical review letters, 99:068002, 2007.
  * [24] Tim Arndt, Antje Brucks, Julio M Ottino, and Richard M Lueptow.  Creeping granular motion under variable gravity levels.  Physical Review E, 74:031307, 2006.
  * [25] Georg Koval, Jean-Noël Roux, Alain Corfdir, and Franifmmode Chevoir.  Annular shear of cohesionless granular materials: From the inertial to quasistatic regime.  Phys. Rev. E, 79:021306, 2009.
  * [26] P. A. Cundall and O. D. L. Strack.  Modeling of Microscopic Mechanisms in Granular Materials.  In J. T. Jenkins and M. Satake, editors, Mechanics of Granular Materials: New Models and Constitutive Relations, pages 137–149, Amsterdam, 1983\. Elsevier.
  * [27] O. R. Walton.  Numerical simulation of inclined chute flows of monodisperse, inelastic, frictional spheres.  Mechanics of Materials, 16:239–247, 1993.
  * [28] S. Luding, M. Lätzel, and H. J. Herrmann.  From discrete element simulations towards a continuum description of particulate solids.  In A. Levy and H. Kalman, editors, Handbook of Conveying and Handling of Particulate Solids, pages 39–44, Amsterdam, The Netherlands, 2001\. Elsevier.
  * [29] Stefan Luding.  Introduction to discrete element methods: basic of contact force models and how to perform the micro-macro transition to continuum theory.  European Journal of Environmental and Civil Engineering, 12:785–826, 2008.
  * [30] S. Luding.  Cohesive, frictional powders: contact models for tension.  Granular Matter, 10:235–246, 2008.
  * [31] D. Fenistein and M. van Hecke.  Kinematics – Wide shear zones in granular bulk flow.  Nature, 425:256, 2003.
  * [32] D. Fenistein, J. W. van de Meent, and M. van Hecke.  Universal and wide shear zones in granular bulk flow.  Phys. Rev. Lett., 92:094301, 2004.
  * [33] S. Luding.  The effect of friction on wide shear bands.  Particulate Science and Technology, 26, 2008.
  * [34] J. A. Dijksman and M. van Hecke.  Granular flows in split-bottom geometries.  Soft Matter, 6:2901–2907, 2010.
  * [35] Stefan Luding.  Constitutive relations for the shear band evolution in granular matter under large strain.  Particuology, 6:501–505, 2008.
  * [36] Abhinendra Singh, Vanessa Magnanimo, Kuniyasu Saitoh, and Stefan Luding.  Effect of cohesion on shear banding in quasistatic granular materials.  Phys. Rev. E, 90:022202, 2014.
  * [37] K. Bagi.  Microstructural stress tensor of granular assemblies with volume forces.  J. Appl. Mech., 66:934–936, 1999.
  * [38] M. Lätzel, S. Luding, and H. J. Herrmann.  Macroscopic material properties from quasi-static, microscopic simulations of a two-dimensional shear-cell.  Granular Matter, 2:123–135, 2000.
  * [39] T. Weinhart, R. Hartkamp, A. R. Thornton, and S. Luding.  Coarse-grained local and objective continuum description of three-dimensional granular flows down an inclined surface.  Physics of Fluids, 25, 2013.
  * [40] O. I. Imole, N. Kumar, V. Magnanimo, and S. Luding.  Hydrostatic and shear behavior of frictionless granular assemblies under different deformation conditions.  KONA, 30:84–108, 2013.
  * [41] N Kumar, S Luding, and V Magnanimo.  Macroscopic model with anisotropy based on micro-macro informations.  Acta Mechanica, 225(8):2319–2343, 2014.
  * [42] M. Depken, W. van Saarloos, and M. van Hecke.  Continuum approach to wide shear zones in quasistatic granular matter.  Phys. Rev. E, 73:031302, 2006.
  * [43] Olukayode I Imole, Mateusz Wojtkowski, Vanessa Magnanimo, and Stefan Luding.  Micro-macro correlations and anisotropy in granular assemblies under uniaxial loading and unloading.  Physical Review E, 89:042210, 2014.
  * [44] Frédéric da Cruz, Sacha Emam, Michaël Prochnow, Jean-Noël Roux, and François Chevoir.  Rheophysics of dense granular materials: Discrete simulation of plane shear flows.  Physical Review E, 72:021309, 2005.
  * [45] P. Jop, Y. Forterre, and O. Pouliquen.  A constitutive law for dense granular flows.  Nature, 441:727–730, 2006.
  * [46] Yoël Forterre and Olivier Pouliquen.  Flows of dense granular media.  Annu. Rev. Fluid Mech., 40:1–24, 2008.
  * [47] Stuart B Savage.  The mechanics of rapid granular flows.  Advances in applied mechanics, 24:289–366, 1984.
  * [48] Geert H Wortel, Joshua A Dijksman, and Martin van Hecke.  Rheology of weakly vibrated granular media.  Physical Review E, 89:012202, 2014.
  * [49] David Muir Wood.  Soil behaviour and critical state soil mechanics.  Cambridge university press, 1990.
  * [50] Balázs Szabó, János Török, Ellák Somfai, Sandra Wegner, Ralf Stannarius, Axel Böse, Georg Rose, Frank Angenstein, and Tamás Börzsönyi.  Evolution of shear zones in granular materials.  Physical Review E, 90:032205, 2014.
  * [51] A. Singh, V. Magnanimo, and S. Luding.  Effect of friction and cohesion on anisotropy in quasi-static granular materials under shear.  In A. Yu and S. Luding, editors, Powders and Grains, AIP Conf. Proc., volume 1542, pages 682–685, 2013.
  * [52] Tamás Unger.  Collective rheology in quasi static shear flow of granular media, September 2010.
  * [53] Stéphane Dorbolo, T Scheller, F Ludewig, G Lumay, and N Vandewalle.  Influence of a reduced gravity on the volume fraction of a monolayer of spherical grains.  Physical Review E, 84:041305, 2011.
  * [54] F. Göncü, O. Duran, and S. Luding.  Constitutive relations for the isotropic deformation of frictionless packings of polydisperse spheres.  C. R. Mécanique, 338:570–586, 2010.
  * [55] H. P. Zhang and H. A. Makse.  Jamming transition in emulsions and granular materials.  Phys. Rev. E, 72:011301, July 2005.
  * [56] C. T. Veje, D. W. Howell, and R. P. Behringer.  Kinematics of a 2D granular Couette experiment at the transition to shearing.  Phys. Rev. E, 59:739–745, 1999.
  * [57] T. S. Majmudar and R. P. Behringer.  Contact force measurements and stress-induced anisotropy in granular materials.  Nature, 435:1079–1082, 2005.
  * [58] Emilien Azéma and Farhang Radjaï.  Internal structure of inertial granular flows.  Physical Review Letters, 112:078001, 2014.
  * [59] A Singh, V Magnanimo, and S Luding.  Effect of friction on the force distribution in sheared granular materials.  In Proceedings of NUMGE2014, pages 409–414. CRC Press/Balkema, NL, 2014.
  * [60] Ken Kamrin and Georg Koval.  Nonlocal Constitutive Relation for Steady Granular Flow.  Phys. Rev. Lett., 108:178301, 2012.
  * [61] Abhinendra Singh.  Micro-macro and rheology in sheared granular matter.  Universiteit Twente, Enschede, 2014.
  * [62] Joshua Albert Dijksman et al.  Granular media: flow & agitations.  Granular and Disordered Media, Leiden Institute of Physics, Faculty of Science, Leiden University, 2009.
  * [63] Eric I Corwin.  Granular flow in a rapidly rotated system with fixed walls.  Physical Review E, 77:031308, 2008.
  * [64] David L. Henann and Ken Kamrin.  A predictive, size-dependent continuum model for dense granular flows.  Proceedings of the National Academy of Sciences, 110:6730–6735, 2013.
  * [65] Ken Kamrin and Georg Koval.  Effect of particle surface friction on nonlocal constitutive behavior of flowing granular media.  Computational Particle Mechanics, 1:169–176, 2014.
  * [66] N. Murdoch, B. Rozitis, S.F. Green, T.-L. Lophem, P. Michel, and W. Losert.  Granular shear flow in varying gravitational environments.  Granular Matter, 15:129–137, 2013.
  * [67] N. Murdoch, B. Rozitis, K. Nordstrom, S. F. Green, P. Michel, T.-L. de Lophem, and W. Losert.  Granular convection in microgravity.  Phys. Rev. Lett., 110:018307, 2013.