2.1 Constitutive laws: $\unicode[STIX]{x1D707}(I)$ -rheology

Analogous to Coulomb’s friction law, using dimensional arguments, experiments and numerical simulations, Iordanoff & Khonsari (Reference Iordanoff and Khonsari2004), Midi (Reference Midi2004) and Da Cruz et al. (Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005) demonstrated that the shear stress $\unicode[STIX]{x1D70F}$ is proportional to the normal stress $P$ for two-dimensional (2-D) dense granular flows of rigid particles so that $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}(I)P$ where the local friction coefficient $\unicode[STIX]{x1D707}$ is only a function of a dimensionless parameter called the inertial number $I=\dot{\unicode[STIX]{x1D6FE}}d/\sqrt{P/\unicode[STIX]{x1D70C}}$ . Here, $\dot{\unicode[STIX]{x1D6FE}}$ is the local shear rate which is related to the macroscopic time scale of the granular flow and $d\sqrt{\unicode[STIX]{x1D70C}/P}$ is the microscopic time scale corresponding to any local rearrangement of grains of diameter $d$ and density $\unicode[STIX]{x1D70C}$ subjected to a local normal stress $P$ . Note that $I$ is also the square root of the Savage or the Coulomb number as used in Savage (Reference Savage1984) or Ancey, Coussot & Evesque (Reference Ancey, Coussot and Evesque1999), respectively. In general, the dimensionless local friction coefficient is given by (Midi Reference Midi2004; Jop et al. Reference Jop, Forterre and Pouliquen2006)

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}(I)=\unicode[STIX]{x1D707}_{s}+{\displaystyle \frac{\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{s}}{\left(1+{\displaystyle \frac{I_{0}}{I}}\right)}}, & \displaystyle\end{eqnarray}$$

whereby $\unicode[STIX]{x1D707}$ saturates towards two fundamental constants for a dense granular medium $\unicode[STIX]{x1D707}_{s}$ or $\unicode[STIX]{x1D707}_{d}$ depending respectively on the inertial number $I\ll 1$ (quasi-static regime) or $I\gg 1$ (kinetic or gaseous regime). Jop et al. (Reference Jop, Forterre and Pouliquen2006) proposed a 3-D generalisation of this scalar constitutive relation for a granular material by decomposing the Cauchy stress tensor into an isotropic contribution from the local pressure $p$ and a traceless deviatoric stress tensor $\unicode[STIX]{x1D70F}_{ij}$ while assuming that $\unicode[STIX]{x1D70F}_{ij}$ is aligned with the strain rate tensor $\dot{\unicode[STIX]{x1D6FE}}_{ij}=1/2(\unicode[STIX]{x2202}_{i}u_{j}+\unicode[STIX]{x2202}_{j}u_{i})$ (where $u_{i}$ represent components of the velocity field). So, if $\boldsymbol{x}$ is the position vector and $t$ represents time, then the Cauchy stress tensor

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}(\boldsymbol{x},t)=-p(\boldsymbol{x},t)\unicode[STIX]{x1D6FF}_{ij}+\unicode[STIX]{x1D70F}_{ij}(\boldsymbol{x},t), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta and $\unicode[STIX]{x1D70F}_{ij}(\boldsymbol{x},t)=\unicode[STIX]{x1D702}(\boldsymbol{x},t)\dot{\unicode[STIX]{x1D6FE}}_{ij}(\boldsymbol{x},t)$ and the local granular liquid viscosity

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=\frac{\unicode[STIX]{x1D707}p}{|\dot{\unicode[STIX]{x1D6FE}}|}, & \displaystyle\end{eqnarray}$$

is, thereby, a nonlinear function of the local pressure $p$ and the local second invariant of the strain rate tensor $|\dot{\unicode[STIX]{x1D6FE}}|=\sqrt{1/2\dot{\unicode[STIX]{x1D6FE}}_{ij}\dot{\unicode[STIX]{x1D6FE}}_{ij}}$ via the local friction coefficient given by (2.1) and the inertial number for the $3$ -D case

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle I=\frac{|\dot{\unicode[STIX]{x1D6FE}}|d}{\sqrt{p/\unicode[STIX]{x1D70C}}}. & \displaystyle\end{eqnarray}$$

In addition, the solid volume fraction $\unicode[STIX]{x1D719}$ is also a linear function of $I$ (see Andreotti et al. Reference Andreotti, Forterre and Pouliquen2011, p. 238) but, in general, it varies very little and so, for the sake of simplicity, it is taken to be a constant in the following.

2.2 Stokes’ first problem and its steady-state solution

Consider an infinite rigid flat plate placed at rest ( $t=0$ ) on top of a semi-infinite dense granular medium. As illustrated in figure 1, the plate is set in motion impulsively at $t>0$ by applying a tangential shear stress $\unicode[STIX]{x1D70F}_{w}$ (along the $x$ -direction) in the presence of a normal stress $p_{w}$ (along the $y$ -direction). It is then natural to restrict the analysis to two dimensions. In fact, the absence of any horizontal length scale implies that the flow properties should depend only on $y$ and $t$ . Incompressibility and the initial condition then imply that the vertical velocity is uniformly zero for all $t\geqslant 0$ . And so the only non-zero components of the strain rate tensor are $\dot{\unicode[STIX]{x1D6FE}}_{xy}=\dot{\unicode[STIX]{x1D6FE}}_{yx}=(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y)/2$ , where $u$ is the $x$ -component of the velocity field. Therefore, the $2$ -D shear-stress tensor is completely determined by a single scalar shear field $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}p$ so that the $x$ and $y$ momentum equations become, respectively,

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}}=-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\unicode[STIX]{x1D707}p), & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y}}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}g, & \displaystyle\end{eqnarray}$$

where $g$ is the acceleration due to gravity. Initially the granular medium is at rest ( $u(y,t=0)=0$ ) but a static granular material can support a wide variety of shear-stress and pressure distributions so long as the yield stress is not exceeded so that $\unicode[STIX]{x1D70F}\leqslant \unicode[STIX]{x1D707}_{s}p$ . Within the context of the $\unicode[STIX]{x1D707}(I)$ -rheology, the only possible configuration where the granular liquid is at rest corresponds to the case $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}_{s}p$ . For all other static states, the constitutive law $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}p$ is violated, so the $\unicode[STIX]{x1D707}(I)$ -rheology is no longer applicable. Therefore, in the following, the initial conditions correspond to a specific static state wherein the shear stress equals the yield stress throughout the granular medium

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}(y,t=0)=\unicode[STIX]{x1D707}_{s}p, & \displaystyle\end{eqnarray}$$

with $p$ the hydrostatic pressure. Whereas the boundary conditions for all $t>0$ are

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle p(y=0,t)=p_{w}, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}(y=0,t)={\displaystyle \frac{\unicode[STIX]{x1D70F}_{w}}{p_{w}}}=\unicode[STIX]{x1D707}_{w}, & \displaystyle\end{eqnarray}$$

