2.1. Flow in the standpipe connecting two hoppers

As shown in figure 2, the DEM study begins with simulations of a hopper-fed vertical pipe with a converging bottom outlet, with each component being constructed from pairs of planar walls. The upper feed hopper has two flat walls of length $L_{in}=121d$ opening with angles $\theta _{in}$ relative to the horizontal axis and joining the top of the pipe at $y=0$. The lower outlet is also composed of two flat walls converging with angles $\theta _{out}$ and has a fixed drop-height $L_h=50d$. This outlet attaches to the bottom of the parallel chute walls which are located at $x=0$ and $x=W$ for $y\in [0,L_p]$, where length $L_p$ is the pipe length. Unlike the hopper walls, the walls of the chute are geometrically roughened by affixing particles which remain static throughout the simulation. These frozen wall particles, which are otherwise identical to the flowing particles, do not overlap with one another and have centres approximately uniformly randomly distributed in $x\in [-d/2,d/2]$ and $x\in [W-d/2,W+d/2]$. Periodic walls are placed in the third direction $z\in [0,L_z]$ where $L_z=10d$ is chosen to minimise computational expense whilst allowing the flow to be fully three-dimensional. Provided that sufficient inflow is available at the top, the primary parameters controlling the mass flux out of the system are $\theta _{out}$ and $W$, which in turn set the outlet opening width $\varDelta$.

Figure 2. DEM particles coloured by vertical velocity in the hopper-controlled set-up after a steady circulating flow has been established. Here, the frozen particles which make up the wall have been rendered translucent black, as have the hopper walls. This example is for a chute width $W=30d$, depth $L_z=10d$ and with a short pipe length $L_p=200d$, chosen for illustrative purposes. The upper inflow hopper opens with angle $\theta _{in}=65^\circ$ whereas the lower outflow hopper has $\theta _{out}=80^\circ$ giving an outlet opening $\varDelta =12.5d$. It should be noted that the velocity scale is limited in panel (b) to better contrast the spread of velocities within the chute. (a) Full simulation, (b) close-up view of the steady flow region and (c) close-up view of the outflow region.

To initialise the flow in the hopper system, the full geometry, including the pipe walls, is first constructed from impenetrable planar panels. The lower hopper is also blocked with an additional wall spanning the outlet. Particles are inserted randomly into this enclosed domain and allowed to settle, ensuring that a sufficient portion is filled. The pipe walls are then formed by freezing in place the particles within the thin wall regions, as illustrated by the black grains in figure 2. Subsequently removing the wall blocking the outlet then starts the gravity-driven flow. This set-up would naturally drain the initial finite mass of the system, leading to related transient effects. Instead, two methods of replenishing material in the upper feed hopper were explored. Firstly, fresh material was added by randomly placing new grains above the top free surface whilst material leaving the outlet was deleted. This worked well, but required tuning the rate of input to match the outlet flux. To overcome this, the domain was instead made periodic in the vertical direction such that material leaving the bottom outlet re-entered the domain just above the surface of the upper hopper material. As shown in figure 2, this set-up reaches a steady recirculating state that is approximately symmetric about $x=W/2$.

Despite this symmetry in the cross-pipe direction, there is potential for much non-uniformity and unsteadiness in the downstream direction, as detailed in figure 3. For very narrow openings there is the possibility of jamming the outlet, either irreversibly or quasi-periodically, as has been widely studied in many granular systems (To, Lai & Pak Reference To, Lai and Pak2001; Zuriguel et al. Reference Zuriguel, Garcimartín, Maza, Pugnaloni and Pastor2005). Here, a snapshot is plotted in panel $(a)$ of figure 3 of an outlet which is narrow but which does not jam irreversibly. In this case the velocity field and flow rate is clearly unsteady and non-uniform throughout the system. Panel $(c)$ in the same figure demonstrates that, in the opposite extreme, the fast flow rates caused by wide outlets may be steady in the laboratory frame but non-uniform in the flow direction. In this case the flow accelerates and detaches from the pipe walls generating a dense fast core surrounded by a very dilute gas. In the present work, the focus will be on the fully steady regions of intermediate flows, for example the one observed in panel $(b)$ of figure 3, which is equivalent to the set-up detailed in figure 2. Experimentation with the opening angle reveals that such flows exist in the approximate range $77^\circ \lesssim \theta _{out}\lesssim 83^\circ$, hence the remainder of this paper will focus on these regimes.

Figure 3. Snapshots of the fully developed flow in hopper-controlled systems with substantially different outlet opening angles. Note that the plotting axes have been rotated by $90^\circ$ to aid illustration. A complementary movie animation is provided to illustrate the unsteady pulsing observed, in particular for the slow and fast flows shown in panels (a) and (c), respectively.

For set-ups tuned to generate an approximately steady uniform region, the noteworthy downstream trends which are summarised in figure 4 for a range of examples. As shown in panel $(c)$ of figure 4, the max flux

(2.1)\begin{equation} Q_M = \rho_* L_z \int_0^W v(x)\phi(x) \,\textrm{d}\kern0.7pt x , \end{equation}

