Skip to main content Accessibility help
×
Home
Hostname: page-component-768ffcd9cc-b9rrs Total loading time: 1.766 Render date: 2022-12-04T16:43:12.917Z Has data issue: true Feature Flags: { "useRatesEcommerce": false } hasContentIssue true

Viscometric flow of dense granular materials under controlled pressure and shear stress

Published online by Cambridge University Press:  20 November 2020

Ishan Srivastava*
Affiliation:
Sandia National Laboratories, Albuquerque, NM87185, USA Center for Computational Sciences and Engineering, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Leonardo E. Silbert
Affiliation:
School of Math, Science, and Engineering, Central New Mexico Community College, Albuquerque, NM87106, USA
Gary S. Grest
Affiliation:
Sandia National Laboratories, Albuquerque, NM87185, USA
Jeremy B. Lechman*
Affiliation:
Sandia National Laboratories, Albuquerque, NM87185, USA
*
Email addresses for correspondence: isriva@lbl.gov, jblechm@sandia.gov
Email addresses for correspondence: isriva@lbl.gov, jblechm@sandia.gov
Rights & Permissions[Opens in a new window]

Abstract

This study examines the flow of dense granular materials under external shear stress and pressure using discrete element method simulations. In this method, the material is allowed to strain along all periodic directions and adapt its solid volume fraction in response to an imbalance between the internal state of stress and the external applied stress. By systematically varying the external shear stress and pressure, the steady rheological response is simulated for: (1) rate-independent quasi-static flow; and (2) rate-dependent inertial flow. The simulated flow is viscometric with non-negligible first and second normal stress differences. While both normal stress differences are negative in inertial flows, the first normal stress difference switches from negative to slightly positive, and second normal stress difference tends to zero in quasi-static flows. The first normal stress difference emerges from a lack of coaxiality between a second-rank contact fabric tensor and strain rate tensor in the flow plane, while the second normal stress difference is linked to an excess of contacts in the shear plane compared with the vorticity direction. A general rheological model of second order (in terms of strain rate tensor) is proposed to describe the two types of flow, and the model is calibrated for various values of interparticle friction from simulations on nearly monodisperse spheres. The model incorporates normal stress differences in both regimes of flow and provides a complete viscometric description of steady dense granular flows.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2020. Published by Cambridge University Press

1. Introduction

Granular flows exhibit several intriguing phenomena that distinguish them from Newtonian fluids, such as the presence of pressure-dependent arrest and flow onset (yield) criteria leading to rate-independent and rate-dependent flows, and a dilute gas-like flow dominated by inelastic particle collisions. A convenient classification defines three distinct types of granular flows (Forterre & Pouliquen Reference Forterre and Pouliquen2008): (1) quasi-static flows; (2) dense inertial flows; and (3) gas-like collisional flows. The behaviour of the three types of granular flows is quite diverse and several constitutive models have been proposed for their description; however, a general constitutive model applicable across all flow types remains elusive. In this work we focus our attention on the rate-independent and rate-dependent flows, where the particle contact lifetimes are relatively long, inertia is important and the material is predominantly dense.

The rate-independent flow regime has been described by various constitutive models, largely inspired by the principles of solid mechanics and plasticity. Beginning with the incipient failure hypothesis of the Coulomb yield criterion (Sokolovskii Reference Sokolovskii1965), further observations of critical state deformation in soils led to the development of several rigid-plastic and elasto-plastic models based on critical state theory and associated plasticity (Schofield & Wroth Reference Schofield and Wroth1968). Recent developments have attempted to include material anisotropy in granular plasticity by introducing state-dependence of material stress, often characterized via material texture or fabric (Sun & Sundaresan Reference Sun and Sundaresan2011; Li & Dafalias Reference Li and Dafalias2012; Gao et al. Reference Gao, Zhao, Li and Dafalias2014). The rate-independent granular plasticity has also been characterized by double shearing models (Spencer Reference Spencer1964) that relax the assumption of homogeneous deformations to explain shear banding along slip planes in granular materials. Further advances on these models have introduced the concepts of dilatancy (Mehrabadi & Cowin Reference Mehrabadi and Cowin1978), work hardening (Anand & Gu Reference Anand and Gu2000) and more recently fabric anisotropy (Nemat-Nasser Reference Nemat-Nasser2000; Zhu, Mehrabadi & Massoudi Reference Zhu, Mehrabadi and Massoudi2006). The reader is referred to a recent review of various constitutive models of rate-independent regime in granular flows (Radjai, Roux & Daouadji Reference Radjai, Roux and Daouadji2017).

Rate-dependent granular flows were first observed to exhibit quadratic scaling of shear and normal stress with strain rate at a constant volume (Bagnold Reference Bagnold1954), which was later verified in several experiments and simulations (Pouliquen Reference Pouliquen1999; Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Da Cruz et al. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005; Lois, Lemaître & Carlson Reference Lois, Lemaître and Carlson2005). Recently, a Bingham-type $\mu (I)$ rheological model for granular materials at moderate shear rates has been proposed (Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006), which attempts to connect the rate-independent and rate-dependent granular flow regimes by introducing pressure as a control variable instead of volume, although the two can be interchanged based on the one-to-one relationship between $\mu$ and $I$ at moderate shear rates. Here $\mu$ is the dimensionless shear stress ratio, and $I$ is the dimensionless inertial number (described later in the text). A detailed discussion of such viscoplastic models can be found in a recent review (Goddard Reference Goddard2014).

Although these rheological models have successfully predicted granular flow profiles in a remarkable number of geometries (MiDi Reference MiDi2004), several rheological effects remain unexplained, such as surface curvature in free-surface flows (Couturier et al. Reference Couturier, Boyer, Pouliquen and Guazzelli2011; McElwaine, Takagi & Huppert Reference McElwaine, Takagi and Huppert2012), negative rod climbing effects (Boyer, Pouliquen & Guazzelli Reference Boyer, Pouliquen and Guazzelli2011b), anomalous stress profile in Couette flows (Mehandia, Gutam & Nott Reference Mehandia, Gutam and Nott2012), and the observation of shear-free sheets in split-bottom Couette flows (Depken et al. Reference Depken, Lechman, Van Hecke, Van Saarloos and Grest2007). Many of these effects arise from a lack of coaxiality between principal directions of stress and strain rate tensors in viscometric flows (Alam & Luding Reference Alam and Luding2003, Reference Alam, Luding, Garcia-Rojo, Herrmann and McNamara2005; Rycroft, Kamrin & Bazant Reference Rycroft, Kamrin and Bazant2009; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013; Saha & Alam Reference Saha and Alam2016; Bhateja & Khakhar Reference Bhateja and Khakhar2018; Seto & Giusteri Reference Seto and Giusteri2018), which is the operating assumption in several constitutive models. As such, there is a need for higher-order constitutive models that incorporate these effects to provide better predictions of granular rheology. Furthermore, microstructural origins of these rheological effects, especially in the dense and quasi-static flow regimes, remain unclear.

In this paper, we describe fully stress-controlled discrete element method (DEM) simulations in both rate-independent and rate-dependent regimes. This simulation method enables the evolution of all strain degrees of freedom of a fully periodic representative volume element of granular material in response to external applied shear stress and pressure. The novelty of this simulation method is fourfold: (1) it naturally captures the pressure-dependent flow-onset (yield) and flow-arrest phenomena (Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019); (2) by prescribing shear stress rather than shear rate, this method can seamlessly traverse across rate-dependent and rate-independent regimes of granular flow; (3) by prescribing pressure rather than solid volume fraction, shear-induced dilation of granular materials is fully captured; and (4) the fully periodic nature of the simulations is devoid of any boundary effects, and thus represents the true bulk response of granular materials to applied stresses. The stress-controlled method is used to simulate shear flows of nearly monodisperse spheres. A second-order rheological model that does not assume coaxiality of stress and strain rate tensors is proposed, and the model is calibrated for various values of interparticle friction from the simulation data. For viscometric flows, the second-order rheological effects result in non-negligible first and normal stress differences. We provide microstructural insights into the origin of normal stress differences by analysing a second-rank contact fabric tensor.

The paper is organized as follows. Section 2 on model and methods describes the second-order rheological model, introduces the stress-controlled DEM simulation method, and provides details on the sphere–sphere contact mechanics model and general simulation parameters. Section 3 provides evidence that the steady flow generated in these simulations by applying shear stress and pressure is viscometric in nature. Section 4 describes the calibration of the rheological model based on the viscometric flow data from DEM simulations. Section 5 describes the normal stress differences measured in these simulations and their microstructural origins.

2. Model and methods

2.1. Rheological model

In anticipation of the results presented below, we introduce a purely dissipative rheological framework proposed by Goddard (Reference Goddard1984), which was utilized to formulate constitutive laws for rate-independent and rate-dependent flow in granular materials (Goddard Reference Goddard1986, Reference Goddard2014). In this framework, stress emerges from dissipation through macroscopic bulk deformation, which dominates over grain-scale inertial relaxation. This is similar to the dense flow of granular materials at low inertial numbers, which is of interest here. Furthermore, elastic effects are ignored, and non-hydrostatic stress components emerge entirely from granular flow, and are zero when there is no flow. In this framework, the Cauchy stress tensor $\boldsymbol {\sigma }$ is given by

(2.1)\begin{equation} \boldsymbol{\sigma} = \boldsymbol{\mathcal{\eta}}\{\mathcal{H}\}:\boldsymbol{D}, \end{equation}

where $\boldsymbol {D}$ is the symmetric strain rate tensor, and $\boldsymbol {\mathcal {\eta }}\{\mathcal {H}\}$ is a positive-definite fourth-rank tensor adhering to the constraints of a purely dissipative material, i.e. $\boldsymbol {D}:\boldsymbol {\mathcal {\eta }}\{\mathcal {H}\}:\boldsymbol {D}>0$, and is dependent upon the history $\mathcal {H}$ of flow, which can be conveniently represented through a deformation gradient $\boldsymbol {F}$ relative to a reference state. For the specific case of $|\boldsymbol {D}|\to 0$ corresponding to rate-independent plastic deformation (here, $|\boldsymbol {D}|=\sqrt {\frac {1}{2}\boldsymbol {D}:\boldsymbol {D}}$), the stress is given as

(2.2)\begin{equation} \boldsymbol{\sigma} = \frac{\boldsymbol{\mu}_{0}\{\mathcal{H}\}:\boldsymbol{D}}{|\boldsymbol{D}|}, \end{equation}

where $\boldsymbol {\mu }_{0}\{\mathcal {H}\}$ is a fourth-rank yield modulus. Therefore, the total stress can be partitioned into its rate-independent and rate-dependent components as

(2.3)\begin{equation} \boldsymbol{\sigma} = \frac{\boldsymbol{\mu}_{0}\{\mathcal{H}\}:\boldsymbol{D}}{|\boldsymbol{D}|} + \boldsymbol{\mathcal{\eta}}_{0}\{\mathcal{H}\}:\boldsymbol{D}, \end{equation}

where $\boldsymbol {\mathcal {\eta }}_{0}\{\mathcal {H}\}$ is a fourth-rank viscosity tensor.

We adapt this rheological framework to model granular rheology through the following assumptions that will be demonstrated to hold true in the present simulations: (1) the flow is homogeneous with a constant stretch history (Noll Reference Noll1962); and (2) the flow is planar and isochoric, i.e. $\boldsymbol {D}$ is characterized by two dominant eigenvalues and $\mathrm {tr}(\boldsymbol {D})=0$. This dependence is introduced in a frame-indifferent manner to produce a second-order rheological model that well-describes non-isotropic flow effects observed in our simulations. With these simplifications, $\boldsymbol {\sigma }$ can be represented as

(2.4)\begin{align} \boldsymbol{\sigma} &= p\boldsymbol{I}+\eta_{1}\boldsymbol{D}+\eta_{2} \left[\boldsymbol{D}^{2}-\frac{\mathrm{tr}\left(\boldsymbol{D}^{2}\right)}{3}\boldsymbol{I}\right]+ \eta_{3}\left[\dot{\boldsymbol{D}} - \boldsymbol{W}\boldsymbol{D}+\boldsymbol{D}\boldsymbol{W}\right] \nonumber\\ &\quad +\kappa_{1}\frac{\boldsymbol{D}}{|\boldsymbol{D}|}+\kappa_{2} \left[\frac{\boldsymbol{D}^{2}}{|\boldsymbol{D}|^{2}}- \frac{\mathrm{tr}\left(\boldsymbol{D}^{2}\right)}{3|\boldsymbol{D}|^{2}}\boldsymbol{I}\right], \end{align}

where

(2.5)\begin{equation} \dot{\boldsymbol{D}}=\frac{\partial \boldsymbol{D}}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{D} \end{equation}