along with the condition that the grains sufficiently far from the plate remain static, i.e. $u(y=\infty ,t)=0$ , and so $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{s}$ at $y=\infty$ . In contrast to the classical Stokes’ problem (Stokes Reference Stokes1851) for a Newtonian fluid, firstly, the frictional force (right-hand side of equation (2.5)) is not only non-uniform due to the hydrostatic pressure $p(y)=p_{w}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}gy$ but also nonlinear since $\unicode[STIX]{x1D707}$ depends on both $p(y)$ and $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y$ via (2.1). And secondly, there exists a non-trivial steady-state solution wherein the shear stress is constant throughout the medium such that $\unicode[STIX]{x1D707}p=\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D707}_{w}p_{w}$ , as already shown in Andreotti et al. (Reference Andreotti, Forterre and Pouliquen2011, pp. $254$ – $256$ ) and Cawthorn (Reference Cawthorn2011, pp. $46$ – $50$ ). Since $\unicode[STIX]{x1D707}\in [\unicode[STIX]{x1D707}_{s},\unicode[STIX]{x1D707}_{d}]$ and $p(y)$ increases linearly with the depth $y$ , it follows that

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}={\displaystyle \frac{\unicode[STIX]{x1D707}_{w}}{1+(\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}gy/p_{w})}}, & \displaystyle\end{eqnarray}$$

for all $y\leqslant \unicode[STIX]{x1D6FF}_{s}$ and $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{s}$ otherwise; here, the critical depth $\unicode[STIX]{x1D6FF}_{s}$ is given by

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}_{s}=\unicode[STIX]{x1D6FD}_{w}\left({\displaystyle \frac{p_{w}}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}g}}\right), & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FD}_{w}=(\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}_{s}-1)$ . This denotes the depth beyond which the granular medium does not flow $y\geqslant \unicode[STIX]{x1D6FF}_{s}$ . Note that the term $\unicode[STIX]{x1D6FD}_{w}$ represents the surplus in wall shear stress over the yield criterion at the wall. The steady-state solution (2.10) can be used to obtain the corresponding velocity profile. As previously shown by Cawthorn (Reference Cawthorn2011, pp. $46$ – $50$ ), the resulting relation between the wall velocity (if the no-slip condition is allowed) and the shear layer thickness compares qualitatively well with the molecular dynamic simulations of Thompson & Grest (Reference Thompson and Grest1991).

2.3 Stokes’ first problem: shear-stress diffusion equation

Most numerical studies obtain the velocity field by solving the momentum equations that account for the constitutive law (2.2). Note that the latter is coupled with the expressions for the local friction coefficient (2.1) and the inertial number $I$ . However, it is possible to write a single equation for the shear stress in the case of Stokes’ first problem within the context of the $\unicode[STIX]{x1D707}(I)$ -rheology. Since $I=-(1/2)(\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y)\sqrt{\unicode[STIX]{x1D70C}d^{2}/p}$ , the local friction coefficient (2.1) can be rewritten as

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}}=-2I_{0}\sqrt{{\displaystyle \frac{p}{\unicode[STIX]{x1D70C}d^{2}}}}\left({\displaystyle \frac{\unicode[STIX]{x1D707}-\unicode[STIX]{x1D707}_{s}}{\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}}}\right), & \displaystyle\end{eqnarray}$$

which when differentiated with respect to time and using $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}p$ yields

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}y}}=-{\displaystyle \frac{2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}I_{0}}{d}}\left(p\sqrt{{\displaystyle \frac{p}{\unicode[STIX]{x1D70C}}}}\right){\displaystyle \frac{1}{\left(\unicode[STIX]{x1D707}_{d}p-\unicode[STIX]{x1D70F}\right)^{2}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}{\unicode[STIX]{x2202}t}}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{s}>0$ . By allowing $\unicode[STIX]{x2202}^{2}u/\unicode[STIX]{x2202}y\unicode[STIX]{x2202}t=\unicode[STIX]{x2202}^{2}u/\unicode[STIX]{x2202}t\unicode[STIX]{x2202}y$ and then introducing (2.5) in the above expression, this leads to a single nonlinear shear-stress diffusion equation

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}{\unicode[STIX]{x2202}t}}=\left({\displaystyle \frac{d}{2\unicode[STIX]{x1D719}\unicode[STIX]{x0394}\unicode[STIX]{x1D707}I_{0}\sqrt{\unicode[STIX]{x1D70C}}}}\right){\displaystyle \frac{\left(\unicode[STIX]{x1D707}_{d}p-\unicode[STIX]{x1D70F}\right)^{2}}{p\sqrt{p}}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D70F}}{\unicode[STIX]{x2202}y^{2}}}. & \displaystyle\end{eqnarray}$$

Finally, by taking $\unicode[STIX]{x1D708}_{g}=(d/2\unicode[STIX]{x1D719}\unicode[STIX]{x0394}\unicode[STIX]{x1D707}I_{0})\sqrt{p_{w}/\unicode[STIX]{x1D70C}}$ as a proper diffusion coefficient and the steady-state shear-layer thickness $\unicode[STIX]{x1D6FF}_{s}$ as the characteristic length scale, the non-dimensional time and space coordinates are $\tilde{t}=\unicode[STIX]{x1D708}_{g}t/\unicode[STIX]{x1D6FF}_{s}^{2}$ and ${\tilde{y}}=y/\unicode[STIX]{x1D6FF}_{s}$ , respectively. Thus, in terms of the normalised pressure $\tilde{p}=p(y)/p_{w}$ and shear stress $\tilde{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D70F}(y,t)/\unicode[STIX]{x1D70F}_{w}$ , the above equation becomes

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70F}}}{\unicode[STIX]{x2202}\tilde{t}}}={\displaystyle \frac{(\unicode[STIX]{x1D707}_{d}\tilde{p}-\unicode[STIX]{x1D707}_{w}\tilde{\unicode[STIX]{x1D70F}})^{2}}{\tilde{p}\sqrt{\tilde{p}}}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D70F}}}{\unicode[STIX]{x2202}{\tilde{y}}^{2}}}, & \displaystyle\end{eqnarray}$$

with boundary conditions $\tilde{\unicode[STIX]{x1D70F}}(0,\tilde{t})=\tilde{\unicode[STIX]{x1D70F}}(1,\tilde{t})=1$ and an initial condition $\tilde{\unicode[STIX]{x1D70F}}({\tilde{y}},0)=\unicode[STIX]{x1D707}_{s}\tilde{p}/\unicode[STIX]{x1D707}_{w}$ . The steady-state solution for the non-dimensional shear stress is $\tilde{\unicode[STIX]{x1D70F}}=1$ . Unlike the steady-state local friction coefficient $\unicode[STIX]{x1D707}$ , the steady-state shear stress is a continuous and infinitely differentiable function for all $y\geqslant 0$ . So, it is expected that $\tilde{\unicode[STIX]{x1D70F}}$ remains smooth also for the unsteady case. In fact, the term $(\unicode[STIX]{x1D707}_{d}\tilde{p}-\unicode[STIX]{x1D707}_{w}\tilde{\unicode[STIX]{x1D70F}})$ is positive–definite. Hence, it is quite straightforward to homogenise the boundary conditions and numerically solve the above equation using a second-order centred finite difference scheme for spatial derivatives and a second-order Crank–Nicolson one for temporal integration. The updated $\tilde{\unicode[STIX]{x1D70F}}$ is obtained by an iterative Richardson minimal residual process.