is strongly dependent on the outlet opening angle $\theta _{out}$ and is effectively uniform throughout the system. However, as shown in panels $(a)$, $(b)$ and $(d)$, this constant flux of grains is achieved by a subtly varying flow in which the velocity, volume fraction and pressure slowly vary down the chute. These variations appear to reach a limiting asymptotic state in the range of $y$ that is bounded by the vertical dashed lines. The lack of significant pressure variation in the steady region of the pipe makes this flow distinct from the linear lithostatic variation observed in inclined plane flow. Known as the ‘Janssen effect’ (Janssen Reference Janssen1895), this signifies that the wall stress accommodates the weight of flowing material. Here, the maximal values of $\phi$ and $v$, which are located close to the chute centre $x=W/2$, are also approximately constant with $\mathrm {max}(\phi )$ close to random close packing $\phi _{rcp}$. It is also these regions for which the temporal fluctuations close to the inflow, indicated by the dashed range in panel $(a)$, are almost negligible. A closer study of these regions will be the focus of the next section.

Figure 4. Variation of the flow along the vertical pipe for three different outlet openings $\theta _{out}$. Here, the remaining parameters are the same as in figure 2 except that a longer pipe $L_p=450d$ is used. The maximum vertical velocity $(a)$, the maximum solid volume fraction $(b)$, the mass flux $(c)$ and the mean value of the pressure $(d)$ are plotted as functions of the vertical coordinate. Dotted curves in panel $(a)$ indicate the range whereas the solid curves are the time-averaged values. For clarity, only the averaged values are given in the other plots. Vertical black dashed lines indicate the proposed steady region which is transferred to the rigid periodic cells. In panel $(b)$ the random close packing $\phi _{rcp}$ fraction is indicated by a red horizontal dashed line.

2.2. Approximation by vertical flow in a periodic cell

Here the full hopper-fed flux-controlled set-up is reduced to one-dimensional flow by matching with smaller periodic cells. To achieve this conversion the section of pipe in which the flow is approximately uniform downstream is copied and placed into a domain in which the vertical $y$ direction is periodic. In essence this allows for the continuation of the flow development that would take place in a very long pipe. Indeed, the trends already observed in the full system, that the velocity increases and that particles migrate to the centre, are smoothly continued in these simulations. This allows for the long-time asymptotic steady state to be approximately realised, one which can be compared with time-independent solutions of the proposed equations. However, as will be seen during the model development, such steady solutions are sensitive to both the mean solid volume fraction $\phi _0$ and the width of the chute $W$. Because the walls of the pipe are constructed from frozen particles, which are irregularly spread over a range of $x$ at each wall, both the width and total volume that the flowing particles occupy are subject to uncertainties.

An alternative idealisation of vertical chute flow is made here which allows for precise control of the parameters $\phi _0$ and $W$ and also guarantees that steady flow is accompanied by no slip at the sidewalls. This is achieved through a doubly wide domain containing two equally partitioned counter-flowing regions. This geometry is perfectly cubic with dimensions $x\in [0,2W]$, $y\in [0,L_y]$ and $z\in [0,L_z]$, and each boundary face is periodic with its opposite. A counter-flow, which represents two simultaneous vertical chute flows, is then driven by an asymmetric gravitational body force such that

(2.2)\begin{equation} \boldsymbol{g} = \begin{cases} g \boldsymbol{\hat{y}} & \text{for } 0\leq x < W , \\ -g \boldsymbol{\hat{y}} & \text{for } W\leq x < 2W , \end{cases} \end{equation}

where $\boldsymbol {\hat {y}}$ is the unit vector in the vertical direction and $g$ is the gravitational acceleration magnitude, as shown in figure 5. Because of the exact asymmetry in $x$ of this arrangement, equilibrium states have equal mass in each half of the domain and no net flux in $y$. This therefore forces the average vertical velocity to be zero at the outer edges $x=0,2W$ and at the centre $x=W$, which demarcates the two flows. This counter-flow arrangement was recently employed by Kim & Kamrin (Reference Kim and Kamrin2020) for hard granular material and previously by Chaudhuri et al. (Reference Chaudhuri, Mansard, Colin and Bocquet2012) for soft jammed particles.

Figure 5. A schematic diagram of the doubly wide fully periodic cell with counter-flow. Here, the vertical periodic boundary is denoted as dot-dashed lines and the horizontal as dashed lines. The regions of downward- and upward-pointing gravity vectors are shaded to distinguish them.

After a sufficiently long time, the periodic flows reach time-independent steady states. The averaged flow fields, as they vary across the chute, are plotted in figure 6 alongside the starting flows from the full hopper-fed set-up. The symmetry about $x=W/2$ is immediately clear as is the similarity in the spatial variation of the fields for each case. The principal evolution which happens between the hopper-fed initial condition and the long-time periodic states is the slow migration of particles towards the centre. This redistribution of mass is accompanied with a slow overall acceleration, leading to faster velocities for all $x$. Whilst this is a subtle effect, the precise value of the volume fraction, particularly close to the walls, clearly plays a key role in setting the velocity magnitude. For the fully periodic cell, perfect matching with the bumpy-walled cases is not straightforward. In figure 6$(b)$ three realisations are presented, one with a mean volume fraction $\phi _0$ equal to that estimated from the coarse-grained hopper fields and another two differing above and below by $1\,\%$ increments. The deviations observed given this small relative change is another clear indication of the importance of $\phi _0$ to vertical chute flow.

Figure 6. Progressive idealisation of the geometry. A comparison between flows within the pipe for the full set-up and with equivalent long-time flows in periodic cells. Diamonds are the hopper-fed data, circles the rigid-walled periodic cells and the solid lines are the doubly wide fully periodic data. Panel $(a)$ is the solid volume fraction, varying across the chute, whereas panel $(b)$ is the velocity. All cases are identical apart from the mass flux, which varies due to different outlet opening angles $\theta _{out}$. Error bars in panel $(b)$ indicate the change in maximal velocity when the mean solid volume fraction is changed by $\pm 1\,\%$ in the fully periodic cells. Horizontal dashed lines in $(a)$ are the critical packing $\phi _c$, in black, and random close packing $\phi _{rcp}$ in red.

