Dilatancy and its coupling to the kinematics in sheared granular media

Abstract

Models for slow flow of dense granular materials often treat the medium as incompressible, thereby neglecting the role of Reynolds dilatancy. However, recent particle simulations have demonstrated the presence of a significant coupling between the volume fraction and velocity fields. The model of Dsouza & Nott (J. Fluid Mech., vol. 888, 2020, R3) incorporates dilatancy and captures the coupling, but it has thus far lacked experimental validation. In this paper, we provide the first experimental demonstration of dilatancy and its coupling to the kinematics in a two-dimensional cylindrical Couette cell. We find a shear layer near the inner cylinder within which there is significant dilation. Within the shear layer, the azimuthal velocity decays roughly exponentially and the volume fraction rises with radial distance from the inner cylinder. The predictions of the model of Dsouza & Nott (2020) are in good agreement with the experimental data for a variety of roughness features of the outer cylinder. Moreover, by comparing the steady states resulting from different initial volume fraction profiles (but having the same average), we show the inter-dependence of the velocity and volume fraction fields, as predicted by the model. Our results establish the importance of shear dilatancy even in systems of constant volume.

1. Introduction

The formulation of reliable and robust continuum models for the flow of granular materials has been an important endeavour, owing to the widespread occurrence of granular flows in industrial processes and natural phenomena. While particle cohesion and complex shape are complicating factors in many practical systems, the simpler case of cohesionless granular media composed of roughly isotropic particles is still sufficiently prevalent that modelling them is informative and, to a reasonable extent, generalisable. Early models have focused on two limiting regimes of slow and rapid flow; the flow regime is determined by the Savage number (Savage & Hutter 1989) or, equivalently, the inertia number $I$,

(1.1)\begin{equation} {\mathit{Sa}} \equiv \frac{\rho_p d_p^2 \dot{\gamma}^2}{N} = I^2, \end{equation}

which characterises the contribution of particle inertia to the stress. Here, $\rho _p$ and $d_p$ are the density and mean diameter of the particles, $\dot {\gamma }$ is the nominal shear rate, and $N$ a stress scale (typically the confining stress). Slow flow corresponds to the limit $Sa \ll 1$, where particle inertia is of no consequence, and rapid flow to $Sa \sim 1$, where particle inertia is dominant. Models for slow flow are based on plasticity theories for soils (Rao & Nott 2008), and models for rapid flows are based on the kinetic theory of dense gases extended to account for inelasticity of particle collisions (Lun et al. 1984). The phenomenological $\mu (I)$ model (Jop, Forterre & Pouliquen 2006), an extension of a simple plasticity model wherein the friction coefficient is assumed to depend on $I$, has been used to model the intermediate regime between slow and rapid flows. The recent paper of Berzi (2024) shows that the $\mu (I)$ model can be derived from kinetic theory if spatial correlation of particle velocity fluctuations is accounted for.

In this paper, we restrict attention to the slow-flow regime. Classical plasticity models used for this regime suffer from two major deficiencies. The first is kinematic indeterminacy, meaning that the deformation rate (and therefore the velocity) cannot be uniquely determined if the stress field is specified. This is a consequence of the models being designed to incorporate rate-independence of the stress, an experimentally observed feature of slow flow. The second deficiency is their inability to incorporate shear dilatancy or volume deformation caused by shear deformation; this feature, first observed by Reynolds (1885), is peculiar to granular materials and has no analogue in fluids. Dilatancy alters the particle volume fraction field $\phi$, which in turn affects the stress and therefore the velocity field. Determining the change in $\phi$ due to dilatancy is therefore important in practical problems such as determining the drag on an intruder (Schröter et al. 2007).

A few models have successfully resolved the kinematic indeterminacy of classical plasticity models. Mohan, Nott & Rao (1999) and Mohan, Rao & Nott (2002) treated granular media as Cosserat continua, in which stress symmetry is not presumed, thereby requiring that the angular momentum balance be enforced along with the balances of linear momentum and mass. However, asymmetry of the stress has thus far not been experimentally verified. Another class of models introduce a scalar order parameter that is governed by a diffusion equation (Aranson & Tsimring 2002; Bouzid et al. 2013; Henann & Kamrin 2013) in the constitutive relation for the stress. These models are phenomenological, with the order parameter chosen by analogy with other physical systems; moreover, the boundary conditions for the order parameter are speculative. Despite their different physical origins, all the models mentioned above share the common feature of introducing a mesoscopic length scale $\ell$, whose effect is to allow the material to deform to a distance O($\ell$) from the point where the yield condition ceases to be satisfied. For this reason, they are called non-local models, though they are not non-local in the formal, mathematical sense (Eringen 1983). Although these models repair kinematic indeterminacy, and their predictions of the velocity profile have been found to be in good agreement with experimental data (Tang et al. 2018; Fazelpour, Tang & Daniels 2022; Fazelpour & Daniels 2023), they all assume the material to be incompressible. As a result, they do not capture shear dilatancy and thereby its effect on the kinematics.

