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

###### 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 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 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 , where , , , and 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 , where , , , and 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 is switched to the sliding force , being the particle friction coefficient when exceeds the critical value, i.e. (with ) [30].

To study the effect of particle softness on macroscopic behaviors, we explore a range of normal contact stiffness , from . While changing , the time step for numerical integration is adjusted to be times the contact duration to ensure accurate dynamic integration [30]. When is changed, and are changed as well, to ensure that the coefficient of restitution remains unchanged.

### 2.2 Split-bottom ring shear cell

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 , , and , respectively, where the concentric cylinders rotate relative to each other at a rate 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 polydispersed spherical particles with density up to height . The average size of particles is mm, and the width of the homogeneous size-distribution (with ) is . The cylindrical walls and the bottom are roughened due to some (about 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 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 and three regimes can be observed as reported in [32, 34]. We keep , 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 . 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 direction. The averaging is thus performed over toroidal volumes over many snapshots in time, leading to fields as function of the radial and vertical positions. Here, the averaging is performed with spacings and 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

(1) |

where Greek letters represent components radial distance , azimuthal angle , and height , while represents the gradient in cylindrical coordinate [35].

The stress tensor is given by

(2) |

with particles , contacts , mass , velocity , force and branch vector , while Greek letters represent components , , and [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.

#### 2.3.2 Isotropic and deviatoric parts

Any given tensor can be uniquely decomposed into isotropic and deviatoric parts as

(4) |

with and the traceless deviator . The latter contains information about the eigensystem of , that is identical to the eigensystem of itself.

Let us assume , , and are the eigenvalues of sorted such that . Based on the normalization, we use following definition to quantify the anisotropy of any tensor using a single scalar quantity:

(5) |

the deviators , , and refer to strain rate , stress and fabric respectively. Other definitions of are reported in [41].

quantifies the local strain rate magnitude, which is very close to the form defined in cylindrical coordinates [42] as tested in [35]. The pressure is the isotropic hydrostatic stress, while quantifies objectively the shear stress. The deviatoric stress ratio, , quantifies the “stress anisotropy”or the macroscopic friction coefficient. The volumetric fabric represents the contact number density, while the deviatoric fabric quantifies the anisotropy of the contact network (as studied in detail in [43]).

0.005 | 100 | 25, 20, 10 | (, , ) | ||||

0.01 | 100 | 10.9, 7.8, 2.7 | (, , ) | ||||

0.01 | 100 | 10.7, 7.5, 2.7 | (, , ) | ||||

0.01 | 100 | 10.3, 7.4,2.6 | (, , ) | ||||

5 | 0.01 | 500 | 10.6, 7.5, 2.1 | (, , ) | |||

0.01 | 100 | 9.7, 7.0, 2.6 | (, , ) | ||||

20 | 0.01 | 400 | 10, 7.1, 2.7 | (, , ) | |||

0.01 | 100 | 8.7, 6.6, 2.5 | (, , ) | ||||

50 | 0.01 | 1000 | 10.1, 7.1, 2.6 | (, , ) | |||

0.01 | 100 | 9.9, 7.0, 2.6 | (, , ) | ||||

10 | 0.01 | 1000 | 9.1, 8.1, 2.6 | (, , ) | |||

10 | 0.01 | 10000 | 10.7, 7.3, 2.8 | (, , ) | |||

10 | 0.1 | 100 | 1.1, 0.7, 0.23 | (, , ) | |||

10 | 0.5 | 100 | 0.2, 0.15, 0.05 | (, , ) | |||

10 | 1.0 | 100 | 0.1, 0.07, 0.02 | (, , ) | |||

10 | 2.0 | 100 | 0.02 0.03, 0.008 | (, , ) |

0.005 | 100 | (, , ) | ||
---|---|---|---|---|

0.01 | 100 | (, , ) | ||

0.01 | 100 | (, , ) | ||

0.01 | 100 | (, , ) | ||

5 | 0.01 | 500 | (, , ) | |

0.01 | 100 | (, , ) | ||

20 | 0.01 | 400 | (, , ) | |

0.01 | 100 | (, , ) | ||

50 | 0.01 | 1000 | (, , ) | |

0.01 | 100 | (, , ) | ||

10 | 0.01 | 1000 | (, , ) | |

10 | 0.01 | 10000 | (, , ) | |

10 | 0.1 | 100 | (, , ) | |

10 | 0.5 | 100 | (, , ) | |

10 | 1.0 | 100 | (, , ) | |

10 | 2.0 | 100 | 0.3, 0.16, 1.51 |

### 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

(6) |

where is the mass of a particle with mean diameter. and can be related to contact time . Next, two time scales associated with external forces, i.e. the gravity and external rotation rate, can be introduced as

(7) |

respectively, where is the time taken by a particle with zero initial velocity to fall a distance .

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, , , , and , are global. On the other hand, two local time scales related to the local shear rate and pressure , as introduced in [18, 44] are:

(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, and 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],

(9) |

provides an estimate of the local rapidity of the flow. For , the flow is quasi-static, where particles interact via enduring contacts and inertial effects are negligible. For , the flow is in the dense inertial regime, and for , the flow is in the rapid, collisional gas like state. A variant of inertial number can be defined by using only the kinetic pressure instead of the total pressure. Other expressions such as Savage or Coulomb number is simply the square of the inertial number [47].

A dimensionless parameter global compressibility can be introduced as the ratio between and :

(10) |

providing a global measure of compressibility of the bulk material. A high 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 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:

(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 and the gravity , the inertial number is well below for all values of . For lower values of gravity , the rotation rate is chosen to be , such that stays in the same range. From this table we infer that the systems for wide a range of and can be safely assumed to be in the rate-independent quasi-static state. However, we observe that with increase in rate of rotation , 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 (where 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 , we perform these tests for that is the smallest 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 () 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 ( or more) which is in agreement with Ries et al[22]. Consequently we consider that the condition to approach the steady state is . 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 .

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 . 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 plotted against the local pressure . We observe that for a given pressure, is higher for larger , however for (with ), becomes almost independent of the local strain rate. This means that is almost constant for all data points with strain rate . In other words, if the dimensionless shear length [35] exceeds unity (where is the time over which averaging of the data is performed), layers can be assumed to be sheared by almost a particle diameter. For , 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 . Our previous works [35, 51] also showed that for rotation rate , 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,

(12) |

unless specified otherwise.

### 3.2 Friction coefficient

We will now turn our attention towards the effect of gravity and softness on the macroscopic friction coefficient 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 and external gravity , respectively. In these plots, for the sake of clarity, both shear stress and pressure are plotted only at the center position of the shear band ( 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 reached in the simulation with particular and (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 decreases with increase in either particle softness or the external gravity. We also see that with increasing softness or gravity, the relation between and becomes non-linear.

#### 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.

(13) |

where is the global friction coefficient which depends neither on the shear rate nor on the pressure.

Figure 4 displays the global friction coefficient plotted against gravity for different values of the normal stiffness , as given in the legend. We observe that decreases with increasing gravity, while it increases with increasing . Figure 4 shows the global friction coefficient with different values of the normal stiffness and gravity , where all results of are collapsed if we plot them against ( is given by (10)).

In figure 4, the solid line is given by

(14) |

where is the global friction coefficient in the rigid particle limit and the exponent and characteristic global compressibility are given by and , 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.

#### 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 () 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:

(15) |

where 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 collapse if we introduce the local compressibility,

(16) |

defined as square of the ratio between two time scales, and . can also be interpreted as non-dimensional average overlap (scaled with mean particle diameter). In this figure, we scanned through a wide range of by systematically varying and . We observe that for , is almost constant, while for higher values, decreases with up to .

This dependence can be written in the form,

(17) |

where 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 and the characteristic local compressibility is . As one extreme of , for the average overlap is almost relative to the mean particle diameter, i.e., the soft particle limit. The upper bound of 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 is almost double as large as for .

### 3.3 Local volume fraction

In figure 6, the local volume fraction is plotted against the local compressibility, . Because of slow quasi-static flows, no strong dilation is observed, i.e., no strong dependence of on local shear rate is observed. We observe that the packing is rather loose for lower and tends to a critical value , as observed in [53]. The data can be fitted well by the functional form

(18) |

with ( can be further expressed in terms of volumetric fabric as reported in [40, 54]). Interestingly, no significant difference in volume fraction is observed for , while for within the fluctuations, increases almost linearly with (the curvature is due to the logarithmic axis). The relation between and 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.

### 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 (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, , plotted against the local compressibility , where 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),

(19) |

where is the anisotropy of contact network in the rigid limit, the exponent is found to be , and . The decrease in with increasing can be explained in terms of the increasing volume fraction with increase in . 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 ).

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 against for different values of , where a linear correlation can be inferred as,

(20) |

where is the friction coefficient in the (extrapolated) limit of the isotropic contact network and 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 .

#### 3.4.2 Shape factor

In this section, we compare the shape factor () for stress and fabric tensors, where , and 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 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 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 , while in the range , 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 and enhanced 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]. tends to positive values for larger , 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 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 , 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 , by nearly an order of magnitude.

## 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 for a fixed gravity and contact stiffness. In the following, we will explore the evolution of the local macroscopic friction coefficient and deviatoric fabric with and combine it with the dependence of both and deviatoric fabric on , to propose the complete rheological law. It is important to mention that in this section the cutoff on local strain rate is lowered to , 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 is plotted against inertial number 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]:

(21) |

with , as given in (17). We observe that data from simulation using a single gravity and contact stiffness does not give a wide variation in and , and fits well the data. A Taylor expansion (in the range ) to above equation is , 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 data collapse on a unique curve. On the other hand, for lower values of , deviations from this relation are observed, depending on the external rotation rate. The friction coefficient in slow flows (steady state) becomes smaller than , i.e. in our system the granular material is able to flow below . 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:

(22) |

where is a constant and is the characteristic inertial number when . 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 and do not show a clear dependence on pressure.

### 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 on . In figure 11, we plot the local as obtained by simulations with different rates of rotation against We observe that like , varies strongly against and its dependence on can be represented as:

(23) |

with being the fabric anisotropy in the quasi-static state, is the limiting fabric anisotropy, and is the typical inertial number, which is an order of magnitude different from . Green, red and black lines show the fit to above relation at pressure levels 100, 200 and 400 , 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 , alone is not able to describe , with the effect of pressure being prominent in case of slow flows i.e., low . 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 , 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 the packing becomes very loose (). 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 become comparable.

100 | 0.1 0.0005 | 0.17 0.0005 | 0.012 |

200 | 0.095 0.0008 | 0.17 0.0001 | 0.011 |

400 | 0.085 0.0001 | 0.17 0.0004 | 0.009 |

Starting from both variations of macroscopic friction and fabric anisotropy as a function of inertial number , 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.

## 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 ) 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 (as defined in (16)). 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 signify rigid particles, while a high 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 (). We show that the deviatoric fabric can also be expressed as a power law of 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 and particle stiffness allows to mimic variable gravity by tweaking/tuning the particle stiffness.

##### Inertial flows

Further, to study the rheology of the system for gravity , 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 and contact anisotropy show a similar increasing behavior with inertial number . This shows that macroscopic friction and contact anisotropy are correlated also in the dense inertial regime. The increase of with 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 values corresponding to a completely different range of . We find that in our system the material is able to flow even below , 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 (), where the effect of softness on the macroscopic behavior is more direct than for large . 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.

## Appendix B Pressure dependence of correlation

In this section, we test the correlation between and 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).

## 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 Congress1991, 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 . 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.