is the material derivative of $\boldsymbol {D}$, $\boldsymbol {v}$ is the material velocity, and $\mathring {\boldsymbol {D}}=\dot {\boldsymbol {D}}-\boldsymbol {W}\boldsymbol {D}+\boldsymbol {D}\boldsymbol {W}$ represents the frame-indifferent corotational derivative of $\boldsymbol {D}$ (Bird & Hassager Reference Bird and Hassager1987). The isotropic pressure is defined as $p=\frac {1}{3}\mathrm {tr}(\boldsymbol {\sigma })$, and $\boldsymbol {I}$ is the unit tensor. The deviatoric stress, $\boldsymbol {\sigma } - p\boldsymbol {I}$, depends on $\boldsymbol {D}=1/2(\boldsymbol {\nabla } \boldsymbol {v} + \boldsymbol {\nabla } \boldsymbol {v}^{\textrm {T}})$ and a vorticity tensor $\boldsymbol {W}=1/2(\boldsymbol {\nabla } \boldsymbol {v} - \boldsymbol {\nabla }\boldsymbol {v}^{\textrm {T}})$. Here, $\dot {\gamma }=|\boldsymbol {D}|$ is the magnitude of the strain-rate tensor. The second, third and fourth terms in (2.4) represent rate-dependent contributions to the total stress that are characterized by the flow functions $\eta _{1}(\dot {\gamma },p)$ and $\eta _{2}(\dot {\gamma },p)$ and $\eta _{3}(\dot {\gamma },p)$, and are similar in form to a second-order description of non-Newtonian fluids using Rivlin–Erickson tensors (Rivlin Reference Rivlin1955). The fifth and sixth terms in (2.4) represent rate-independent contributions to the total stress that are characterized by plastic yield-like functions $\kappa _{1}(p)$ and $\kappa _{2}(p)$, which generally depend on the flow history. The pressure dependence of the flow functions is similar in spirit to the implicit constitutive theory of Rajagopal (Reference Rajagopal2006). In this work we focus on simple shear flows, but in general, the coefficients $\eta _{1}$, $\eta _{2}$ and $\eta _{3}$ depend on $\mathrm {tr}(\boldsymbol {D}^2)$ and $\mathrm {tr}(\mathring {\boldsymbol {D}}^2)$, which is important when modelling non-viscometric flows (Giusteri & Seto Reference Giusteri and Seto2018). Similarly, the coefficients $\kappa _{1}$ and $\kappa _{2}$ depend on $\mathrm {tr}(\boldsymbol {D}^2)/|\boldsymbol {D}|^{2}$, which can be calibrated from anisotropic models of granular yield criterion. Such anisotropy was demonstrated in simulations (Thornton & Zhang Reference Thornton and Zhang2010; Li & Dafalias Reference Li and Dafalias2012), resulting in deviations from the Drucker–Prager-like isotropic yield criterion that is implicit in the $\mu (I)$ rheology (Jop et al. Reference Jop, Forterre and Pouliquen2006). Furthermore, the rheological model can be extended to multiaxial flows that are observed in practice (Cortet et al. Reference Cortet, Bonamy, Daviaud, Dauchot, Dubrulle and Renouf2009) by introducing additional dependence of flow functions on $\mathrm {tr}(\boldsymbol {D}^3)$ and $\mathrm {tr}(\mathring {\boldsymbol {D}}^3)$ (Wang Reference Wang1965; Larson Reference Larson1985).

In this paper, we will consider steady homogeneous planar shear flow of granular materials resulting from a constant applied external shear stress and pressure, in which the memory of the flow has decayed and the deformation history is unimportant. In such steady homogeneous flows $\dot {\boldsymbol {D}}=0$, indicating that the eigenvectors of $\boldsymbol {D}$ are uniform in space and time, and local material rotation arises entirely from flow vorticity (Schunk & Scriven Reference Schunk and Scriven1990; Giusteri & Seto Reference Giusteri and Seto2018). Consider a uniform velocity gradient $\boldsymbol {\nabla } \boldsymbol {v}$ with the following viscometric form:

(2.6)\begin{equation} \boldsymbol{\nabla} \boldsymbol{v} = \left[\begin{array}{@{}ccc@{}} 0 & 2\dot{\gamma} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array}\right],\end{equation}

for flow along $x$ direction, velocity gradient along $y$ direction and vorticity along $z$ direction, and where $\mathrm {tr}(\boldsymbol {D}^3)=0$. In such viscometric flows, $\eta _{1}$, $\eta _{2}$ and $\eta _{3}$ represent the standard viscometric flow functions for non-Newtonian fluids (Coleman, Markovitz & Noll Reference Coleman, Markovitz and Noll1966) corresponding to shear stress, second normal stress difference and first normal difference, respectively. Similarly, $\kappa _{1}$ and $\kappa _{2}$ represent the analogous rate-independent flow functions. Correspondingly, for such viscometric flows, the stress tensor takes the following general form:

(2.7)\begin{equation} \boldsymbol{\sigma} = \left[\begin{array}{@{}ccc@{}} \sigma_{xx} & \sigma_{xy} & 0 \\ \sigma_{xy} & \sigma_{yy} & 0 \\ 0 & 0 & \sigma_{zz} \end{array}\right],\end{equation}

where $\sigma _{xx}\neq \sigma _{yy}\neq \sigma _{zz}$. Previous simulations on sheared granular flows have proposed a similar form for the stress tensor, such as in granular flows down an incline (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013), free surface flows (McElwaine et al. Reference McElwaine, Takagi and Huppert2012) and in the shear-free sheets model (Depken, Van Saarloos & Van Hecke Reference Depken, Van Saarloos and Van Hecke2006) that proposed $\sigma _{xx}=\sigma _{yy}\neq \sigma _{zz}$ for quasi-static granular flows in a split-bottom Couette cell (Depken et al. Reference Depken, Lechman, Van Hecke, Van Saarloos and Grest2007).

The rheological model reduces to the well known $\mu (I)$ relationship (Jop et al. Reference Jop, Forterre and Pouliquen2006) for sheared granular flows when the second-order coefficients $\eta _{2,3}=0$ and $\kappa _{2}=0$. In this case, the stress tensor is assumed to be coaxial with the strain rate tensor, and the two are related to each other by a scalar relationship,

(2.8)\begin{equation} \boldsymbol{\sigma}=p\boldsymbol{I}+\mu(I)p\frac{\boldsymbol{D}}{|\boldsymbol{D}|}, \end{equation}

where the stress ratio $\mu =|\boldsymbol {\sigma }-p\boldsymbol {I}|/p$ and the inertial number $I=|\boldsymbol {D}|a/(p/\rho )^{0.5}$, for an average particle diameter $a$ and material density $\rho$. The $\mu (I)$ function is related to the flow coefficients of the rheological model in (2.4) through

(2.9)\begin{equation} \mu(I) = \frac{1}{p}\left(\eta_{1}|\boldsymbol{D}|+\kappa_{1}\right). \end{equation}

2.2. Constant stress simulations

Steady sheared flows can be simulated by applying a constant strain rate or a constant stress on the granular material. Previous simulations on granular flows have imposed a constant strain rate either through a solid wall-driven flow (Da Cruz et al. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005; Koval et al. Reference Koval, Roux, Corfdir and Chevoir2009; Kamrin & Koval Reference Kamrin and Koval2014; Salerno et al. Reference Salerno, Bolintineanu, Grest, Lechman, Plimpton, Srivastava and Silbert2018) or by shearing the periodic simulation domain (Campbell Reference Campbell2002, Reference Campbell2005; Otsuki & Hayakawa Reference Otsuki and Hayakawa2011; Sun & Sundaresan Reference Sun and Sundaresan2011; Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). Wall-driven granular flows often result in flow localization near the walls (Shojaaee et al. Reference Shojaaee, Roux, Chevoir and Wolf2012b), which requires careful calibration of wall properties to extract the bulk rheological properties (Shojaaee et al. Reference Shojaaee, Brendel, Török and Wolf2012a; Schuhmacher, Radjai & Roux Reference Schuhmacher, Radjai and Roux2017). While a constant strain rate at the periodic boundaries can produce a viscometric flow field without walls (Campbell Reference Campbell2002, Reference Campbell2005; Peyneau & Roux Reference Peyneau and Roux2008; Sun & Sundaresan Reference Sun and Sundaresan2011), it often results in large shear stress fluctuations (Peyneau & Roux Reference Peyneau and Roux2008), especially in the quasi-static flow regime, which makes it challenging to calibrate the rate-independent part of granular rheology. Additionally, it was recently demonstrated that near the critical yield stress, granular flows are highly intermittent with a stochastic flow-arrest transition behaviour (Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). As such, a constant stress boundary condition is able to provide an accurate prediction of the rheology near the yield stress (Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). In this work, we simulate granular flows by applying a constant shear stress at the periodic boundaries, in which material is allowed to flow or not depending on the magnitude of applied stress. We will demonstrate that this boundary condition results in a well-defined viscometric flow.

Granular flows can also be simulated either at constant volume (isochoric) (Campbell Reference Campbell2002; Otsuki & Hayakawa Reference Otsuki and Hayakawa2011; Sun & Sundaresan Reference Sun and Sundaresan2011) or by imposing a constant normal stress (Campbell Reference Campbell2005; Sun & Sundaresan Reference Sun and Sundaresan2011; de Coulomb et al. Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017; Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). Granular materials dilate upon shearing, resulting in significant differences in the rheology between the two conditions (Campbell Reference Campbell2005). When the applied normal stress is constant, the material can dilate or compact upon shearing (depending on the initial condition) towards a ‘critical state’ solid volume fraction in the quasi-static regime (Schofield & Wroth Reference Schofield and Wroth1968; Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). Furthermore, granular materials exhibit shear-induced dilation in the inertial regime. Isochoric granular flows are not commonly observed in practice, and various experiments often naturally correspond to a constant normal stress condition, such as in free surface flows (Jop et al. Reference Jop, Forterre and Pouliquen2006; McElwaine et al. Reference McElwaine, Takagi and Huppert2012) or flows in Couette cells (Lu, Brodsky & Kavehpour Reference Lu, Brodsky and Kavehpour2007; Dijksman et al. Reference Dijksman, Wortel, Van Dellen, Dauchot and Van Hecke2011). Furthermore, a constant volume condition precludes the possibility of simulating granular flows near the yield stress in the quasi-static regime. If the solid volume fraction is set lower than the critical solid volume fraction at the onset of flow, then a $\mu (I)$ frictional rheology cannot be extracted as the shear stress will go to zero (rather than its yield threshold value) as the strain rate goes to zero. Similarly, if the volume fraction is set greater than its critical value, the flow is prohibited for any applied stress in the limit of rigid grains. In this work, we simulate granular flows at a constant applied pressure where the material is free to adapt its volume. A constant pressure condition is different from the case where all the normal stress components are specified equal to each other, as simulated previously in Peyneau & Roux (Reference Peyneau and Roux2008). This allows an efficient estimation of normal stress differences that will be described later in the text.

To simulate the evolution of a granular system under constant external stress and pressure, we utilize a modularly invariant adaptation (Shinoda, Shiga & Mikami Reference Shinoda, Shiga and Mikami2004) of the Parrinello–Rahman method (Parrinello & Rahman Reference Parrinello and Rahman1981) for molecular dynamics. This method was originally introduced to simulate the bulk properties of molecular systems in an isoenthalpic–isotension ensemble, including any phase transitions induced by the applied external stress (Parrinello & Rahman Reference Parrinello and Rahman1981). Such stress-controlled simulation methods adapted from molecular dynamics were previously implemented to study jamming (Smith et al. Reference Smith, Srivastava, Fisher and Alam2014) and creep (Srivastava & Fisher Reference Srivastava and Fisher2017) in granular packings, and recently to analyse flow-arrest transition in granular flows (Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). However, this is the first study that utilizes these methods to simulate steady frictional granular flows under external shear stress in order to extract their constitutive rheological behaviour. Simplified versions of such methods were also previously reported in simulations of non-equilibrium simple shear flows of Lennard-Jones fluids at a constant pressure and temperature to estimate their viscosity (Evans & Ely Reference Evans and Ely1986; Hood, Evans & Morriss Reference Hood, Evans and Morriss1987), and recently for simulating the rheology of colloidal suspensions (Wang & Brady Reference Wang and Brady2015). The simulation framework described here can robustly simulate more complex flows beyond simple shear.

In the present simulations, a collection of particles contained within a three-dimensional triclinic periodic cell is allowed to evolve under the application of a constant external stress tensor $\boldsymbol {\sigma }_{{ext}}$, which is constrained by (i) $(1/3)\sum\sigma _{{ext},ii}=p_{{ext}}$; (ii) $\sigma _{{ext},ij}=\tau _{{ext}}$ for $i,j=1,2$ and $2,1$; and (iii) $\sigma _{{ext},ij}=0$ for all other Einstein indices $i\neq j$, as shown in the schematic in figure 1. Because the traction at the boundaries of the periodic cell is prescribed, the periodic cell itself is allowed to dilate (or compact) and deform its shape in all possible ways, thus simulating the true bulk response of the granular material under external stress and pressure. The triclinic periodic cell is represented by a cell matrix ${\boldsymbol{\mathsf{H}}}$ which is a concatenation of the three lattice cell vectors that define the periodicity of the system. The cell matrix is constrained to be upper-triangular and the internal stress tensor is symmetrized to prevent any spurious cell rotations, which was a problem in the original Parrinello–Rahman method. This was achieved differently using a positive-definite metric tensor in another variant of this method reported previously by Souza & Martins (Reference Souza and Martins1997). Upon the application of $\boldsymbol {\sigma }_{{ext}}$, the equations of motion for $N$ particle positions and momenta $\{\boldsymbol {r}_{i},\boldsymbol {p}_{i}\}$, and the periodic cell matrix and its associated momentum tensor $\{{\boldsymbol{\mathsf{H}}},{\boldsymbol{\mathsf{P}}}_{g}\}$ are given by