A coupling between the $\phi$ and $\boldsymbol {u}$ fields was observed in the particle dynamics simulations of Krishnaraj & Nott (2016) and Dsouza & Nott (2021), who show that dilatancy in conjunction with gravity drives large-scale secondary flows in cylindrical Couette and split-bottom shear cells. However, very few experimental studies have measured the variation of $\phi$ in three-dimensional (3-D) sheared granular media in the slow-flow regime, as non-invasive imaging of opaque media requires relatively specialised and expensive facilities such as X-ray computed tomography (X-ray CT) and magnetic resonance imaging (MRI). The only study that we are aware of that measured the density field is that of Mueth et al. (2000), who imaged a cylindrical Couette shear cell using a specialised high-speed MRI facility. While they do find dilation in the shear layer, they report data only very close to the rotating inner cylinder.

Recently, Dsouza & Nott (2020) proposed a model that addresses both the shortcomings of classical plasticity, namely kinematic indeterminacy and absence of dilatancy, by a systematic non-local extension of the critical state plasticity theory. The model poses the flow rule, and the relation between volume fraction and the critical state pressure in a non-local sense without introducing additional variables or equations. Its predictions for plane shear flow were found to be in good agreement with 3-D discrete particle simulations. However, the model predictions have thus far not been validated by experiments.

In this paper, we provide the first experimental measurement of the coupling between the volume fraction and velocity fields in the slow-flow regime, and thereby a verification of the model of Dsouza & Nott (2020), by simultaneously measuring the area fraction and velocity fields in a two-dimensional (2-D) cylindrical Couette shear cell. By imaging a horizontal layer of circular disks sheared between two concentric cylinders, as developed by Tang et al. (2018), Fazelpour et al. (2022) and Fazelpour & Daniels (2023), we measure directly the variation of the azimuthal velocity $u_\theta$ and $\phi$ with radial distance $r$. While measurements in a 2-D monolayer are not expected to correspond quantitatively with those in three dimensions, it has been shown in several studies that the qualitative features of the kinematics and even the statistics of particle velocity (Ananda, Moka & Nott 2008) in two and three dimensions are very similar. By determining the steady-state profiles $u_\theta (r)$ and $\phi (r)$ for different initial $\phi$ profiles (but having the same spatial average), we demonstrate the coupling between the two fields. We compare our experimental data with the predictions of the non-local model of Dsouza & Nott (2020) and find good agreement.

2. Experimental method

Our experiments are conducted in the 2-D cylindrical Couette rheometer shown in figure 1 described by Fazelpour & Daniels (2023). The granular material is confined between coaxial inner and outer cylinders of radii $r_{i}=16$ cm and $r_o=27$ cm, respectively. The inner cylinder has semi-circular cavities of diameter $d_p$ so that disks settle into them and are carried by the cylinder; this generates sufficient traction to shear the assembly. The inner cylinder is rotated by a motor at a constant angular speed $\varOmega$ such that its azimuthal velocity of $u_{w} \equiv r_i \varOmega = 1.1\ {\rm cm}\ {\rm s}^{-1}$. Each experiment was conducted using an outer cylinder with a specified pattern of roughness (figure 1): two with cavities, two with convex bumps (protrusions) and one made of compliant leaf springs. The dimensions of the roughness features of all the outer cylinders are given in the caption of figure 1.

Figure 1. Top view of 2-D cylindrical Couette device. The inner cylinder of radius $r_{i}=16$ cm is rotated at constant angular speed $\varOmega$ and the outer cylinder of radius $r_{o}=27$ cm is stationary. The inner cylinder is machined to have semi-circular cavities of diameter $d_p$. The outer cylinder is changeable, and five different roughness designs were used: semi-circular bumps of diameter $d_p$, hereafter called small bumps (SB); bumps of diameter $8 d_p$ and chord length $5.4 d_p$, hereafter called large bumps (LB); semi-circular cavities of diameter $d_p$, hereafter called small cavities (SC); cavities of diameter $8 d_p$ and chord length $5.4 d_p$, hereafter called large cavities (LC); constructed with leaf springs (LS) (see Fazelpour & Daniels 2023) with the same roughness dimensions as the LB cylinder.

