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

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