2.4 Well-posedness of the shear-stress diffusion equation

Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015) demonstrated that the $\unicode[STIX]{x1D707}(I)$ -rheology is well-posed for intermediate values of inertial number $I$ , but that it is ill-posed for both high and low inertial numbers. In the present case, $I\gg 1$ in the neighbourhood of the wall when either the applied wall shear $\unicode[STIX]{x1D70F}_{w}$ is close to the critical wall shear $\unicode[STIX]{x1D70F}_{d}=\unicode[STIX]{x1D707}_{d}p_{w}$ or when wall speed $u_{w}$ is sufficiently large. In addition, there is always a zone where $I\ll 1$ (or $\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D707}_{s}$ ) since the medium is slowly moving or stationary when $y\sim \unicode[STIX]{x1D6FF}_{s}$ . As already shown by Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015) this should provoke Hadamard instability (Joseph & Saut Reference Joseph and Saut1990) whereby infinitesimally small short-wave perturbations are amplified indefinitely. Thus, numerical solutions may not converge as the grid is refined and so, before proceeding any further, it is important to verify that the shear-stress diffusion (2.15) is indeed well-posed for all values of the inertial number $I$ .

In fact, this nonlinear diffusion equation (2.15) is one-dimensional and when it is linearised about an arbitrary base state, as in Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015, see p. 799), the resulting dispersion equation of the normal mode analysis in the high-wavenumber limit should be

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}=-\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D709}_{y}^{2}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}$ is the complex frequency, $\unicode[STIX]{x1D709}_{y}$ is the wavenumber in the ${\tilde{y}}$ -direction and $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}(\tilde{p}^{0},\tilde{\unicode[STIX]{x1D70F}}^{0})$ is some function of the base state pressure $\tilde{p}^{0}$ and shear stress $\tilde{\unicode[STIX]{x1D70F}}^{0}$ . Since the real part of the complex frequency $\unicode[STIX]{x1D706}$ provides the perturbation growth rate, it is straightforward to see that (2.15) is stable for short waves along the $y$ -coordinate.

Note that this is a one-dimensional analysis since it takes into account only plane waves along the ${\tilde{y}}$ -axis in order to analyse the well-posedness of (2.15) for one-dimensional, time-dependent computations. For two-dimensional codes that consider the granular Stokes’ problem as a test case, the situation is more complex. In this case, as previously shown by Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015), the equations are still ill-posed since oblique two-dimensional short waves are unstable in the region $y\sim \unicode[STIX]{x1D6FF}_{s}$ (or close to the wall for $\unicode[STIX]{x1D707}_{w}\sim \unicode[STIX]{x1D707}_{d}$ ).

3.1 Self-similarity at $\tilde{t}\ll 1$

When $\tilde{t}\ll 1$ , the non-dimensional shear layer thickness $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{s}$ is small, as observed in figure 2 (bottom). By taking the non-dimensional local pressure $\tilde{p}=1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}$ and local shear stress $\tilde{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D707}(1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}})/\unicode[STIX]{x1D707}_{w}$ , the diffusion equation (2.15) becomes

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}\tilde{t}}}=(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707})^{2}\left({\displaystyle \frac{2\unicode[STIX]{x1D6FD}_{w}}{\sqrt{1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}{\tilde{y}}}}+\sqrt{1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}{\tilde{y}}^{2}}}\right), & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D707}({\tilde{y}},0)=\unicode[STIX]{x1D707}_{s}$ , $\unicode[STIX]{x1D707}(0,\tilde{t})=\unicode[STIX]{x1D707}_{w}$ and $\unicode[STIX]{x1D707}(1,\tilde{t})=\unicode[STIX]{x1D707}_{s}$ . Note that, in general, $\unicode[STIX]{x1D6FD}_{w}<1$ and hence, for $\tilde{t}\ll 1$ , the highest-order derivative of $\unicode[STIX]{x1D707}$ should dominate if ${\tilde{y}}\leqslant \unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{s}$ so that $\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D707}_{w}-(\unicode[STIX]{x1D707}_{w}-\unicode[STIX]{x1D707}_{s})y/\unicode[STIX]{x1D6FF}$ . And in the outer region, $\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D707}_{s}$ . As the spatial variations of the local friction coefficient $\unicode[STIX]{x1D707}$ are stronger inside the shear layer (i.e. when $\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}\ll 1$ ), it is reasonable to simplify (3.1) to

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}\tilde{t}}}=\unicode[STIX]{x0394}\unicode[STIX]{x1D707}^{2}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}{\tilde{y}}^{2}}}, & \displaystyle\end{eqnarray}$$

at the leading order with the same boundary conditions as before. By taking $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D707}-(\unicode[STIX]{x1D707}_{w}-(\unicode[STIX]{x1D707}_{w}-\unicode[STIX]{x1D707}_{s}){\tilde{y}})$ , the above equation admits a self-similar solution for $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ with $\unicode[STIX]{x1D702}={\tilde{y}}/(2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}})$ satisfying the initial and the boundary conditions for all $\tilde{t}\geqslant 0$ . Thereby, the local friction coefficient is deduced to be

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D707}_{s}+(\unicode[STIX]{x1D707}_{w}-\unicode[STIX]{x1D707}_{s})(1-{\tilde{y}})\text{ erfc }\left({\displaystyle \frac{{\tilde{y}}}{2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}}}\right), & \displaystyle\end{eqnarray}$$

which implies that the shear layer thickness should grow as $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{s}\sim 4.66\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}$ .

The expression (3.3) is presented in figure 3 (top) where the time evolution of the local friction coefficient is displayed as function of ${\tilde{y}}$ at different times $\tilde{t}=10^{-4},10^{-3},10^{-2},10^{-1},10^{0},10^{1}$ for the same values of $(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{w})/\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ as in figure 2. When compared with numerical solutions of $\unicode[STIX]{x1D707}$ (and also, $\unicode[STIX]{x1D6FF}$ ) as seen in figure 3 (top) (and figure 2, respectively), these approximations are very satisfactory for all $\tilde{t}\ll 1$ . Indeed, the expression (3.3) for the local friction coefficient, and especially the estimations of the time evolution of the shear layer thickness $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{s}\sim 4.66\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}$ , are reasonably good even when $\tilde{t}$ is of order $1$ .

Figure 3.

