A note on Stokes' problem in dense granular media using the $μ(I)$--rheology
J. John Soundar Jerome Bastien Di Pierro


Optical conductivity of the Weyl semimetal NbP
    Artem V. Pronin
    March 26, 2018
==============================================




The classical Stokes' problem describing the fluid motion due to a steadily moving infinite wall is revisited in the context of dense granular flows of mono-dispersed beads using the recently proposed μ(I)–rheology. In Newtonian fluids, molecular diffusion brings about a self-similar velocity profile and the boundary layer in which the fluid motion takes place increases indefinitely with time t as √(ν t), where ν is the kinematic viscosity. For a dense granular visco-plastic liquid, it is shown that the local shear stress, when properly rescaled, exhibits self-similar behaviour at short-time scales and it then rapidly evolves towards a steady-state solution. The resulting shear layer increases in thickness as √(ν_g t) analogous to a Newtonian fluid where ν_g is an equivalent granular kinematic viscosity depending not only on the intrinsic properties of the granular media such as grain diameter d, density ρ and friction coefficients but also on the applied pressure p_w at the moving wall and the solid fraction ϕ (constant). In addition, the μ(I)–rheology indicates that this growth continues until reaching the steady-state boundary layer thickness δ_s = β_w (p_w/ϕρ g ), independent of the grain size, at about a finite time proportional to β_w^2 (p_w/ρ g d)^3/2√(d/g), where g is the acceleration due to gravity and β_w = (τ_w - τ_s)/τ_s is the relative surplus of the steady-state wall shear-stress τ_w over the critical wall shear stress τ_s (yield stress) that is needed to bring the granular media into motion. For the case of Stokes' first problem when the wall shear stress τ_w is imposed externally, the μ(I)–rheology suggests that the wall velocity simply grows as √(t) before saturating to a constant value whereby the internal resistance of the granular media balances out the applied stresses. In contrast, for the case with an externally imposed wall speed u_w, the dense granular media near the wall initially maintains a shear stress very close to τ_d which is the maximum internal resistance via grain-grain contact friction within the context of the μ(I)–rheology. Then the wall shear stress τ_w decreases as 1/√(t) until ultimately saturating to a constant value so that it gives precisely the same steady state solution as for the imposed shear stress case. Thereby, the steady-state wall velocity, wall shear stress and the applied wall pressure are related as u_w ∼  (g δ_s^2/ν_g) f( β_w ) where f( β_w ) is either 𝒪(1) if τ_w ∼τ_s or logarithmically large as τ_w approaches τ_d. 









§ INTRODUCTION



Flowing matter containing a dense collection of grains like sand, gravel, cereals, sugar, etc. is ubiquitous in nature as well as in many industrial processes. Such granular media exist in various states at any given common flow situation depending strongly on the energy supplied by external deformation and/or shear stresses <cit.>. And so, they show a very rich phenomenology <cit.> : a gaseous regime wherein the flow is very rapid and dilute, and the particles interact by collision <cit.> and a quasi-static regime in which the material deformation is extremely slow wherein frictional contacts between particles dominate the rheology as often in soil mechanics <cit.>. Indeed, there exists an intermediate regime in the presence of both collisions and friction that result in huge dissipation. Here a dense granular media behaves like a viscoplastic liquid <cit.>. A decade ago, generalising the scalar rheology of <cit.>, <cit.> proposed the so-called μ(I)–rheology to describe such a dense granular liquid state.  It has since been well-exploited often via direct numerical simulations to study and model many a common flow configurations <cit.>.

However, recent works by <cit.> illustrate that the governing equations of the μ(I)–rheology can exhibit ill-posed behaviour in the parameter range corresponding to quasi-gaseous and quasi-static regimes, respectively. <cit.> showed that ill-posed problems suffer from the so-called Hadamard instability and so, they characterized ill-posedness through a stability analysis that identifies exponential temporal growth of short-wavelength perturbations. And as a consequence, grid-dependent numerical results may not converge as the spatial refinement is enhanced for these cases <cit.>. In particular, <cit.> demonstrated both theoretically and numerically the governing equations of the μ(I)–rheology are Hadamard unstable even for the simple case of Bagnold flow. More recently, <cit.> also observed it in their numerical simulations for the case of granular column collapse on inclined channels. Nonetheless novel attempts to regularise the governing equations via a proper functional form of μ(I), atleast in the quasi-static regime, have already been proposed by <cit.>. They successfully simulated granular roll-waves in two dimensions and it now remains to see if their regularisation is valid also in direct computations of other unsteady granular flows.



In this context, this work aims to determine, both numerically and theoretically, the time evolution characteristics of an unsteady, non-uniform velocity and shear stress fields arising in the μ(I)–rheology for a canonical flow situation, namely, the so-called Stokes' first problem <cit.> of the fluid motion that is brought about by impulsively starting an infinite wall. Unlike the classical case, the granular media is placed underneath the plate (see figure <ref>). It is the simplest unsteady parallel flow in which some important features of fluid flows such as transverse momentum transfer and the resulting boundary layer development due to direct balance between local fluid acceleration and the friction forces can be treated. It is also known as the dragged-plate problem in <cit.> where only a steady-state analytical solution can so far been found. The objective of the present work is to treat the transient solution and its characteristics.

Note that there has been a steady interest in Stokes's first problem for non-Newtonian fluids, in particular, visco-elastic fluids <cit.>. Similarly, Stokes' second problem <cit.> that considers the time evolution of the velocity field due to an horizontally oscillating infinite flat plate is recently studied for viscoelastic fluids by <cit.>. Its applications include high-frequency microfluidics <cit.> for viscoelastic materials and rheometers for viscoplastic <cit.> and power law fluids <cit.>. Recent literature also considers a third type of Stokes' problem wherein a transient velocity field is set-up by suddenly applying a body force to the fluid that is initially at rest. In fact, for a granular media, <cit.> used this configuration to numerically validate their proposed μ(I)–rheology. This has later been referred to as Stokes' third problem by <cit.> for a Herschel–Bulkley material and interestingly, for the numerical resolution, the authors resorted to a Stefan problem, with a moving interface (boundary condition) that separates the sheared and unsheared regions.

The constitutive laws for many non-Newtonian fluids are often non-linear but they can be simplified in the case of Stokes' problems. Whereas the yield stress of the granular material varies in space since the μ(I)–rheology proposes a constitutive law for dense granular flows wherein the media behaves like a visco-plastic liquid with the local viscosity non-linearly related to the local strain–rate as well as the local pressure. Moreover, care should be taken to express a well-posed initial-value problem using the μ(I)–rheology to avoid Hadamard instability. Finally, it is only recently dense granular flows have been successfully studied using a continuum model and so, the governing equations have so far been unexplored even for Stokes' first problem. In addition, apart from the simple case of steady Bagnold flow over an inclined plane, the constant shear flow case, the steady state solution of the Stokes' first-problem <cit.> there exist only few analytical results describing unsteady dense granular flows (see the notable recent exception of <cit.> who gave entrainment rates in the case of transient heap flows from the depth-integrated layer dynamics assuming a local μ(I)–rheology). Therefore, this brief note is aimed at bringing out the key features of this canonical problem as predicted by the μ(I)–rheology. 