(2.10a)\begin{gather} \dot{\boldsymbol{r}}_{i}=\frac{\boldsymbol{p}_{i}}{m_{i}}+\frac{{\boldsymbol{\mathsf{P}}}_{g}}{W_{g}}\boldsymbol{r}_{i}, \end{gather}
(2.10b)\begin{gather}\dot{\boldsymbol{p}}_{i}=\boldsymbol{f}_{i}-\frac{{\boldsymbol{\mathsf{P}}}_{g}}{W_{g}}\boldsymbol{p}_{i}-\frac{1}{3N} \frac{\mathrm{Tr\left[{\boldsymbol{\mathsf{P}}}_{g}\right]}}{W_{g}}\boldsymbol{p}_{i}, \end{gather}
(2.10c)\begin{gather}\dot{\boldsymbol{\mathsf{H}}}=\frac{{\boldsymbol{\mathsf{P}}}_{g}}{W_{g}}{\boldsymbol{\mathsf{H}}}, \end{gather}
(2.10d)\begin{gather}\dot{{\boldsymbol{\mathsf{P}}}_{g}}=V\left(\boldsymbol{\sigma}_{{int}}- \boldsymbol{I}p_{{ext}}\right)-{\boldsymbol{\mathsf{H}}}\boldsymbol{\varSigma}{\boldsymbol{\mathsf{H}}}^{\textrm{T}}, \end{gather}

where $\boldsymbol {f}_{i}$ is the net force on a particle $i$, $V$ is the variable volume of the periodic cell, $\boldsymbol {I}$ is the identity tensor and $W_{g}$ is a ‘fictitious’ mass associated with the inertia of the periodic cell. The stress quantities $\boldsymbol {\sigma }_{{int}}$ and $\boldsymbol {\varSigma }$ are defined below.

Figure 1. Schematic of the simulation method: from left to right, the three images represent the configurations of a granular system at three consecutive simulation times during steady flow, while subjected to an external pressure $p_{{ext}}$ and shear stress $\tau _{{ext}}$. The triclinic periodic cell boundaries (in black) at three times are, respectively, represented by matrices $\boldsymbol {H}_{0}$, $\boldsymbol {H}_{1}$ and $\boldsymbol {H}_{2}$. The triclinic periodic cell volume is almost equal at all three times in steady flow. The dotted lines in the global coordinate system represent directions into the plane.

In the original Parrinello–Rahman method for molecular systems, the fictitious mass is suggested to be set as $W_{g}=Nk_{B}T/\omega _{g}^{2}$ for an efficient sampling of the isoenthalpic–isotension ensemble (Martyna et al. Reference Martyna, Tuckerman, Tobias and Klein1996). Here, $k_{B}$ is the Boltzmann constant, $T$ is intended temperature of the ensemble and $\omega _{g}$ is the characteristic phonon frequency of the system. Such suggestions do not apply to the athermal flow simulations considered here. Analogously, the fictitious mass in the present case can be set as $W_{g}=Nk_{n}a^{2}/\omega _{g}^{2}$, where $k_{n}$ is the elastic constant associated with particle contacts (see § 2.4), $a$ is the mean particle diameter and $k_{n}a^{2}$ set the energy scale of system. The choice of $\omega _{g}$ controls the magnitude of stress fluctuations during steady granular flow, but it does not affect the rheology of flow within some upper and lower bounds of $\omega _{g}$, as was established by testing various values of $\omega _{g}$. A convenient value is $\omega _{g}=2.2\sqrt {m/k_{n}}$, where $m$ is the mean particle mass. Smaller values of $\omega _{g}$ resulted in larger stress fluctuations, whereas larger values of $\omega _{g}$ took longer simulation times to achieve steady flow. Similar analyses of the effect on $\omega _{g}$ on stress-controlled simulations were previously presented for non-equilibrium flow of Lennard-Jones fluids (Evans & Ely Reference Evans and Ely1986; Hood et al. Reference Hood, Evans and Morriss1987). A comprehensive numerical analysis of the effect of $\omega _{g}$ on stress-controlled simulations of granular flows is a part of our ongoing work.

The first two terms of the right-hand side of (2.10d), respectively, represent the imbalance between bulk internal stress of the granular system $\boldsymbol {\sigma }_{{int}}$ and external applied stress, which drives the motion of the periodic cell. The components of the bulk internal stress $\boldsymbol {\sigma }_{{int}}$ are calculated as (Walton & Braun Reference Walton and Braun1986; Da Cruz et al. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005)

(2.11)\begin{equation} \sigma_{\alpha\beta,{int}}=\frac{1}{V}\sum_{i}\left[\sum_{\,j \neq i} \frac{1}{2}r_{\alpha,ij} \,f_{\beta,ij} + m_{i}v^{'}_{\alpha,i}v^{'}_{\beta,i}\right], \end{equation}

where $\boldsymbol {r}_{ij}$ and $\boldsymbol {f}_{ij}$ are the branch vector and the force between two contacting particles $i$ and $j$. The fluctuating velocity $\boldsymbol {v}^{'}_{i}$ of particle $i$ is defined as the difference between velocity $\boldsymbol {v}_{i}$ of particle $i$ and mean shearing field velocity, such that $\boldsymbol {v}^{'}_{i}=\boldsymbol {v}_{i}-(\boldsymbol {\nabla } \boldsymbol {v})\boldsymbol {x}_{i}$, where $\boldsymbol {\nabla } \boldsymbol {v}$ is bulk velocity gradient, and $\boldsymbol {x}_{i}$ is the position of particle $i$. Hereafter, the subscript ‘${int}$’ will be dropped while referring to the internal state of the stress of the granular system. In (2.10d), the tensor $\boldsymbol {\varSigma }$ is defined as (Shinoda et al. Reference Shinoda, Shiga and Mikami2004)

(2.12)\begin{equation} \boldsymbol{\varSigma}=\boldsymbol{H}_{0}^{-1}\left(\boldsymbol{\sigma}_{{ext}}-\boldsymbol{I}p_{{ext}}\right) \boldsymbol{H}_{0}^{\textrm{T}-1}, \end{equation}

where $\boldsymbol {H}_{0}$ is some reference state of the periodic cell, and $J^{-1}{\boldsymbol{\mathsf{H}}}\boldsymbol {\varSigma }{\boldsymbol{\mathsf{H}}}^{\textrm {T}}$ represents the ‘true’ measure of the external deviatoric stress, which is defined with respect to the reference state (Souza & Martins Reference Souza and Martins1997). Here $J=\mathrm {det}[\boldsymbol {F}]$ is the Jacobian of the deformation gradient $\boldsymbol {F}$, which is defined as

(2.13)\begin{equation} \boldsymbol{F}=\boldsymbol{H}\boldsymbol{H}_{0}^{-1}. \end{equation}

It is evident from (2.10d) that a difference between the internal stress and external applied stress drives the perpetual motion of the periodic cell during steady flow. In the case where the internal and external stress balance each other, the motion of the cell eventually stops because the external stress is not sufficient to continually drive the motion of the cell, thus enabling the precise identification of the yield stress of the granular system (Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). As a result, this implementation of a constant external stress on the granular system prescribes the second Piola–Kirchoff measure of the external stress, or equivalently the thermodynamic tension (Souza & Martins Reference Souza and Martins1997). In the present simulations, the reference state is updated to the current state at the end of every time step of integration of the equations of motion, in order to minimize the deviation of internal strain energy from work done by the external stress. All the simulations are performed using the large-scale molecular dynamics software LAMMPS (Plimpton Reference Plimpton1995).

2.3. Bulk deformation

Upon applying an external pressure $p_{{ext}}$ and shear stress $\tau _{{ext}}$ to a granular system, all the components of the macroscopic internal stress tensor $\boldsymbol {\sigma }$ evolve independently with time. Correspondingly, the triclinic periodic cell, represented by the matrix $\boldsymbol {H}$, also evolves with time from bulk volumetric and shear deformation. Figure 1 shows a schematic of the evolution of deformation of a triclinic periodic cell in steady flow as it is subjected to a constant external shear stress and pressure. The states of the triclinic periodic cell $\boldsymbol {H}$ are stored at every simulation time step (such as $\boldsymbol {H}_{0}$, $\boldsymbol {H}_{1}$ and $\boldsymbol {H}_{2}$ shown in figure 1), and are used to compute the bulk velocity gradient in the periodic system, as described below.

Consider the position $\boldsymbol {r}(t)$ of a particle at a simulation time $t$ within the periodic cell, defined with respect to an origin (typically, one of the corners of the periodic cell). Its reduced coordinates $\boldsymbol {s}(t)$ can be defined by

(2.14)\begin{equation} \boldsymbol{r}(t)=\boldsymbol{H}(t)\boldsymbol{s}(t), \end{equation}

such that $0<\boldsymbol {s}(t)<1$. The periodic tiling of the space by the triclinic cell $\boldsymbol {H}(t)$ implies that a spatial coordinate $\boldsymbol {r}^{'}(t)=\boldsymbol {H}(t)[\boldsymbol {s}(t)+\boldsymbol {\varDelta }]$ represents the periodic image of $\boldsymbol {r}$, where $\boldsymbol {\varDelta }$ is a vector of integers. The velocity $\boldsymbol {v}(t)=\dot {\boldsymbol {r}}(t)$ of the particle is defined such that

(2.15)\begin{equation} \boldsymbol{v}(t)=\dot{\boldsymbol{H}}(t)\boldsymbol{s}(t) + \boldsymbol{H}(t)\dot{\boldsymbol{s}}(t), \end{equation}

where the first term represents the contribution from the bulk periodic cell deformation and the second term represents the fluctuating non-affine velocity. Consequently, a bulk velocity gradient can be defined as $\boldsymbol {\nabla }\boldsymbol {v}(t)=\nabla _{\boldsymbol {r}}(\dot {\boldsymbol {H}}(t)\boldsymbol {s}(t))$. Upon substituting (2.14) we get

(2.16)\begin{equation} \boldsymbol{\nabla} \boldsymbol{v}(t) = \dot{\boldsymbol{H}}(t)\boldsymbol{H}^{-1}(t). \end{equation}

2.4. Contact mechanics

In the present simulations, frictional spherical particles interact only upon contact. The contact forces are modelled using a spring and a dashpot along with a static yield criterion to model contact friction. This model was first developed by Cundall & Strack (Reference Cundall and Strack1979), and since has been tested and implemented in various granular flow simulations (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Campbell Reference Campbell2005; Rycroft et al. Reference Rycroft, Kamrin and Bazant2009; Sun & Sundaresan Reference Sun and Sundaresan2011). Two contacting particles $\{i,j\}$ of diameters $\{a_{i},a_{j}\}$, masses $\{m_{i},m_{j}\}$, at positions $\{\boldsymbol {r}_i,\boldsymbol {r}_j\}$ with velocities $\{\boldsymbol {v}_i,\boldsymbol {v}_j\}$ and angular velocities $\{\boldsymbol {\omega }_i,\boldsymbol {\omega }_j\}$ are considered to be in contact if $\delta _{ij}=\frac {1}{2}(a_{i}+a_{j})-|\boldsymbol {r}_{ij}|>0$, where $\boldsymbol {r}_{ij}=\boldsymbol {r}_{i}-\boldsymbol {r}_{i}$ is the vector connecting their centroids; these quantities are tracked at every time step as they evolve from particle collisions or affine particle motion caused by triclinic cell deformation, as described in (2.10a). The contact normal force $\boldsymbol {f}_{nij}$ and tangential force $\boldsymbol {f}_{tij}$ on particle $i$ are given by

(2.17a)\begin{gather} \boldsymbol{f}_{nij}=k_{n} \delta_{ij} \boldsymbol{n}_{ij} - \gamma_{n} m_{e} \boldsymbol{v}_{nij}, \end{gather}
(2.17b)\begin{gather}\boldsymbol{f}_{tij}=-k_{t} \boldsymbol{u}_{tij} - \gamma_{t} m_{e} \boldsymbol{v}_{tij}, \end{gather}

where $k_{n,t}$ and $\gamma _{n,t}$ are contact stiffness and damping constants, and $m_{e}=m_{i}m_{j}/(m_{i}+m_{j})$ is the effective mass. The corresponding force on particle $j$ is given by Newton's third law such that $\boldsymbol {f}_{ji}=\boldsymbol {f}_{ij}$. The unit normal along the axis of contact is given by $\boldsymbol {n}_{ij}=\boldsymbol {r}_{ij}/|\boldsymbol {r}_{ij}|$, and $\boldsymbol {v}_{nij}$ and $\boldsymbol {v}_{tij}$ are, respectively, the normal and tangential components of the relative velocity $\boldsymbol {v}_{ij}=\boldsymbol {v}_{i}-\boldsymbol {v}_{j}$ given by

(2.18a)\begin{gather} \boldsymbol{v}_{nij}=\left(\boldsymbol{v}_{ij}\boldsymbol{\cdot}\boldsymbol{n}_{ij}\right)\boldsymbol{n}_{ij}, \end{gather}
(2.18b)\begin{gather}\boldsymbol{v}_{tij}=\boldsymbol{v}_{ij}-\boldsymbol{v}_{nij}-\tfrac{1}{2}\left(\boldsymbol{\omega}_{i}+\boldsymbol{\omega}_{j}\right)\times\boldsymbol{r}_{ij}. \end{gather}

An elastic displacement $\boldsymbol {u}_{tij}$ representing shear in the tangential direction is tracked during the lifetime of a contact, and it evolves according to the following ordinary differential equation:

(2.19)\begin{equation} \frac{\mathrm{d} \boldsymbol{u}_{tij}}{\mathrm{d} t}=\boldsymbol{v}_{tij}- \frac{\left(\boldsymbol{u}_{tij}\boldsymbol{\cdot}\boldsymbol{v}_{ij}\right) \boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij}|^{2}}, \end{equation}

with $\boldsymbol {u}_{tij}=0$ at the initiation of the contact.