Comparison between self-similar approximation (top) and the long-time diffusion approximation (bottom) with the direct numerical solution (continuous lines) for various normalised wall friction coefficients $(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{w})/(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{s})$ (a) $0.9$ , (b) $0.5$ and (c) $0.1$ .

3.2 Diffusion at $\tilde{t}\gg 1$

As soon as $\tilde{t}\sim O(1)$ the non-dimensional shear layer thickness is no longer small, and hence a singular perturbation of (3.1) cannot be obtained with the present scaling for the $y$ -coordinate. However, by using the non-dimensional pressure $\tilde{p}$ as an equivalent normalised spatial variable ${\hat{y}}=1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}$ , it is possible to show that

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70F}}}{\unicode[STIX]{x2202}\tilde{t}}}=(\unicode[STIX]{x1D6FD}_{w}(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}))^{2}\sqrt{{\hat{y}}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D70F}}}{\unicode[STIX]{x2202}{\hat{y}}^{2}}}, & \displaystyle\end{eqnarray}$$

which is singular if $(\unicode[STIX]{x1D6FD}_{w}\unicode[STIX]{x0394}\unicode[STIX]{x1D707})^{2}$ tends to zero. This is often true since $(\unicode[STIX]{x1D6FD}_{w}(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}))^{2}$ is of the order of $(\unicode[STIX]{x1D6FD}_{w}\unicode[STIX]{x0394}\unicode[STIX]{x1D707})^{2}$ and ${\hat{y}}\sim O(1)$ . Thus, taking $\unicode[STIX]{x1D6FD}_{w}(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707})\sim \unicode[STIX]{x1D6FD}_{w}\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ at the leading order in (3.4), it becomes linear and admits a WKB approximation in the ${\hat{y}}$ -coordinate. Thereby, the local shear stress $\tilde{\unicode[STIX]{x1D70F}}$ can be shown to be

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D70F}}\sim 1+{\hat{y}}^{-1/8}\mathop{\sum }_{m=1}^{\infty }\left[\unicode[STIX]{x1D6EC}_{m}\exp (-\unicode[STIX]{x1D706}_{m}^{2}\tilde{t})\sin \left(m\unicode[STIX]{x03C0}{\displaystyle \frac{{\hat{y}}^{3/4}-1}{(\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}_{s})^{3/4}-1}}\right)\right], & \displaystyle\end{eqnarray}$$

where

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6EC}_{m}=2\int _{1}^{\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}_{s}}{\hat{y}}^{1/8}\left({\displaystyle \frac{\unicode[STIX]{x1D707}_{s}}{\unicode[STIX]{x1D707}_{w}}}{\hat{y}}-1\right)\sin \left(m\unicode[STIX]{x03C0}{\displaystyle \frac{{\hat{y}}^{3/4}-1}{(\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}_{s})^{3/4}-1}}\right)\,\text{d}{\hat{y}}, & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}_{m}={\displaystyle \frac{3}{4}}{\displaystyle \frac{m\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FD}_{w}\unicode[STIX]{x0394}\unicode[STIX]{x1D707}}{(\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}_{s})^{3/4}-1}}. & \displaystyle\end{eqnarray}$$

Note that it is possible to solve the linearised version of (3.4) directly by separation of variables as well. In that case, the general solution is an infinite sum of Bessel functions in the ${\hat{y}}$ -direction. But it is advantageous to work with the approximate solution (3.5) when $m$ is large. Nevertheless, the comparison between the WKB approximation (3.5) (dot-dash lines in figure 3) and the numerical solution for the local friction coefficient is good. In particular, the agreement is excellent when the applied friction coefficient $\unicode[STIX]{x1D707}_{w}$ is close to yield friction coefficient $\unicode[STIX]{x1D707}_{s}$ .

3.3 Velocity field development

It is now possible to compute the temporal evolution of the velocity field using local shear stress $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D707}p$ . By rewriting (2.5) in terms of $\tilde{t}$ and ${\tilde{y}}$ , a natural normalisation for the velocity $u$ can be shown to be

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle u=\tilde{u} \left({\displaystyle \frac{g\unicode[STIX]{x1D6FF}_{s}^{2}}{\unicode[STIX]{x1D6FD}_{w}\unicode[STIX]{x1D708}_{g}}}\right), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{s}$ is the steady-state shear layer thickness (2.11), $\unicode[STIX]{x1D6FD}_{w}=(\unicode[STIX]{x1D707}_{w}-\unicode[STIX]{x1D707}_{s})/\unicode[STIX]{x1D707}_{s}$ and $\unicode[STIX]{x1D708}_{g}=(d/\unicode[STIX]{x1D719}\unicode[STIX]{x0394}\unicode[STIX]{x1D707}I_{0})\sqrt{p_{w}/\unicode[STIX]{x1D70C}}$ is the diffusion coefficient which appears in the granular Stokes’ equation (2.14). Using this normalisation and by taking

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D707}}={\displaystyle \frac{\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D707}}}, & \displaystyle\end{eqnarray}$$

the equation for the shear rate (2.12) becomes

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{u} }{\unicode[STIX]{x2202}{\tilde{y}}}}=-{\displaystyle \frac{\sqrt{\tilde{p}}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D707}}}\left({\displaystyle \frac{1}{\tilde{\unicode[STIX]{x1D707}}}}-1\right), & \displaystyle\end{eqnarray}$$

which can then be integrated to study the velocity field. In the following, the no-slip condition is assumed so that the wall speed and the velocity of the granular medium right next to it ( ${\tilde{y}}=0$ ) are equal.

Figure 4.

Normalised wall velocity as a function of the non-dimensional time $\tilde{t}$ for various applied wall shear stresses: $(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{w})/(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{s})$ (a) $0.9$ , (b) $0.5$ and (c) $0.1$ . The different symbols denote calculations from the numerical solution using (2.15) (filled circles), the self-similar solution (3.3) (▫) and the long-time approximate solution (3.5) ( $+$ ).