The granular material is composed of polyurethane disks (Vishay Precision Group, PhotoStress material PSM-4) of Young's modulus $E=4$ MPa and density ${\rho =1.06\ {\rm g}\ {\rm m}^{-3}}$ (Fazelpour & Daniels 2023). The disks are bidisperse in size, with equal numbers of disks of diameter $0.9d_p$ and $1.1d_p$, where the mean diameter $d_p = 1$ cm. The disks are of thickness 6.35 mm and rest on a flat acrylic base while being sheared between the cylinders. They are placed by hand randomly and distributed as uniformly as possible within the Couette gap. For the 2-D monolayer of disks, the dimensionless density is characterised by the area fraction; hereafter, we use $\phi$ to represent the area fraction of the disks. The roughness of a boundary typically determines the shear rate and velocity slip adjacent to it (Ananda et al. 2008; Dsouza & Nott 2020). While they are small near the outer cylinder, due to the roughly exponentially decay of $u_{\theta }$ with distance from the inner cylinder, there is a measurable difference in their values for the different outer cylinders. The boundary roughness also plays a role in another way: the initial volume fraction profiles $\phi (r)$ are not identical for the different boundaries. We shall see below that these two effects yield measurably different steady-state profiles of $u_{\theta }$ and $\phi (r)$ in the shear layer near the inner cylinder.

Because the net area available for the disks varies slightly for the different outer cylinders, the number of particles is adjusted to maintain an average area fraction of $\bar {\phi } = 0.76$. For the outer cylinder composed of LS, the area fraction cannot be held precisely constant due to their compliance; the number of disks used is equal to that for the outer cylinder with LB, which have the same shape.

We track the motion of all particles by imaging the Couette cell by a video camera mounted above the apparatus. The movies are taken at $30\ {\rm frames}\ {\rm s}^{-1}$ for a duration of 10 min (7 revolutions of the inner cylinder). Each frame in a movie is processed to identify the circular disks and determine their centres $\boldsymbol {x}_i$. The velocities of the disks are computed from their displacement between two successive frames, $\boldsymbol {u}_i = \Delta \boldsymbol {x}_i/\Delta t$. As the flow is axisymmetric and steady, continuum estimates of $u_\theta$ and $\phi$ at radial position $r$ are determined by averaging the particle velocities and their areas within an annular bin of width $\Delta r = 1 d_p$ (centred at $r$) and azimuthal span $0 \leqslant \theta \leqslant {\rm \pi}/2$, and further averaging over a time interval of 200 s (6000 frames). The profiles of $u_\theta$ and $\phi$ are obtained by moving the annular bin by small radial increments to generate data at 48 radial positions. We note that such fine spatial resolution is achieved because we average over large intervals of time and $\theta$, which may not be possible when the flow is unsteady or not fully developed.

The stress field $\boldsymbol {\sigma }(r)$ was measured for the same set-up and all the outer cylinders by Fazelpour & Daniels (2023). As discussed below, we use their measurements of the pressure to determine one of the parameters in the continuum model.

3. The non-local model of Dsouza & Nott

The model of Dsouza & Nott (2020) is based on the idea that the deformation rate $\boldsymbol{\mathsf{D}}$ at a point depends on the stress $\boldsymbol {\sigma }$ not just at that point, but in a mesoscopic volume surrounding it, as a result of the stress and strain rate being correlated over a mesoscale (such as through force chains). Mathematically, this idea translates to a non-local flow rule (the relation between the strain rate and the stress) and similarly between the volume fraction and the critical state pressure. A mesoscopic length $\ell$, a material parameter, characterises the extent of non-locality. The constitutive relation for $\boldsymbol {\sigma }$ resulting from these arguments is

(3.1a)$$\begin{gather} \boldsymbol{\sigma} ={-}p \boldsymbol{\delta} + \frac{2\mu}{\dot{\gamma}}(p_c \boldsymbol{\mathsf{D}}'-\ell^2 \varPi \nabla^2 \boldsymbol{\mathsf{D}}'), \end{gather}$$

(3.1b)$$\begin{gather}p = p_c \left(1- \frac{\mu _b}{\dot{\gamma}} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}\right) + \ell^2 \varPi \frac{\mu _b}{\dot{\gamma}}\nabla ^2 \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}, \end{gather}$$

(3.1c)$$\begin{gather}p_c =\varPi - \ell^2 \frac{{\rm d}\varPi}{{\rm d}\phi}\nabla ^2 \phi, \end{gather}$$

where $\boldsymbol{\mathsf{D}}' \equiv \boldsymbol{\mathsf{D}} - \frac {1}{3}\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} \, \boldsymbol {\delta }$ is the deviatoric part of $\boldsymbol{\mathsf{D}}$ and $\dot {\gamma } \equiv (2\boldsymbol{\mathsf{D}}' : \boldsymbol{\mathsf{D}}')^{1/2}$ is its scalar invariant, $p$ is the pressure and $p_c$ the pressure at the critical state. The terms multiplied by $\ell ^2$ are the non-local contributions to the constitutive relation. The material properties are the friction coefficients $\mu$ and $\mu _b$ for shear and volume deformation, respectively, the mesoscopic length $\ell$ and the local form of the critical state pressure $\varPi (\phi )$.