Tangential friction between two contacting particles is modelled by a static yield criterion $|\,\boldsymbol {f}_{tij}|<\mu _{s}|\,\boldsymbol {f}_{nij}|$, which is always satisfied by limiting the tangential shear displacement $\boldsymbol {u}_{tij}$. The particle coefficient of sliding friction $\mu _{s}$ is a measure of its surface roughness, and significantly impacts the rheology of granular flow. The normal and tangential viscous damping at a contact are controlled by the coefficients of restitution $e_{n,t}=\exp (-\gamma _{n,t}t_{c}/2)$, where $t_{c}={\rm \pi} (2k_{n}/m_{e}-\gamma _{n}^{2}/4)^{-1/2}$ is the collision time between two contacting particles (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001).

2.5. Simulation details

The contact stiffness between particles $k_{n}$ and $k_{t}$ is set equal to each other. The normal damping constant $\gamma _{n}=0.5/t_{c}$ and the tangential damping constant $\gamma _{t}=0.5\gamma _{n}$. Initially, dilute configurations of granular systems at a solid volume fraction $\phi =0.05$ are subjected to a constant external shear stress and hydrostatic pressure. We simulate granular flow at three external pressures $p_{{ext}}a/k_{n}=10^{-4},10^{-5},10^{-6}$, all in the limit of the rigid particle regime where the rheology is unaffected by the applied pressure and particle stiffness (Da Cruz et al. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005; de Coulomb et al. Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017). The external shear stress $\tau _{{ext}}$ is varied from $\tau _{{ext}}/p_{{ext}}=0.0$ to $\tau _{{ext}}/p_{{ext}}=1.2$ to simulate flows at various shear rates, and three different realizations are simulated for each shear rate. Each simulation consists of $N\sim 10^{4}$ frictional particles whose diameters are uniformly distributed between $0.9a$ and $1.1a$. Several particle coefficients of sliding friction ranging from $\mu _{s}=0.0$ to $\mu _{s}=0.3$ are analysed to study the effect of friction on stress-controlled granular rheology. Contact mechanics between two particles is resolved by setting the simulation time step to $0.02t_{c}$. In the results presented below, time is scaled by $t_{c}$, length is scaled by $a$, energy is scaled by $k_{n}a^{2}$ and stress is scaled by $k_{n}/a$.

3. Evolution towards viscometric flow

When the external shear stress $\tau _{{ext}}$ and pressure $p_{{ext}}$ are switched on at $t=0$, a dilute assembly of particles at an initial solid volume fraction $\phi =0.05$ responds with rapid volumetric compaction, as shown by the evolution of $\phi$ in figure 2(a) for a particular case of interparticle friction $\mu _{s}=0.3$, $p_{{ext}}=10^{-4}$ and $\tau _{{ext}}=5\times 10^{-5}$. To estimate the total deformation accumulated by the material beyond isotropic compaction, we calculate the deformation gradient $\boldsymbol {F}(t)=\boldsymbol {H}(t)\boldsymbol {H}_{0}^{-1}$ as defined in (2.13), where $\boldsymbol {H}_{0}$ is the periodic cell at $t=0$. The rapid volumetric compaction at early times is seen by an equivalent decrease in $F_{ii}$ in figure 2(b) for $i=x,y,z$. The shear component $F_{xy}$ exhibits a super-linear increase at early times as a result shear deformation at low solid volume fractions in the absence of any significant resistance to the applied shear. At long times, $F_{xy}$ increases linearly with time, while $F_{ii}$ is constant and the other two shear components are negligible, thus indicating the achievement of steady incompressible viscometric flow, i.e. $\boldsymbol {F}(t)=\boldsymbol {I}+t\boldsymbol {M}$ (Coleman et al. Reference Coleman, Markovitz and Noll1966), where $\boldsymbol {M}$ is a constant tensor, and which is a special case of motion with constant stretch history (Noll Reference Noll1962) where the deviatoric stress depends on the form of $\boldsymbol {M}$ (Coleman et al. Reference Coleman, Markovitz and Noll1966). Several important and well-studied flows such as Couette flow, Poiseuille flow, simple shearing flow and some specific cases of torsional flows can be categorized as viscometric flows (Coleman et al. Reference Coleman, Markovitz and Noll1966). Although our focus here is on steady flows, the simulation method and rheological analysis described above provide the capability to calibrate a general history-dependent rheological model defined in (2.4) for transient granular flows under constant or time-varying applied stresses.

Figure 2. Evolution with time $t$ of (a) solid volume fraction $\phi$, (b) components of the deformation gradient tensor $F_{ij}$, (c) $d_{i}/d_1$, where $d_{i}$ are the eigenvalues of $\boldsymbol {D}$, (d) vorticity parameter $\beta$ for a particular case of interparticle friction $\mu _{s}=0.3$, applied pressure $p_{{ext}}=10^{-4}$ and applied shear stress $\tau _{{ext}}=5\times 10^{-5}$.

Further insight into the nature of viscometric flow is given by the eigenvalue decomposition of the symmetric tensor $\boldsymbol {D}(t)$. We use the convention that the three orthonormal eigenvectors of $\boldsymbol {D}$: $\hat {\boldsymbol {d}_{1}}$, $\hat {\boldsymbol {d}_{2}}$ and $\hat {\boldsymbol {d}_{3}}$ are ordered in decreasing order of signed eigenvalues $d_1$, $d_2$ and $d_3$. Figure 2(c) shows the evolution of $d_2/d_1$ and $d_3/d_1$ as a function of time. At early times, the sum of eigenvalues is positive, which corresponds with rapid volumetric compaction as described above. The long-time steady-state flow is characterized by $d_3=-d_1$ and $d_2=0$, which is a signature of planar flow, where the flow plane is spanned by $\hat {\boldsymbol {d}_{1}}$ and $\hat {\boldsymbol {d}_{3}}$. To further ascertain the nature of planar flow, we calculate a vorticity parameter $\beta$ defined as

(3.1)\begin{equation} \beta = \frac{1}{\dot{\gamma}}\frac{\boldsymbol{W}:\boldsymbol{G}}{\boldsymbol{G}:\boldsymbol{G}}, \end{equation}

where $\boldsymbol {G}=\hat {\boldsymbol {d}_{3}}\hat {\boldsymbol {d}_{1}}-\hat {\boldsymbol {d}_{1}}\hat {\boldsymbol {d}_{3}}$. Figure 2(d) shows the evolution of $\beta$ with time. When the system transition into steady state flow at long times, $\beta =1$, indicating simple shear deformation in the flow plane, thus confirming the viscometric nature of flow. During the transient evolution at early times, $0<\beta <1$, indicating a complex flow behaviour that is a mixture of vorticity-free elongational flow ($\beta =0$) and simple shear flow ($\beta =1$) (Wagner & Mckinley Reference Wagner and Mckinley2016; Giusteri & Seto Reference Giusteri and Seto2018). However, the flow is homogeneous at all times within the periodic cell during steady state, with no spatial gradients of the local strain rate.

We emphasize that the steady homogeneous shear flow states achieved in the present simulations considerably simplify the rheological model in (2.4). Because the eigenvectors of $\boldsymbol {D}$ are uniform in space and time, the material derivative $\dot {\boldsymbol {D}}=0$, and the local material rotation is equivalent to flow vorticity. In this particular case of steady homogeneous flow with constant stretch history, the rheology is equally well-represented by the following form of (2.4) (Larson Reference Larson1985; Brunn & Asoud Reference Brunn and Asoud2003; Giusteri & Seto Reference Giusteri and Seto2018):

(3.2)\begin{align} \boldsymbol{\sigma} &= p\boldsymbol{I}+\eta_{1}\boldsymbol{D}+\eta_{2} \left[\boldsymbol{D}^{2}-\frac{\mathrm{tr}\left(\boldsymbol{D}^{2}\right)}{3}\boldsymbol{I}\right]+ \eta_{3}\left[\boldsymbol{D}\boldsymbol{W}-\boldsymbol{W}\boldsymbol{D}\right] \nonumber\\ &\quad +\kappa_{1}\frac{\boldsymbol{D}}{|\boldsymbol{D}|}+\kappa_{2} \left[\frac{\boldsymbol{D}^{2}}{|\boldsymbol{D}|^{2}}- \frac{\mathrm{tr}\left(\boldsymbol{D}^{2}\right)}{3|\boldsymbol{D}|^{2}}\boldsymbol{I}\right]. \end{align}

We emphasize that in unsteady or inhomogeneous flows where the material rate of rotation can differ from flow vorticity, several criteria for classifying local flow kinematics have been prescribed (Schunk & Scriven Reference Schunk and Scriven1990; Thompson & Mendes Reference Thompson and Mendes2005), and they can be incorporated in the current rheological model.

In steady state, the bulk rheological quantities fluctuate around their mean values, as seen in figure 2(ad). In order to achieve robust statistics, every simulation is run for at least $10^{7}$ time steps to guarantee the achievement of steady state flow. This is especially necessary near the critical yield stress, where steady state equilibration takes a long time. Upon achieving steady state, each simulation is continued to run for at least another $10^{6}$ time steps, during which the all data of interest are recorded at every $10$ time steps and averaged to estimate their steady mean value. The statistical uncertainty associated with mean estimation is measured by its standard error using a block averaging method (Flyvbjerg & Petersen Reference Flyvbjerg and Petersen1989). This method not only provides robust estimates of uncertainty around a mean value, but also indicates if the data has any long-time correlations, which would necessitate longer simulation runs for meaningful averaging.

4. Model calibration

In this section, friction-dependent functional forms of all flow coefficients in (3.2) will be described. The material constants associated with these flow coefficients are extracted from the DEM simulation data by utilizing the fact that the four tensors $\boldsymbol {I}$, $\boldsymbol {D}$, $(\boldsymbol {D}^{2}-({\mathrm {tr}(\boldsymbol {D}^{2})}/{3})\boldsymbol {I})$ and $(\boldsymbol {D}\boldsymbol {W}-\boldsymbol {W}\boldsymbol {D})$ are orthogonal to each other in viscometric flows.

4.1. Flow functions: $\eta _{1}$ and $\kappa _{1}$

The two flow coefficients $\eta _{1}$ and $\kappa _{1}$ have a first-order contribution (in terms of $\boldsymbol {D}$) to the total stress $\boldsymbol {\sigma }$, and they provide a measure of the shear stress in viscometric flow. These coefficients are estimated by

(4.1)\begin{equation} \eta_{1}\dot{\gamma}+\kappa_{1} = \frac{1}{2\dot{\gamma}}\boldsymbol{\sigma}:\boldsymbol{D}, \end{equation}

where $\tau =({1}/{2\dot {\gamma }})\boldsymbol {\sigma }:\boldsymbol {D}$ is the total flow-induced shear stress in the system.

Previous research has shown that shear flow of granular materials can be well-described by a local rheological relationship between a stress ratio $\mu$ and an inertial number $I$ (Jop et al. Reference Jop, Forterre and Pouliquen2006). In the present model, the stress ratio (hereby written with a subscript $1$) is $\mu _{1}=(\eta _{1}\dot {\gamma }+\kappa _{1})/p$, where $\eta _{1}\dot {\gamma }/p$ is the rate-dependent contribution and $\kappa _{1}/p$ is the rate-independent contribution. As such, $\eta _{1}$ represents the effective shear viscosity and $\kappa _{1}/p$ represents the yield coefficient as $I\to 0$, which is associated with a critical volume fraction discussed in § 4.4. Figure 3(a) shows the variation of $\mu _{1}$ with $I$ for five interparticle frictions $\mu _{s}$ at three $p_{{ext}}$. All the curves at various pressures collapse onto a master curve for each $\mu _{s}$, which can be approximated by a power law for dense granular flows described in several previous studies (DeGiuli et al. Reference DeGiuli, Düring, Lerner and Wyart2015; DeGiuli, McElwaine & Wyart Reference DeGiuli, McElwaine and Wyart2016; de Coulomb et al. Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017; Salerno et al. Reference Salerno, Bolintineanu, Grest, Lechman, Plimpton, Srivastava and Silbert2018),

(4.2)\begin{equation} \mu_{1}=\mu_{1}^{0}+A_{1}I^{\alpha_1}, \end{equation}

where, $\mu _{1}^{0}$, $A_1$ and $\alpha _1$ are fitting parameters. In the quasi-static, rate-independent regime where $I\to 0$, the stress ratio reaches a constant value $\mu _{1}\to \mu _{1}^{0}$, which is equivalent to $\kappa _{1}/p$ in the rheological model. In this regime shear stress saturates towards a threshold value, while the pressure is well-controlled at its prescribed value, thus indicating the approach towards a yield stress. Although the rheology is unaffected by the applied pressure, as also observed in de Coulomb et al. (Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017), lower values of inertial numbers are achieved when the confining pressure is low, as seen by the variation of $\mu _{1}-\mu _{1}^{0}$ with $I$ for three pressures and two $\mu _{s}$ in figure 3(b). Previous pressure and shear rate controlled simulations had demonstrated that the transition from quasi-static to inertial flow regimes occurs at lower inertial number for lower confining pressure (de Coulomb et al. Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017), thus confirming the current observations. However, the present simulations produce highly stochastic flows in the quasi-static regime, which often arrest in the vicinity of the static yield coefficient (Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019). Therefore, the rheology at low inertial numbers is not well-resolved for low pressures, especially for intermediate interparticle friction, as seen in figure 3(a). Recent experiments (Perrin et al. Reference Perrin, Clavaud, Wyart, Metzger and Forterre2019) and simulations (Degiuli & Wyart Reference Degiuli and Wyart2017) have indicated that the local rheology of frictional granular materials possibly exhibits hysteresis at very low inertial numbers, which would also prohibit very slow flows in the present stress-controlled simulations.