Figure 4 compares the normalised wall velocity $\tilde{u} _{w}=\tilde{u} ({\tilde{y}}=0,\tilde{t})$ when obtained from the numerical solution of (2.15) (filled circles), the self-similar solution (3.3) (▫) and the long-time approximate solution (3.5) ( $+$ ) for a given boundary condition at the wall. Each panel shows the temporal evolution of $\tilde{u} _{w}$ when the wall friction coefficient is taken as $(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{w})/(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{s})$ (a) $0.9$ , (b) $0.5$ and (c) $0.1$ . In all cases, the wall velocity grows monotonically by following a power law in $\tilde{t}$ as long as $\tilde{t}\leqslant 1$ and then it saturates at the steady-state value $\tilde{u} _{w}^{\infty }=\tilde{u} ({\tilde{y}}=0,\tilde{t}=\infty )$ . Firstly, the self-similar solution (3.3) gives good qualitative agreement with the numerical results, and also, it matches well with the power law $\tilde{u} _{w}(\tilde{t})\sim \sqrt{\tilde{t}}$ . The long-time WKB approximation (3.5) matches very well with the numerical solution at all times $\tilde{t}>1$ . Secondly, it is observed that the steady-state wall velocity $\tilde{u} _{w}^{\infty }$ depends largely on the imposed wall shear via $\tilde{\unicode[STIX]{x1D707}}_{w}=(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{w})/\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ as already shown by previous works (Andreotti et al. Reference Andreotti, Forterre and Pouliquen2011; Cawthorn Reference Cawthorn2011). In fact, for the steady-state solution, Cawthorn (Reference Cawthorn2011, p. 49) had already given an asymptotic solution for the wall speed $\tilde{u} _{w}^{\infty }$ when the applied shear stress is just above the yield shear, i.e. when $\unicode[STIX]{x1D707}({\tilde{y}}=0,\tilde{t})$ approaches $\unicode[STIX]{x1D707}_{s}$ from above (or $\tilde{\unicode[STIX]{x1D707}}_{w}\sim 1$ ). However, it is also possible to obtain expressions for $\tilde{u} _{w}^{\infty }$ for a wide range of $\tilde{\unicode[STIX]{x1D707}}_{w}$ via the self-similar solution (3.3) since, as suggested by figure 4, this gives a good approximation to the steady-state wall speed. In the limit when $\tilde{t}\gg 1$ , the approximate solution (3.3) becomes a function only of ${\tilde{y}}$ , and so, in terms of the normalised friction coefficient (3.9), it is given by $\tilde{\unicode[STIX]{x1D707}}\sim \tilde{\unicode[STIX]{x1D707}}_{w}+(1-\tilde{\unicode[STIX]{x1D707}}_{w}){\tilde{y}}$ . Using this expression in (3.10), it reads

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{w}^{\infty }\sim \int _{0}^{1}{\displaystyle \frac{\sqrt{1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D707}}}\left({\displaystyle \frac{1}{\tilde{\unicode[STIX]{x1D707}}_{w}+(1-\tilde{\unicode[STIX]{x1D707}}_{w}){\tilde{y}}}}-1\right)\,\text{d}{\tilde{y}}, & \displaystyle\end{eqnarray}$$

and since $\unicode[STIX]{x1D6FD}_{w}=(\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}_{s}-1)<1$ , it could be further developed to obtain a simple expression for the steady-state wall velocity

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{w}^{\infty }\sim -{\displaystyle \frac{1}{\unicode[STIX]{x0394}\unicode[STIX]{x1D707}}}\left[1+{\displaystyle \frac{\log \tilde{\unicode[STIX]{x1D707}}_{w}}{(1-\tilde{\unicode[STIX]{x1D707}}_{w})}}\right]+O(\unicode[STIX]{x1D6FD}_{w}). & \displaystyle\end{eqnarray}$$

Noting that $\tilde{\unicode[STIX]{x1D707}}_{w}=1-\unicode[STIX]{x1D707}_{s}\unicode[STIX]{x1D6FD}_{w}/\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ , it is straightforward to see that the first term in the above expression cancels out when $\tilde{\unicode[STIX]{x1D707}}_{w}\sim 1$ (or $\unicode[STIX]{x1D6FD}_{w}=(\unicode[STIX]{x1D707}_{w}/\unicode[STIX]{x1D707}_{s}-1)\ll 1$ ), and thereby, it gives $\tilde{u} _{w}^{\infty }\sim O(\unicode[STIX]{x1D6FD}_{w})$ . In this case, as already deduced by Cawthorn (Reference Cawthorn2011, p. 49), the above integral leads to

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{w}^{\infty }\sim {\displaystyle \frac{\unicode[STIX]{x1D707}_{s}}{2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}^{2}}}\left(\unicode[STIX]{x1D6FD}_{w}-{\displaystyle \frac{5}{2}}\unicode[STIX]{x1D6FD}_{w}^{2}+O(\unicode[STIX]{x1D6FD}_{w}^{3})\right). & \displaystyle\end{eqnarray}$$

For the case when $\tilde{\unicode[STIX]{x1D707}}_{w}$ tends to zero (or $\unicode[STIX]{x1D707}\rightarrow \unicode[STIX]{x1D707}_{d}$ ), the integral (3.11) will exhibit a logarithmic singularity at ${\tilde{y}}=0$ , and so $\tilde{u} _{w}^{\infty }\sim -\log \tilde{\unicode[STIX]{x1D707}}_{w}/\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ . For the sake of completeness, it is pointed out that, even when $\unicode[STIX]{x1D6FD}_{w}$ is not smaller than one, an expression similar to (3.12) could be developed at $\tilde{\unicode[STIX]{x1D707}}_{w}\ll 1$ by exploiting the logarithmic singularity in the integral (3.11).

Figure 5.

Comparison between direct computations and approximate expressions of the steady-state wall speed $\tilde{u} _{w}^{\infty }=\tilde{u} ({\tilde{y}}=0,\tilde{t}\gg 1)$ as a function for the entire range of normalised wall friction coefficient $\tilde{\unicode[STIX]{x1D707}}_{w}=(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{w})/(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707}_{s})$ . Black circles denote the numerical solution using (2.15). The continuous line (blue) and the dashed line (red) are obtained using the expressions (3.12) and (3.13), respectively. The latter gives a good match when the applied shear stress is just above the yield shear stress $\tilde{\unicode[STIX]{x1D707}}_{w}\sim 1$ (or $\unicode[STIX]{x1D707}_{w}\sim \unicode[STIX]{x1D707}_{s}$ ), whereas the former captures the trend for all values of $\tilde{\unicode[STIX]{x1D707}}_{w}\in [0,1]$ .

Figure 6.

Velocity profile at different times $\tilde{t}=10^{-3},10^{-2},10^{-1},10^{0},10^{1}$ in the shear layer for various normalised wall friction coefficients $\tilde{\unicode[STIX]{x1D707}}_{w}=0.9$ (a), $\tilde{\unicode[STIX]{x1D707}}_{w}=0.5$ (b) and $\tilde{\unicode[STIX]{x1D707}}_{w}=0.1$ (c). The thick dashed line represents the steady-state solution. These data are obtained by directly integrating the numerical solution of the non-dimensional shear-stress equation (2.15). The data collapse features in the second and third rows imply that the $\tilde{u} ({\tilde{y}},\tilde{t})\sim \tilde{u} _{w}^{\infty }\sqrt{\tilde{t}}f({\tilde{y}}/2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}})$ , where $\tilde{u} _{w}^{\infty }$ is given by the expression (3.12).