We now apply the non-local model (3.1) to shear in a 2-D cylindrical Couette cell, for which the deviatoric deformation rate tensor is $\boldsymbol{\mathsf{D}}' \equiv \boldsymbol{\mathsf{D}} - \frac {1}{2}\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} \, \boldsymbol {\delta }$ and $\phi$ represents the area fraction. The governing equations are the balances of mass, and the $r$ and $\theta$ components of momentum. For steady axisymmetric flow, the mass balance is trivially satisfied, and the $r$ and $\theta$ momentum balances reduce to

(3.2a,b)\begin{equation} \partial \sigma_{rr}/\partial r = 0, \quad (1/r^2) \partial/\partial r (r^2 \sigma_{r \theta}) +\mu_{base}\, \phi \rho_p g = 0. \end{equation}

The second term in (3.2b) is the frictional resistance offered by the base. On substituting the constitutive relation for $\boldsymbol {\sigma }$ (3.1), the momentum balances (3.2) take the form

(3.3a)$$\begin{gather} \frac{{\rm d}}{{\rm d}r} \left[ \varPi - l^2\, \frac{{\rm d}\varPi}{{\rm d} \phi} \left( \frac{{\rm d}^2 \phi}{{\rm d}r^2} + \frac{1}{r} \frac{{\rm d} \phi}{{\rm d} r} \right) \right] = 0, \end{gather}$$

(3.3b)$$\begin{gather} \frac{1}{r^2} \frac{{\rm d}}{{\rm d} r} \left[ r^2 \frac{\mu p_c}{\dot{\gamma}} \left(\frac{{\rm d}u_{\theta}}{{\rm d}r} - \frac{u_{\theta}}{r} \right)- r^2 \ell^2 \, \frac{\mu \varPi}{\dot{\gamma}} \left(\frac{{\rm d}^3 u_{\theta}}{{\rm d}r^3} + \frac{1}{r^2} \frac{{\rm d}u_{\theta}}{{\rm d}r} - \frac{u_{\theta}}{r^3} \right) \right] - \mu_{base}\, \phi \rho_pg = 0, \end{gather}$$

where $\dot {\gamma } = |{\rm d}u_{\theta }/{\rm d}r - u_{\theta }/r|$.

To close (3.3), four boundary conditions for $u_{\theta }$ and three for $\phi$ are required. Unlike fluids, granular media usually slip at rigid boundaries (Natarajan, Hunt & Taylor 1995; Ananda et al. 2008), and accounting for slip is often important. A boundary condition that accommodates slip was proposed by Mohan et al. (2002) for their Cosserat continuum model in which the slip velocity is $K d_p \boldsymbol {n} \times \boldsymbol {\omega }$, where $K$ is the slip coefficient and $\boldsymbol {\omega }$ is the rate of material spin at the boundary of unit normal $\boldsymbol {n}$. For a classical continuum, $\boldsymbol {\omega }$ is equal to half the vorticity $\boldsymbol {w}$; the two fields can differ in Cosserat continua and the difference is proportional to the asymmetry of the Cauchy stress (Dahler 1959; Lun 1991; Mohan et al. 2002). As there is little experimental evidence for asymmetry of the Cauchy stress in granular media, we proceed with the assumption that the granular medium is a classical continuum, and hence $\boldsymbol {\omega } = \boldsymbol {w}$. Using this condition and the wall friction boundary condition (Rao & Nott 2008), we have

(3.4a,b)\begin{equation} u_{\theta} - u_{w} = K d_p \left( \frac{{\rm d}u_{\theta}}{{\rm d} r} + \frac{u_{\theta}}{r} \right), \quad \frac{-\sigma_{r \theta}}{\sigma_{rr}} = \mu_w \quad \mbox{at}\ r = r_i, \end{equation}

where $\mu _w$ is the friction coefficient between the inner wall and the granular material. At steady state, the shear stress $\sigma _{r\theta }$ decreases with $r$ while the normal stress $\sigma _{rr}$ is constant, whence $|\sigma _{r\theta }/\sigma _{rr}|$ decreases with $r$. For large enough Couette gap, non-local effects cannot transmit shear to the outer wall, whence the shear rate there must vanish (Mohan et al. 2002; Dsouza & Nott 2020). Together with the slip condition (3.4a), we get the boundary conditions

(3.5a,b)\begin{equation} u_{\theta}=0, \quad \frac{{\mathrm d} u_{\theta}}{{\mathrm d} r} = 0 \quad \mbox{at}\ r = r_o. \end{equation}