Figure 3. (a) Stress ratio $\mu _{1}$ as a function of inertial number $I$ for five interparticle frictions $\mu _{s}$ (see legend) at three applied pressures: $p_{{ext}}=10^{-4},10^{-5},10^{-6}$. The vertical and horizontal error bars represent the standard error in the calculation of $\mu _1$ and $I$, respectively. The black dashed lines represent fits for each $\mu _{s}$ given in (4.2). (b) Variation of $\mu _{1}-\mu _{1}^{0}$ with $I$ at three applied pressures (see legend) for particles with $\mu _{s}=0.0$ (red) and $\mu _{s}=0.3$ (black). The dotted lines represent power-law fits from (4.2). (c) Variation of $\mu _{1}^{0}$ and (d) $\alpha _{1}$ with $\mu _{s}$. The open symbols in panel (c) and panel (d) indicate the values for zero friction.

The quasi-static stress ratio increases with friction from $\mu _{1}^{0}=0.09$ for frictionless particles to $\mu _{1}^{0}=0.33$ for particles with high friction, as shown in figure 3(c). The non-zero value of $\mu _{1}^{0}$ for frictionless particles is consistent with previous simulations (Peyneau & Roux Reference Peyneau and Roux2008) and experiments (Clavaud et al. Reference Clavaud, Bérut, Metzger and Forterre2017; Perrin et al. Reference Perrin, Clavaud, Wyart, Metzger and Forterre2019) that demonstrated a non-zero internal friction angle for frictionless granular material. The value of $\mu _{1}^{0}$ at high friction is consistent with previous simulations and experiments (Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011a; Salerno et al. Reference Salerno, Bolintineanu, Grest, Lechman, Plimpton, Srivastava and Silbert2018; Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019), and is also similar to the critical stress ratio from the critical state theory (Schofield & Wroth Reference Schofield and Wroth1968). The power-law exponent varies monotonically between $\alpha _1=0.37$ for frictionless particles to $\alpha _1=0.7$ for particles with high friction, as seen in figure 3(d). Although the exponent for frictionless particles corresponds well with prior theoretical predictions (DeGiuli et al. Reference DeGiuli, Düring, Lerner and Wyart2015, Reference DeGiuli, McElwaine and Wyart2016), the exponent at high friction is smaller than theoretical predictions of $\alpha _1=1.0$ (DeGiuli et al. Reference DeGiuli, Düring, Lerner and Wyart2015, Reference DeGiuli, McElwaine and Wyart2016). This could be attributed to a lack of data at low inertial numbers and the associated sensitivity of power-law fitting.

4.2. Flow functions: $\eta _{2}$ and $\kappa _{2}$

In addition to the shear stress contribution to the total internal stress, there are non-negligible second-order contributions that are typically observed in the flow of non-Newtonian fluids. In a viscometric description of such fluids, these effects are characterized by normal stress difference functions (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). In the present rheological model, $\eta _{2}$ and $\kappa _{2}$ represent one set of such rate-dependent and rate-independent contributions. These coefficients are estimated by

(4.3)\begin{equation} \eta_{2}\dot{\gamma}^{2}+\kappa_{2} = \frac{3}{2\dot{\gamma}^{2}}\boldsymbol{\sigma}: \left(\boldsymbol{D}^{2}-\frac{\mathrm{tr}\left(\boldsymbol{D}^{2}\right)}{3}\boldsymbol{I}\right), \end{equation}

and they represent the difference between mean normal stress in the flow plane and normal stress in the vorticity direction. A second stress ratio similar to $\mu _{1}$ is defined as $\mu _{2}=(\eta _{2}\dot {\gamma }^{2}+\kappa _{2})/p$, where $\eta _{2}\dot {\gamma }^{2}/p$ is the rate-dependent contribution and $\kappa _{2}/p$ is the rate-independent contribution. As such, $\eta _{2}$ represents a normal viscosity and $\kappa _{2}/p$ represents the threshold value as $I\to 0$. Figure 4(a) shows the variation of $\mu _{2}$ with the square of inertial number $I$ for five $\mu _{s}$ and three $p_{{ext}}$. All the curves at various pressures collapse onto a master curve for each $\mu _{s}$, which can be approximated by a power law,

(4.4)\begin{equation} \mu_{2}=\mu_{2}^{0}+A_{2}\left(I^{2}\right)^{\alpha_2}, \end{equation}

where, $\mu _{2}^{0}$, $A_2$ and $\alpha _2$ are fitting parameters. In a manner similar to $\mu _{1}$, the quasi-static values of $\mu _{2}$ at low inertial numbers are accessed for low confining pressures, as shown by the variation of $\mu _{2}-\mu _{2}^{0}$ with $I^{2}$ in figure 4(b) for two different $\mu _{s}$ that appear to collapse onto a single curve. However, the data at low inertial numbers is also noisy, resulting from the stochastic nature of slow granular flows, and increased noise in the measured data.

Figure 4. (a) Second stress ratio $\mu _{2}$ as a function of inertial number $I^2$ for five interparticle frictions $\mu _{s}$ (see legend) at three applied pressures: $p_{{ext}}=10^{-4},10^{-5},10^{-6}$. The vertical and horizontal error bars represent the standard error in the calculation of $\mu _2$ and $I^2$, respectively. The black dashed lines represent fits for each $\mu _{s}$ given in (4.4). (b) Variation of $\mu _{2}-\mu _{2}^{0}$ with $I^2$ at three applied pressures (see legend) for particles with $\mu _{s}=0.0$ (red) and $\mu _{s}=0.3$ (black). The dotted lines represent power-law fits from (4.4). (c) Variation of $\mu _{2}^{0}$ and (d) $\alpha _{2}$ with $\mu _{s}$. The open symbols in panel (c) and panel (d) indicate the values for zero friction.

In the quasi-static regime, $\mu _{2}$ tends towards a constant value $\mu _{2}\to \mu _{2}^{0}$, which is equivalent to $\kappa _{2}/p$ in the rheological model. Its value varies monotonically from $\mu _{2}^{0}=0.01$ for frictionless particles to $\mu _{2}^{0}=0.1$ for particles with high friction, as shown in figure 4(c). The non-zero value of $\mu _{2}^{0}$ for particles with high friction indicates that normal stress effects are present even in the quasi-static regime of flow, thus indicating a mild anisotropic nature of the yield surface that is commonly assumed to be isotropic (in the Drucker–Prager sense) within the $\mu (I)$ rheology (Jop et al. Reference Jop, Forterre and Pouliquen2006), but has been shown to be anisotropic in recent simulations (Thornton & Zhang Reference Thornton and Zhang2010; Li & Dafalias Reference Li and Dafalias2012). The power-law exponent varies monotonically between $\alpha _2=0.28$ for frictionless particles to $\alpha _2=0.44$ for particles with high friction, as shown in figure 4(d).

4.3. Flow function: $\eta _{3}$

An additional second-order contribution to the total stress emerges through the rate-dependent flow coefficient $\eta _{3}$, which is estimated by

(4.5)\begin{equation} \eta_{3}\dot{\gamma}^{2} = \frac{1}{8\dot{\gamma}^{2}}\boldsymbol{\sigma}: \left(\boldsymbol{D}\boldsymbol{W}-\boldsymbol{W}\boldsymbol{D}\right), \end{equation}

and for viscometric flows it represents the difference between the two normal stresses in the flow plane. A third stress ratio $\mu _{3}$ is defined as $\mu _{3}=\eta _{3}\dot {\gamma }^{2}/p$, where $\eta _{3}$ represents an additional normal viscosity. Figure 5(a) shows the variation of $\mu _{3}$ as a decreasing function of $I^{2}$ for five $\mu _{s}$ at three $p_{{ext}}$. All the curves at various pressures collapse onto a master curve for each $\mu _{s}$, which can be approximated by the following power law:

(4.6)\begin{equation} \mu_{3}=-A_{3}\left(I^{2}\right)^{\alpha_3}, \end{equation}

where $A_{3}$ and $\alpha _3$ are fitting parameters.

Figure 5. (a) Third stress ratio $\mu _{3}$ as a function of inertial number $I^2$ for five interparticle frictions $\mu _{s}$ (see legend) at three applied pressures: $p_{{ext}}=10^{-4},10^{-5},10^{-6}$. The vertical and horizontal error bars represent the standard error in the calculation of $\mu _3$ and $I^2$, respectively. The black dashed lines represent fits for each $\mu _{s}$ given in (4.6). (b) Variation of $-\mu _{3}$ with $I^2$ at three applied pressures (see legend) for particles with $\mu _{s}=0.0$ (red) and $\mu _{s}=0.3$ (black). The dotted lines represent power-law fits from (4.6). (c) Variation of $\alpha _{3}$ with $\mu _{s}$. The open symbol in panel (c) indicates the value for zero friction.

The stress ratio $\mu _{3}$ exhibits a transition from negative values at high inertial numbers to small positive values in the quasi-static regime for all $\mu _{s}$, as seen in figure 5(a). Although the small positive value of $\mu _{3}$ in the quasi-static regime is intriguing, its existence is debated (see § 5.2 below) and this effect is not included in our rheological model. As such, a simple power law in (4.6) well-predicts the variation of $\mu _{3}$ with $I$, as also seen by the variation of $-\mu _{3}$ with $I^{2}$ in figure 5(b). The power-law exponent $\alpha _3$ varies slightly between $0.85$ and $0.75$ from low to high friction, as shown in figure 5(c).

For steady homogeneous simple shear flows simulated in this work, the constitutive model for viscometric flows in (3.2) reduces to the following relationships between the components of symmetric stress tensor $\boldsymbol {\sigma }$ and the three stress ratios, in the case of a shear flow along the $x$ direction and flow gradient along the $y$ direction:

(4.7a)\begin{gather} \sigma_{xx}=p\left(1+\frac{\mu_{2}}{3}-2\mu_{3}\right), \end{gather}
(4.7b)\begin{gather}\sigma_{yy}=p\left(1+\frac{\mu_{2}}{3}+2\mu_{3}\right), \end{gather}
(4.7c)\begin{gather}\sigma_{zz}=p\left(1-\frac{2\mu_{2}}{3}\right), \end{gather}
(4.7d)\begin{gather}\sigma_{xy}=p\mu_{1}, \end{gather}
(4.7e)\begin{gather}\sigma_{yz}=\sigma_{xz}=0. \end{gather}

4.4. Granular flow-induced dilation

The solid volume fraction $\phi$ of granular materials is highly sensitive to pressure and the rate of shear flow. These materials compact (jam) under the action of external pressure. However, under the action of external shear stress they dilate in order to flow, and the extent of dilation is higher for faster flows. In the present simulations, $\phi$ is not prescribed, and the system is allowed to freely attain its steady state solid volume fraction in response to the external stress and pressure. As such, we extract a dilatancy law relating the steady-state $\phi$ with the inertial number $I$ of the flow. Figure 6(a) shows the variation of $\phi$ with $I$ for five $\mu _{s}$ at three $p_{{ext}}$. All the curves at various pressures collapse onto a master curve for each $\mu _{s}$, which can be approximated by a power law as described in several previous studies (DeGiuli et al. Reference DeGiuli, Düring, Lerner and Wyart2015, Reference DeGiuli, McElwaine and Wyart2016; de Coulomb et al. Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017),

(4.8)\begin{equation} \phi=\phi^{0}-A_{4}I^{\alpha_4}, \end{equation}

where $\phi ^{0}$, $A_{4}$ and $\alpha _4$ are fitting parameters. The applied pressure moderately affects the volume fraction $\phi$, with lower $\phi$ at lower pressures, as seen in figure 6(b). It has been previously demonstrated that for sufficiently rigid particles (or equivalently, low enough applied pressures) in the hard particle limit, the effect of pressure is negligible on the volume fraction of granular material at onset of flow (de Coulomb et al. Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017).

Figure 6. (a) Solid volume fraction $\phi$ as a function of inertial number $I$ for five interparticle frictions $\mu _{s}$ (see legend) at three applied pressures: $p_{{ext}}=10^{-4},10^{-5},10^{-6}$. The vertical and horizontal error bars represent the standard error in the calculation of $\phi$ and $I$, respectively. The black dashed lines represent fits for each $\mu _{s}$ given in (4.8). The black crosses represent the data from Peyneau & Roux (Reference Peyneau and Roux2008). (b) Variation of $\phi ^{0}-\phi$ with $I$ at three applied pressures (see legend) for particles with $\mu _{s}=0.0$ (red) and $\mu _{s}=0.3$ (black). The dotted lines represent power-law fits from (4.8). (c) Variation of $\phi ^{0}$ and (d) $\alpha _{4}$ with $\mu _{s}$. The open symbols in panel (c) and panel (d) indicate the values for zero friction.