3.1. Vertical flow between rough parallel walls

Given the steady one-dimensional approximation, momentum conservation in $x$ reduces to

(3.1)\begin{equation} \frac{\partial p}{\partial x}=0 , \end{equation}

and in $y$

(3.2)\begin{equation} \frac{\partial\tau}{\partial x}={-}\rho_* \phi g , \end{equation}

where $p$ is the constant pressure and $\tau =\tau _{xy}$ is the only non-zero component of shear stress. These equations are closed here using the $\mu (I),\varPhi (I)$-rheology (Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006) which states that the bulk friction

(3.3)\begin{equation} \frac{|\tau|}{p} = \mu(I) , \end{equation}

and solid volume fraction

(3.4)\begin{equation} \phi=\varPhi(I) , \end{equation}

are functions solely of the inertial number

(3.5)\begin{equation} I = \frac{\textrm{d}\dot{\gamma}}{\sqrt{p/\rho_*}} , \end{equation}

where $\dot {\gamma }=|d_x v|$ is the strain rate in this geometry.

To complete the problem statement, boundary conditions and functional forms are needed. Firstly, the pipe walls are assumed to be sufficiently frictional to ensure no slip and are impenetrable such that mass is conserved. This implies

(3.6)\begin{equation} v = 0 \quad \text{at }x=0,W, \end{equation}

and that the mean solid volume fraction $\phi _0$ is a constant such that

(3.7)\begin{equation} \frac{1}{W}\int_0^W \phi(x) \,\textrm{d}\kern0.7pt x = \phi_0 . \end{equation}

Symmetry about $x=W/2$ suggests that the weight of material will be equally supported by the stresses on the walls so

(3.8)\begin{equation} \tau ={\pm}\rho_*\phi_0g \frac{W}{2} \quad \text{at }x=0,W, \end{equation}

and thus $\tau =0$ at $x=W/2$. Finally, the analysis will be restricted to the linear functions

(3.9)$$\begin{gather} \varPhi(I) = \phi_c - a I , \end{gather}$$

(3.10)$$\begin{gather}\mu(I) = \mu_c + b I , \end{gather}$$

with parameters given in table 2, which are good approximations when $I$ is not very large (see da Cruz et al. (Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005) and Hung, Stark & Capart (Reference Hung, Stark and Capart2016) for related works using linear forms).

Table 2. Parameters for the linear $\mu (I)$ and $\varPhi$ relations to match with the DEM simulations. See Appendix B for details of the rheological data and fitting.

Given that the inertial number (3.5) is strictly positive, the $\mu (I)$ function (3.10) has a minimum value $\mu _c$ when $I=0$. In turn this means that a static yield stress exists and so this model cannot accommodate the vanishing shear stress at $x=W/2$. Instead this relation implies that there are two regions of yielding material close to the walls and a region of solid-like material in the centre of the pipe. Denoting the width of the yielding regions as $\delta$, the remaining analysis will only consider solutions within the range $0\leq x\leq \delta$, which will be termed the shear zone. The central core will be assumed to behave as a plug flow with $d_x v=0$ and $\phi =\phi _P$, where the mean plug density $\phi _P$ is a free parameter. Finally, the flow in the other yielding region $(W-\delta )\leq x\leq W$ can be readily reconstructed by symmetry.

An important feature of the $\varPhi (I)$ compressible rheology is that at steady state $I$ and $\phi$ are interchangeable variables. Specifically, it is useful in the following to define the inverse function

(3.11)\begin{equation} \varPsi(\phi) \equiv \varPhi^{{-}1}(\phi) = \frac{\phi_c-\phi}{a} , \end{equation}

which is equivalent to the inertial number $I$ at steady state given the functional form (3.9). Substituting this into the $\mu (I)$ relation gives the shear stress

(3.12)\begin{equation} |\tau(x)| = \mu(\varPsi(\phi(x)))p , \end{equation}

which converts the $y$-component of momentum balance (3.2) into an autonomous ordinary differential equation (ODE) for $\phi$ as

(3.13)\begin{equation} \frac{\partial\phi}{\partial x} ={-}\frac{\rho_*g}{p}\frac{\phi}{\partial_\varPsi\mu(\varPsi(\phi))\partial_\phi\varPsi(\phi)} = \frac{a \rho_* g}{b p}\phi . \end{equation}

General solutions for $\phi$ in the shear zone are then simply exponentials, which match the edge of the yielding region ($\phi =\phi _c$ at $x=\delta$) with the form

(3.14)\begin{equation} \phi(x) = \phi_c \exp\left(\frac{x-\delta}{l_p} \right) , \end{equation}

where some of the constants have been grouped into a useful length scale

(3.15)\begin{equation} l_p = \frac{b p}{a \rho_* g}. \end{equation}

Given this solution, it is straightforward to derive an ODE for the velocity by equating the two definitions of the steady inertial number (3.5) and (3.11) which gives

(3.16)\begin{equation} \frac{\partial v}{\partial x} = \sqrt{\frac{p}{\rho_* d^2}}\varPsi(\phi(x)) = \sqrt{\frac{p}{\rho_* d^2}} \frac{\phi_c}{a} \left[1- \exp\left(\frac{x-\delta}{l_p} \right)\right] , \end{equation}