The value of the slip coefficient $K$ in (3.4a) will depend on the nature of the boundary; in studies where the boundaries are topographically roughened, such as by gluing coarse sandpaper to the wall, the slip velocity has been found to be small, implying $K = 0$. For a given wall, $K$ must be determined by fitting with experimental data for the velocity. Jing et al. (2016) correlated the slip velocity in an inclined chute with the roughness parameters of the base using discrete element method simulations; generalising such a study to determine $K$ as a function of the wall roughness features and $\phi$ would be a useful exercise.

For the area fraction $\phi$, one obvious condition is that the mean solids fraction $\bar {\phi }$, or alternatively the normal stress at one of the boundaries, must be specified. Apart from that, we do not have other well-established boundary conditions for $\phi$ at the walls. We follow Dsouza & Nott (2020) and assume that $\phi$ at the two boundaries are known; we discuss this choice in § 5. As a result, we have the following conditions for $\phi$:

(3.6a–c)\begin{equation} \int_{r_i}^{r_o} \phi 2 {\rm \pi}r \,{\mathrm d} r = \bar{\phi} \, {\rm \pi}(r^2_o - r^2_i), \quad \phi (r_i) = \phi_i, \quad \phi (r_o) = \phi_o. \end{equation}

Equations (3.3)–(3.6) constitute a well-posed boundary value problem. Note that the momentum balances (3.3a) and (3.3b) are uncoupled; solution of the former yields the area fraction profile $\phi (r)$ and the latter then yields the velocity profile $u_\theta (r)$. The equations are solved by the finite difference method using second-order discretisation. The resulting nonlinear algebraic equations for $u_\theta$ and $\phi$ at the $N$ discretised values of $r$ are solved using the nonlinear equation solver fsolve in MATLAB. The results reported in § 4 are for $N = 100$; the solutions remain unchanged (to within a relative tolerance of $10^{-3}$) for larger $N$.

4. Results and discussion

For numerical solution of the governing equations (3.3)–(3.6), the values of parameters that characterise the granular material and the boundaries are needed. To solve (3.4a), we need the form of the function $\varPi (\phi )$ in the expression for the critical state pressure (3.1c); as done by Dsouza & Nott (2020), we assume the form proposed by Nott & Jackson (1992),

