1 Introduction
The development of a robust continuum-mechanical description of granular flows has been a subject of active interest, due to its utility in modelling industrial processes and natural phenomena. In the regime of dense, slow flow, granular flows exhibit distinctive features, such as independence (or very weak dependence) of the stress on the magnitude of the strain rate (Rao & Nott Reference Rao and Nott2008), the formation of thin shear layers (Natarajan, Hunt & Taylor Reference Natarajan, Hunt and Taylor1995; Losert et al. Reference Losert, Bocquet, Lubensky and Gollub2000; Mueth et al. Reference Mueth, Debregeas, Karczmar, Eng, Nagel and Jaeger2000; Ananda, Moka & Nott Reference Ananda, Moka and Nott2008), and dilation or compaction of the medium upon shear deformation. While rate independence of the stress is incorporated in classical plasticity theories (Schofield & Wroth Reference Schofield and Wroth1968; Jackson Reference Jackson and Meyer1983; Rao & Nott Reference Rao and Nott2008), and considerable attention has been devoted in recent years to the kinematics (Mohan, Rao & Nott Reference Mohan, Rao and Nott2002; Pouliquen & Forterre Reference Pouliquen and Forterre2009; Bouzid et al. Reference Bouzid, Trulsson, Claudin, Clément and Andreotti2013; Henann & Kamrin Reference Henann and Kamrin2013), there is no satisfactory continuum description of dilatancy, or density change driven by shear. Dilatancy in granular materials was first observed by Reynolds (Reference Reynolds1885), and its occurrence within shear layers has been well documented (Desrues et al. Reference Desrues, Chambon, Mokni and Mazerolle1996; Mueth et al. Reference Mueth, Debregeas, Karczmar, Eng, Nagel and Jaeger2000; Kabla & Senden Reference Kabla and Senden2009). It was recently shown by Krishnaraj & Nott (Reference Krishnaraj and Nott2016) that dilatancy can drive a large-scale secondary flow which leaves a striking rheological signature (Mehandia, Gutam & Nott Reference Mehandia, Gutam and Nott2012; Gutam, Mehandia & Nott Reference Gutam, Mehandia and Nott2013).
Though rate independence of stress is incorporated in classical plasticity theories, such as the critical state theory (Schofield & Wroth Reference Schofield and Wroth1968; Jackson Reference Jackson and Meyer1983; Rao & Nott Reference Rao and Nott2008), the theory suffers from kinematic indeterminacy in problems of sustained flow: if the deformation rate tensor 
$\unicode[STIX]{x1D63F}\equiv \frac{1}{2}[\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}]$ (where 
$\boldsymbol{u}$ is the velocity field) satisfies the quasistatic balance of momentum, so will any multiple of 
$\unicode[STIX]{x1D63F}$ – as a result, 
$\boldsymbol{u}$ cannot be uniquely determined. In the past two decades, several models have been proposed to regularize classical plasticity and repair this shortcoming (Mohan et al. Reference Mohan, Rao and Nott2002; Pouliquen & Forterre Reference Pouliquen and Forterre2009; Bouzid et al. Reference Bouzid, Trulsson, Claudin, Clément and Andreotti2013; Henann & Kamrin Reference Henann and Kamrin2013). While these models are able to predict deformation in regions where the yield condition is apparently not satisfied, as observed in experiments, they all assume the granular medium to be incompressible, thereby ignoring density variation. The critical state theory allows variation of density with pressure (as in a compressible fluid), but not due to shear, which can occur even at constant pressure.
In this paper, we construct a non-local model that makes quasistatic flow kinematically determinate, and incorporates dilatancy. Our model is based on the simple idea that plastic deformation and dilation at any location are caused not just by yielding at that point, but in a mesoscopic region around it, thereby leading to a non-local constitutive relation for the stress. It is a systematic non-local extension of critical state plasticity, and introduces no additional variables or conservation equations. The model predictions of the velocity and density fields in simple shear in the absence and presence of gravity compare favourably with the results of particle dynamics simulations.
2 Model description
2.1 Local constitutive model
We briefly recall the local model, namely critical state plasticity (Schofield & Wroth Reference Schofield and Wroth1968; Jackson Reference Jackson and Meyer1983; Rao & Nott Reference Rao and Nott2008), before proceeding to its non-local extension. The model comprises a yield condition and a flow rule, which are generally written as
$$\begin{eqnarray}F(\unicode[STIX]{x1D748})\equiv \unicode[STIX]{x1D70F}(\unicode[STIX]{x1D748}^{\prime })-Y(p,\unicode[STIX]{x1D719})=0,\quad \unicode[STIX]{x1D60B}_{ij}=\dot{\unicode[STIX]{x1D706}}\,\frac{\displaystyle \unicode[STIX]{x2202}F}{\displaystyle \unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ji}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D748}$ is the stress tensor, 
$\unicode[STIX]{x1D748}^{\prime }\equiv \unicode[STIX]{x1D748}+p\,\unicode[STIX]{x1D739}$ its deviatoric part (
$\unicode[STIX]{x1D739}$ being the identity tensor), 
$p=-\frac{1}{3}\text{tr}(\unicode[STIX]{x1D748})$ is the pressure, 
$\unicode[STIX]{x1D719}$ is the volume fraction of particles, and 
$\dot{\unicode[STIX]{x1D706}}$ is the scalar ‘fluidity’ field. The functional forms of 
$\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D748}^{\prime })$ and the yield function 
$Y(p,\unicode[STIX]{x1D719})$ are specified by the choice of the yield condition. We use the extended von Mises condition (Rao & Nott Reference Rao and Nott2008; Nott Reference Nott2009), as it is simple and possesses all the features required for our description:
$$\begin{eqnarray}\unicode[STIX]{x1D70F}=(\unicode[STIX]{x1D748}^{\prime }\boldsymbol{ : }\unicode[STIX]{x1D748}^{\prime })^{1/2},\quad Y=\sqrt{2}\unicode[STIX]{x1D707}p_{c}\widehat{Y}(p/p_{c}),\end{eqnarray}$$ where the dimensionless function 
$\widehat{Y}=n(p/p_{c})-(n-1)(p/p_{c})^{n/(n-1)}$ specifies the shape of the yield surface (Prakash & Rao Reference Prakash and Rao1988). Here 
$p_{c}=\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D719})$ is the pressure at the critical state, a state of isochoric deformation. The material properties in (2.2) are the friction coefficient 
$\unicode[STIX]{x1D707}$, the exponent 
$n$ (typically slightly greater than unity) that determines the convexity the yield surface, and the function 
$\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D719})$. It is then straightforward to derive the expressions for the deviatoric and volumetric deformation rates (Srivastava & Sundaresan Reference Srivastava and Sundaresan2003; Nott Reference Nott2009)
$$\begin{eqnarray}\unicode[STIX]{x1D63F}^{\prime }=\frac{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}{\displaystyle 2\unicode[STIX]{x1D707}p_{c}\widehat{Y}}\,\unicode[STIX]{x1D748}^{\prime },\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\dot{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D707}p_{c}\frac{\displaystyle \unicode[STIX]{x2202}\widehat{Y}}{\displaystyle \unicode[STIX]{x2202}p},\end{eqnarray}$$ where 
$\dot{\unicode[STIX]{x1D6FE}}\equiv \left(2\unicode[STIX]{x1D63F}^{\prime }\boldsymbol{ : }\unicode[STIX]{x1D63F}^{\prime }\right)^{1/2}=\sqrt{2}\,\dot{\unicode[STIX]{x1D706}}$. Equation (2.3) can be readily inverted to a stress-explicit form (Nott Reference Nott2009):
$$\begin{eqnarray}\unicode[STIX]{x1D748}=-p\,\unicode[STIX]{x1D739}+\frac{\displaystyle 2\unicode[STIX]{x1D707}p_{c}\widehat{Y}}{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}\,\unicode[STIX]{x1D63F}^{\prime },\quad p=p_{c}\left(1-\frac{1}{n\unicode[STIX]{x1D707}\dot{\unicode[STIX]{x1D6FE}}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}\right)^{n-1}.\end{eqnarray}$$ We note that 
$p_{c}$ provides the stress scale in the critical state theory, and is based on observations in experiments and particle dynamics simulations of a one-to-one relation between the pressure and 
$\unicode[STIX]{x1D719}$ at critical state. The micromechanical origin of the stress scale in an assembly of nearly rigid particles is not clear, and calls for further study.
2.2 The non-local extension
We now describe how non-locality is introduced into the constitutive relation in a systematic and physically reasonable manner. Our principal argument is that the flow rule is non-local, in that the deformation rate at a point depends on the stress not just at that point, but in a small volume surrounding it. This is in consonance with a significant body of experimental evidence that fluctuations in the velocity and strain rate in dense granular materials are correlated over a mesoscale of 
${\sim}10$ grain diameters (Moka & Nott Reference Moka and Nott2005; Reddy, Forterre & Pouliquen Reference Reddy, Forterre and Pouliquen2011), whereby plastic yield at one location leads to deformation at a distance from it. Force chains provide one such mechanism – a break in a chain at one location will be felt at points which are within the persistence length of the chain (Majmudar & Behringer Reference Majmudar and Behringer2005). To incorporate this physical picture, our proposal is that the flow rule (2.3) must be modified to hold only in the following volume-averaged sense:
$$\begin{eqnarray}\unicode[STIX]{x1D63F}^{\prime }=\int _{\boldsymbol{ y}}^{}\left[\frac{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}{\displaystyle 2\unicode[STIX]{x1D707}p_{c}\widehat{Y}}\,\unicode[STIX]{x1D748}^{\prime }\right]g(\boldsymbol{y}-\boldsymbol{x})\,\text{d}\boldsymbol{y},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\int _{\boldsymbol{ y}}^{}\left[\dot{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D707}p_{c}\frac{\displaystyle \unicode[STIX]{x2202}\widehat{Y}}{\displaystyle \unicode[STIX]{x2202}p}\right]g(\boldsymbol{y}-\boldsymbol{x})\,\text{d}\boldsymbol{y},\end{eqnarray}$$ where the terms within square brackets are functions of 
$\boldsymbol{y}$, and the integrals are over all space. Here 
$g(\unicode[STIX]{x1D743})$ is a smoothly varying, isotropic weight function that decays fast enough as 
$|\unicode[STIX]{x1D743}|\rightarrow \infty$ so that all its moments are finite. For the purpose of our analysis, we do not need the detailed form of 
$g(\unicode[STIX]{x1D743})$, but only require that its zeroth and second moments satisfy
$$\begin{eqnarray}\int g(\unicode[STIX]{x1D743})\,\text{d}\unicode[STIX]{x1D743}=1,\quad \int g(\unicode[STIX]{x1D743})|\unicode[STIX]{x1D743}|^{2}\,\text{d}\unicode[STIX]{x1D743}=\ell ^{2}.\end{eqnarray}$$ The length 
$\ell$ is the effective radius of the averaging volume, and is a measure of non-locality. Equation (2.5) is similar in spirit to the proposal of Pouliquen & Forterre (Reference Pouliquen and Forterre2009), but quite different in the details; a more detailed discussion is given in § 4. Though (2.5) is derived for the extended von Mises yield condition and the associated flow rule, it can be generalized to any isotropic yield condition and flow rule.
As already discussed, deformation is accompanied by dilation (or compaction), and hence the volume fraction 
$\unicode[STIX]{x1D719}$ too must be non-local. To derive the corresponding relation for 
$\unicode[STIX]{x1D719}$, we first recognize that the local relation is 
$\unicode[STIX]{x1D719}=\unicode[STIX]{x1D6F1}^{-1}(p_{c})$. Following the same arguments that lead to (2.5), we then have
$$\begin{eqnarray}\unicode[STIX]{x1D719}=\int _{y}[\unicode[STIX]{x1D6F1}^{-1}(p_{c})]g(\boldsymbol{y}-\boldsymbol{x})\,\text{d}\boldsymbol{y}.\end{eqnarray}$$ Equations (2.5) and (2.7) are reminiscent of the memory integral formulation for the stress in polymer solutions (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1977), but here information on the state of the material is communicated in space rather than time. As in polymeric fluids, it is convenient to write (2.5) and (2.7) in differential form to solve boundary value problems. To do so, we expand the terms within square brackets in (2.5) and (2.7) in Taylor series about 
$\boldsymbol{y}=\boldsymbol{x}$, and use the moments of 
$g(\unicode[STIX]{x1D743})$ (2.6) to get
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}^{\prime }=\frac{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}{\displaystyle 2\unicode[STIX]{x1D707}p_{c}\widehat{Y}}\,\unicode[STIX]{x1D748}^{\prime }+\ell ^{2}\,\unicode[STIX]{x1D6FB}^{2}\left[\frac{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}{\displaystyle 2\unicode[STIX]{x1D707}p_{c}\widehat{Y}}\,\unicode[STIX]{x1D748}^{\prime }\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\dot{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D707}p_{c}\frac{\displaystyle \unicode[STIX]{x2202}\widehat{Y}}{\displaystyle \unicode[STIX]{x2202}p}+\ell ^{2}\,\unicode[STIX]{x1D6FB}^{2}\left[\dot{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D707}p_{c}\frac{\displaystyle \unicode[STIX]{x2202}\widehat{Y}}{\displaystyle \unicode[STIX]{x2202}p}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}=\unicode[STIX]{x1D6F1}^{-1}(p_{c})+\ell ^{2}\,\unicode[STIX]{x1D6FB}^{2}[\unicode[STIX]{x1D6F1}^{-1}(p_{c})], & \displaystyle\end{eqnarray}$$
$O(\ell ^{4})$ and higher. To invert (2.8) to a stress-explicit form, we write the stress as a perturbation expansion in 
$\ell ^{2}$, 
$\unicode[STIX]{x1D748}^{\prime }=\unicode[STIX]{x1D748}^{\prime (0)}+\ell ^{2}\unicode[STIX]{x1D748}^{\prime (1)}$, 
$p=p^{(0)}+\ell ^{2}p^{(1)}$, and 
$p_{c}=p_{c}^{(0)}+\ell ^{2}p_{c}^{(1)}$; it suffices to determine them to 
$O(\ell ^{2})$. The leading terms are the local quantities, given by (2.4). The 
$O(\ell ^{2})$ corrections are then obtained by simply substituting the leading order expressions for the terms within the square brackets in (2.8a–c), which are the left-hand sides of the respective equations. We then arrive at the constitutive relation 
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D748}=-p\,\unicode[STIX]{x1D739}+\frac{\displaystyle 2\unicode[STIX]{x1D707}}{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}(p_{c}\,\unicode[STIX]{x1D63F}^{\prime }-\ell ^{2}\,\unicode[STIX]{x1D6F1}\,\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}^{\prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p_{c}=\unicode[STIX]{x1D6F1}-\ell ^{2}\,\frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}\unicode[STIX]{x1D719}}\,\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p=p_{c}\left(1-\frac{\unicode[STIX]{x1D707}_{b}}{\dot{\unicode[STIX]{x1D6FE}}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}\right)+\ell ^{2}\,\unicode[STIX]{x1D6F1}\,\frac{\unicode[STIX]{x1D707}_{b}}{\dot{\unicode[STIX]{x1D6FE}}}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$O(\ell ^{2})$ terms represent the non-local corrections, and 
$\unicode[STIX]{x1D707}_{b}\equiv (n-1)/(n\unicode[STIX]{x1D707})$ is the bulk plastic modulus. We have made the simplifying assumptions 
$\widehat{Y}\approx 1$ in (2.9a) and 
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}/\dot{\unicode[STIX]{x1D6FE}}\ll 1$ in (2.9c), both of which may be easily relaxed. Equation (2.9) is the main result of this paper, wherein it is evident that non-locality arises in the form of the Laplacian of the deviatoric and volumetric strain rates, and the volume fraction. We first illustrate the key features of the model for a simple flow, before making a critical comparison with previously proposed higher gradient models.3 Model predictions for simple shear
Schematic of plane Couette flow. The top plate is translated with constant velocity 
$u_{w}$ and the bottom plate is held stationary. In the DEM simulations, the granular material was a mixture of spheres of diameters 
$0.9\,d_{p}$, 
$d_{p}$ and 
$1.1\,d_{p}$, of number fractions 0.3, 0.4 and 0.3, respectively, and the walls were constructed of spheres of diameter 
$d_{w}$ in a close-packed square array. Periodic boundary conditions in the 
$x$ and 
$z$ directions were enforced with unit cell dimensions 
$L=60\,d_{p}$, 
$D=40\,d_{p}$.
To evaluate the predictions of the model, we consider shear between plane parallel plates with and without gravity, where the bottom plate is held stationary and top plate moved at constant velocity 
$u_{w}$ in the 
$x$ direction (figure 1). For steady fully developed flow, the relevant momentum balances to be satisfied are 
$\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{yy}/\unicode[STIX]{x2202}y+\unicode[STIX]{x1D70C}_{p}\unicode[STIX]{x1D719}g=0$ and 
$\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{yx}/\unicode[STIX]{x2202}y=0$, where 
$\unicode[STIX]{x1D70C}_{p}$ is the intrinsic density of the particles. Substituting the constitutive equation for the stress (2.9) in the balances, we obtain the equations that govern 
$u_{x}(y)$ and 
$\unicode[STIX]{x1D719}(y)$,
$$\begin{eqnarray}\frac{\displaystyle \text{d}}{\displaystyle \text{d}y}\left(\frac{\displaystyle \unicode[STIX]{x1D707}\,p_{c}}{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}\frac{\displaystyle \text{d}u_{x}}{\displaystyle \text{d}y}-\ell ^{2}\,\frac{\displaystyle \unicode[STIX]{x1D707}\,\unicode[STIX]{x1D6F1}}{\displaystyle \dot{\unicode[STIX]{x1D6FE}}}\frac{\displaystyle \text{d}^{3}u_{x}}{\displaystyle \text{d}y^{3}}\right)=0,\quad \frac{\displaystyle \text{d}p_{c}}{\displaystyle \text{d}y}=\unicode[STIX]{x1D70C}_{p}\unicode[STIX]{x1D719}g.\end{eqnarray}$$ The equations are of order four in 
$u_{x}$ and three in 
$\unicode[STIX]{x1D719}$, requiring as many boundary conditions. Boundary conditions are an aspect that has been quite neglected in the study in granular flows – for instance, the majority of studies assume no slip at solid walls, despite ample evidence to the contrary (Natarajan et al. Reference Natarajan, Hunt and Taylor1995; Ananda et al. Reference Ananda, Moka and Nott2008). We impose the slip boundary condition proposed by Mohan et al. (Reference Mohan, Rao and Nott2002), wherein the slip velocity is proportional to the mean particle spin, which in a classical continuum is half the vorticity. In addition, we employ the wall friction boundary condition that is often used in granular plasticity (Rao & Nott Reference Rao and Nott2008). For the problem under consideration, these conditions take the form
$$\begin{eqnarray}u_{x}-u_{w}=n_{y}\,Kd_{p}\frac{\displaystyle \text{d}u_{x}}{\displaystyle \text{d}y},\quad \frac{\unicode[STIX]{x1D70E}_{yx}}{\unicode[STIX]{x1D70E}_{yy}}=\unicode[STIX]{x1D707}_{w},\end{eqnarray}$$ where 
$u_{w}$ is the wall velocity, 
$n_{y}$ is the 
$y$-component of unit normal at the wall pointing into the granular medium, 
$K$ is the slip coefficient, 
$d_{p}$ the mean particle diameter, and 
$\unicode[STIX]{x1D707}_{w}$ is the wall friction coefficient. The above conditions only apply when the material adjacent to the wall is in plastic deformation; in the case of plane shear under gravity, the material does not deform at the bottom wall, whence we use the following boundary conditions (Mohan et al. Reference Mohan, Rao and Nott2002) in place of (3.2):
$$\begin{eqnarray}u_{x}=0,\quad \frac{\displaystyle \text{d}u_{x}}{\displaystyle \text{d}y}=0.\end{eqnarray}$$ While the above kinematic and stress boundary conditions are supported by plausible physical arguments and some experimental validation (Mohan et al. Reference Mohan, Rao and Nott2002; Ananda et al. Reference Ananda, Moka and Nott2008), very little experimental data exist to guide the choice of the boundary conditions for 
$\unicode[STIX]{x1D719}$. One clear condition is that either the mean volume fraction or the normal stress at a boundary must be specified; the latter is typically how experiments are performed. We assume for simplicity that 
$\unicode[STIX]{x1D719}$ is specified at the boundary where the medium is shearing, and 
$\text{d}\unicode[STIX]{x1D719}/\text{d}y=0$ at the boundary where the medium is not shearing, whence the boundary conditions for 
$\unicode[STIX]{x1D719}$ are
$$\begin{eqnarray}\frac{1}{W}\int _{0}^{W}\unicode[STIX]{x1D719}(y)\,\text{d}y=\overline{\unicode[STIX]{x1D719}}\quad \text{or}\quad \unicode[STIX]{x1D70E}_{yy}(0)=N,\quad \unicode[STIX]{x1D719}(0)=\unicode[STIX]{x1D719}_{0},\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D719}(W)=\unicode[STIX]{x1D719}_{1}\quad \text{or}\quad \frac{\displaystyle \text{d}\unicode[STIX]{x1D719}}{\displaystyle \text{d}y}(W)=0.\end{eqnarray}$$ Equations (3.1)–(3.4) form a well-posed boundary value problem, and their numerical solution is readily obtained if the form of 
$\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D719})$ is known. Data for the critical state pressure are scarce, especially at low loads. We use the simple form proposed by Nott & Jackson (Reference Nott and Jackson1992),
$$\begin{eqnarray}\unicode[STIX]{x1D6F1}=\unicode[STIX]{x1D6FC}\frac{(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{min})^{2}}{(\unicode[STIX]{x1D719}_{max}-\unicode[STIX]{x1D719})^{5}}\quad \unicode[STIX]{x1D719}\geqslant \unicode[STIX]{x1D719}_{min},\quad \unicode[STIX]{x1D6F1}=0\quad \unicode[STIX]{x1D719}<\unicode[STIX]{x1D719}_{min},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D719}_{min}$ and 
$\unicode[STIX]{x1D719}_{max}$ are the volume fractions corresponding to loose and dense random packing, respectively, and 
$\unicode[STIX]{x1D6FC}$ is a material constant. The form for 
$\unicode[STIX]{x1D6F1}$ in (3.5) has the desired features that it vanishes as 
$\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{min}$, and diverges as 
$\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{max}$. We set it to 
$4\times 10^{-7}\unicode[STIX]{x1D70C}_{p}gd_{p}$, which yields a volume fraction of approximately 0.617 for a pressure of 
$15\,\unicode[STIX]{x1D70C}_{p}gd_{p}$, as observed in our previous study (Krishnaraj & Nott Reference Krishnaraj and Nott2016).
The model predictions are compared with particle dynamics simulations using the discrete element method (DEM), a widely used computational tool for granular materials, wherein the motion of each particle is tracked in time, with particle interactions assumed to be elastoplastic and frictional. The interaction parameters used for the DEM simulations are the same as in Krishnaraj & Nott (Reference Krishnaraj and Nott2016). The simulations were performed with spheres of slight size dispersity to avoid crystalline order, and the walls were constructed of particles of diameter 
$d_{w}$ in a close-packed square array (see figure 1). To assess the effect of the geometric wall roughness on the velocity and volume fraction profiles, 
$d_{w}$ was varied. In all the simulations, the inertia number 
$I\equiv (\unicode[STIX]{x1D70C}_{p}d_{p}^{2}\dot{\unicode[STIX]{x1D6FE}}^{2}/p)^{1/2}$, a dimensionless measure of particle inertia, was less than 
$10^{-3}$, to ensure that the flow is in the quasistatic regime. The profiles of 
$u_{x}$ and 
$\unicode[STIX]{x1D719}$ were obtained by averaging spatially over bins of width 
$\unicode[STIX]{x0394}y=2\,d_{p}$ and over time corresponding to 10 strain units at steady state, and assigning the values for each bin to its centre.
Profiles of the volume fraction (a) and velocity (b) for plane shear in the absence of gravity for a Couette gap 
$W=40d_{p}$ and mean volume fraction 
$\overline{\unicode[STIX]{x1D719}}=0.585$. The dashed lines are the results of DEM simulations using wall particles (see figure 1) of two sizes. The solid lines are the model predictions with 
$K=1.65$ and 
$\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}=0.33$ (fitting 
$d_{w}=d_{p}$) and 
$K=0.95$ and 
$\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}=0.1$ (fitting 
$d_{w}=\frac{1}{2}d_{p}$). In panel (b) the (red) dash-dot line is the model prediction for 
$d_{w}=\frac{1}{2}d_{p}$ when 
$\unicode[STIX]{x1D719}(y)$ is assumed to be constant, and the (black) dotted line is the prediction for the ‘fully rough’ wall 
$\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}=1$.
Figure 2 shows the velocity and volume fraction profiles in plane shear in the absence of gravity, compared with DEM simulations for the same. The model predictions are obtained by solving (3.1) with boundary conditions (3.2) at both walls, and (3.4a,c,d). The value of 
$\unicode[STIX]{x1D719}_{0}=\unicode[STIX]{x1D719}_{1}$ is set to that observed in the simulations, the slip coefficient 
$K$ is adjusted to match the observed slip at the walls, and the ratio 
$\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}$ chosen to achieve a fit for shear in the absence of gravity. Apart from these, no attempt is made to fit the parameters; importantly, the mesoscopic length 
$\ell$ was set to 
$10\,d_{p}$, as in the Cosserat plasticity model of Mohan et al. (Reference Mohan, Nott and Rao1999, Reference Mohan, Rao and Nott2002). The expected symmetry of 
$\unicode[STIX]{x1D719}$ and antisymmetry of 
$u_{x}$ (about 
$y=W/2$) are seen in the simulations. The agreement between the model predictions and the DEM simulations is remarkably good, but it is not an aspect we wish to emphasize – the model predictions are sensitive to the form of 
$\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D719})$, and our choice is based on scanty data (Nott & Jackson Reference Nott and Jackson1992). (They are, however, independent of the value of 
$\unicode[STIX]{x1D6FC}$, as it factors out of the momentum balances.) Our principal aim here is to display the qualitative features of the model, and emphasize the strong coupling between the velocity and volume fraction fields – a greater variation of 
$\unicode[STIX]{x1D719}$ leads to a sharper variation of the velocity (and thence the shear rate). If we assume the volume fraction to be constant (as in all prior studies), the inhomogeneity in the shear rate is significantly lower (red dash-dot curve in figure 2b). The lower values for 
$K$ and 
$\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}$ for the smoother walls are in agreement with intuitive expectations and experimental data (Nedderman & Laohakul Reference Nedderman and Laohakul1980; Ananda et al. Reference Ananda, Moka and Nott2008) that the velocity slip will increase, and wall friction decrease, as the geometric roughness of the wall decreases.
When 
$\unicode[STIX]{x1D719}$ does not vary with 
$y$, equations (3.1) and (3.2) yield the analytical solution 
$u_{x}/u_{w}=\frac{1}{2}+(1/2\unicode[STIX]{x1D6FD})\sinh [\unicode[STIX]{x1D701}(W/2-y)]$, where 
$\unicode[STIX]{x1D701}\equiv (1-\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707})^{1/2}/\ell$ and 
$\unicode[STIX]{x1D6FD}\equiv Kd_{p}\unicode[STIX]{x1D701}\cosh (W\unicode[STIX]{x1D701}/2)+\sinh (W\unicode[STIX]{x1D701}/2)$, which in the limit of large 
$W$ tends to the familiar exponential decay 
$\frac{1}{2}[1+\text{e}^{-\unicode[STIX]{x1D701}y}/(1+Kd_{p}\unicode[STIX]{x1D701})]$. The decay length of the velocity is 
$\ell \unicode[STIX]{x1D707}^{1/2}/(\unicode[STIX]{x1D707}-\unicode[STIX]{x1D707}_{w})^{1/2}$, while in the Cosserat model of Mohan et al. (Reference Mohan, Rao and Nott2002) it was shown to be 
$\ell \unicode[STIX]{x1D707}_{w}/(\unicode[STIX]{x1D707}^{2}-\unicode[STIX]{x1D707}_{w}^{2})^{1/2}$. In both the models, inhomogeneity in the shear rate is caused by the difference between the coefficients of material and wall friction.
Next we consider shear in the presence of gravity, where the normal stress at the top wall is fixed at 
$N=1\,\unicode[STIX]{x1D70C}_{p}gW$. As 
$\unicode[STIX]{x1D70E}_{yy}$ rises with distance from the top wall while 
$\unicode[STIX]{x1D70E}_{yx}$ remains constant, it is expected that the medium will not shear beyond a certain distance from the top wall. However, we are unable to specify a priori the boundary between the plastically deforming and undeforming regions – determining this boundary is indeed an important open problem in granular mechanics. A limiting case is to assume that plastic deformation occurs everywhere, as all previous theoretical studies have done (Mohan et al. Reference Mohan, Rao and Nott2002; Pouliquen & Forterre Reference Pouliquen and Forterre2009; Henann & Kamrin Reference Henann and Kamrin2013): we therefore impose boundary conditions (3.2) and (3.4b,c) at the top wall, and (3.3) and (3.4e) at the bottom wall. The results of the DEM simulation for a Couette gap of 
$W=23\,d_{p}$ are shown in figure 3 along with model predictions. In the latter, the values of 
$\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}$ and 
$K$ used are the fits obtained for shear in the absence of gravity (see caption of figure 2). While the model qualitatively captures the features observed in the simulations, namely dilation near the top wall, a constant 
$\unicode[STIX]{x1D719}$ asymptote at large 
$y$, and a roughly exponential decay in the velocity, there is a significant quantitative discrepancy in the velocity profile. In the DEM simulation, a significant fraction of the Couette gap does not shear; thus, not unexpectedly, the assumption of plastic deformation everywhere is not a good one. If we impose the lower boundary conditions at 
$y=\ell$ instead, the agreement is seen to be much better. The inset in figure 3(b) shows the velocity profiles obtained from DEM simulations for 
$W=23$ and 
$40\,d_{p}$ (for the same imposed normal stress 
$N$) – apart from the slight difference in the wall slip, the two velocity profiles are virtually identical. Nevertheless, further studies are needed to formulate a condition based on physical principles that determines the interface between plastically deforming and undeforming regions.
The variation of the volume fraction across the shearing layer in figures 2 and 3 might appear to be small, but it is well known that small changes in 
$\unicode[STIX]{x1D719}$ can cause large changes in the kinematics and stress – indeed, 
$\unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D719})$ varies sharply with 
$\unicode[STIX]{x1D719}$ and diverges as 
$\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{max}$. We have considered only steady plane shear flow to illustrate the key features of our model, and demonstrate the strong coupling between the density and velocity fields. We plan to conduct shortly a more detailed study that considers the unsteady development of the dilation in plane shear, oscillatory shear (Sun & Sundaresan Reference Sun and Sundaresan2011), and the dilation-driven secondary flow in a cylindrical Couette device (Krishnaraj & Nott Reference Krishnaraj and Nott2016).
Profiles of the volume fraction (a) and velocity (b) for plane shear in the presence of gravity for a Couette gap 
$W=23d_{p}$ and normal stress on the top plate 
$N=1\,\unicode[STIX]{x1D70C}_{p}gW$. The dashed lines are the results of DEM simulations with 
$d_{w}=d_{p}$, and the solid lines are the model predictions. The (red) dash-dot lines are the model predictions when the lower boundary condition is imposed at 
$y=\ell$. The inset in (b) shows the velocity profile obtained from simulations for 
$W=23$ and 
$40\,d_{p}$ for an imposed normal stress of 
$N=23\,\unicode[STIX]{x1D70C}_{p}gd_{p}$.
4 Discussion
We now discuss our model with reference to other higher gradient or non-local models proposed in recent years (Mohan et al. Reference Mohan, Rao and Nott2002; Pouliquen & Forterre Reference Pouliquen and Forterre2009; Bouzid et al. Reference Bouzid, Trulsson, Claudin, Clément and Andreotti2013; Henann & Kamrin Reference Henann and Kamrin2013). With the exception of the model of Pouliquen & Forterre (Reference Pouliquen and Forterre2009), the others are not non-local in the conventional sense (Eringen Reference Eringen1983), but we refer to them as such as they all allow plastic deformation even when the yield condition is locally not satisfied. The major difference between these models and ours is, of course, that they assume the granular medium to be incompressible, while we allow density variation and capture dilation at constant pressure. But even apart from dilation, the physical arguments invoked in these models differ in significant ways from ours. Mohan et al. (Reference Mohan, Nott and Rao1999, Reference Mohan, Rao and Nott2002) treated the granular medium as a Cosserat continuum, allowing asymmetry of the Cauchy stress and the presence of a couple stress. While there is some evidence of a kinematic Cosserat effect (Ananda et al. Reference Ananda, Moka and Nott2008), no definitive measurement of the asymmetry of the Cauchy stress or the couple stress is available. In the model presented here (2.9), the granular medium is a classical continuum in which the stress is symmetric. However, we have usefully adopted the stress and slip boundary condition proposed by Mohan et al. Pouliquen & Forterre (Reference Pouliquen and Forterre2009) proposed a non-local model for unidirectional simple shear, having the form of an integral equation for the deformation rate, as we do in (2.5). Their model is based on probabilistic arguments of plastic ‘jumps’ caused by stress fluctuations and their spatial correlation, for which they use arbitrary (but plausible) forms. Equation (2.5) in this paper is the continuum analogue of this idea in tensorial form, and requires no phenomenological functions or assumptions.
Two recent higher gradient models have been proposed as extensions of the 
$\unicode[STIX]{x1D707}(I)$ model (GDR MiDi 2004; Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006). Henann & Kamrin (Reference Henann and Kamrin2013, Reference Henann and Kamrin2014) coupled the equations of motion (with (2.4a) for the stress) with a diffusion equation for the fluidity 
$\dot{\unicode[STIX]{x1D706}}$ that was proposed earlier by Bocquet, Colin & Ajdari (Reference Bocquet, Colin and Ajdari2009), based on several assumptions and simplifications of the statistics and spatial interaction of plastic events. The appropriateness of a hydrodynamic equation governing 
$\dot{\unicode[STIX]{x1D706}}$ is open to question, as it is neither a conserved quantity nor a broken symmetry order parameter. Treating the model as a phenomenological one, it is not clear what the appropriate boundary conditions for fluidity are – the authors assume zero flux at solid boundaries, with the justification that it corresponds to ‘minimal microscopic wall effects’, but it is clear from our DEM simulations that the effects of the wall are far from minimal (figure 2), and that the velocity and solids fraction profiles are sensitive to the topological and frictional roughness of the boundaries; a boundary condition that is independent of these features is therefore unrealistic. Bouzid et al. (Reference Bouzid, Trulsson, Claudin, Clément and Andreotti2013, Reference Bouzid, Izzet, Trulsson, Clément, Claudin and Andreotti2015) argued that the yield function 
$Y$ in (2.1) has a non-local correction of inertial origin; by perturbing about the state 
$I=0$, and assuming that the correction is finite as 
$I\rightarrow 0$, they reasoned that 
$Y=\unicode[STIX]{x1D707}(I)\,p\,(1-\ell ^{2}\unicode[STIX]{x1D6FB}^{2}I/I)$ is the corrected yield function to be used in (2.4). Conceptually it is not clear how yielding can be modulated by inertia – indeed, this is contrary to its definition as the condition to be satisfied for plastic deformation to occur – and why the perturbation should be singular at 
$I=0$. Our model introduces no additional variables and equations, and no phenomenological assumptions other than the flow rule being non-local (2.5). It is a much simpler embodiment of the argument of Pouliquen & Forterre (Reference Pouliquen and Forterre2009) and Bocquet et al. (Reference Bocquet, Colin and Ajdari2009) that plastic events are spatially communicated. It requires boundary conditions only on the primary hydrodynamic variables 
$\boldsymbol{u}$ and 
$\unicode[STIX]{x1D719}$, which are directly measurable in experiments and particle simulations. We note, however, that our relation for the deviatoric stress 
$\unicode[STIX]{x1D748}^{\prime }$ (2.9a) coincides with that of Bouzid et al. (Reference Bouzid, Trulsson, Claudin, Clément and Andreotti2013) for incompressible simple shear.
A brief note on the Hadamard well-posedness of our model is in order. It was shown by Barker et al. (Reference Barker, Schaeer, Bohorquez and Gray2015) that the incompressible 
$\unicode[STIX]{x1D707}(I)$ model is Hadamard ill-posed. The inclusion of compressibility partially alleviates the problem (Heyman et al. Reference Heyman, Delannay, Tabuteau and Valance2017), but the model remains ill-posed in the limit 
$I\rightarrow 0$. Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) showed that if the yield condition is allowed to depend on 
$I$, the model is well-posed under certain conditions – however, as discussed above, dependence of the yield condition on 
$I$ is conceptually problematic. Goddard & Lee (Reference Goddard and Lee2017) showed that introducing a (rate-dependent) correction to the stress of the form 
$-\unicode[STIX]{x1D712}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D735}\boldsymbol{u}$, where 
$\unicode[STIX]{x1D712}$ is a constant, repairs the ill-posedness of the incompressible 
$\unicode[STIX]{x1D707}(I)$ model. The non-local term in (2.9a) is the generalization of this correction for compressible flows that also preserves rate independence of the stress. It is straightforward to show that our model (2.9) is unconditionally well-posed, with or without compressibility. Finally, we note that while our model does not ensure positive dissipation locally, the dissipation is positive in the non-local sense, i.e. 
$(\unicode[STIX]{x1D748}^{\prime (0)}+\ell ^{2}\unicode[STIX]{x1D748}^{\prime (1)})\boldsymbol{ : }(\unicode[STIX]{x1D63F}^{\prime }+\ell ^{2}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}^{\prime })>0$ to 
$O(\ell ^{2})$, in accord with the assumptions the model is based on.
5 Conclusion
We have derived a non-local model for slow granular flows with the only assumption that the flow rule is valid not locally, but in a volume-averaged sense. No other phenomenological assumption is invoked. This simple and physically intuitive idea leads to a constitutive relation for the stress that captures shear localization and dilatancy, and requires boundary conditions only on the primary hydrodynamic variables. For exposition of the key features of the model, we have compared its predictions with DEM simulations for steady simple shear with and without gravity; more complex flows, such as the emergence of a dilation-driven secondary flow (Krishnaraj & Nott Reference Krishnaraj and Nott2016), will be considered in a subsequent paper. Our study underscores the importance of formulating physically realistic boundary conditions, an aspect that has received little attention in previous studies, and the importance of measuring the density and velocity fields in experiments and particle simulations. The model we have proposed is for slow flow, or negligibly small inertia number 
$I$; however, it is amenable to extension to finite 
$I$ in the same spirit as the 
$\unicode[STIX]{x1D707}(I)$ model, by treating 
$\unicode[STIX]{x1D707}$ and 
$\unicode[STIX]{x1D707}_{b}$ as functions of 
$I$, and 
$\unicode[STIX]{x1D6F1}$ as a function of 
$\unicode[STIX]{x1D719}$ and 
$I$ in (2.9).
Acknowledgements
We acknowledge useful discussions with K. K. Rao, O. Pouliquen, J. Goddard, M. Bouzid and R. Mari. This work was supported by the Science and Engineering Research Board under grant EMR/2016/002817. A part of the work was conducted while P.R.N. was visiting the Kavli Institute for Theoretical Physics, supported by the National Science Foundation under grant no. NSF PHY-1748958.
Declaration of interests
The authors report no conflict of interest.