where the assumption that $d_x v\geq 0$ has been made. As with $\phi (x)$ this equation can be solved exactly, with the boundary condition coming from the no-slip walls (3.6), to give

(3.17)\begin{equation} v(x) = \sqrt{\frac{p}{\rho_* d^2}} \frac{\phi_c}{a} t(x-l_p\exp(-\delta/l_p)[\exp(x/l_p)-1]). \end{equation}

The final task is to ascertain the pressure $p$ and shear zone width $\delta$ from the remaining constraints. Given the assumption that $\phi =\phi _P$ inside the central plug, conservation (3.7) implies

(3.18)\begin{equation} \int_0^\delta \phi(x)\,\textrm{d}\kern0.7pt x+\phi_P\left(\frac{W}{2}-\delta\right) = \phi_0 \frac{W}{2} , \end{equation}

which leads to an equation relating $p$ (through $l_p$) and $\delta$ as

(3.19)\begin{equation} l_p\phi_c [1-\exp(-\delta/l_p)] = \phi_0 \frac{W}{2} - \phi_P\left(\frac{W}{2}-\delta\right) , \end{equation}

making use of the exact solution (3.14) for $\phi (x)$ in the shear zone. Similarly, the shear stress at the wall (3.8) can be equated to (3.12) evaluated at $\phi (x=0)\equiv \phi ^W$ to give

(3.20)\begin{equation} \mu(\varPsi(\phi^W))p=\rho_*\phi_0g \frac{W}{2} , \end{equation}

which in turn gives another equation

(3.21)\begin{equation} \frac{a}{b}\mu_c + \phi_c [1-\exp(-\delta/l_p)]=\phi_0 \frac{W}{2l_p} , \end{equation}

linking $l_p$ and $\delta$. Rearranging (3.21) gives an expression for the shear zone width

(3.22)\begin{equation} \delta ={-}l_p\log\left[1+\frac{a\mu_c}{b\phi_c}-\frac{\phi_0}{\phi_c}\frac{W}{2l_p}\right] . \end{equation}

Substituting this expression into (3.19) gives an equation of the form

(3.23)\begin{equation} Al_p-l_p\log\left(B-\frac{C}{l_p}\right)+D =0 , \end{equation}

where the parameter groupings are

(3.24a–d)\begin{equation} A= \frac{a\mu_c}{b\phi_P} ,\quad B=1+\frac{a\mu_c}{b\phi_c} ,\quad C= \frac{\phi_0}{\phi_c}\frac{W}{2} ,\quad \text{and} \quad D={-}\frac{W}{2}. \end{equation}

Finally, this equation has solution

(3.25)\begin{align} l_p = \frac{CD}{BD-C\mathcal{W}\left(\frac{D}{C}\exp[A+BD/C]\right)} = \frac{W}{2} \frac{\phi_0}{\phi_c+a\mu_c/b+\phi_0\mathcal{W}(-\phi_c\exp[X]/\phi_0)} , \end{align}

where $\mathcal {W}$ is the Lambert-W function and $X=A+BD/C$ is a dimensionless variable depending on $\phi _0$, $\phi _P$ and the rheological parameters.

When combined with the exact solutions for the volume fraction (3.14) and the velocity (3.17), the scaled pressure (3.25) completes the solutions for flow between parallel vertical walls. In order to remove the remaining free parameter, the plug density $\phi _P$, the approximation $\phi _P=(\phi _c+\phi _{RCP})/2$ is made throughout this paper based on the data in figure 6 showing a parabolic $\phi (x)$ bounded by the critical and random close packed values. Because this range is quite narrow, the solutions are not very sensitive to this choice.

Figure 7 shows that these solutions provide a good approximation to many aspects of the DEM simulation results in the fully periodic cells. Because the exact solutions are based on the linearised $\mu (I)$ and $\varPhi (I)$ relations, there is clearly some nonlinear spatial variation which is missing. However, as highlighted in panels $(a,c)$, the DEM values of the volume fraction closely straddle the exact solutions provided that the mean packing fraction is not too low. Due to the $\varPhi (I)$ relation, linking the volume fraction to the strain rate, this straddling conspires to ensure realistic values for the total strain rate in the shear zone and hence good values for the maximum velocity, as shown in panels $(b,d)$.

Figure 7. Long-time steady solid volume fractions and vertical velocities from the fully periodic DEM simulations (open squares) and the equivalent exact solutions (solid lines). Panels $(a)$ and $(b)$ are for variation of the chute width, given $\phi _0=0.59$, whereas panels $(c)$ and $(d)$ detail differing mean solid volume fractions for a fixed chute width $W=110d$.

For small values of the mean volume fraction ($\phi _0\lesssim 0.59$), the discrepancies between the exact solutions and the DEM results become noticeable, especially close to the walls of the chute. This is to be expected as the volume fraction at the wall is low and hence the inertial number is high. In this limit the linear approximations do not provide a good fit to the rheological data, as detailed in Appendix B. In spite of these differences, the simplified linear theory is shown here to provide very good predictions for the magnitude of the flow velocity. In every case in figure 7 the plug velocity is within $5\,\%$ of the prediction, except for the extremely dense case $\phi _0=0.61$ which has a relative error of $21\,\%$. As this outlier value was ruled out by the hopper-fed simulations, the conclusion is that the present model may be employed directly as a practical tool for estimating flow rates.

3.2. Scaling laws