(4.1)\begin{equation} \varPi(\phi) = \left\{ \begin{array}{@{}ll@{}} \alpha \dfrac{(\phi - \phi_{min})^{n_1}}{(\phi_{max} - \phi)^{n_2}}, & \phi \geqslant \phi_{min},\\ 0, & \phi < \phi_{min}, \end{array} \right. \end{equation}

where $\phi _{min} = 0.65$ and $\phi _{max} = 0.82$ are the area fractions for loose and dense random packing, respectively. For a 2-D layer of particles, we expect the exponents $n_1$ and $n_2$ to differ from the values used for 3-D assemblies by Nott & Jackson (1992) and Dsouza & Nott (2020). We determine them by minimising the error between the model predictions and the experimental data for $\phi (r)$ for the outer cylinder with LB (see figure 1). The constant $\alpha$ is estimated from the measurement by Fazelpour & Daniels (2023) of the pressure $p \approx 10^3$ Pa for $\bar {\phi } = 0.76$, which yields the value $\alpha = 30$ Pa. To solve (3.4b), we need values for the three friction coefficients, the slip coefficient $K$ and the mesoscopic length $\ell$. We retain the value $\ell = 10$ used in all our previous studies (Mohan et al. 2002; Ananda et al. 2008; Dsouza & Nott 2020). By integrating (3.3b) once from $r_i$ to $r$ and using the boundary condition (3.4b), we see that the friction coefficients occur only as the ratios $\mu _w/\mu$ and $\mu _{base}/\mu$. These ratios and $K$ are determined by minimising the error between the model predictions and the experimental data for $u_\theta (r)$ for the LB outer cylinder. Details of parameter estimation by minimising the error are given in the Appendix. The best-fit set of model parameters are $n_1 = 2$, $n_2 = 1$, $\mu _w/\mu = 1$, $\mu _{base}/\mu = 3.76$ and $K = 0.63$. The parameters are now fixed at these values and used for comparison against all the other experimental data sets.

The experimental data for the LB outer cylinder are shown in figure 2, along with the model predictions using the best-fit parameter values. A notable feature of the experimental data is that the initial area fraction profile (i.e. before shearing has commenced) has significant spatial fluctuations about the mean $\bar {\phi } = 0.76$. They arise because the disks are placed in the Couette gap by hand, whence some non-uniformity is unavoidable. Moreover, $\phi$ is lower near the inner and outer cylinders due to packing constraints near rigid boundaries. As the material is sheared, we observe the development of pronounced dilation (reduction of $\phi$) in the shear layer near the inner cylinder; concurrently, the fluctuations in $\phi$ become smoother and ultimately vanish at steady state. Conservation of mass leads to some compaction at larger radial positions, but the fluctuations in $\phi$ in the outer half of the Couette gap remain even after prolonged shear. This is because the velocity and shear rate decay roughly exponentially with distance from the inner cylinder (figure 2(b) and inset of figure 4b), thereby causing little rearrangement of the packing away from the inner cylinder.

Figure 2. Experimental data and model predictions for the outer cylinder with LB, see figure 1. ($a$) Data for the area fraction $\phi (r)$ before shear is commenced and at steady state compared with the model prediction using the best-fit parameters values $n_1 = 2$ and $n_2 = 1$ and $\alpha = 30$ Pa. ($b$) Data for the velocity $u_\theta (r)$ at steady state compared with the model predictions using the best-fit parameters values $\mu _w/\mu = 0.99$, $\mu _{base}/\mu = 3.76$ and $K = 0.63$. The best-fit parameters are obtained by minimising the error between the model predictions and the data (see the Appendix).

A brief discussion on how dilatancy arises in the model is pertinent. At the steady state, which in the present problem is the critical state, the momentum balance in the $r$ direction (3.3a) requires the critical state pressure $p_c$ to be constant across the Couette gap. The non-local relation between $p_c$ and $\phi$ results in the spatial variation in $\phi$ across the Couette gap such that $\phi$ is lower in the shear layer, i.e. dilatancy. In the absence of the non-local contribution to $p_c$, namely the second term on the right-hand side of (3.1c), $\phi$ is constant and there is no dilatancy. The physical mechanism for the evolution of $\phi$ from the initial state, by depletion of particles from the shear layer and accumulation further away, becomes clear when we consider the unsteady development of the flow fields in § 4.1.

As described in § 2, experiments were conducted with outer cylinders of different roughness patterns. Figure 3 compares the area fraction and velocity profiles and the corresponding model predictions (solid lines) for all the cylinders. We observe that the profiles of $\phi$ and $u_\theta$ in all the cases are quite similar: there is substantial dilation near the inner cylinder (figure 3a), where the velocity gradient is high (figure 3b), and the velocity profiles decay roughly exponentially.

Figure 3. Data for the five different outer cylinders (see figure 1) and the corresponding model predictions. ($a$) Area fraction profiles, with the mean area fraction being $\bar {\phi } = 0.76$ in all the cases. ($b$) Azimuthal velocity profiles scaled by the velocity of the inner cylinder. The solid lines are the model predictions and the symbols with error bars are experimental data.

Nevertheless, there are slight variations between the $\phi (r)$ and $u_\theta (r)$ profiles for outer cylinders of different roughness, possibly due to the differences in the initial area fraction profile. Since the shear rate is very small (${\sim }10^{-5}\ {\rm s}^{-1}$) in the outer portion of the Couette gap, it is likely that the small differences in the mean area fraction within the shear layer remain even at steady state. Serendipitously, these differences turn out to be useful in demonstrating a key aspect of the model, namely the coupling between the packing fraction and velocity fields. Figure 4 compares the steady-state profiles for the outer cylinder with LS and LB. Despite the error bars, a small but systematic difference in the steady-state area fraction profiles is evident. For the cylinder with LS, $\phi$ is lower near the inner cylinder and higher at larger $r$. Figure 4($b$) shows that this results in a measurable difference in the $u_\theta$ profiles, which are in agreement with the model predictions. To further establish the coupling, we have shown in figure 4($b$) the velocity profile obtained by assuming the area fraction to be constant, i.e. $\phi (r) = \bar {\phi }$; it is clear the velocity profile differs substantially from those corresponding to the dilated steady states. Thus, the data show clearly the inter-dependence of the packing fraction and velocity profiles, which was demonstrated earlier in Dsouza & Nott (2020) using particle simulations.

Figure 4. Experimental data (symbols) for ($a$) area fraction and ($b$) azimuthal velocity $u_\theta$ for the outer cylinders with LS and LB compared with the model predictions (solid lines). The mean area fraction in both the cases is $\bar {\phi } = 0.76$. The blue dashed line in panel ($b$) is the model prediction obtained by assuming the area fraction to be constant, i.e. $\phi (r) = \bar {\phi }$. The inset in panel ($b$) shows the velocity on a logarithmic scale; data for $y > 0.5$ are not shown, as $u_\theta /u_w$ approaches the measurement resolution (${\approx }10^{-3}$).

4.1. Unsteady evolution of the packing fraction and velocity fields

Thus far, we have shown the predictions of the model for the steady state and compared them with the experimental data. As already noted, for the steady state, momentum balances in the $r$ and $\theta$ directions are uncoupled, allowing the determination of $\phi$ first and then $u_\theta$. We now consider the unsteady evolution of the flow fields from the initial un-sheared state, where the packing fraction and velocity fields evolve together, and the mechanism in the model that captures dilatancy becomes much clearer.

We expect the presence of a radial velocity $u_r$ and therefore a non-zero bulk deformation rate $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = \partial u_r/\partial r + u_r/r$. Discarding the inertial terms in the unsteady momentum balances, as our interest is in the slow-flow regime, the equations of motion for axisymmetric unsteady flow are

(4.2a)\begin{align} &\frac{\partial \phi}{\partial t} + \frac{1}{r} \frac{\partial}{\partial r} (r \phi u_r)=0, \end{align}

(4.2b)\begin{align} &- \frac{\partial p}{\partial r} + \left(\frac{2}{r} + \frac{\partial}{\partial r} \right) \left[ \frac{\mu p_c}{\dot{\gamma}} \left(\frac{\partial u_r}{\partial r} - \frac{u_r}{r} \right) - \ell^2 \, \frac{\mu \varPi}{\dot{\gamma}} \left(\frac{\partial^3 u_r}{\partial r^3} + \frac{1}{r^2} \frac{\partial u_r}{\partial r} - \frac{u_r}{r^3} \right) \right] \nonumber\\ &\quad- \mu_{base}\, \phi \rho_pg \frac{u_r}{(u_r^2 + u_\theta^2)^{1/2}} = 0, \end{align}

(4.2c)\begin{align} &\frac{1}{r^2} \frac{\partial}{\partial r} \left[ r^2 \frac{\mu p_c}{\dot{\gamma}} \left(\frac{\partial u_{\theta}}{\partial r} - \frac{u_{\theta}}{r} \right) - r^2 \ell^2 \, \frac{\mu \varPi}{\dot{\gamma}} \left(\frac{\partial^3 u_{\theta}}{\partial r^3} + \frac{1}{r^2} \frac{\partial u_{\theta}}{\partial r} - \frac{u_{\theta}}{r^3} \right) \right] \nonumber\\ &\quad - \mu_{base}\, \phi \rho_pg \frac{u_\theta}{(u_r^2 + u_\theta^2)^{1/2}} = 0, \end{align}

with the pressure $p$ given by (3.1b) and $\dot {\gamma } = [ (\partial u_{\theta }/\partial r - u_{\theta }/r)^2 + (\partial u_r/\partial r - u_r/r)^2 ]^{1/2}$.

Since the momentum balances are quasistatic, solution of (4.2) requires an initial condition only for $\phi$, for which we use the initial distribution $\phi _0(r)$,

(4.3)\begin{equation} \phi(0,r) = \phi_0(r). \end{equation}

Additionally, we require four boundary conditions for $u_r$. Two obvious conditions are no penetration at the walls. For the two additional conditions, we propose that the ratio of dilation rate normal to a wall to the tangential shear rate is proportional to the departure from the critical state,

(4.4)\begin{equation} (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u})_n /\dot{\gamma}_t = \beta (\phi - \phi_c), \end{equation}