The expressions (3.12) and (3.13) for $\tilde{u} _{w}^{\infty }$ can now be verified by plotting the steady-state wall speed with respect to the normalised wall friction coefficient $\tilde{\unicode[STIX]{x1D707}}_{w}$ as in figure 5. Here, the exact wall speed (open circles) is computed by substituting the steady-state solution (2.10) in (3.10) and integrating it numerically. The asymptotic results (3.12) and (3.13) are displayed as continuous (blue) and dashed (red) lines, respectively. The normalised wall velocity varies slowly with the normalised friction coefficient $\tilde{\unicode[STIX]{x1D707}}_{w}$ as long as $1-\tilde{\unicode[STIX]{x1D707}}_{w}$ is small (or $\unicode[STIX]{x1D707}_{w}\sim \unicode[STIX]{x1D707}_{s}$ ). Thus, for a given normal stress $p_{w}$ at the wall, $\tilde{u} _{w}$ varies linearly with the applied shear stress $\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D707}_{w}p_{w}$ if the latter is sufficiently close to the yield shear stress $\unicode[STIX]{x1D70F}_{s}=\unicode[STIX]{x1D707}_{s}p_{w}$ . Then the wall velocity increases rapidly with $1-\tilde{\unicode[STIX]{x1D707}}_{w}$ and it becomes logarithmically large as $\unicode[STIX]{x1D707}_{w}$ approaches $\unicode[STIX]{x1D707}_{d}$ . This is not surprising since when $\unicode[STIX]{x1D70F}_{w}$ approaches $\unicode[STIX]{x1D707}_{d}p_{w}$ , the inertial number satisfies $I\gg 1$ in such a manner that frictional grain–grain contacts become less dominant compared to grain–grain collisions, as internal grain rearrangements are very frequent compared to the local deformation rate (Andreotti et al. Reference Andreotti, Forterre and Pouliquen2011). Thus, a highly agitated flow can occur near the wall. Indeed, beyond this critical value, there is no longer an equilibrium between the applied shear stress and the internal resistance via frictional contacts as previously pointed out by Cawthorn (Reference Cawthorn2011). A more relevant description is given by models inspired by the kinetic theory of gases (Jenkins & Savage Reference Jenkins and Savage1983; Goldhirsch Reference Goldhirsch2003). Furthermore, figure 5 indicates that the expression (3.13) provides a very good approximation to the steady-state wall speed when the applied shear stress is just above the yield stress $\unicode[STIX]{x1D70F}_{s}=\unicode[STIX]{x1D707}_{s}p_{w}$ and the expression (3.12) remarkably captures the wall speed variation for all applied shear stresses $\unicode[STIX]{x1D70F}_{w}=\unicode[STIX]{x1D707}_{w}p_{w}\in [\unicode[STIX]{x1D707}_{s}p_{w},\unicode[STIX]{x1D707}_{d}p_{w}]$ .

The above result that the wall speed $\tilde{u} _{w}(\tilde{t})\propto \sqrt{\tilde{t}}$ along with the fact that $\tilde{u} _{w}(\tilde{t}\gg 1)$ follows a universal trend (3.12) suggests that the velocity profile can be approximately deduced from the short-time asymptotic solution (3.3). This hypothesis is explored in figure 6, which shows three different normalisations of the velocity field as a function of time. Each column corresponds to specific wall boundary conditions corresponding to $\tilde{\unicode[STIX]{x1D707}}_{w}=0.9$ (a), $\tilde{\unicode[STIX]{x1D707}}_{w}=0.5$ (b) and $\tilde{\unicode[STIX]{x1D707}}_{w}=0.1$ (c). The steady-state solution is displayed as a thick dashed line only in the first row but it is left out in the rest of the graphs for the sake of clarity. In this row, all figures present the time development of the velocity $\tilde{u} ({\tilde{y}},\tilde{t})$ as computed by numerically integrating (3.10) at each $\tilde{t}=10^{-3},10^{-2},10^{-1},10^{0},10^{1}$ . As seen before in figure 4, it can be readily observed that the wall velocity $\tilde{u} ({\tilde{y}}=0,\tilde{t})$ increases in time and attains a steady-state value which in turn depends on the applied wall shear stress. The numerical solution at $\tilde{t}=10^{1}$ is already superimposed on the steady-state solution (dashed line) for all the cases of $\tilde{\unicode[STIX]{x1D707}}_{w}$ shown here. The second and the third rows in figure 6 display the same data when the velocity field is normalised with the wall speed $\tilde{u} _{w}(\tilde{t})=\tilde{u} ({\tilde{y}}=0,\tilde{t})$ and $\tilde{u} _{w}^{\infty }\sqrt{\tilde{t}}$ , respectively, as a function of ${\tilde{y}}/2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}$ . Note that $\tilde{u} _{w}^{\infty }$ is taken from the expression (3.12). Irrespective of the velocity normalisation, all the velocity profiles collapse onto a unique curve except for the cases when $\tilde{t}>O(1)$ , as expected from the previous results. However, the collapse is only marginally good when $\tilde{\unicode[STIX]{x1D707}}_{w}=0.1$ when $\tilde{t}\sim O(1)$ or greater. This implies that this observation may apply at shorter and shorter times as $\unicode[STIX]{x1D707}_{w}$ tends towards $\unicode[STIX]{x1D707}_{d}$ . Nonetheless, as long as $\tilde{t}\leqslant O(1)$ , the velocity field should be given by $\tilde{u} ({\tilde{y}},\tilde{t})\sim \tilde{u} _{w}^{\infty }\sqrt{\tilde{t}}f({\tilde{y}}/2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}})$ , where $\tilde{u} _{w}^{\infty }$ is given by the expression (3.12).

3.4 Stokes’ first problem with imposed wall velocity

So far in this article Stokes’ first problem has been considered for the case when the wall shear stress $\unicode[STIX]{x1D70F}_{w}$ is imposed externally. Therefore, the wall velocity developed with time and, as the internal resistance of the granular medium balances out the applied stresses, it saturated to a constant value $u_{w}$ . In contrast, it should be possible to set the granular medium in motion by imposing the wall speed $u_{w}$ . Here, the resulting shear stress experienced by the wall should vary temporally as the internal resistance of the granular medium develops with time. However, it should also ultimately saturate to a constant value $\unicode[STIX]{x1D70F}_{w}$ so that it gives precisely the same steady-state solution as for the imposed shear-stress case. In this subsection, a brief note on this variant of Stokes’ first problem is presented.

As already seen in figure 5, in the steady-state solution, for each wall friction coefficient there exists one and only one wall velocity. Therefore, it is reasonable to leave the normalised variables of the previous sections as such. Now by using the normalised local friction coefficient $\tilde{\unicode[STIX]{x1D707}}=(\unicode[STIX]{x1D707}_{d}-\unicode[STIX]{x1D707})/\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ , (3.1) can be rewritten as

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}\tilde{t}}}=\unicode[STIX]{x0394}\unicode[STIX]{x1D707}^{2}\tilde{\unicode[STIX]{x1D707}}^{2}\left({\displaystyle \frac{2\unicode[STIX]{x1D6FD}_{w}}{\sqrt{1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}}}}{\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}{\tilde{y}}}}+\sqrt{1+\unicode[STIX]{x1D6FD}_{w}{\tilde{y}}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}{\tilde{y}}^{2}}}\right), & \displaystyle\end{eqnarray}$$