One of the most useful outcomes of the simplified analysis presented above is the establishment of clear relationships between the input parameters describing the material and geometry and the resultant flow fields. Unlike non-local theories, which include other length scales, this local theory is formulated solely in terms of the pipe width $W$ and grain diameter $d$. One surprising outcome of the analysis is that the spatial structure of the solutions scales with $W$ only and that the values of $W$ and $d$ control the velocity magnitude but not the values of $\phi$.

These features can be best seen by re-writing $l_p$ from (3.25) as

(3.26)\begin{equation} l_ p = W F(\phi_0;\phi_c,\mu_c,a,b) , \end{equation}

i.e. the width $W$ multiplied by a non-dimensional function $F$ of the mean concentration $\phi _0$ and the rheological constants. Given this relation, it is also clear from (3.22) that the shear zone width $\delta$ can be written in an equivalent form

(3.27)\begin{equation} \delta = W G(\phi_0;\phi_c,\mu_c,a,b). \end{equation}

Defining $\hat {x}=x/W$ gives the solid volume fraction

(3.28)\begin{equation} \phi(\hat{x}) = \phi_c \exp\left(\frac{\hat{x}-G}{F}\right) , \end{equation}

from (3.14) and the velocity as

(3.29)\begin{equation} v(\hat{x}) =\phi_c \frac{W^{3/2}\sqrt{g}}{d}\sqrt{\frac{F}{ab}} (\hat{x}-F\exp({-}G/F)[\exp(\hat{x}/F)-1]) , \end{equation}

from its exact solution (3.17).

Another important quantity to consider is the mass flux of material

(3.30)\begin{equation} Q_M = 2\rho_*L_z \int_0^\delta v(x)\phi(x)\, \textrm{d}\kern0.7pt x + \rho_*L_z\phi_P v_P(W-2\delta) , \end{equation}

where the velocity in the plug

(3.31)\begin{equation} v_P = \phi_c \frac{W^{3/2}\sqrt{g}}{d}\sqrt{\frac{F}{ab}}\left(G-F + F \exp\left(-\frac{G}{F}\right)\right) , \end{equation}

comes from evaluating (3.29) at the inner edge of the shear zone $\hat {x}=G$. As with the other derivations in the parallel plate geometry, the integral in (3.30) can be found exactly. The key features of the resultant solution can be most conveniently expressed in the compact form

(3.32)\begin{equation} Q_M = \rho_* W^{5/2} L_z\frac{\phi_c}{d}\sqrt{\frac{g}{ab}} H(\phi_0,\phi_P;\phi_c,\mu_c,a,b) , \end{equation}

where $H$ is a grouping of the terms involving $F$ and $G$. Incidentally, the dependence on $W$ and $d$ in these scaling laws is mirrored by the incompressible $\mu (I)$-rheology (see § 4.1 and Cawthorn Reference Cawthorn2010) and ultimately comes from the dimensional arguments inherent in the formulation of the inertial number. However, the key roles played by $\phi _0$ and $\phi _p$ are naturally absent from incompressible theories.

The power and veracity of these scaling relations is demonstrated in figure 8. Despite the simplified nature of the exact solutions, these plots clearly show that many important global features of the flow are accurately reproduced. Panels (a–c), for which the chute width is varied at fixed $\phi _0$, show a particularly good match between the DEM simulations and the linear relations for pressure (3.26) and shear band with (3.27) as well as the $5/2$ power law for mass flux (3.32). Varying the mean solid fraction, as is done in panels (d–e) of figure 8, reinforces the viewpoint that flow in the vertical chute is very sensitive to the mass flux. Both the DEM results and the theory predict rising pressures and narrower shear zones as the flux diminishes due at higher packings.

Figure 8. Scalings and trends for the fully periodic vertical flow. The pressure, shear zone width and mass flux are plotted in (a–c) with varying chute widths and in (d–f) for different mean solid volume fractions. Circles are the averaged DEM data whereas solid lines are the predictions of the exact solution.

It is also worth noting that, due to the definition (3.15) of $l_p$, the pressure

(3.33)\begin{equation} p = \frac{ag\rho_*}{b} W F , \end{equation}

scales with $W$, $g$ and $\rho _*$ identically to the shear stress at the walls (3.8). However, unlike the wall traction, the dependence of $p$ on $\phi _0$ is not linear due to the product-log term in the denominator of (3.25). As such, this finding is contradictory to the assumption of a constant wall friction coefficient $\mu _W$ as the model presented here instead has

(3.34)\begin{equation} \mu_W = \frac{\tau_W}{p} = \frac{\phi_0 b}{2aF(\phi_0;\phi_c,\mu_c,a,b)} , \end{equation}

where $\tau _W$ is the shear stress at $x=0$.

3.3. Vertical flow in rough cylindrical pipes

Many practical apparatus involve grains flowing in pipes with a circular cross-section, rather than the elongated rectangular channels detailed in the previous sections. Here, a cylindrical vertical pipe with gravity-aligned axis $z$ and diameter $W$ is considered, as shown in figure 9. Like the rigid parallel-plate geometry, the pipe walls in the DEM simulations are constructed from fixed particles because this geometry does not allow for a doubly periodic arrangement. As will be seen, these walls are sufficiently frictional to ensure no slip.

Figure 9. DEM particles coloured by vertical velocity at steady state in the periodic cylindrical flow geometry. Particles in the slice $-{\rm \pi} /3<\theta <{\rm \pi} /3$ have been removed for illustration purposes.

As in the parallel plate case, steady-state flow is described by one non-zero velocity component $u_z$ which, along with the volume fraction $\phi$, is a function of the radial coordinate $r$ only. Due to the natural azimuthal symmetry, the non-trivial momentum balances are in $r$