The article is set as follows. Firstly, the governing equations are shown to result in a single non-linear shear stress diffusion equation. Its numerical solution is then computed for the case when the wall shear stress is imposed while letting the wall velocity to develop with time. Some approximate unsteady solutions are obtained and compared with computations. Finally, a brief note on the Stokes' first problem with imposed wall speed is given. 




§ GOVERNING EQUATIONS





 §.§ Constitutive laws : μ(I)–rheology


Analogous to Coulomb's friction law, using dimensional arguments, experiments and numerical simulations, <cit.> demonstrated that the shear stress τ  is proportional to the normal stress P for 2D dense granular flows of rigid particles so that τ = μ (I) where the local friction coefficient μ is only a function of a dimensionless parameter called the Inertial number I =γ̇d/√(P/ρ). Here, γ̇ is the local shear rate which is related to the macroscopic timescale of the granular flow and d√(ρ/) is the microscopic timescale corresponding to any local rearrangement of grains of diameter d and density ρ 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 <cit.> or <cit.>, respectively. In general, the dimensionless local friction coefficient is given by <cit.>


    μ (I) = μ_s + μ_d -μ_s/(1 + I_0/I),


whereby μ saturates towards two fundamental constants for a dense granular media μ_s or μ_d depending respectively on the inertial number I ≪ 1 (quasi-static regime) or I ≫ 1 (kinetic or gaseous regime). <cit.> proposed a 3D 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 τ_ij while assuming that τ_ij is aligned with the strain rate tensor γ̇_ij = 1/2( ∂_i u_j + ∂_j u_i ) (where u_i represent components of the velocity field). So, if x⃗ is the position vector and t represents time, then Cauchy stress tensor


    σ_ij(x⃗, t) =  -p(x⃗, t) δ_ij + τ_ij(x⃗, t),


where δ_ij is the kronecker delta and τ_ij(x⃗, t) =  η(x⃗, t) γ̇_ij(x⃗, t) with the local granular liquid viscosity


    η = μ p/|γ̇|,


is, thereby, a non-linear function of the local pressure p and the local second invariant of the strain rate tensor |γ̇| = √(1/2γ̇_ijγ̇_ij) via the local friction coefficient given by <ref> and the Inertial number for the 3D case


    I =|γ̇| d/√(p/ρ).