The quasi-static solid volume fraction $\phi ^{0}$ varies significantly with $\mu _{s}$ ranging from $\phi ^{0}=0.63$ for frictionless particles to $\phi ^{0}=0.59$ for particles with high friction, as shown in figure 6(c). Such a dependence of $\phi ^{0}$ on friction has been previously demonstrated in two-dimensional (Da Cruz et al. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005) and three-dimensional simulations (Sun & Sundaresan Reference Sun and Sundaresan2011), and confirmed in recent experiments (Tapia, Pouliquen & Guazzelli Reference Tapia, Pouliquen and Guazzelli2019). The similarity between $\phi ^{0}=0.63$ for frictionless particles and the solid volume fraction of random close packing of monodisperse spheres indicates that frictionless particles do not dilate at the onset of flow, which is consistent with prior simulations (Peyneau & Roux Reference Peyneau and Roux2008) and experiments (Clavaud et al. Reference Clavaud, Bérut, Metzger and Forterre2017). For particles with high friction, $\phi ^{0}$ is consistent with the critical state solid volume fraction from the critical state theory (Schofield & Wroth Reference Schofield and Wroth1968), and previous simulations (Sun & Sundaresan Reference Sun and Sundaresan2011; de Coulomb et al. Reference de Coulomb, Bouzid, Claudin, Clément and Andreotti2017; Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2019) and experiments (Boyer et al. Reference Boyer, Guazzelli and Pouliquen2011a; Tapia et al. Reference Tapia, Pouliquen and Guazzelli2019).

The power-law exponent varies between $\alpha _4=0.82$ for frictionless particles and $\alpha _4=0.92$ for high friction particles, as shown in figure 6(d). The value of this exponent at high friction is similar to previous theoretical predictions of a unity exponent (DeGiuli et al. Reference DeGiuli, Düring, Lerner and Wyart2015, Reference DeGiuli, McElwaine and Wyart2016). For frictionless particles, our prediction of $\alpha _4$ does not correspond well with the theoretical prediction (DeGiuli et al. Reference DeGiuli, Düring, Lerner and Wyart2015, Reference DeGiuli, McElwaine and Wyart2016) of $\alpha _4=0.35$, and a previous simulation study (Peyneau & Roux Reference Peyneau and Roux2008) that demonstrated $\alpha _4=0.39$. However, as shown in figure 6(a), our data corresponds well with the simulations of Peyneau & Roux (Reference Peyneau and Roux2008) at low and moderate inertial numbers, but deviates slightly at higher inertial numbers, resulting in large changes to the power-law exponent.

5. Normal stress differences and their microstructural origins

The non-negligible second-order contributions to stress in viscometric granular flows indicate the presence of normal stress differences. Previous research on sheared granular and suspension flows has demonstrated the existence of normal stress differences (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Alam & Luding Reference Alam, Luding, Garcia-Rojo, Herrmann and McNamara2005; Rycroft et al. Reference Rycroft, Kamrin and Bazant2009; Boyer et al. Reference Boyer, Pouliquen and Guazzelli2011b; Couturier et al. Reference Couturier, Boyer, Pouliquen and Guazzelli2011; Sun & Sundaresan Reference Sun and Sundaresan2011; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013; Saha & Alam Reference Saha and Alam2016; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018; Seto & Giusteri Reference Seto and Giusteri2018), and these differences have been attributed to flow-induced fluctuating velocity effects in dilute granular flows (Saha & Alam Reference Saha and Alam2016) and microstructural effects in dense suspension flows (Seto & Giusteri Reference Seto and Giusteri2018). Particularly, normal stress differences can arise either from: (1) a misalignment of $\boldsymbol {\sigma }$ and $\boldsymbol {D}$ in the flow plane, known as the first normal stress difference; or (2) from the anisotropy of normal stress between the flow plane and the vorticity direction, known as the second normal stress difference.

In this section, we describe normal stress differences and their microstructural origins in dense viscometric granular flows. The microstructure of a granular material is quantified through a second-rank contact fabric tensor $\boldsymbol {A}$, which provides a convenient description of the directional distribution of the particle contact network and inherent structural anisotropy (Oda Reference Oda1982; Kanatani Reference Kanatani1984). The orientational distribution $P(\boldsymbol {n})$ of contact normal unit vectors $\boldsymbol {n}$ can be expanded to the second order in Fourier series as (Rothenburg & Bathurst Reference Rothenburg and Bathurst1989; Bathurst & Rothenburg Reference Bathurst and Rothenburg1990)

(5.1)\begin{equation} P(\boldsymbol{n})=\frac{1}{4{\rm \pi}}\left[1+\boldsymbol{A}: \left(\boldsymbol{n}\otimes\boldsymbol{n}\right)\right], \end{equation}

where $\boldsymbol {A}$ is trace free and symmetric. For dense granular materials where internal stress $\boldsymbol {\sigma }$ is dominated by particle contacts, it can be expressed as the following integral in the orientational space $\varOmega$ (Rothenburg & Bathurst Reference Rothenburg and Bathurst1989; Bathurst & Rothenburg Reference Bathurst and Rothenburg1990; Srivastava et al. Reference Srivastava, Lechman, Grest and Silbert2020):

(5.2)\begin{equation} \sigma_{ij}=\frac{N_{c}\langle l_{0} \rangle \langle \,f_{i} \rangle}{V} \int_{\varOmega} P(\boldsymbol{n}) n_{j}\,\mathrm{d}\boldsymbol{n}, \end{equation}

where $\langle l_{0} \rangle$ and $\langle\, f_{i} \rangle$ are the average magnitudes of the branch vector and the $i$th component of the normal force between two contacting particles, respectively. This representation of the stress tensor ensures that $\boldsymbol {\sigma }$ and $\boldsymbol {A}$ have the same structure.

5.1. Second normal stress difference

Significant normal stress anisotropy emerges from the difference between the mean normal stress in the flow plane and normal stress in the vorticity direction, and is represented by the viscometric flow function $N_{0}/\tau =(2\sigma _{zz}-\sigma _{xx}-\sigma _{yy})/2\tau$ (Seto & Giusteri Reference Seto and Giusteri2018), where $x$ is the flow direction, $y$ is the flow gradient direction and $z$ is the vorticity direction, and the stresses are defined positive in the compressive sense, since the forces are all repulsive. In the present simulations, $N_{0}/\tau$ is computed by

(5.3)\begin{equation} \frac{N_{0}}{\tau} = \frac{-3\boldsymbol{\sigma}:\left(\boldsymbol{D}^{2}-\frac{\mathrm{tr} \left(\boldsymbol{D}^{2}\right)}{3}\boldsymbol{I}\right)}{\dot{\gamma}\boldsymbol{\sigma}:\boldsymbol{D}}, \end{equation}

which is equivalent to $N_{0}/\tau =-\mu _{2}/\mu _{1}$. In figure 7(b), $N_{0}/\tau$ is plotted as a function of the distance to quasi-static solid volume fraction $\phi -\phi _{0}$ for five $\mu _{s}$. The negative value of $N_{0}$ implies that there is more normal stress in the flow plane than in the vorticity direction, as seen in figure 7(a,b), and the ratio of the two normal stresses is consistent with previous simulations on dry granular flows (c.f. figure 7c) (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013). In present simulations, an imbalance between the external applied pressure and the internal pressure drives isochoric periodic cell deformation, as shown by the equal values of $F_{ii}(t)$ in figure 2(b). In another scenario where each $\sigma _{ii}$ is individually balanced, we observed a rapid compaction of the cell in the vorticity direction leading to simulation instability arising from the second normal stress difference. The magnitude of second normal stress difference is larger for frictional particles than for frictionless particles; however, even frictionless particles exhibit non-zero second normal stress difference during flow at finite inertial numbers, as seen in figure 7(a). As the solid volume fraction increases towards quasi-static $\phi _{0}$, the anisotropy consistently decreases for all $\mu _{s}$. For frictionless particles, $N_{0}$ appears to tend to zero in the quasi-static limit corresponding to $\phi ^{0}=0.63$, which is similar to the random close packing density for monodisperse spheres. However, the out of flow plane stress anisotropy is demonstrably non-zero for frictional particles even in the quasi-static limit, as also observed previously by Seto & Giusteri (Reference Seto and Giusteri2018). The notion of non-zero anisotropy in the quasi-static regime is also consistent with recent observations of an anisotropic yield surface in frictional granular materials (Thornton & Zhang Reference Thornton and Zhang2010; Li & Dafalias Reference Li and Dafalias2012).

Figure 7. Variation of scaled second normal stress difference $N_{0}/\tau$ with inertial number (a) $I$ and (b) distance to quasi-static solid volume fraction $\phi -\phi _{0}$, for five interparticle frictions $\mu _{s}$ (see legend) at applied pressure $p_{{ext}}=10^{-4}$. (c) Variation of $2\sigma _{zz}/(\sigma _{xx}+\sigma _{yy})$ with $I$ for five interparticle frictions. Variation of scaled first normal stress difference $N_{1}/\tau$ with inertial number (d) $I$ and (e) distance to quasi-static solid volume fraction $\phi -\phi _{0}$, for five interparticle frictions $\mu _{s}$ (see legend) at applied pressure $p_{{ext}}=10^{-4}$. (f) Variation of $\sigma _{yy}/\sigma _{xx}$ with $I$ for five interparticle frictions.

An implication of these findings is that the flow of frictional granular materials is not codirectional, i.e. the hypothesis $\boldsymbol {\sigma }\propto \boldsymbol {D}$, which has been assumed within the $\mu (I)$ rheological model is not accurate. Here, $N_{0}/\tau$ increases with $I$ for all $\mu _{s}$, as seen in figure 7(a), and remains measurably non-zero for frictional particles even at low $I$. Prior simulations on quasi-static simple shear granular flows (Sun & Sundaresan Reference Sun and Sundaresan2011), granular flows down an incline (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013), gravity-driven granular flows through an orifice (Rycroft et al. Reference Rycroft, Kamrin and Bazant2009) and granular flows in a split-bottom Couette cell (Depken et al. Reference Depken, Lechman, Van Hecke, Van Saarloos and Grest2007) have questioned the codirectionality hypothesis. Two previously proposed theoretical models – double shearing (Spencer Reference Spencer1964) and shear-free sheets (Depken et al. Reference Depken, Lechman, Van Hecke, Van Saarloos and Grest2007) – have also incorporated these effects for quasi-static and dense granular flows.

The second normal stress difference results from an excess of contacts oriented in the flow plane rather than those oriented in the vorticity direction. Figure 8(b) shows the variation of $N_{0}^{a}/Z_{2}$ with $\phi -\phi _{0}$ for various interparticle friction. Here, $N_{0}^{a}=(A_{zz} - {(A_{xx}+A_{yy})}/{2})$ is the contact fabric ‘second normal difference’, which represents the anisotropy in average orientation of contacts between the flow plane and the vorticity direction. The rattler-free coordination number is computed as $Z_{2}=2N_{c}/(N-N_{r})$, where $N_{r}$ is the number of rattler particles with fewer than two contacts, and $N_{c}$ is the total number of contacts with non-zero normal force belonging to non-rattler particles (Sun & Sundaresan Reference Sun and Sundaresan2011). At high $I$ corresponding to low $\phi$, a higher fraction of contacts are oriented in the flow plane, which results in a large normal stress difference, as shown figure 8(a). Furthermore, all the data collapses onto a single curve for all interparticle friction. Upon approach to the quasi-static regime at high $\phi$, the contact distribution becomes more isotropic, resulting in reduced normal stress difference. For frictionless particles, the orientational distribution of contacts expectedly becomes isotropic in the quasi-static regime at random close packing volume fraction, as seen by $N_{0}^{a} \to 0$.

Figure 8. Variation of contact fabric ‘second normal difference’ $N_{0}^{a}$ scaled by rattler-free coordination $Z_{2}$ with (a) inertial number $I$ and (b) distance to quasi-static solid volume fraction $\phi -\phi _{0}$ for the five $\mu _{s}$ (see legend) at applied pressure $p_{{ext}}=10^{-4}$. (c) A schematic depicting the misalignment angle $\theta _{c}$ between the principal directions of $\boldsymbol {D}$ and $\boldsymbol {A}$ in the flow plane (shown in red). Variation of $\theta _{c}$ with (d) $N_{1}/\tau$ and (e) inertial number $I$ for the five $\mu _{s}$ (see legend) at applied pressure $p_{{ext}}=10^{-4}$. The vertical and horizontal error bars in panel (a,b) and panel (d,e) represent the standard error in the calculations.

5.2. First normal stress difference

The first normal stress difference, which characterizes the anisotropy between $\boldsymbol {\sigma }$ and $\boldsymbol {D}$ in the flow plane, is represented by the viscometric flow function $N_{1}/\tau =(\sigma _{yy}-\sigma _{xx})/\tau$ (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). In the present simulations, this is estimated from

(5.4)\begin{equation} \frac{N_{1}}{\tau} = \frac{\boldsymbol{\sigma}:\left(\boldsymbol{D}\boldsymbol{W}-\boldsymbol{W} \boldsymbol{D}\right)}{\dot{\gamma}\boldsymbol{\sigma}:\boldsymbol{D}}, \end{equation}

which is equivalent to $N_{1}/\tau = 4\mu _{3}/\mu _{1}$. The variation of $N_{1}/\tau$ with $\phi -\phi _{0}$ for five $\mu _{s}$ is displayed in figure 7(e). At low $\phi$, $N_{1}$ is negative for all $\mu _{s}$, and its magnitude increases with increasing $\mu _{s}$ for a given distance from the quasi-static solid volume fraction $\phi -\phi ^{0}$. At low $\phi$, the ratio of $N_{0}/N_{1}$ is approximately 3–4 for all interparticle friction, consistent with previous findings (Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014). The stress anisotropy in the flow plane increases with inertial number, as shown in figure 7(d), and also by the ratio of the two normal stresses in the flow plane, as shown in figure 7(f). The value of this ratio is consistent with previous simulations (c.f. figure 7f) on granular flows down an incline (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013).