(3.35)\begin{equation} \frac{\partial p}{\partial r} = 0 , \end{equation}

so the pressure is again a constant, and the $z$-component gives

(3.36)\begin{equation} \frac{1}{r}\frac{\partial(r\tau)}{\partial r} = \rho_*\phi g , \end{equation}

where $\tau =\tau _{rz}$. Conservation of the initial mass of material, with a presumed uniform density $\phi _0$, implies

(3.37)\begin{equation} \int_0^{2{\rm \pi}}\int_0^{W/2} \phi(r) r \,\textrm{d}r \,\textrm{d}\theta = 2{\rm \pi}\int_{W/2-\delta}^{W/2} \phi(r) r \,\textrm{d}r + {\rm \pi}\phi_P\left(\frac{W}{2}-\delta\right)^2 = {\rm \pi}\phi_0\frac{W^2}{4} , \end{equation}

and that at steady state the traction on the walls of pipe balances the weight of this material such that integrating along the wall

(3.38)\begin{equation} \int_0^{2{\rm \pi}} \tau_W \frac{W}{2} \,\textrm{d}\theta = {\rm \pi}\frac{W^2}{4} \rho_* \phi_0 g , \end{equation}

gives the wall stress boundary condition

(3.39)\begin{equation} \tau_W = \frac{W}{4}\rho_*\phi_0 g . \end{equation}

The equations can be processed, as in § 3.1, by substituting $\tau =\mu (\varPsi (\phi ))p$ into the vertical momentum balance (3.36) to recover an ODE for $\phi (r)$. For cylindrical pipes this takes the form

(3.40)\begin{equation} \frac{\partial\phi}{\partial r} ={-}\frac{\phi}{l_p}+\frac{a\mu_c/b+\phi_c-\phi}{r} , \end{equation}

where the pressure-dependent length scale $l_p=pb/(a\rho _*g)$ is the same as the parallel plate case (3.15). The condition of no yielding on the edge of the shear zone is given by $\phi (W/2-\delta )=\phi _c$ in this geometry, so that the volume fraction has the exact solution

(3.41)\begin{align} \phi(r) &= \frac{l_p}{r}(a\mu_c/b+\phi_c)+\frac{\phi_c(W/2-\delta)-l_p(a\mu_c/b+\phi_c)}{r}\nonumber\\ &\quad \times \exp\left(-\frac{r-W/2+\delta}{l_p}\right). \end{align}

As the assumptions of the $\mu (I)$,$\varPhi (I)$-rheology limit the deformation to be within the annulus $(W/2-\delta )< r\leq W/2$, the singular behaviour of $\phi$ at $r=0$ is not of concern unless $\delta =W/2$. However, this artefact of the chosen coordinate system does lead to greater complexity in the subsequent derivation of the full solution. For example, analogously to (3.16), the velocity may be recovered from

(3.42)\begin{equation} \frac{\partial v}{\partial r} = \mathrm{sign}(d_r v) \sqrt{\frac{p}{\rho_*d^2a^2}}(\phi_c-\phi(r)) , \end{equation}

which, given $\mathrm {sign}(d_r v)=-1$, results in solutions of the form

(3.43)\begin{align} v(r)& = a_1\sqrt{l_p}\left(r+(\delta-W/2+a_2l_p) \exp\left[\frac{W/2-\delta}{l_p}\right]\mathrm{Ei}({-}r/l_p)-a_2l_p\log(r/l_p) \right)\nonumber\\ &\quad +c_1 , \end{align}

where

(3.44)\begin{equation} \mathrm{Ei}(X)={-}\int_{{-}X}^\infty t^{{-}1}\,\textrm{e}^{{-}t} \,\textrm{d}t , \end{equation}

is the exponential integral, $a_{1,2}$ are non-dimensional parameter groupings and $c_1$ is a constant chosen to satisfy no slip $u(r=W/2)=0$ at the walls.

Given such divergent general solutions for $\phi$ and $v$, it is perhaps not surprising that exact expressions for the pressure and shear-band width are in general not forthcoming. Deriving simultaneous equations, as in the parallel plate case, based on conservation (3.37) and the traction condition (3.39) are instead most conveniently solved numerically, here using MATLAB's ‘vpasolve’ function. The resultant flow field solutions (3.41) and (3.43), given numerically estimated $p$ and $\delta$, are plotted for a range of pipe widths in panel $(a)$ of figure 10 and the dependence on the mean solid volume fraction $\phi _0$ is explored in panel $(b)$. Interestingly, these solutions provide a better fit to the DEM data than in the parallel walls case at small values of the mean packing fraction $\phi _0$. Part of this success may be due to the circular geometry weighting more highly the packings at larger $r$, due to the $r\,dr$ component of the elemental volume in (3.37). This means that lower overall packings come with relatively larger fractions close to the pipe walls compared with the parallel case so stay closer to the range of the linear rheology.

Figure 10. Solid volume fraction and vertical velocity profiles for a range of flows in cylindrical pipes. DEM simulation results, averaged in $z$ and $\theta$, are plotted as open symbols whereas the predictions of the exact solutions are plotted as solid curves. Panels $(a)$ and $(b)$ show the results for $\phi _0=0.565$, with different pipe widths, and panels $(c)$ and $(d)$ have results with fixed $W=50d$ and a variety of mean solid volume fractions.