where $(\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u})_n \equiv \boldsymbol {n} \boldsymbol {\cdot } \boldsymbol{\mathsf{D}} \boldsymbol {\cdot } \boldsymbol {n}$ and $\dot {\gamma }_t \equiv |\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol{\mathsf{D}}' \boldsymbol {\cdot } \boldsymbol {t}|$, $\boldsymbol {n}$ and $\boldsymbol {t}$ are the unit normal and unit tangent in the flow direction, $\phi _c$ is the packing fraction at the critical state (which here is the steady state), and $\beta$ is a constant. As $\dot {\gamma }_t$ vanishes at the outer cylinder (see (3.5)), (4.4) implies that the dilation rate also vanishes. The boundary conditions for $u_r$ then are

(4.5a,b)$$\begin{gather} u_r = 0, \quad \frac{2}{|\partial u_{\theta}/\partial r - u_{\theta}/r|} \frac{\partial u_r}{\partial r} = \beta (\phi - \phi_i) \quad \mbox{at}\ r = r_i, \end{gather}$$

(4.5c,d)$$\begin{gather} u_r = 0, \quad \frac{\partial u_r}{\partial r} = 0 \quad \mbox{at}\ r = r_o. \end{gather}$$

For the reasons mentioned at the end of this section, we do not attempt to determine $\beta$ by fitting with experimental data; however, we find that $\beta = 1$ gives a reasonable match with experiments for the time scale to approach steady state.