When the flow becomes dense, $N_{1}$ increases towards zero and becomes slightly positive for highly dense flows in the quasi-static regime. The change of sign of $N_{1}$ at high $\phi$ has been previously observed in simulations (Alam & Luding Reference Alam, Luding, Garcia-Rojo, Herrmann and McNamara2005; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013; Seto & Giusteri Reference Seto and Giusteri2018) and experiments (Couturier et al. Reference Couturier, Boyer, Pouliquen and Guazzelli2011), but its existence is debated, and has been attributed to interparticle friction (Dbouk, Lobry & Lemaire Reference Dbouk, Lobry and Lemaire2013) and boundary wall effects in experiments (Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014). In the present simulations, we observe slightly positive $N_{1}$ for all values of $\mu _{s}$ at large $\phi$, and there are no boundary effects in these bulk simulations. Recently it was demonstrated that finite particle stiffness – which is often used as numerical regularization in hard particle simulations – causes $N_{1}$ to become positive at large $\phi$ in simulations on inertialess frictional suspensions (Seto & Giusteri Reference Seto and Giusteri2018). However, our granular simulations do not provide any conclusive evidence of vanishing positive $N_{1}$ as a result of increasing particle stiffness. A careful analysis about this effect constitutes a part of our future work.

The first normal stress difference is related to the angular misalignment $\theta _{c}$ between the principal directions of $\boldsymbol {D}$ ($\hat {\boldsymbol {d}_{1}}$ and $\hat {\boldsymbol {d}_{3}}$) and $\boldsymbol {A}$ ($\hat {\boldsymbol {a}}_{1}$ and $\hat {\boldsymbol {a}}_{3}$) in the flow plane, as described in the schematic in figure 8(c). Here, $\hat {\boldsymbol {d}_{1}}$ and $\hat {\boldsymbol {d}_{3}}$ represent the compression and expansion directions of shear flow, respectively, and $\theta _{c}$ represents the angle between $\hat {\boldsymbol {d}_{1}}$ and the major principal direction $\hat {\boldsymbol {a}}_{1}$ of $\boldsymbol {A}$. The misalignment angle $\theta _{c}$ and $N_{1}/\tau$ are strongly correlated, as depicted in figure 8(d) for various $\mu _{s}$ that largely collapse onto a single curve, indicating a one-to-one correspondence between stress anisotropy and the misalignment (Seto & Giusteri Reference Seto and Giusteri2018). The misalignment between $\boldsymbol {A}$ and $\boldsymbol {D}$ results in excess stress along the flow direction as compared with the gradient direction, which sets the negative sign of first normal stress difference, similar to previous observations in dry granular flows (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001; Weinhart et al. Reference Weinhart, Hartkamp, Thornton and Luding2013). Such microstructural origins of $N_{1}$ arising from a misalignment between the projected contact vectors and principal flow direction in the flow plane were previously demonstrated for inertialess frictional suspensions (Seto & Giusteri Reference Seto and Giusteri2018), and in the case of dilute granular flows through a similar misalignment between fluctuating velocity moment tensor and principal flow direction in the flow plane (Saha & Alam Reference Saha and Alam2016). Remarkably, $\theta _{c}\to 0$ as $N_{1}\to 0$, indicating a vanishing misalignment between $\boldsymbol {D}$ and $\boldsymbol {A}$ at high solid volume fractions. This is also seen by the variation of $\theta _{c}$ with $I$ in figure 8(e), where the data for all $\mu _{s}$ collapse onto a single curve. A small positive first normal stress difference exists at high solid volume fractions in the vicinity of yield stress despite a near-complete alignment of $\boldsymbol {D}$ and $\boldsymbol {A}$ in the flow plane, thus indicating that either a different underlying physical phenomenon is responsible for positive $N_{1}$, or it is possibly a consequence of finite system size.

The observations of microstructure-induced normal stress differences in dense granular flows, especially the collapse of $N_{0}^{a}/Z_{2}$ and $\theta _{c}$ with $I$ in figure 8(a,e), indicate that a fabric tensor is an appropriate internal state variable that can be used to construct a rheological model with evolution equations for the microstructure, as was done for quasi-static granular flows (Sun & Sundaresan Reference Sun and Sundaresan2011; Parra & Kamrin Reference Parra and Kamrin2019). Such an approach has been previously used in modelling suspension rheology (Goddard Reference Goddard2006; Stickel, Phillips & Powell Reference Stickel, Phillips and Powell2006), and a general framework adaptable to dry granular flows has been provided by Goddard (Reference Goddard2014).

6. Conclusions

In this paper, we described a discrete element method to simulate dense granular flows under external applied stress in a fully periodic representative volume element. Rather than prescribing solid volume fraction and/or strain rate, this method enables independent evolution of a solid volume fraction and three-dimensional strain rate tensor in response to an imbalance between internal state of stress and external applied stress. Using this method, bulk viscometric granular flows were simulated under external pressure and shear stress, which was devoid of any boundary effects, and thus closely represented the boundary conditions often found in practice.

We developed a second-order rheological model to relate the internal Cauchy stress $\boldsymbol {\sigma }$ with the strain rate tensor $\boldsymbol {D}$ for various interparticle friction. The model considers both rate-dependent and rate-independent contributions to the total stress, where the latter is often described using models of granular plasticity. The rheological model well-predicts the $\mu (I)$ rheology of granular materials. Additionally, it also predicts normal stress differences in steady viscometric granular flows, which have often been observed in simulations and experiments, but have not been well-characterized. A major implication of this model is that it does not impose coaxiality between $\boldsymbol {\sigma }$ and $\boldsymbol {D}$ in dense granular flows, which is often assumed in several other constitutive models.

A major focus of this work has been to highlight the role of interparticle friction on viscometric granular rheology in the dense flowing regime, particularly on the two normal stress differences. We found that friction not only increases the quasi-static shear stress ratio, but also the quasi-static value of the second normal stress difference, thus indicating the presence of an anisotropic yield stress, whereas frictionless particles do not exhibit such anisotropy in the quasi-static regime at solid volume fraction similar to the random close packing of monodisperse spheres. At higher flow rates in the inertial regime, friction consistently increases the magnitude of both normal stress differences, indicating an increasing departure from the coaxiality of $\boldsymbol {\sigma }$ and $\boldsymbol {D}$. Although the second normal stress difference is always negative, the first normal stress difference changes sign from negative to positive at high solid volume fractions in the quasi-static regime. Further microstructural investigations highlighted that negative first normal stress difference results from a misalignment between $\boldsymbol {D}$ and a second-rank contact fabric tensor $\boldsymbol {A}$ in the flow plane, which describes the orientational distribution of sphere–sphere contacts in granular flows. Furthermore, the magnitude of misalignment increases with the inertial number similarly for all interparticle frictions. The second normal stress difference results from an excess of contacts oriented in the flow plane rather than in the vorticity direction, which is also observed from the anisotropy in the normal components of the fabric tensor. Upon appropriate normalization with friction-dependent coordination number, the fabric tensor anisotropy was shown to collapse onto a single curve for all interparticle frictions.

These results demonstrate the importance of developing rheological models beyond simple scalar models to predict granular rheology in even simple shear flows, and certainly for complex and heterogeneous flow fields that are observed in practice. The breakdown of coaxiality of stress and strain rate tensors highlights the need for an anisotropic rheological model that includes a contact fabric tensor as an internal variable. A general form of such anisotropic models $\boldsymbol {\sigma }=\mathcal {F}(\boldsymbol {D},\boldsymbol {A})$ was recently proposed for granular materials and suspensions (Goddard Reference Goddard2006; Stickel et al. Reference Stickel, Phillips and Powell2006; Goddard Reference Goddard2014), and was calibrated for rate-independent granular flows (Sun & Sundaresan Reference Sun and Sundaresan2011; Parra & Kamrin Reference Parra and Kamrin2019). The calibration of these models in rate-dependent flows is required, along with other non-viscometric flows such as uniaxial or triaxial compression, as well as transient evolving inertial granular flows. These topics are currently a subject of our ongoing study. Lastly, the constitutive model described here could also be extended to include other important granular flow phenomena such as hysteresis (Degiuli & Wyart Reference Degiuli and Wyart2017; Perrin et al. Reference Perrin, Clavaud, Wyart, Metzger and Forterre2019) and non-locality (Henann & Kamrin Reference Henann and Kamrin2013) at low inertial numbers.

Acknowledgements

The authors acknowledge helpful discussions with D. Henann and K. Kamrin. This work was performed at the Center for Integrated Nanotechnologies, a US Department of Energy and Office of Basic Energy Sciences user facility. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the US Department of Energy's National Nuclear Security Administration under Contract no. DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the US Department of Energy or the US Government.

Declaration of interest

The authors report no conflict of interest.

Appendix. Fitting parameters of the rheological model

In this Appendix, Table 1 provides the fitting parameters $A_{1}$, $A_{2}$, $A_{3}$ and $A_{4}$ of the rheological model, defined in (4.2), (4.4), (4.6) and (4.8), respectively, as a function of interparticle friction $\mu _{s}$.

Table 1. Fitting parameters corresponding to (4.2), (4.4), (4.6) and (4.8), as a function of interparticle friction $\mu _{s}$.

References

REFERENCES