which for the case of Stokes’ first problem with imposed wall velocity should satisfy the initial condition $\tilde{\unicode[STIX]{x1D707}}({\tilde{y}},\tilde{t}=0)=1$ (or $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{s}$ ) along with the boundary conditions

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{w}=\int _{0}^{1}{\displaystyle \frac{\sqrt{\tilde{p}}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D707}}}\left({\displaystyle \frac{1}{\tilde{\unicode[STIX]{x1D707}}}}-1\right)\,\text{d}{\tilde{y}}, & \displaystyle\end{eqnarray}$$

and $\tilde{\unicode[STIX]{x1D707}}({\tilde{y}}=1,\tilde{t})=1$ . Here, the initial and lower boundary conditions are chosen so as to satisfy $\tilde{u} =0$ , which is possible only when $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{s}$ in the $\unicode[STIX]{x1D707}(I)$ -rheology. The integral condition imposes the wall velocity on the choice of the vertical distribution of the normalised local friction coefficient $\tilde{\unicode[STIX]{x1D707}}$ . In particular, note that the parameter $\unicode[STIX]{x1D6FD}_{w}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}g\unicode[STIX]{x1D6FF}_{s}/p_{w}$ as in (2.11) since the steady-state solution (2.10) for Stokes’ first problem with applied shear stress should also apply to this case where the wall speed is externally imposed. Thus $\unicode[STIX]{x1D6FD}_{w}=\unicode[STIX]{x1D707}_{w}^{\infty }/\unicode[STIX]{x1D707}_{s}-1$ , where $\unicode[STIX]{x1D707}_{w}^{\infty }$ is the steady-state wall friction coefficient that is needed to sustain the applied wall speed $u_{w}$ .

Figure 7.

Numerical results for Stokes’ first problem with imposed wall velocity showing the evolution of the velocity $\tilde{u} ({\tilde{y}},\tilde{t})$ (top) and the normalised friction coefficient $\tilde{\unicode[STIX]{x1D707}}({\tilde{y}},\tilde{t})$ (bottom) at different times $\tilde{t}=10^{-3},10^{-2},10^{-1},10^{0},10^{1}$ for various imposed wall speeds $\tilde{u} _{w}$ . These data are obtained by directly integrating the numerical solution of the non-dimensional shear-stress equation (3.14) with the wall boundary condition (3.15). The thick dashed line represents the steady-state solution. The imposed wall speed $\tilde{u} _{w}$ was chosen to match with the steady-state wall speed in figure 6.

Figure 8.

Temporal evolution of the shear layer thickness $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{s}$ for the profiles presented in figure 7. Different symbols correspond to various imposed wall velocities (▫, $\tilde{u} _{w}=0.198$ ; ○, $\tilde{u} _{w}=1.334$ ; ♢, $\tilde{u} _{w}=4.629$ ). As already seen in figure 2 (bottom) for the case when the wall shear stress is imposed, the shear layer grows as $\sqrt{\tilde{t}}$ until approximately $\tilde{t}\sim O(1)$ .

Figure 7 presents the velocity $\tilde{u} ({\tilde{y}},\tilde{t})$ (top) and the normalised friction coefficient $\tilde{\unicode[STIX]{x1D707}}({\tilde{y}},\tilde{t})$ (bottom) profiles that are obtained by numerically solving (3.14) for $\tilde{\unicode[STIX]{x1D707}}$ satisfying the imposed wall velocity condition as given by (3.15). Each continuous line represents a different time as indicated in the figure $(\tilde{t}=10^{-3},10^{-2},10^{-1},10^{0},10^{1})$ and the thick dashed line represents the steady-state solution. Note that particular wall speeds ((a)  $\tilde{u} _{w}=0.198$ , (b) $\tilde{u} _{w}=1.334$ and (c) $\tilde{u} _{w}=4.629$ ) were chosen so that the resulting steady-state velocity and local friction coefficient profiles are the same as that obtained in figure 6 for the case when wall shear stress is imposed. As expected, in all figures, the velocity profile (top) and the related boundary layer develops with time in such a way that the velocity at the wall is equal to the applied wall speed $\tilde{u} _{w}$ at all times, and the size of the mobile layer increases until approximately $\tilde{t}\sim O(1)$ . This is true for the normalised wall friction coefficient (bottom) as well. Also, both $\tilde{u} ({\tilde{y}},\tilde{t})$ and $\tilde{\unicode[STIX]{x1D707}}({\tilde{y}},\tilde{t})$ match with their respective steady-state profiles (thick dashed line) for sufficiently large time $\tilde{t}>O(1)$ .

With these data, it is then possible to calculate the size of the mobile layer, say $\unicode[STIX]{x1D6FF}(\tilde{t})$ , where the local shear stress is just above the threshold shear $\unicode[STIX]{x1D70F}_{s}=\unicode[STIX]{x1D707}_{s}p$ . Figure 8 illustrates that $\unicode[STIX]{x1D6FF}(\tilde{t})$ increases with time as $\sqrt{\tilde{t}}$ until $\tilde{t}\sim O(1)$ , after which it saturates to the steady-state limit. Even though this is similar to what was observed previously for the case of Stokes’ problem when the wall shear is imposed, as seen in figure 2 (bottom), the evolution of the wall friction coefficient $\unicode[STIX]{x1D707}_{w}(\tilde{t})$ takes place in two stages (see figure 9). It is observed that $\unicode[STIX]{x1D707}_{w}(\tilde{t})$ is close to $\unicode[STIX]{x1D707}_{d}$ for sufficiently small times, say up to when $\tilde{t}\leqslant \tilde{t}_{\unicode[STIX]{x1D707}d}$ . Then it decreases as $1/\sqrt{\tilde{t}}$ until it attains the steady state at $\tilde{t}\sim \tilde{t}_{\unicode[STIX]{x1D707}w}$ . Thus, the corresponding wall shear stress $\unicode[STIX]{x1D70F}_{w}(\tilde{t})=\unicode[STIX]{x1D707}_{w}(\tilde{t})p_{w}$ initially remains sufficiently close to the critical shear stress $\unicode[STIX]{x1D707}_{d}p_{w}$ before decreasing monotonically towards the wall shear stress that is needed to sustain the applied wall velocity $\tilde{u} _{w}$ . It suggests that the time at which these two stages occur should depend on $\tilde{u} _{w}$ and these time scales are different from $\unicode[STIX]{x1D6FF}_{s}^{2}/\unicode[STIX]{x1D708}_{g}$ .

Figure 9.

