## 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 Guazzelli2011*b*), 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

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

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

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

where

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:

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:

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,

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.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 Wolf2012*b*), which requires careful calibration of wall properties to extract the bulk rheological properties (Shojaaee *et al.* Reference Shojaaee, Brendel, Török and Wolf2012*a*; 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

*a*)\begin{gather} \dot{\boldsymbol{r}}_{i}=\frac{\boldsymbol{p}_{i}}{m_{i}}+\frac{{\boldsymbol{\mathsf{P}}}_{g}}{W_{g}}\boldsymbol{r}_{i}, \end{gather}

*b*)\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}

*c*)\begin{gather}\dot{\boldsymbol{\mathsf{H}}}=\frac{{\boldsymbol{\mathsf{P}}}_{g}}{W_{g}}{\boldsymbol{\mathsf{H}}}, \end{gather}

*d*)\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.

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.10*d*), 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)

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.10*d*), the tensor $\boldsymbol {\varSigma }$ is defined as (Shinoda *et al.* Reference Shinoda, Shiga and Mikami2004)

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

It is evident from (2.10*d*) 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

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

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.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.10*a*). The contact normal force $\boldsymbol {f}_{nij}$ and tangential force $\boldsymbol {f}_{tij}$ on particle $i$ are given by

*a*)\begin{gather} \boldsymbol{f}_{nij}=k_{n} \delta_{ij} \boldsymbol{n}_{ij} - \gamma_{n} m_{e} \boldsymbol{v}_{nij}, \end{gather}

*b*)\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

*a*)\begin{gather} \boldsymbol{v}_{nij}=\left(\boldsymbol{v}_{ij}\boldsymbol{\cdot}\boldsymbol{n}_{ij}\right)\boldsymbol{n}_{ij}, \end{gather}

*b*)\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:

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.

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

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

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(*a*–*d*). 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

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),

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.

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 Pouliquen2011*a*; 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

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,

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.

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

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:

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

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:

*a*)\begin{gather} \sigma_{xx}=p\left(1+\frac{\mu_{2}}{3}-2\mu_{3}\right), \end{gather}

*b*)\begin{gather}\sigma_{yy}=p\left(1+\frac{\mu_{2}}{3}+2\mu_{3}\right), \end{gather}

*c*)\begin{gather}\sigma_{zz}=p\left(1-\frac{2\mu_{2}}{3}\right), \end{gather}

*d*)\begin{gather}\sigma_{xy}=p\mu_{1}, \end{gather}

*e*)\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),

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

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 Pouliquen2011*a*; 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 Guazzelli2011*b*; 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)

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

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

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 7*c*) (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).

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

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

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 7*f*) 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.