Alam, M. & Luding, S. 2003 First normal stress difference and crystallization in a dense sheared granular fluid. Phys. Fluids 15 (8), 22982312.CrossRefGoogle Scholar
Alam, M. & Luding, S. 2005 Non-newtonian granular fluid: simulation and theory. In Powders and Grains (ed. Garcia-Rojo, R., Herrmann, H. J. & McNamara, S.), pp. 11411144. Balkema.Google Scholar
Anand, L. & Gu, C. 2000 Granular materials: constitutive equations and strain localization. J. Mech. Phys. Solids 48 (8), 17011733.CrossRefGoogle Scholar
Bagnold, R. A. 1954 Experiments on a gravity-free dispersion of large solid spheres in a newtonian fluid under shear. Proc. R. Soc. Lond. A 225 (1160), 4963.Google Scholar
Bathurst, R. J. & Rothenburg, L. 1990 Observations on stress-force-fabric relationships in idealized granular materials. Mech. Mater. 9 (1), 6580.CrossRefGoogle Scholar
Bhateja, A. & Khakhar, D. V. 2018 Rheology of dense granular flows in two dimensions: comparison of fully two-dimensional flows to unidirectional shear flow. Phys. Rev. Fluids 3 (6), 062301.CrossRefGoogle Scholar
Bird, R. B. & Hassager, O. 1987 Dynamics of Polymeric Liquids: Fluid mechanics, vol. 1. Wiley.Google Scholar
Boyer, F., Guazzelli, É. & Pouliquen, O. 2011 a Unifying suspension and granular rheology. Phys. Rev. Lett. 107 (18), 188301.CrossRefGoogle ScholarPubMed
Boyer, F., Pouliquen, O. & Guazzelli, É. 2011 b Dense suspensions in rotating-rod flows: normal stresses and particle migration. J. Fluid Mech. 686, 525.CrossRefGoogle Scholar
Brunn, P. O. & Asoud, H. 2003 An explicit constitutive equation of a simple fluid in motions with constant stretch history. J. Non-Newtonian Fluid Mech. 112 (2-3), 129139.CrossRefGoogle Scholar
Campbell, C. S. 2002 Granular shear flows at the elastic limit. J. Fluid Mech. 465, 261291.CrossRefGoogle Scholar
Campbell, C. S. 2005 Stress-controlled elastic granular shear flows. J. Fluid Mech. 539 (1), 273297.CrossRefGoogle Scholar
Clavaud, C., Bérut, A., Metzger, B. & Forterre, Y. 2017 Revealing the frictional transition in shear-thickening suspensions. Proc. Natl Acad. Sci. USA 114, 51475152.CrossRefGoogle ScholarPubMed
Coleman, B. D, Markovitz, H. & Noll, W. 1966 Viscometric Flows of Non-Newtonian Fluids: Theory and Experiment. Springer-Verlag.CrossRefGoogle Scholar
Cortet, P.-P., Bonamy, D., Daviaud, F., Dauchot, O., Dubrulle, B. & Renouf, M. 2009 Relevance of visco-plastic theory in a multi-directional inhomogeneous granular flow. Europhys. Lett. 88 (1), 14001.CrossRefGoogle Scholar
de Coulomb, A. F., Bouzid, M., Claudin, P., Clément, E. & Andreotti, B. 2017 Rheology of granular flows across the transition from soft to rigid particles. Phys. Rev. Fluids 2 (10), 102301.CrossRefGoogle Scholar
Couturier, É., Boyer, F., Pouliquen, O. & Guazzelli, É. 2011 Suspensions in a tilted trough: second normal stress difference. J. Fluid Mech. 686, 2639.CrossRefGoogle Scholar
Cundall, P. A. & Strack, O. D. L. 1979 A discrete numerical model for granular assemblies. Géotechnique 29 (1), 4765.CrossRefGoogle Scholar
Da Cruz, F., Emam, S., Prochnow, M., Roux, J. N. & Chevoir, F. 2005 Rheophysics of dense granular materials: discrete simulation of plane shear flows. Phys. Rev. E 72 (2), 021309.CrossRefGoogle ScholarPubMed
Dbouk, T., Lobry, L. & Lemaire, E. 2013 Normal stresses in concentrated non-brownian suspensions. J. Fluid Mech. 715, 239272.CrossRefGoogle Scholar
DeGiuli, E., Düring, G., Lerner, E. & Wyart, M. 2015 Unified theory of inertial granular flows and non-Brownian suspensions. Phys. Rev. E 91, 062206.CrossRefGoogle ScholarPubMed
DeGiuli, E., McElwaine, J. N. & Wyart, M. 2016 Phase diagram for inertial granular flows. Phys. Rev. E 94, 012904.CrossRefGoogle ScholarPubMed
Degiuli, E. & Wyart, M. 2017 Friction law and hysteresis in granular materials. Proc. Natl Acad. Sci. USA 114 (35), 92849289.CrossRefGoogle ScholarPubMed
Depken, M., Lechman, J. B., Van Hecke, M., Van Saarloos, W. & Grest, G. S. 2007 Stresses in smooth flows of dense granular media. Europhys. Lett. 78 (5), 58001.CrossRefGoogle Scholar
Depken, M., Van Saarloos, W. & Van Hecke, M. 2006 Continuum approach to wide shear zones in quasistatic granular matter. Phys. Rev. E 73 (3), 031302.CrossRefGoogle ScholarPubMed
Dijksman, J. A., Wortel, G. H., Van Dellen, L. T. H., Dauchot, O. & Van Hecke, M. 2011 Jamming, yielding, and rheology of weakly vibrated granular media. Phys. Rev. Lett. 107 (10), 108303.CrossRefGoogle ScholarPubMed
Evans, D. J. & Ely, J. F. 1986 Viscous flow in the stress ensemble. Mol. Phys. 59, 10431048.CrossRefGoogle Scholar
Flyvbjerg, H. & Petersen, H. G. 1989 Error estimates on averages of correlated data. J. Chem. Phys. 91 (1), 461466.CrossRefGoogle Scholar
Forterre, Y. & Pouliquen, O. 2008 Flows of dense granular media. Annu. Rev. Fluid Mech. 40 (1), 124.CrossRefGoogle Scholar
Gallier, S., Lemaire, E., Peters, F. & Lobry, L. 2014 Rheology of sheared suspensions of rough frictional particles. J. Fluid Mech. 757, 514549.CrossRefGoogle Scholar
Gao, Z., Zhao, J., Li, X. & Dafalias, Y. F. 2014 A critical state sand plasticity model accounting for fabric evolution. Intl J. Numer. Anal. Meth. Geomech. 38 (4), 370390.CrossRefGoogle Scholar
Giusteri, G. G. & Seto, R. 2018 A theoretical framework for steady-state rheometry in generic flow conditions. J. Rheol. 62 (3), 713723.CrossRefGoogle Scholar
Goddard, J. D. 1984 Dissipative materials as models of thixotropy and plasticity. J. Non-Newtonian Fluid Mech. 14, 141160.CrossRefGoogle Scholar
Goddard, J. D. 1986 Dissipative materials as constitutive models for granular media. Acta Mechanica 63 (1-4), 313.CrossRefGoogle Scholar
Goddard, J. D. 2006 A dissipative anisotropic fluid model for non-colloidal particle dispersions. J. Fluid Mech. 568, 117.CrossRefGoogle Scholar
Goddard, J. D. 2014 Continuum modeling of granular media. Appl. Mech. Rev. 66 (5), 050801.CrossRefGoogle Scholar
Guazzelli, É. & Pouliquen, O. 2018 Rheology of dense granular suspensions. J. Fluid Mech. 852, P1.CrossRefGoogle Scholar
Henann, D. L. & Kamrin, K. 2013 A predictive, size-dependent continuum model for dense granular flows. Proc. Natl Acad. Sci. USA 110, 67306735.CrossRefGoogle ScholarPubMed
Hood, L. M., Evans, D. J. & Morriss, G. P. 1987 Time correlation functions in the stress ensemble. Mol. Phys. 62, 419428.CrossRefGoogle Scholar
Jop, P., Forterre, Y. & Pouliquen, O. 2006 A constitutive law for dense granular flows. Nature 441 (7094), 727730.CrossRefGoogle ScholarPubMed
Kamrin, K. & Koval, G. 2014 Effect of particle surface friction on nonlocal constitutive behavior of flowing granular media. Comput. Part. Mech. 1 (2), 169176.CrossRefGoogle Scholar
Kanatani, K.-I. 1984 Distribution of directional data and fabric tensors. Intl J. Engng Sci. 22 (2), 149164.Google Scholar
Koval, G., Roux, J.-N., Corfdir, A. & Chevoir, F. 2009 Annular shear of cohesionless granular materials: from the inertial to quasistatic regime. Phys. Rev. E 79 (2), 021306.CrossRefGoogle ScholarPubMed
Larson, R. G. 1985 Flows of constant stretch history for polymeric materials with power-law distributions of relaxation times. Rheol. Acta 24 (5), 443449.CrossRefGoogle Scholar
Li, X. S. & Dafalias, Y. F. 2012 Anisotropic critical state theory: role of fabric. J. Engng Mech. ASCE 138 (3), 263275.CrossRefGoogle Scholar
Lois, G., Lemaître, A. & Carlson, J. M. 2005 Numerical tests of constitutive laws for dense granular flows. Phys. Rev. E 72 (5), 051303.CrossRefGoogle ScholarPubMed
Lu, K., Brodsky, E. E. & Kavehpour, H. P. 2007 Shear-weakening of the transitional regime for granular flow. J. Fluid Mech. 587, 347372.CrossRefGoogle Scholar
Martyna, G. J., Tuckerman, M. E., Tobias, D. J. & Klein, M. L. 1996 Explicit reversible integrators for extended systems dynamics. Mol. Phys. 87, 11171157.CrossRefGoogle Scholar
McElwaine, J. N., Takagi, D. & Huppert, H. E. 2012 Surface curvature of steady granular flows. Granul. Matt. 14 (2), 229234.CrossRefGoogle Scholar
Mehandia, V., Gutam, K. J. & Nott, P. R. 2012 Anomalous stress profile in a sheared granular column. Phys. Rev. Lett. 109 (12), 128002.CrossRefGoogle Scholar
Mehrabadi, M. M. & Cowin, S. C. 1978 Initial planar deformation of dilatant granular materials. J. Mech. Phys. Solids 26 (4), 269284.CrossRefGoogle Scholar
MiDi, G. D. R. 2004 On dense granular flows. Eur. Phys. J. E 14 (4), 341365.CrossRefGoogle Scholar
Nemat-Nasser, S. 2000 A micromechanically-based constitutive model for frictional deformation of granular materials. J. Mech. Phys. Solids 48 (6-7), 15411563.CrossRefGoogle Scholar
Noll, W. 1962 Motions with constant stretch history. Arch. Rat. Mech. Anal. 11 (1), 97105.CrossRefGoogle Scholar
Oda, M. 1982 Fabric tensor for discontinuous geological materials. Soils Found. 22 (4), 96108.CrossRefGoogle Scholar
Otsuki, M. & Hayakawa, H. 2011 Critical scaling near jamming transition for frictional granular particles. Phys. Rev. E 83 (5), 051301.CrossRefGoogle ScholarPubMed
Parra, E. R. & Kamrin, K. 2019 Capturing transient granular rheology with extended fabric tensor relations. Granul. Matt. 21 (89), 17.Google Scholar
Parrinello, M. & Rahman, A. 1981 Polymorphic transitions in single crystals: a new molecular dynamics method. J. Appl. Phys. 52 (12), 71827190.CrossRefGoogle Scholar
Perrin, H., Clavaud, C., Wyart, M., Metzger, B. & Forterre, Y. 2019 Interparticle friction leads to nonmonotonic flow curves and hysteresis in viscous suspensions. Phys. Rev. X 9, 031027.Google Scholar
Peyneau, P. E. & Roux, J. N. 2008 Frictionless bead packs have macroscopic friction, but no dilatancy. Phys. Rev. E 78 (1), 011307.CrossRefGoogle ScholarPubMed
Plimpton, S. 1995 Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117 (1), 119.CrossRefGoogle Scholar
Pouliquen, O. 1999 Scaling laws in granular flows down rough inclined planes. Phys. Fluids 11 (3), 542548.CrossRefGoogle Scholar
Radjai, F., Roux, J. N. & Daouadji, A. 2017 Modeling granular materials: century-long research across scales. J. Engng Mech. ASCE 143 (4), 04017002.CrossRefGoogle Scholar
Rajagopal, K. R. 2006 On implicit constitutive theories for fluids. J. Fluid Mech. 550, 243249.CrossRefGoogle Scholar
Rivlin, R. S. 1955 Stress-deformation relations for isotropic materials. Arch. Rat. Mech. Anal. 4, 323425.Google Scholar
Rothenburg, L. & Bathurst, R. J. 1989 Analytical study of induced anisotropy in idealized granular materials. Géotechnique 39 (4), 601614.CrossRefGoogle Scholar
Rycroft, C. H., Kamrin, K. & Bazant, M. Z. 2009 Assessing continuum postulates in simulations of granular flow. J. Mech. Phys. Solids 57 (5), 828839.CrossRefGoogle Scholar
Saha, S. & Alam, M. 2016 Normal stress differences, their origin and constitutive relations for a sheared granular fluid. J. Fluid Mech. 795, 549580.CrossRefGoogle Scholar
Salerno, K. M., Bolintineanu, D. S., Grest, G. S., Lechman, J. B., Plimpton, S. J., Srivastava, I. & Silbert, L. E. 2018 Effect of shape and friction on the packing and flow of granular materials. Phys. Rev. E 98 (5), 050901.CrossRefGoogle Scholar
Schofield, A. & Wroth, P. 1968 Critical State Soil Mechanics, vol. 310. McGraw-Hill.Google Scholar
Schuhmacher, P., Radjai, F. & Roux, S. 2017 Wall roughness and nonlinear velocity profiles in granular shear flows. In EPJ Web of Conferences, vol. 140, p. 03090. EDP Sciences.CrossRefGoogle Scholar
Schunk, P. R. & Scriven, L. E. 1990 Constitutive equation for modeling mixed extension and shear in polymer solution processing. J. Rheol. 34 (7), 10851119.CrossRefGoogle Scholar
Seto, R. & Giusteri, G. G. 2018 Normal stress differences in dense suspensions. J. Fluid Mech. 857, 200215.CrossRefGoogle Scholar
Shinoda, W., Shiga, M. & Mikami, M. 2004 Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Phys. Rev. B 69 (13), 134103.CrossRefGoogle Scholar
Shojaaee, Z., Brendel, L., Török, J. & Wolf, D. E. 2012 a Shear flow of dense granular materials near smooth walls. II. Block formation and suppression of slip by rolling friction. Phys. Rev. E 86, 011302.CrossRefGoogle ScholarPubMed
Shojaaee, Z., Roux, J. N., Chevoir, F. & Wolf, D. E. 2012 b Shear flow of dense granular materials near smooth walls. I. Shear localization and constitutive laws in the boundary region. Phys. Rev. E 86, 011301.CrossRefGoogle Scholar
Silbert, L. E., Ertaş, D., Grest, G. S., Halsey, T. C., Levine, D. & Plimpton, S. J. 2001 Granular flow down an inclined plane: bagnold scaling and rheology. Phys. Rev. E 64 (5), 051302.CrossRefGoogle ScholarPubMed
Smith, K. C., Srivastava, I., Fisher, T. S. & Alam, M. 2014 Variable-cell method for stress-controlled jamming of athermal, frictionless grains. Phys. Rev. E 89 (4), 042203.CrossRefGoogle ScholarPubMed
Sokolovskii, V. V. 1965 Statics of Granular Media. Pergamon Press.Google Scholar
Souza, I. & Martins, J. L. 1997 Metric tensor as the dynamical variable for variable-cell-shape molecular dynamics. Phys. Rev. B 55 (14), 8733.CrossRefGoogle Scholar
Spencer, A. J. M. 1964 A theory of the kinematics of ideal soils under plane strain conditions. J. Mech. Phys. Solids 12 (5), 337351.CrossRefGoogle Scholar
Srivastava, I. & Fisher, T. S. 2017 Slow creep in soft granular packings. Soft Matt. 13 (18), 34113421.CrossRefGoogle ScholarPubMed
Srivastava, I., Lechman, J. B., Grest, G. S. & Silbert, L. E. 2020 Evolution of internal granular structure at the flow-arrest transition. Granul. Matt. 22 (41), 18.CrossRefGoogle Scholar
Srivastava, I., Silbert, L. E., Grest, G. S. & Lechman, J. B. 2019 Flow-arrest transitions in frictional granular matter. Phys. Rev. Lett. 122 (4), 048003.CrossRefGoogle ScholarPubMed
Stickel, J. J., Phillips, R. J. & Powell, R. L. 2006 A constitutive model for microstructure and total stress in particulate suspensions. J. Rheol. 50 (4), 379413.CrossRefGoogle Scholar
Sun, J. & Sundaresan, S. 2011 A constitutive model with microstructure evolution for flow of rate-independent granular materials. J. Fluid Mech. 682, 590616.CrossRefGoogle Scholar
Tapia, F., Pouliquen, O. & Guazzelli, É. 2019 Influence of surface roughness on the rheology of immersed and dry frictional spheres. Phys. Rev. Fluid 4, 104302.CrossRefGoogle Scholar
Thompson, R. L. & Mendes, P. R. S. 2005 Persistence of straining and flow classification. Intl J. Engng Sci. 43 (1–2), 79105.CrossRefGoogle Scholar
Thornton, C. & Zhang, L. 2010 On the evolution of stress and microstructure during general 3D deviatoric straining of granular media. Géotechnique 60 (5), 333341.CrossRefGoogle Scholar
Wagner, C. E. & Mckinley, G. H. 2016 The importance of flow history in mixed shear and extensional flows. J. Non-Newtonian Fluid Mech. 233, 133145.CrossRefGoogle Scholar
Walton, O. R. & Braun, R. L. 1986 Viscosity, granular–temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J. Rheol. 30 (5), 949980.CrossRef