In addition, the solid volume fraction ϕ is also a linear function of I <cit.> but, in general, it varies very little and so, for the sake of simplicity, it is taken to be a constant in the following.





 §.§ 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 media. As illustrated in figure <ref>, the plate is set to motion impulsively at t > 0 by applying a tangential shear stress τ_w (along x-direction) in the presence of a normal stress p_w (along 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 ≥ 0. And so the only non-zero components of the strain rate tensor are γ̇_xy = γ̇_yx = 1/2( ∂ u / ∂ y ), where u is the x-component of the velocity field. Therefore, the 2D shear stress tensor is completely determined by a single scalar shear field τ =   μ p so that the x and y momentum equations become, respectively,


    ϕρ∂ u/∂ t = -∂/∂ y( μ p),



    ∂ p/∂ y =  ϕρ,


where  is the acceleration due to gravity. Initially the granular media 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 in not exceeded so that τ≤μ_s p. Within the context of the μ(I)–rheology, the only possible configuration were the granular liquid is at rest corresponds to the case τ = μ_s p. For all other static states, the constitutive law τ = μ p is violated and so, the μ(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 media


    τ (y, t = 0) = μ_s p


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


    p(y = 0, t)    =    p_w  
    μ(y = 0, t)     =    τ_w/p_w = μ_w ,


along with the condition that the grains sufficiently far from the plate remain static i.e. u(y = ∞, t) = 0 and so, μ = μ_s at y = ∞. In contrast to the classical Stokes' problem <cit.> for a Newtonian fluid, firstly, the frictional force (R.H.S. of eqn. <ref>) is not only non-uniform due to the hydrostatic pressure p(y) = p_w + ϕρ y but also, non-linear since μ depends on both p(y) and ∂ u / ∂ y via eqn. <ref>. And secondly, there exists a non-trivial steady-state solution wherein the shear stress is constant throughout the media such that μ p = τ_w = μ_w p_w as already shown in <cit.> and <cit.>. Since μ∈ [μ_s, μ_d] and p(y) increases linearly with the depth y, it follows that 


    μ = μ_w/1 + (ϕρ y/p_w),


for all y ≤δ_s and μ = μ_s otherwise; here the critical depth δ_s is given by


    δ_s = β_w (p_w/ϕρ),


with β_w = (μ_w/μ_s - 1). It denotes the depth beyond which the granular media does not flow y ≥δ_s. Note that the term β_w represents the surplus in wall shear stress over the yield criterion at the wall. The steady-state solution <ref> can be used to obtain the corresponding velocity profile. As previously shown by <cit.>, the resulting relation between the wall velocity (if no-slip condition is allowed) and the shear layer thickness compares qualitatively well with the molecular dynamic simulations of <cit.>. 




 §.§ 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 <ref>. Note that the latter is coupled with the expressions for the local friction coefficient (<ref>) and the inertial number I. However, it is possible to write a single equation for the shear stress in the case of the Stokes' first problem with in the context of the μ(I)–rheology. Since I=-1/2∂ u/∂ y√(ρ d^2/p), the local friction coefficient <ref> can be rewritten as


    ∂ u/∂ y = - 2 I_0 √(p/ρ d^2)( μ - μ_s/μ_d - μ),


which when differentiated w.r.t. time and using τ = μ p yields

    ∂^2 u/∂ t ∂ y = - 2 Δμ I_0/d(p √(p/ρ)) 1/( μ_d p - τ)^2∂τ/∂ t,


where Δμ = μ_d - μ_s > 0. By allowing ∂^2 u/∂ y ∂ t = ∂^2 u/∂ t ∂ y and then introducing <ref> in the above expression, it leads to a single non-linear shear-stress diffusion equation


    ∂τ/∂t = ( d/2 ϕΔμ I_0 √(ρ))(μ_d p - τ)^2/p √(p)∂^2 τ/∂y^2.


Finally, by taking ν_g = (d/2ϕΔμ I_0) √(p_w/ρ) as a proper diffusion coefficient and the steady–state shear layer thickness δ_s as the characteristic length scale, the non-dimensional time and space coordinates are  t̃ = ν_g t/δ_s^2 and ỹ = y/δ_s, respectively. Thus, in terms of the normalised pressure p̃ = p(y)/p_w and shear stress τ̃ = τ(y, t)/τ_w, the above equation becomes


    ∂τ̃/∂t̃ = (μ_d p̃ - μ_w τ̃)^2/p̃√(p̃)∂^2 τ̃/∂ỹ^2,


with boundary conditions τ̃(0, t̃) = τ̃(1, t̃) = 1 and an initial condition τ̃(ỹ, 0) = μ_s p̃/μ_w. The steady-state solution for the non-dimensional shear stress is τ̃ = 1. Unlike the steady-state local friction coefficient μ, the steady-state shear stress is a continuous and infinitely differentiable function for all y ≥ 0. So, it is expected that τ̃ remains smooth also for the unsteady case. In fact, the term (μ_d p̃ - μ_w τ̃) is positive–definite. Hence, it is quite straight-forward 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-Nicholson one for temporal integration. The updated τ̃ is obtained by an iterative Richardson Minimal Residual process.



 §.§ Well-posedness of the shear stress diffusion equation


<cit.> demonstrated that the μ(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 ≫ 1 in the neighbourhood of the wall either when the applied wall-shear τ_w is close to the critical wall shear τ_d = μ_d p_w or when wall speed u_w is sufficiently large. In addition, there is always a zone where I ≪ 1 (or μ∼μ_s) since the media is slowly-moving or stationary when y ∼δ_s. As already shown by <cit.> this should provoke Hadamard instability <cit.> 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 equation <ref> is indeed well-posed for all values of inertial number I. 




In fact, this non-linear diffusion equation <ref> is one-dimensional and when it is linearized about an arbitrary base state as in <cit.>, the resulting dispersion equation of the normal mode analysis in the high-wavenumber limit should be

    λ = - α^2 ξ_y^2,

where λ is the complex frequency, ξ_y is the wave number in the ỹ-direction and  α = α(p̃^0, τ̃^0) is some function of the base state pressure p̃^0 and shear stress τ̃^0. Since the real-part of the complex frequency λ provides the perturbation growth rate, it is straight-forward to see that the equation <ref> is stable for short-waves along the y-coordinate.

Note that this is an one-dimensional analysis since it takes into account only plane waves along the ỹ-axis in order to analyse the well-posedness of <ref> for one-dimensional, time dependent computations. For two-dimensional codes that consider granular Stokes problem as a test case, the situation remains more complex. In this case, as previously shown by <cit.>, the equations are still ill-posed since oblique two-dimensional short waves are unstable in the region y ∼δ_s (or close to the wall for μ_w ∼μ_d).



§ UNSTEADY SOLUTIONS AND THEIR CHARACTERISTICS


Figure <ref> presents the results of such numerical solutions (continuous lines) for three typical values of applied wall shear stress when μ_s = tan 21^∘ and μ_d = tan 33^∘ (typical values for spherical mono–dispersed glass beads as in <cit.>). Each graphic (top) depicts the normalised shear stress τ̃ = μp̃/μ_w variation in the y-direction at various time t̃ = 10^-4, 10^-3, 10^-2, 10^-1, 10^0, 10^1 for typical values of the normalised wall friction coefficient (μ_d - μ_w)/Δμ (where, Δμ = μ_d - μ_s). In all cases, the static initial condition τ̃(0, ỹ) = μ_s p̃/μ_w (dashed line), wherein the local shear stress is taken to be the yield criterion μ_s p, evolves continuously towards the steady–state solution τ̃ = μp̃/μ_w = 1. The spatial variation of τ̃ shows that there exists a layer in which τ̃ is greater than the yield stress and so, the granular media should flow in this region. If the size of this shear layer, say δ(t̃), is defined as the region where τ̃ = 0.999 μ_s p̃/μ_w, figure <ref> (bottom) clearly illustrates that δ(t̃) increases with time as √(t̃) until t̃∼ (1), after which it saturates to the steady-state limit. Therefore, it is expected from these results that, for any general μ_d and μ_s, approximate solutions to <ref> can be obtained at both t̃≪ 1 and t̃≫ 1 by properly linearising it.





 §.§ Self-similarity at t̃≪ 1


When t̃≪ 1, the non-dimensional shear layer thickness δ/δ_s is small as observed in figure <ref> (bottom). By taking the non-dimensional local pressure p̃ = 1 + β_w ỹ and local shear stress τ̃ = μ (1 + β_w ỹ)/μ_w, the diffusion equation <ref> becomes


    ∂μ/∂t̃ = (μ_d - μ)^2 ( 2 β_w/√(1 + β_w ỹ)∂μ/∂ỹ + √(1 + β_w ỹ)∂^2 μ/∂ỹ^2),


with μ(ỹ, 0) = μ_s, μ(0, t̃) = μ_w and μ(1, t̃) = μ_s. Note that, in general, β_w < 1 and hence, for t̃≪ 1, the highest–order derivative of μ should dominate if ỹ≤δ/δ_s so that μ∼μ_w - (μ_w - μ_s) y/δ. And in the outer region, μ∼μ_s. As the spatial variations of the local friction coefficient μ are stronger inside the shear layer (i.e., when β_w ỹ≪ 1), it is reasonable to simplify <ref> to


    ∂μ/∂t̃ = Δμ^2  ∂^2 μ/∂ỹ^2,


at the leading order with the same boundary conditions as before. By taking ψ = μ - ( μ_w - (μ_w - μ_s)ỹ), the above equation admits a self-similar solution for ψ = ψ(η) with η = ỹ/(2 Δμ√(t̃)) satisfying the initial and the boundary conditions for all t̃≥ 0. Thereby, the local friction coefficient is deduced to be


    μ∼μ_s + (μ_w - μ_s)(1 - ỹ) ( ỹ/2 Δμ√(t̃)),


which implies that the shear layer thickness should grow as δ / δ_s ∼ 4.66 Δμ√(t̃). 

The expression <ref> is presented in figure <ref> (top) where the time evolution of the local friction coefficient is displayed as function of ỹ at different time t̃ = 10^-4, 10^-3, 10^-2, 10^-1, 10^0, 10^1 for the same values of (μ_d - μ_w)/Δμ as in figure <ref>. When compared with numerical solutions of μ (and also, δ) as seen in figure <ref> (top) (and figure <ref>, respectively), these approximations are very satisfactory for all t̃≪ 1. Indeed, the expression <ref> for the local friction coefficient and especially, the estimations of the time evolution of the shear layer thickness δ / δ_s ∼ 4.66 Δμ√(t̃) are reasonably good even when t̃ is of order 1.





 §.§ Diffusion at t̃≫ 1





As soon as t̃∼ (1) the non-dimensional shear layer thickness is no longer small and hence, a singular perturbation of <ref> cannot be obtained with the present scaling for the y-coordinate. However, by using the non-dimensional pressure p̃ as an equivalent normalised spatial variable ŷ = 1 + β_w ỹ, it is possible to show that


    ∂τ̃/∂t̃ = (β_w (μ_d - μ))^2 √(ŷ)∂^2 τ̃/∂ŷ^2,


which is singular if (β_w Δμ)^2 tends to zero. This is often true since (β_w (μ_d - μ))^2 is of the order of (β_w Δμ)^2 and ŷ∼ (1). Thus, taking β_w (μ_d - μ) ∼β_w Δμ at the leading order in <ref>, it becomes linear and admits a WKB–approximation in the ŷ-coordinate. Thereby, the local shear stress τ̃ can be shown to be


    τ̃∼ 1 + ŷ^-1/8∑^∞_m = 1[ Λ_m exp(- λ_m^2 t̃)sin( m πŷ^3/4 - 1/(μ_w/μ_s)^3/4 - 1)],


where


    Λ_m    =    2 ∫^μ_w/μ_s_1ŷ^1/8( μ_s/μ_wŷ - 1 ) sin( m πŷ^3/4 - 1/(μ_w/μ_s)^3/4 - 1) d ŷ, 
    λ_m    =    3/4m πβ_w Δμ/(μ_w/μ_s)^3/4 - 1.


Note that it is possible to solve the linearised version of <ref> directly by separation of variables as well. In that case, the general solution is an infinite sum of Bessel functions in the ŷ–direction. But it is advantageous to work with the approximate solution <ref> when m is large. Nevertheless, the comparison between the WKB–approximation <ref> (dot-dash lines in figure <ref>) and the numerical solution for the local friction coefficient is good. And, in particular, the agreement is excellent when the applied friction coefficient μ_w is close to yield friction coefficient μ_s.





 §.§ Velocity field development


It is now possible to compute the temporal evolution of the velocity field using local shear stress τ = μ p. By rewriting eqn. <ref> in terms of t̃ and ỹ, a natural normalisation for the velocity u can be shown to be



    u = ũ(g δ_s^2/β_w ν_g),


where δ_s is the steady-state shear layer thickness (eqn. <ref>), β_w = (μ_w - μ_s)/μ_s and ν_g  = (d/ϕΔμ I_0) √(p_w/ρ) is the diffusion coefficient which appears in the granular Stokes' equation <ref>. Using this normalisation and by taking 


    μ̃ = μ_d - μ/Δμ,


the equation for the shear rate <ref> becomes


    ∂ũ/∂ỹ = -√(p̃)/Δμ( 1/μ̃ - 1 ),


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


Figure <ref> compares normalised wall velocity ũ_w = ũ(ỹ = 0, t̃) when obtained from the numerical solution of <ref> (filled circles), self-similar solution <ref> (□) and the long–time approximate solution <ref> (+) at a given boundary condition at the wall. Each graphic shows the temporal evolution of ũ_w when the wall friction coefficient is taken as (μ_d - μ_w)/(μ_d - μ_s) (a) 0.9, (b) 0.5 & (c) 0.1. In all cases, the wall velocity grows monotonically by following a power law in t̃ as long as t̃≤ 1 and then it saturates at the steady state value ũ_w^∞ = ũ(ỹ = 0, t̃ = ∞). Firstly, the self-similar solution <ref> gives a good qualitative agreement with numerical results and also, it matches well with the power law ũ_w(t̃) ∼√(t̃). Whereas the long–time WKB approximation <ref> matches very well with the numerical solution at all times t̃ > 1. Secondly, it is observed that the steady-state wall velocity ũ_w^∞ depends largely on the imposed wall-shear via μ̃_w = ( μ_d - μ_w )/Δμ as already shown by previous works <cit.>. In fact, for the steady state solution, <cit.> had already given an asymptotic solution for the wall speed  ũ_w^∞ when the applied shear stress is just above the yield shear i.e, when μ(ỹ = 0, t̃) approaches μ_s from above (or μ̃_w ∼ 1). However, it is also possible to obtain expressions for ũ_w^∞ for a wide range of μ̃_w  via the self-similar solution <ref> since, as suggested by figure <ref>, it gives a good approximation to the steady state wall speed.  In the limit when t̃≫ 1, the approximate solution <ref> becomes a function only of ỹ and so, in terms of the normalised friction coefficient <ref>, it is given by μ̃∼μ̃_w + (1 - μ̃_w)ỹ. Using this expression in <ref>, it reads


    ũ_w^∞∼∫_0^1√(1 + β_w ỹ)/Δμ( 1/μ̃_w + (1 - μ̃_w)ỹ -1 ) dỹ,


and since β_w = (μ_w/μ_s -1 ) < 1, it could be further developed to obtain







a simple expression for the the steady-state wall velocity


    ũ_w^∞∼ -1/Δμ[ 1 + logμ̃_w/(1 - μ̃_w)]+ 𝒪(β_w).


Noting that μ̃_w = 1 - μ_s β_w/ Δμ, it is straightforward to see that the first term in the above expression cancels out when μ̃_w ∼ 1 (or β_w = (μ_w/μ_s -1) ≪ 1) and thereby, it gives ũ_w^∞∼𝒪(β_w). In this case, as already deduced by <cit.>, the above integral leads to











    ũ_w^∞∼μ_s/2 Δμ^2( β_w - 5/2β_w^2 + 𝒪(β_w^3) ).


For the case when μ̃_w tends to zero (or μ→μ_d), the integral <ref> will exhibit a logarithmic singularity at ỹ = 0 and so, ũ_w^∞∼ -logμ̃_w/Δμ. For the sake of completeness, it is pointed out that, even when β_w is not smaller than one, an expression similar to <ref> could be developed at μ̃_w ≪ 1 by exploiting the logarithmic singularity in the integral <ref>.





















The expressions <ref> & <ref> for ũ_w^∞ can now be verified by plotting the steady-state wall speed with respect to the normalised wall friction coefficient μ̃_w as in figure <ref>. Here, the exact wall speed (open circles) is computed by substituting the steady state solution <ref> in equation <ref> and integrating it numerically. The asymptotic results <ref> & <ref> are displayed as continuous (blue) and  dashed (red) lines, respectively. The normalised wall velocity varies slowly with the normalised friction coefficient μ̃_w as long as 1 - μ̃_w is small (or μ_w ∼μ_s). Thus, for a given normal stress  p_w at the wall, ũ_w varies linearly with the applied shear stress τ_w = μ_w p_w if the latter is sufficiently close to the yield shear stress τ_s = μ_s p_w. Then the wall velocity increases rapidly with 1 - μ̃_w and it becomes logarithmically large as μ_w approaches μ_d. This is not surprising since when τ_w approaches μ_d p_w, the Inertial number I ≫ 1 in such a manner that frictional grain-grain contacts become less dominant compared to a grain-grain collisions as internal grain rearrangements are much frequent compared to the local deformation rate <cit.>. 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 <cit.>. A more relevant description is given by models inspired from kinetic theory of gases <cit.>. Furthermore, figure <ref> indicates that the expression <ref> provides a very good approximation to the steady-state wall speed when the applied shear stress is just above the yield stress τ_s = μ_s p_w. And the expression <ref> remarkably captures the wall speed variation for all applied shear stress τ_w = μ_w p_w ∈ [μ_s p_w, μ_d p_w].

The above result that the wall speed ũ_w(t̃) ∝√(t̃) along with the fact that ũ_w(t̃≫ 1) follows an universal trend <ref> suggests that the velocity profile can be approximately deduced from the short-time asymptotic solution <ref>. This hypothesis is explored in figure <ref> which shows three different normalisations of the velocity field as a function of time. Each column corresponds to specific wall boundary conditions corresponding to μ̃_w = 0.9 (left),  μ̃_w = 0.5 (middle) and  μ̃_w = 0.1 (right). 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 ũ(ỹ, t̃) as computed by numerically integrating <ref> at each t̃ = 10^-3, 10^-2, 10^-1, 10^0, 10^1. As seen before in figure <ref>, it can be readily observed that the wall velocity ũ(ỹ=0, t̃) increases in time and attains a steady-state value which in turn depends on the applied wall shear stress. The numerical solution at t̃ = 10^1 is already superposed on the steady state solution (dashed line) for all the cases of μ̃_w shown here. The second and the third row in figure <ref> display the same data when the velocity field is normalised with the wall speed ũ_w(t̃) = ũ(ỹ=0, t̃) and ũ_w^∞√(t̃), respectively, as a function of ỹ/2 Δμ√(t̃). Note that ũ_w^∞ is taken from the expression <ref>. Irrespective of the velocity normalisation, all the velocity profiles collapse on to an unique curve expect for the cases when t̃ > 𝒪(1) as expected from the previous results. However, the collapse is only marginally good when μ̃_w = 0.1 when t̃∼𝒪(1) or greater. This implies that this observation may apply at shorter and shorter times as μ_w tends towards μ_d. Nonetheless, as long as t̃≤𝒪(1), the velocity field should be given by ũ(ỹ, t̃) ∼ũ_w^∞√(t̃) f( ỹ/2 Δμ√(t̃)), where ũ_w^∞ is given by the expression <ref>.










 §.§ Stokes's first problem with imposed wall velocity


So far, in this article Stokes' first problem is considered for the case when the wall shear stress τ_w is imposed externally. Therefore, the wall velocity developed with time and, as the internal resistance of the granular media balances out the applied stresses, it saturated to a constant value u_w. In contrast, it should be possible to set the granular media under 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 media develops with time. However, it should also ultimately saturate to a constant value τ_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 the Stokes' first problem is presented.

As already seen in figure <ref>, 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 μ̃ = (μ_d - μ)/Δμ, the equation <ref> can be rewritten as


    ∂μ̃/∂t̃ = Δμ^2 μ̃^2 ( 2 β_w/√(1 + β_w ỹ)∂μ̃/∂ỹ + √(1 + β_w ỹ)∂^2 μ̃/∂ỹ^2),


which for the case of Stokes's first problem with imposed wall velocity should satisfy the initial condition μ̃(ỹ, t̃ = 0) = 1 (or μ = μ_s) along with the boundary conditions


    ũ_w =  ∫_0^1√(p̃)/Δμ( 1/μ̃ - 1 ) d ỹ,


and μ̃(ỹ=1, t̃) = 1. Here, the initial and lower boundary conditions are chosen so as to satisfy ũ = 0 which is possible only when μ = μ_s in the μ(I)-rheology. The integral condition imposes the wall velocity on the choice of the vertical distribution of the normalised local friction coefficient μ̃. In particular, note that the parameter β_w = ϕρ g δ_s / p_w as in <ref> since the steady-state solution <ref> for the Stokes' first problem with applied shear stress should also apply to this case where the wall speed is externally imposed. Thus, β_w = μ_w^∞/μ_s -1 where μ_w^∞ is the steady state wall friction coefficient that is needed to sustain the applied wall speed u_w.




Figure <ref> presents the velocity ũ(ỹ, t̃) (top) and the normalised friction coefficient μ̃(ỹ, t̃) (bottom) profiles that are obtained by numerically solving eqn. <ref> for μ̃ satisfying the imposed wall velocity condition as given by <ref>. Each continuous line represents different time as indicated in the figure (t̃ = 10^-3, 10^-2, 10^-1, 10^0, 10^1) and thick dashed line represents the steady-state solution. Note that particular wall speeds ((a)ũ_w = 0.198, (b) ũ_w = 1.334 & (c) ũ_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 <ref> 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 ũ_w at all times and the size of the mobile layer increases until about t̃∼𝒪(1). This is true for the normalised wall friction coefficient (bottom) as well. Also, both ũ(ỹ, t̃) and μ̃(ỹ, t̃) match with their respective steady-state profiles (thick dashed line) for sufficiently large time t̃ > 𝒪(1).

With these data, it is then possible to calculate the size of the mobile layer, say δ(t̃), where the local shear stress is just above the threshold shear τ_s = μ_s p. Figure <ref> illustrates that δ(t̃) increases with time as √(t̃) until t̃∼𝒪(1), after which it saturates to the steady-state limit. Even though this is similar to what was observed previously for the case of the Stokes' problem when the wall shear is imposed as seen in figure <ref> (bottom), the evolution of the wall friction coefficient μ_w(t̃) takes place in two stage (see figure <ref>). It is observed that μ_w(t̃) is close to μ_d for sufficiently small time, say up to when t̃≤t̃_μ d. Then it decreases as 1/√(t̃) until it attains the steady-state at t̃∼t̃_μ w.  Thus, the corresponding wall shear stress τ_w(t̃) = μ_w(t̃) p_w initially remains sufficiently close the critical shear stress μ_d p_w before decreasing monotonically towards the wall shear stress that is needed to sustain the applied wall velocity ũ_w. It suggests that the time at which these two stages occur should depend on ũ_w and these timescales are different from δ_s^2/ν_g. 

For t̃_μ d≤t̃≤t̃_μ w, an order of magnitude analysis of eqn. <ref> gives u_w/t∼(μ_w(t) - μ_s)p_w/δ(t). In terms of the non-dimensional time t̃ and wall velocity ũ_w, this can be rewritten as


    μ_w(t̃) - μ_s ∼ũ_w/3.66 Δμ√(t),


since δ∼ 3.66 Δμ√(t̃) as shown by figure <ref>. Since when t̃∼t̃_μ d, μ_w(t̃) - μ_s ≈Δμ and this implies that the wall friction coefficient should be around μ_d until some time


    t̃_μ d∝ũ_w^2.


Thus, according to μ(I)–rheology in order to sustain an applied wall speed ũ_w with no slip condition, the granular media develops a strong shear stress (about μ_d p_w) near the wall for all time t̃≤t̃_μ d. From there onwards, the wall shear stress should decrease monotonically as 1/√(t̃) until when the local friction coefficient (or the local shear stress) distribution has reached the stead-state solution to give μ_w(t̃≫ 1) = μ_w^∞. Here, μ_w^∞ is a function of the applied wall speed and it can be either estimated directly by integrating the steady-state profile <ref> or approximately via the expression <ref>. Using <ref>, it is seen that the wall friction coefficient and hence, the wall shear stress should attain a steady-state value at some time


    t̃_μ w∝(ũ_w/μ_w^∞ - μ_s)^2.


These two timescales can be verified by plotting the data from figure <ref> (left) as a function of properly rescaled time with respect to the relations <ref> & <ref>. This is done in figure <ref> (see plots on the right) where the data collapse indicates a good agreement with the above scaling laws.





Finally, an approximate expression can now be elaborated for μ(ỹ,t̃) due to the fact that the shear layer δ̃ is small up to t̃∼𝒪(1). As done before in section <ref>, the equation <ref> becomes


    ∂μ̃/∂t̃ = ( Δμ^2 μ̃^2 ) ∂^2 μ̃/∂ỹ^2,


at short-time t̃. This equation is singular near the wall region ỹ≪ 1 since μ̃∼ 0 at least until t̃_μ d (see figure <ref>). Nonetheless, the chosen initial condition τ (y, t = 0) = μ_s p implies that there exists a zone where μ∼μ_s (or μ̃∼ 1) away from the wall until the shear layer is completely developed i.e. for all t̃≤𝒪(1). In this outer zone Δμμ̃∼Δμ and hence, the above equation reduces to a simple diffusion equation wherein the outer solution should be μ = μ_s + A (ỹ/2 Δ√(t)). Here A is an arbitrary constant that could be deduced by matching this solution to the inner region where the friction coefficient μ̃∼ 0. Here, μ∼μ_w (t̃) and therefore


    μ (ỹ, t̃) = μ_s + (μ_w(t̃) - μ_s) (ỹ/2 Δμ√(t̃)),


is obtained as an approximate expression for the friction coefficient μ (ỹ, t̃) as long as t̃ is sufficiently small. When this expression is compared with the numerical results (see figure <ref>) a reasonable fit is observed for all t̃≤𝒪(1). As previously shown from computational data in figure <ref>, it is pointed out here that μ_w(t̃) - μ_s ∼Δμ for all t̃≤t̃_μ d and μ_w(t̃) - μ_s ∼ũ_w/Δμ√(t̃) when t̃_μ d≤t̃≤t̃_μ w.



§ CONCLUSION


Using the μ(I)–rheology, the so-called Stokes' first problem on the motion of a granular liquid set-in by an impulsively started flat plate is studied both numerically and theoretically. The problem was first well-posed in terms of a non-linear diffusion equation for the local shear stress with proper initial and boundary conditions in order to avoid Hadamard instability. Numerical solutions are then obtained for both externally imposed wall stress and speed. Approximate solutions at short and long times are also illustrated to capture the main features of the numerical results.

For the case when the dense granular flow is brought about by applying constant shear stress τ_w at t > 0, if τ_w is greater than the yield stress at the wall μ_s p_w (where p_w is the applied pressure at the wall) then the μ(I)–rheology implies that it is diffused into the granular media until the shear stress is uniform throughout the media. Such that, at any time t, the applied shear stress reaches a depth proportional to √(ν_g t) where ν_g = (d/2ϕΔμ I_0) √(p_w/ρ) denotes the diffusion coefficient for the local shear stress. A steady-state, wherein a finite zone of grains (of thickness, say, δ_s) yield, and hence flow, due to the applied shear, is thus shown to occur at about a finite time of the order δ_s^2/ν_g. Here, if β_w = (μ_w - μ_s)/μ_s, the shear layer depth δ_s is β_w p_w / ϕρ g as already obtained by <cit.>. If no-slip condition is allowed, the wall velocity develops in time as u_w^∞√(ν_g t)/ δ_s until about a time t ∼𝒪(δ_s^2/ν_g).

If the dense granular flow is set-up by suddenly imparting a constant speed u_w on the wall at t > 0, the internal resistance of the media develops in a small region close to the wall and later it is diffused into the shear layer. As a result, the μ(I)–rheology suggests that, initially when the wall is set into motion, the shear stress experienced by the granular media in the neighbourhood of the velocity driven wall should be sufficiently close to the critical shear stress τ_d = μ_d p_w until some time t_μ d proportional to ν_g (u_w/g δ_s)^2. Thereafter, the wall shear stress τ_w decreases with time as a power law 1/√(t) before reaching its steady state value necessary to support the externally imposed wall speed u_w. At this stage the shear stress becomes uniform in the bulk of the mobile layer.

In both variants of Stokes' problem in granular media, a properly rescaled friction coefficient (or the shear stress) is illustrated to be approximately self-similar with respect to the variable y/Δμ√(ν_g t). Moreover, the steady-state wall speed, wall shear stress and the applied wall pressure are related by a simple approximate expression <ref> which in terms of dimensional parameters can be given as u_w ∼ (g δ_s^2/ν_g) f( β_w ) where f(β_w) is a function of the surplus steady-state wall shear stress β_w = (τ_w - μ_s p_w)/μ_s p_w such that it is either 𝒪(1) when the wall shear stress is just above the yield stress τ_w ∼τ_s or logarithmically large as τ_w approaches τ_d.

Note that when the local friction coefficient μ approaches μ_d, the local Inertial number which compares the timescale of grain-grain rearrangements with that arising from the macroscopic deformation of the media should be very large (I ≫ 1). Therefore, the aforementioned result that, for t ≤ t_μ d, the local shear stress τ = μ p tends towards τ_d = μ_d p in the immediate vicinity of the region where velocity is imposed, suggests that a highly agitated granular flow could occur in this zone. A viscoplastic description of the μ(I)–rheology might, in fact, be not suitable in such zones  as the local shear stress therein cannot be supported by internal grain-grain frictional resistance alone. Here, a proper model for rapid granular flows <cit.> should be more pertinent. Furthermore, if the applied wall speed u_w is much larger than ν_g / g δ_s^2 so that t_μ d becomes sufficiently large, the resulting unsteady granular flow as computed from the μ(I)–rheology may even be incorrect.

It is expected that this study motivates investigations on the further validity of the μ(I)–rheology for unsteady dense granular flows using simple experiments. These results should also be helpful to better understand shear layers, effective viscosity, drag force and characteristic diffusion timescales in future studies with the μ(I)–rheology. Especially in the context of ill-posedness of the μ(I)–rheology as an initial-value problem, it might be essential to identify what features predicted by this rheology are still meaningful. It will be of some interest to include the spatio-temporal variation of the solid fraction ϕ as well – via, for example, a linear function of the inertial number I as in <cit.>. Non-local effects by which a granular media can yield even if the local shear stress is below the yield criterion are omitted in the present short note for sake of simplicity but they might play an important role under common experimental conditions. Finally, it is pointed out that the Stokes' 2nd problem with an oscillating wall boundary condition remains a very interesting open problem as it might shed light upon how static and dynamic zones can simultaneously appear and move around in unsteady flow fields predicted by μ(I)–rheology of dense granular flows. However, the μ(I)–rheology can neither account for the history of the shear stress in the bulk of the media nor consider other static initial conditions that are different from the yield criterion. It is nonetheless important to study these configurations using the μ(I)–rheology to further advance knowledge about continuum models for unsteady dense granular flows.


The authors acknowledge vital inputs from Simon Dagois Bohy and all correspondences with Delphine Doppler and Pierre Jop. 






48
natexlab#1#1#1#1 #1#1 #1#1#1#1#1#1
  #1#1#1#1 #1#1 #1#1
  #1#1#1#1#1#1#1#1

[Ancey & Bates(2017)]Ancey_JNonNewt2017
Ancey, Christophe & Bates, Belinda M. 2017  Stokes’
  third problem for herschel–bulkley fluids.  Journal of Non-Newtonian
  Fluid Mechanics  243 (Supplement C),  27 – 37.

[Ancey et al.(1999)Ancey, Coussot &
  Evesque]ancey1999theoretical
Ancey, Christophe, Coussot, Philippe & Evesque, Pierre
  1999  A theoretical framework for granular suspensions in a steady
  simple shear flow.  Journal of Rheology  43 (6),
  1673–1699.

[Andreotti et al.(2011)Andreotti, Forterre &
  Pouliquen]Andreotti_Forterre_Pouliquen_2011
Andreotti, Bruno, Forterre, Yoël & Pouliquen, Olivier
  2011 Les milieux granulaires: entre fluide et solide.  EDP
  sciences.

[Aranson & Tsimring(2006)]aranson2006patterns
Aranson, Igor S & Tsimring, Lev S 2006  Patterns and
  collective behavior in granular media: Theoretical concepts.  Reviews of
  modern physics  78 (2),  641.

[Baker et al.(2016)Baker, Barker & Gray]baker2016two
Baker, JL, Barker, T & Gray, JMNT 2016  A
  two-dimensional depth-averaged μ(I)–rheology for dense
  granular avalanches.  Journal of Fluid Mechanics  787,
  367–395.

[Balmforth et al.(2009)Balmforth, Forterre &
  Pouliquen]Balforth_JNonNewt2009
Balmforth, Neil J., Forterre, Y. & Pouliquen, O.
  2009  The viscoplastic stokes layer.  Journal of Non-Newtonian
  Fluid Mechanics  158 (1),  46 – 53, visco-plastic fluids: From
  theory to application.

[Barker & Gray(2017)]barker2017partial
Barker, T & Gray, JMNT 2017  Partial regularisation
  of the incompressible μ(I)-rheology for granular flow.
  Journal of Fluid Mechanics  828,  5–32.

[Barker et al.(2015)Barker, Schaeffer, Bohorquez &
  Gray]barker2015well
Barker, T, Schaeffer, DG, Bohorquez, P & Gray, JMNT
  2015  Well-posed and ill-posed behaviour of the
  μ(I)–rheology for granular flow.  Journal of Fluid
  Mechanics  779,  794–818.

[Barker et al.(2017)Barker, Schaeffer, Shearer &
  Gray]barker2017well
Barker, Thomas, Schaeffer, DG, Shearer, Michael &
  Gray, JMNT 2017 Well-posed continuum equations for granular flow
  with compressibility and μ(I)-rheology.  In Proc. R.
  Soc. A, ,  vol. 473,  p. 20160846. The Royal Society.

[Capart et al.(2015)Capart, Hung & Stark]capart2015depth
Capart, H, Hung, C-Y & Stark, CP 2015
  Depth-integrated equations for entraining granular flows in narrow
  channels.  Journal of Fluid Mechanics  765.

[Cawthorn(2011)]cawthorn2011thesis
Cawthorn, Christopher John 2011  Several applications of a
  model for dense granular flows. PhD thesis, University of Cambridge.

[Chauchat & Médale(2014)]chauchat2014three
Chauchat, Julien & Médale, Marc 2014  A
  three-dimensional numerical model for dense granular flows based on the
  μ(I)–rheology.  Journal of Computational Physics
  256,  696–712.

[Da Cruz et al.(2005)Da Cruz, Emam, Prochnow, Roux &
  Chevoir]da2005rheophysics
Da Cruz, Frédéric, Emam, Sacha, Prochnow,
  Michaël, Roux, Jean-Noël & Chevoir, François
  2005  Rheophysics of dense granular materials: Discrete simulation
  of plane shear flows.  Physical Review E  72 (2),  021309.

[Devakar & Iyengar(2008)]devakar2008stokes
Devakar, M & Iyengar, TKV 2008  Stokes’ problems
  for an incompressible couple stress fluid.  Nonlinear Analysis:
  Modelling and Control  1 (2),  181–190.

[Devakar & Iyengar(2009)]devakar2009stokes
Devakar, M & Iyengar, TKV 2009  Stokes’ first
  problem for a micropolar fluid through state-space approach.  Applied
  Mathematical Modelling  33 (2),  924–936.

[Ekinci et al.(2008)Ekinci, Karabacak &
  Yakhot]YakhotOscillatingFlows_PRL2008
Ekinci, K. L., Karabacak, D. M. & Yakhot, V. 2008
  Universality in oscillating flows.  Phys. Rev. Lett.  101,
  264501.













[Forterre & Pouliquen(2008)]forterre2008flows
Forterre, Yoël & Pouliquen, Olivier 2008  Flows
  of dense granular media.  Annu. Rev. Fluid Mech.  40,
  1–24.

[Goddard & Lee(2017)]goddard2017stability
Goddard, JD & Lee, Jaesung 2017  On the stability of
  the μ(I) rheology for granular flow.  Journal of Fluid
  Mechanics  833,  302–331.

[Goldhirsch(2003)]goldhirsch2003rapid
Goldhirsch, Isaac 2003  Rapid granular flows.  Annual
  review of fluid mechanics  35 (1),  267–293.

[Gray & Edwards(2014)]gray2014depth
Gray, JMNT & Edwards, AN 2014  A depth-averaged
  μ(I)–rheology for shallow granular free-surface flows.
  Journal of Fluid Mechanics  755,  503.

[Gray et al.(2003)Gray, Tai & Noelle]gray2003shock
Gray, JMNT, Tai, Y-C & Noelle, S 2003  Shock
  waves, dead zones and particle-free regions in rapid granular free-surface
  flows.  Journal of Fluid Mechanics  491,  161–181.

[Hutter & Rajagopal(1994)]hutter1994flows
Hutter, K & Rajagopal, KR 1994  On flows of granular
  materials.  Continuum Mechanics and Thermodynamics  6 (2),
  81–139.

[Iordanoff & Khonsari(2004)]iordanoff2004granular
Iordanoff, Ivan & Khonsari, MM 2004  Granular
  lubrication: toward an understanding of the transition between kinetic and
  quasi-fluid regime.  Journal of Tribology  126 (1),
  137–145.

[Jaeger et al.(1996)Jaeger, Nagel &
  Behringer]jaeger1996granular
Jaeger, Heinrich M, Nagel, Sidney R & Behringer, Robert P
  1996  Granular solids, liquids, and gases.  Reviews of modern
  physics  68 (4),  1259.

[Jenkins & Savage(1983)]jenkins1983theory
Jenkins, James T & Savage, Stuart B 1983  A theory
  for the rapid flow of identical, smooth, nearly elastic, spherical
  particles.  Journal of fluid mechanics  130,  187–202.

[Jop et al.(2006)Jop, Forterre &
  Pouliquen]jop2006constitutive
Jop, Pierre, Forterre, Yoël & Pouliquen, Olivier
  2006  A constitutive law for dense granular flows.  Nature
  441 (7094),  727–730.

[Jop et al.(2007)Jop, Forterre &
  Pouliquen]Jop2007Stokes3rdProb
Jop, Pierre, Forterre, Yoël & Pouliquen, Olivier
  2007  Initiation of granular surface flows in a narrow channel.
  Physics of Fluids  19 (8),  088102,  arXiv:
  https://doi.org/10.1063/1.2753111.







[Joseph & Saut(1990)]joseph1990short
Joseph, Daniel D & Saut, Jean Claude 1990  Short-wave
  instabilities and ill-posed initial-value problems.  Theoretical and
  Computational Fluid Dynamics  1 (4),  191–227.

[Kamrin(2010)]kamrin2010nonlinear
Kamrin, Ken 2010  Nonlinear elasto-plastic model for dense
  granular flow.  International Journal of Plasticity  26 (2),
  167–188.

[Lagrée et al.(2011)Lagrée, Staron &
  Popinet]lagree2011granular
Lagrée, P-Y, Staron, Lydie & Popinet, Stéphane
  2011  The granular column collapse as a continuum: validity of a
  two-dimensional navier–stokes model with a μ(I)–rheology.
  Journal of Fluid Mechanics  686,  378–408.

[Liu & Nagel(1998)]liu1998nonlinear
Liu, Andrea J & Nagel, Sidney R 1998  Nonlinear
  dynamics: Jamming is not just cool any more.  Nature
  396 (6706),  21–22.

[Martin et al.(2017)Martin, Ionescu, Mangeney, Bouchut &
  Farin]martin2017continuum
Martin, Nathan, Ionescu, IR, Mangeney, Anne, Bouchut,
  François & Farin, Maxime 2017  Continuum viscoplastic
  simulation of a granular column collapse on large slopes: μ(I)
  rheology and lateral wall effects.  Physics of Fluids  29 (1),
  013301.

[MiDi(2004)]midi2004dense
MiDi, GDR 2004  On dense granular flows.  European
  Physical Journal E–Soft Matter  14 (4).

[Morrison(1956)]morrison1956wave
Morrison, JA 1956  Wave propagation in rods of voigt
  material and visco-elastic materials with three-parameter models.
  Quarterly of Applied Mathematics  14 (2),  153–169.

[Panton(1968)]panton1968transient
Panton, Ronald 1968  The transient for stokes's oscillating
  plate: a solution in terms of tabulated functions.  Journal of Fluid
  Mechanics  31 (4),  819–825.

[Preziosi & Joseph(1987)]preziosi1987stokes
Preziosi, Luigi & Joseph, Daniel D 1987  Stokes'
  first problem for viscoelastic fluids.  Journal of non-newtonian fluid
  mechanics  25 (3),  239–259.

[Pritchard et al.(2011)Pritchard, McArdle &
  Wilson]Pritchard_JNonNewt2011
Pritchard, David, McArdle, Catriona R. & Wilson,
  Stephen K. 2011  The stokes boundary layer for a power-law fluid.
   Journal of Non-Newtonian Fluid Mechanics  166 (12),  745 –
  753.

[Savage(1984)]savage1984mechanics
Savage, Stuart B 1984  The mechanics of rapid granular
  flows.  Advances in applied mechanics  24,  289–366.

[Schlichting(1968)]schlichting1968boundary
Schlichting, Hermann 1968 Boundary-layer theory.
  McGraw-Hill.

[Staron et al.(2012)Staron, Lagrée &
  Popinet]staron2012granular
Staron, Lydie, Lagrée, P-Y & Popinet, Stéphane
  2012  The granular silo as a continuum plastic flow: The hour-glass
  vs the clepsydra.  Physics of Fluids  24 (10),  103301.

[Stokes(1851)]stokes1851effect
Stokes, George Gabriel 1851 On the effect of the internal
  friction of fluids on the motion of pendulums, ,  vol. 9.
  Pitt Press Cambridge.






[Tanner(1962)]tanner1962note
Tanner, Roger I 1962  Note on the rayleigh problem for a
  visco-elastic fluid.  Zeitschrift für Angewandte Mathematik und
  Physik (ZAMP)  13 (6),  573–580.

[Thompson & Grest(1991)]Thompson_PRL1991granular
Thompson, Peter A & Grest, Gary S 1991  Granular
  flow: friction and the dilatancy transition.  Physical Review Letters
  67 (13),  1751.

[Yakhot & Colosqui(2007)]YakhotColosqui_JFM2007
Yakhot, Victor & Colosqui, Carlos 2007  Stokes'
  second flow problem in a high-frequency limit: application to nanomechanical
  resonators.  Journal of Fluid Mechanics  586,  249–258.