Despite the additional complexities arising in the cylindrical geometry, the effect of changing the wall separation distance is found to be almost equivalent to the parallel wall case. As the shear stress at the wall (3.39) scales with $W$, the $\mu (I)$ relation means that the pressure (and hence $l_p$) also scales linearly with pipe diameter. As such, introducing the scaled radial position $\hat {r}=r/W$ converts (3.41) into $\phi (\hat {r})$ which is invariant of the value of $W$. Similarly to (3.29) the same scalings reveal that the velocity (3.43) in the cylindrical pipe also scales with $W^{3/2}$. The only significant difference is that the mass flux $Q_M\propto W^{7/2}$ for the cylindrical case, compared with the $5/2$ power law for parallel walls, because the cross-sectional area scales with $W^2$ rather than with $W$. These relations are plotted alongside equivalent DEM simulations in figure 11. The close fit of these results demonstrates the universality of the scaling laws, across the geometries considered here.

Figure 11. Scalings and trends for cylindrical pipe flow. The pressure, shear zone width and mass flux are plotted in (a–c) with varying chute widths and in (d–f) for different mean solid volume fractions. Circles are the averaged DEM data whereas solid lines are the predictions of the exact solution.

4.1. The incompressible $\mu (I)$-rheology

Neglecting the $\varPhi (I)$ relation (3.4) and fixing the volume fraction to a constant value $\phi =\phi _*$ reduces the modelling to the incompressible $\mu (I)$-rheology, introduced by Jop et al. (Reference Jop, Forterre and Pouliquen2006). For flow in a vertical chute, the task is then to find the constant pressure $p_*$ along with the shear stress $\tau$ and velocity $v$, which are functions of $x$. Solving momentum balance in the vertical direction (3.2), given the wall stress condition (3.8), gives a precisely linear variation of shear stress

(4.1)\begin{equation} \tau(x) = \rho_*\phi_*g\left(\frac{W}{2}-x\right) . \end{equation}

The $\mu (I)$-rheology then sets the width of the shear zone $x=\delta$ because $\tau /p_*=\mu _c$ at

(4.2)\begin{equation} \delta = \frac{W}{2} - \frac{\mu_c p_*}{\rho_*\phi_*g} . \end{equation}

This expression has precisely the same linear scaling in $W$ as predicted by the compressible $\mu (I),\varPhi (I)$ theory (3.27). However, the theory is severely limited beyond this as the two remaining unknown variables $p_*$ and $v(x)$ must be found from just one equation

(4.3)\begin{equation} \frac{\partial v}{\partial x} ={\pm} \sqrt{\frac{p_*}{\rho_*d^2}}\mu^{{-}1}\left(\frac{|\tau(x)|}{p_*}\right) , \end{equation}

along with the no-slip condition (3.6) at the walls. As such, the theory would require either $p_*$ to be specified, or an additional closure to be introduced. Incidentally, taking instead a constant value for $\mu$ leads to a classical Coulomb-type theory. This theory is even more limited as spatial variation of $\tau$ is prohibited at yield so a sliding plug flow would be the only admissible solution.

4.2. Kinetic theory

Savage (Reference Savage1998) modelled vertical chute flow using a theory which combined the collisional kinetic terms of Jenkins & Savage (Reference Jenkins and Savage1983) with the critical state soil mechanics of Schofield & Wroth (Reference Schofield and Wroth1968). The resultant formulation includes equations of state for the stresses which depend on both the solid volume fraction and the granular temperature

(4.4)\begin{equation} T_g = \frac{\langle \boldsymbol{c}\boldsymbol{\cdot}\boldsymbol{c}\rangle }{N_D} , \end{equation}

where $\boldsymbol {c}=\boldsymbol {u}_i-\langle \boldsymbol {u}\rangle$ is the per-particle instantaneous velocity fluctuation vector, $N_D$ is the number of spatial dimensions and angle brackets denote ensemble averaging. In application, this extension necessitates an additional balance law, to complement conservation of mass and linear momentum, which represents the balance of work done by stresses, body forces and the fluctuating collisions. The version of this equation derived by Savage (Reference Savage1998) is analogous to diffusion of the granular temperature with a variable nonlinear source term and diffusion coefficient.

Analysis of the full system of balance equations shows that $T_g$ tends to zero smoothly as the centre of the chute is approached from the walls. In turn the constitutive equations and momentum balance imply that the strain rate tends to zero whilst $\phi$ approaches a maximal value in the same limit. These aspects are all qualitatively reflected in the DEM results of the present paper. However, Savage (Reference Savage1998) goes on to claim that the solutions are not spatially similar as the chute width is changed. Instead, the boundary layers, which are equivalent to shear zones, remain roughly constant such that $\delta \simeq 10d$ irrespective of $W$. A comparison of this relation with the DEM simulations and linear scaling proposed by the present study is given in figure 12. As with the experiments of Pouliquen & Gutfraind (Reference Pouliquen and Gutfraind1996) and DEM simulations of Gutfraind & Pouliquen (Reference Gutfraind and Pouliquen1996), this key difference is most likely due to the very large mean solid volume fractions studied. Indeed, as shown in Appendix B, the rheological behaviour is markedly different to the linear $\mu (I),\varPhi (I)$-rheology when $\phi$ is large and hence $I$ is very small.

Figure 12. A comparison of the shear zone width $\delta$ from the DEM simulations for different chute widths $W$. Here, the data for parallel walls, from figure 8$(b)$, are plotted as red squares and the cylindrical pipe data, as detailed in figure 11$(b)$, are plotted as blue circles. Both axes are scaled logarithmically and the black lines represent different proposed scaling relations from the literature and the present study.