Equations (4.2) with initial condition (4.3) and boundary conditions (3.4)–(3.6) and (4.5) form a well-posed problem. They are solved numerically in the same manner as in § 3. As an illustrative example, we determine the unsteady solution starting from an initial state corresponding to the data for the outer cylinder with large bumps, shown in figure 2($a$). We take $\phi _0(r)$ to be a smoothed fit of the data, shown in the inset of figure 5($a$). The model prediction of the unsteady evolution of the $\phi$ and velocity fields are shown in figure 5. The mechanism for dilatancy now becomes clear – we see that a radial velocity $u_r$ develops in the transient state (figure 5$b$), which convects particles away from the region of high shear rate, thereby reducing $\phi$ near the inner cylinder (figure 5$a$). The radial velocity decays with time, and at steady (critical) state, $\phi (r)$ is determined by the constraint of constant pressure (3.2a). Thus, the non-local contribution to $p_c$ in (3.1c) is essential for the non-uniform distribution of $\phi$. As dilatancy proceeds, we see in figure 5($c$) that the sharpness in the decay of the primary velocity $u_\theta$ increases. Thus, the unsteady evolution illustrates the physical mechanism in the model that leads to dilatancy via convection of particles away from the shear layer, and how that couples with the velocity field.

Figure 5. Unsteady evolution of the packing fraction and velocity profiles with dimensionless time $\hat {t} \equiv t u_w/(r_o - r_i)$. The flow is taken to be axisymmetric. ($a$) Evolution of the $\phi$ profile. The inset shows the initial condition $\phi _0(r)$, which is a smoothed fit of initial state in the experiment for the LB outer cylinder (see figure 2$a$). ($b$) Evolution of the radial velocity $u_r$; note that $u_r$ is everywhere positive, largest near (but not at) the inner cylinder and decays to zero at large $\hat {t}$. ($c$) Evolution of the azimuthal velocity $u_\theta$; the decay with distance from the inner cylinder becomes more rapid with time, due to depletion of particles near the inner cylinder.

We are unable to make a comparison of the transient profiles in figure 5 with experimental data, as the unsteady data are noisy. As explained in § 2, obtaining smooth profiles with reasonable spatial resolution requires averaging over a long period of time, which is not possible in the transient state.

5. Discussion and conclusion

We have presented experimental data for the profiles of the area fraction $\phi$ and the azimuthal velocity $u_\theta$ for the shear of disks in a 2-D cylindrical Couette cell. At steady state, we find the velocity to decay roughly exponentially with radial distance from the rotating inner cylinder, in agreement with the findings of previous studies (Howell, Behringer & Veje 1999; Losert et al. 2000; Mueth et al. 2000). Shear causes pronounced dilation, or reduction in $\phi$, in the shear layer, and mass conservation results in compaction at larger radial positions. The experimental data are in good agreement with the predictions of the non-local model of Dsouza & Nott (2020). To our knowledge, this is the first systematic experimental measurement of dilatancy at constant volume and the coupling between packing fraction and velocity fields for granular materials in the slow-flow regime.

As elaborated in § 1, classical plasticity models, which have been traditionally used for slow granular flow, suffer from two major inadequacies: kinematic indeterminacy and inability to capture shear dilatancy. The former is associated with the mathematical ill-posedness of the models (Barker et al. 2015). While a few models have been proposed to repair kinematic indeterminacy, they treat the medium as incompressible. The non-local model of Dsouza & Nott (2020) treats the medium as compressible and builds in dilatancy. The experimental validation of the model is therefore a useful advance, and highlights the importance of incorporating dilatancy in models for slow granular flow. It is pertinent to note that Andrade et al. (2012) proposed a model for the rate-dependent regime, where they introduced a phenomenological dilatancy function $\beta \equiv \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}/\dot {\gamma }$ which decays with increasing shear strain; to capture the effect of dilatancy, they solved the time-dependent problem. In the model of Dsouza & Nott (2020), dilatancy arises naturally from the momentum balances; no additional function needs to be introduced and the steady-state solution can be obtained directly.

Though the model of Dsouza & Nott (2020) is formulated for the (non-inertial) slow-flow regime, it can be extend to include the effects of particle inertia by using the ansatz of the $\mu (I)$ model (GDR MiDi 2004; Jop et al. 2006), i.e. by considering $\mu$ to be a function of the inertia number $I$. We hope that our study will spur other investigations to probe the effect of inertia on dilatancy and its coupling to kinematics. A shortcoming of the model prediction is that we have used the experimentally measured values of $\phi$ at the two cylinders as the boundary conditions for that variable. Ideally, one would like a boundary condition that specifies $\phi$ (or its gradient) as a function of the normal stress and the kinematic state of the material. Deriving such a boundary condition from a first principles approach is an important exercise and requires further investigation.

We were only able to study the effect of small variations in the $\phi$ profile (keeping its average $\bar {\phi }$ constant) on the velocity profile. While our data do provide evidence of the inter-dependence of the two fields, it would nevertheless be useful to achieve larger variations in $\phi (r)$, such as by altering the roughness of the inner cylinder, in a future investigation. Finally, our study points to the need for precise non-invasive measurement of the volume fraction field in granular flow. Despite X-ray CT and MRI being expensive and not widely available, using them to study a few complex three-dimensional flows would be worthwhile to validate and refine models.