Evolution of the wall friction coefficient $\unicode[STIX]{x1D707}({\tilde{y}}=0,\tilde{t})=\unicode[STIX]{x1D707}_{w}(\tilde{t})$ for various imposed wall velocities (▫, $\tilde{u} _{w}=0.198$ ; ○, $\tilde{u} _{w}=1.334$ ; ♢, $\tilde{u} _{w}=4.629$ ). Figures on the right display the same data as a function of two different normalisations of the time variable based on the scaling laws (3.17) and (3.18).

For $\tilde{t}_{\unicode[STIX]{x1D707}d}\leqslant \tilde{t}\leqslant \tilde{t}_{\unicode[STIX]{x1D707}w}$ , an order of magnitude analysis of (2.5) gives $u_{w}/t\sim (\unicode[STIX]{x1D707}_{w}(t)-\unicode[STIX]{x1D707}_{s})p_{w}/\unicode[STIX]{x1D6FF}(t)$ . In terms of the non-dimensional time $\tilde{t}$ and wall velocity $\tilde{u} _{w}$ , this can be rewritten as

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{w}(\tilde{t})-\unicode[STIX]{x1D707}_{s}\sim {\displaystyle \frac{\tilde{u} _{w}}{3.66\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{t}}}, & \displaystyle\end{eqnarray}$$

since $\unicode[STIX]{x1D6FF}\sim 3.66\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}$ as shown by figure 8. Since when $\tilde{t}\sim \tilde{t}_{\unicode[STIX]{x1D707}d}$ , $\unicode[STIX]{x1D707}_{w}(\tilde{t})-\unicode[STIX]{x1D707}_{s}\approx \unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ , this implies that the wall friction coefficient should be around $\unicode[STIX]{x1D707}_{d}$ until some time

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{t}_{\unicode[STIX]{x1D707}d}\propto \tilde{u} _{w}^{2}. & \displaystyle\end{eqnarray}$$

Thus, according to $\unicode[STIX]{x1D707}(I)$ -rheology, in order to sustain an applied wall speed $\tilde{u} _{w}$ with the no-slip condition, the granular medium develops a strong shear stress (about $\unicode[STIX]{x1D707}_{d}p_{w}$ ) near the wall for all times $\tilde{t}\leqslant \tilde{t}_{\unicode[STIX]{x1D707}d}$ . From this point onwards, the wall shear stress should decrease monotonically as $1/\sqrt{\tilde{t}}$ until the local friction coefficient (or the local shear stress) distribution has reached the steady-state solution to give $\unicode[STIX]{x1D707}_{w}(\tilde{t}\gg 1)=\unicode[STIX]{x1D707}_{w}^{\infty }$ . Here, $\unicode[STIX]{x1D707}_{w}^{\infty }$ is a function of the applied wall speed and it can be either estimated directly by integrating the steady-state profile (2.10) or approximately via the expression (3.12). Using (3.16), it is seen that the wall friction coefficient, and hence the wall shear stress, should attain a steady-state value at some time

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{t}_{\unicode[STIX]{x1D707}w}\propto \left({\displaystyle \frac{\tilde{u} _{w}}{\unicode[STIX]{x1D707}_{w}^{\infty }-\unicode[STIX]{x1D707}_{s}}}\right)^{2}. & \displaystyle\end{eqnarray}$$

These two time scales can be verified by plotting the data from figure 9(a) as a function of properly rescaled time with respect to the relations (3.17) and (3.18). This is done in figure 9 (see plots b,c) where the data collapse indicates a good agreement with the above scaling laws.

Figure 10.

The same data as in figure 7 (bottom) for different times $\tilde{t}=10^{-3},10^{-2},10^{-1},10^{0},10^{1}$ but plotted here against the self-similar variable $\unicode[STIX]{x1D702}={\tilde{y}}/2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}$ . A reasonably good collapse is observed.

Finally, an approximate expression can now be elaborated for $\unicode[STIX]{x1D707}({\tilde{y}},\tilde{t})$ due to the fact that the shear layer $\tilde{\unicode[STIX]{x1D6FF}}$ is small up to $\tilde{t}\sim O(1)$ . As done before in § 3.1, (3.14) becomes

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}\tilde{t}}}=(\unicode[STIX]{x0394}\unicode[STIX]{x1D707}^{2}\tilde{\unicode[STIX]{x1D707}}^{2}){\displaystyle \frac{\unicode[STIX]{x2202}^{2}\tilde{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}{\tilde{y}}^{2}}}, & \displaystyle\end{eqnarray}$$

at short time $\tilde{t}$ . This equation is singular near the wall region ${\tilde{y}}\ll 1$ since $\tilde{\unicode[STIX]{x1D707}}\sim 0$ at least until $\tilde{t}_{\unicode[STIX]{x1D707}d}$ (see figure 9). Nonetheless, the chosen initial condition $\unicode[STIX]{x1D70F}(y,t=0)=\unicode[STIX]{x1D707}_{s}p$ implies that there exists a zone where $\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D707}_{s}$ (or $\tilde{\unicode[STIX]{x1D707}}\sim 1$ ) away from the wall until the shear layer is completely developed, i.e. for all $\tilde{t}\leqslant O(1)$ . In this outer zone $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\tilde{\unicode[STIX]{x1D707}}\sim \unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ and hence the above equation reduces to a simple diffusion equation wherein the outer solution should be $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{s}+A\,\text{erfc}\,({\tilde{y}}/2\unicode[STIX]{x0394}\sqrt{t})$ . Here $A$ is an arbitrary constant that could be deduced by matching this solution to the inner region where the friction coefficient $\tilde{\unicode[STIX]{x1D707}}\sim 0$ . Here, $\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D707}_{w}(\tilde{t})$ and therefore

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}({\tilde{y}},\tilde{t})=\unicode[STIX]{x1D707}_{s}+(\unicode[STIX]{x1D707}_{w}(\tilde{t})-\unicode[STIX]{x1D707}_{s})\,\text{erfc}\left({\displaystyle \frac{{\tilde{y}}}{2\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}}}\right) & \displaystyle\end{eqnarray}$$

is obtained as an approximate expression for the friction coefficient $\unicode[STIX]{x1D707}({\tilde{y}},\tilde{t})$ as long as $\tilde{t}$ is sufficiently small. When this expression is compared with the numerical results (see figure 10) a reasonable fit is observed for all $\tilde{t}\leqslant O(1)$ . As previously shown from computational data in figure 9, it is pointed out here that $\unicode[STIX]{x1D707}_{w}(\tilde{t})-\unicode[STIX]{x1D707}_{s}\sim \unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ for all $\tilde{t}\leqslant \tilde{t}_{\unicode[STIX]{x1D707}d}$ and $\unicode[STIX]{x1D707}_{w}(\tilde{t})-\unicode[STIX]{x1D707}_{s}\sim \tilde{u} _{w}/\unicode[STIX]{x0394}\unicode[STIX]{x1D707}\sqrt{\tilde{t}}$ when $\tilde{t}_{\unicode[STIX]{x1D707}d}\leqslant \tilde{t}\leqslant \tilde{t}_{\unicode[STIX]{x1D707}w}$ .