4.3. Cosserat continuum

In addition to translation, grains may rotate relative to one another. This particle-level spin may thus contribute to a bulk angular velocity vector $\boldsymbol {\omega }$ and a couple-stress tensor $\boldsymbol{\mathsf{M}}$ at the continuum level. Cosserat theory (Mühlhaus & Vardoulakis Reference Mühlhaus and Vardoulakis1987; Harris & Grekova Reference Harris and Grekova2005) then provides a framework in which these quantities are modelled and linked to the existing conservation laws of mass and momentum. Significantly, in these models the vorticity $\boldsymbol {\nabla }\times \boldsymbol {u}$ is not necessarily equal to the bulk angular velocity, so that $\boldsymbol {\omega }$ is a new independent vector, and the rotational motion may be accompanied by asymmetry of the stress tensor.

As illustrated by Mohan, Nott & Rao (Reference Mohan, Nott and Rao1999), application of Cosserat theory to vertical chute flow generates an additional differential equation, coming from angular momentum conservation, for the single non-zero component of the couple-stress tensor $M_{xz}$. This combines with the shear-stress components to modify the Mohr–Coulomb yield condition, thus allowing spatial variation of the flow in $x$. As with the kinetic theory, these modifications imply that all material in the vertical chute is yielding and that no precisely static central plug exists. Despite this, Mohan et al. (Reference Mohan, Nott and Rao1999) performed an asymptotic analysis for wide chutes to find

(4.5)\begin{equation} \lim_{W\to\infty} \delta \to \hat{l}L^{2/3}W^{1/3} , \end{equation}

where $\hat {l}$ is a non-dimensional constant and $L$ is a proposed length scale, over which the couple-stress moments act. This relation is plotted alongside the DEM data and exact linear scaling of the present study in figure 12. Part of the mismatch may be explained by Ananda, Moka & Nott (Reference Ananda, Moka and Nott2008), who have since shown that full numerical solutions approximate this scaling for chute widths much greater than $W\gtrapprox 100d$ and that quite different behaviour is observed for narrow chutes. However, this trend is not observed here.

4.4. Finite-size effect non-local modelling

Recent models have been proposed with the aim of including the effects of both fluctuations and finite-size effects. In particular, the theories of Kamrin & Koval (Reference Kamrin and Koval2012) and Bouzid et al. (Reference Bouzid, Trulsson, Claudin, Clément and Andreotti2013) propose that grain motions may be correlated over so-called ‘cooperative length scales’ and that these connections enable deformations to occur even when local stress fields are not sufficient to overcome the static yield stress. In effect, the stress solutions are those of the ‘local’ equations, but deformations obey dynamic equations modified by the presence, or lack, of fluctuations.

Originally the non-local granular fluidity (NGF) theory of Kamrin & Koval (Reference Kamrin and Koval2012) was formulated phenomenologically in terms of a scalar $\mathcal {G}$ reflecting the ‘granular fluidity’ (labelled $g$ in the literature). The connection of this concept to flow propensity and yield is made clear through the definition

(4.6)\begin{equation} \mathcal{G} \equiv \frac{\dot{\gamma}}{\mu} , \end{equation}

and a more direct link between ‘fluidity’ and velocity fluctuations has been established by Zhang & Kamrin (Reference Zhang and Kamrin2017) who suggest $\mathcal {G}=\sqrt {\langle \boldsymbol {c}\boldsymbol {\cdot }\boldsymbol {c}\rangle }F(\phi )/d$, where $F$ is a non-dimensional function of the volume fraction. This fluidity field is governed by an additional differential equation which is loosely based on the partial-fluidisation theory of Aranson & Tsimring (Reference Aranson and Tsimring2002) and shares structural similarity with the Ginzburg–Landau theory of phase transitions. As in the kinetic theory, at steady state the NGF model predicts smooth diffusion of the fluctuations, and as in Cosserat theory, a proposed length scale mediates the spatial variations.

The NGF theory was applied to steady vertical flow by Kamrin & Koval (Reference Kamrin and Koval2012) who used a pressure-controlled cell in which the fixed vertical direction was periodic and the bounding transverse walls moved dynamically to set a prescribed pressure. Coupled to the incompressible $\mu (I)$-rheology, as detailed above in § 4.1, the non-local theory smoothed the transition between the shear regions and the plug and thus provided a good fit for DEM data. A recent extension to this is the work of Kim & Kamrin (Reference Kim and Kamrin2020) who detail collapses of both $\mu$ and $T_g$ against $\phi$ in the shear-plug transition region for vertical flow and a range of other geometries. Their proposed linear $\phi -\mu$ fit is similar to the laws here, which can be seen by combining (3.10) with (3.9) to give

(4.7)\begin{equation} \varPhi(\mu)=\phi_c-\frac{a}{b}(\mu-\mu_c) . \end{equation}

However, the parameter values are different as the fits in the present work account for the high-$I$ values close to the walls whereas Kim & Kamrin (Reference Kim and Kamrin2020) neglect these regions. This suggests that a nonlinear or piecewise curve could bridge the gap between these fits.

Liu & Henann (Reference Liu and Henann2018) also considered vertical chute flow as a test case for the NGF model. The focus of that paper being the isolation of the transition between fully static material and flowing material. This transition was found to occur when the channel is sufficiently wide, a finite-size effect that would not be present in the $\mu (I)$,$\varPhi (I)$-rheology. It remains to be seen, however, how the non-local length scale affects the fully developed flow solutions and hence $\delta$.