4.1. DEM

We first consider the trends shown by the DEM, and then compare them with the model predictions. For smooth walls, the magnitude of the scaled vertical velocity $-u_y/\sqrt {g W}$, where $g$ is the acceleration due to gravity and $W$ is the half-width of the channel in the $x$-direction, decreases as $x$ increases (see the symbols in figure 3). This is in accord with expectation and experimental observations that the walls retard the flow of the material. For a fixed value of $2\,W/d_p$, the flow is expected to become faster as the bulk solids fraction $\bar {\phi }$ decreases. Hence $-u_y$ increases. Similarly, as $2\,W/d_p$ increases, the confining effect of the walls decreases and hence the flow is faster. Rough walls retard the flow and hence the velocities are smaller in magnitude than for the smooth walls (the insets in figure 3). Even though it is often assumed that rough walls imply no slip, this is not true for the parameter values used here.

Figure 3.

Variation of the scaled vertical velocity $-u_y/\sqrt {g W}$ with $x/W$ for smooth (main figures) and rough (insets) walls, where $g$ is the acceleration due to gravity, $W$ is the half-width of the channel and $x$ is the horizontal coordinate. The first row (a–d) and second row (e–h) correspond to $2\,W/d_p = 40$ and 80, respectively, where $d_p$ is the particle diameter: $\bar {\phi } = 0.62$, (a,e); 0.61, (b, f); 0.60, (c,g); 0.59, (d,h). Here, $\bar {\phi }$ is the bulk solids fraction. The symbols correspond to the DEM results, and the curves to the predictions of the models: Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (brown solid line); Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) (red dashed line); Henann & Kamrin (Reference Henann and Kamrin2013) (green dotted line); Dsouza & Nott (Reference Dsouza and Nott2020) (magenta dash-dotted line – numerical solution, blue – approximate solution).

The thickness $\varDelta$ of the shear layer is chosen as $W - x_*$, where $u_y (x_*) - u_y(W) = 0.95 \, (u_y(0) - u_y(W))$. For smooth walls, the dimensionless thickness $\varDelta /d_p$ is in the range 4–9 (the open black symbols in figure 4a). If $\varDelta$ scales with $d_p$, $\varDelta /d_p$ should be independent of $W$, but it appears to increase with $W$. Computations should be done over a larger range of $W$ values to check whether this is a general trend. The main obstacle here is that the DEM calculations become very time consuming when the number of particles increases. As expected, rough walls lead to thicker shear layers (the filled black symbols in figure 4a). As $\bar {\phi }$ decreases, the thickness increases for both smooth and rough walls.

Figure 4.

Variation of the dimensionless thickness of the shear layer $\varDelta /d_p$ (a,b); (a) DEM results, (b) model results; mass flow rate $\dot {M}/\dot {M}_b$ (c) and slip length $l_u/d_p$ (d) with $\bar {\phi }$: $2\,W/d_p = 30$ ($\vartriangle$); 40 ($\diamond$); 50 ($\circ$); 80 ($\square$). The open and filled symbols correspond to the smooth and rough walls, respectively, and the black symbols to the DEM results. The model predictions are shown by the coloured symbols and the curves joining them: Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (brown solid line); Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) (red dashed line); Henann & Kamrin (Reference Henann and Kamrin2013) (green dotted line); Dsouza & Nott (Reference Dsouza and Nott2020) (magenta – numerical solution).

It is of interest to compare the mass flow rate $\dot {M}$ obtained from the DEM with that predicted by the correlation of Beverloo, Leniger & Van de Velde (Reference Beverloo, Leniger and Van de Velde1961) ($\dot {M}_b$). The latter has been deduced from data for flow through bins with circular exit slots, and can be modified as suggested by Nedderman et al. (Reference Nedderman, Tüzün, Savage and Houlsby1982) for rectangular exit slots. The solids fraction $\phi$ used in the correlation is taken as 0.6, which is within the range of $\bar {\phi }$ values used in the present work. The comparison of $\dot {M}$ with $\dot {M}_b$ is not strictly valid, as the actual bin has impermeable walls in the $z$-direction, whereas the present work is based on the use of periodic b.c.s in this direction. For smooth walls, $\dot {M}/\dot {M}_b$ is in the range 1–18 (figure 4c). Thus the DEM overestimates $\dot {M}$, probably because the bin width $2\,W$ and the slot width $\boldsymbol{\mathsf{D}}$ are equal in the former case (figure 1b,c) but $2\,W > \boldsymbol{\mathsf{D}}$ in the latter case (figure 1a). The discrepancy increases when $\bar {\phi }$ decreases, in keeping with the trends shown by the DEM for the scaled centreline velocity $-u_y(x = 0)/\sqrt {g W}$ (figure 3). (The mass flow rate $\dot {M}$ depends on $u_y$ and $\phi$, but the dependence on $u_y$ is stronger as $\phi$ varies only slightly except close to the wall (figure 5).) For a fixed value of $\bar {\phi }$, $\dot {M}/\dot {M}_b$ increases as $2\,W/d_p$ increases. This is because $\dot {M}_b$ is independent of $W$ but $\dot {M}$ increases as $W$ increases.

Figure 5.

Variation of the solids fraction $\phi$ with $x/W$: $\bar {\phi } = 0.62$ (a,e); 0.61, (b, f); 0.60, (c,g); 0.59, (d,h). Panels (a–d) and (e–h) correspond to $2\,W/d_p = 40$ and 80, respectively. The symbols correspond to the DEM results for smooth ($\circ$) and rough ($\square$) walls, and the curves to the model predictions of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (brown solid line); Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) (red dashed line); Henann & Kamrin (Reference Henann and Kamrin2013) (green dotted line); Dsouza & Nott (Reference Dsouza and Nott2020) (magenta dash dot line – numerical solution, blue – approximate solution). The variation of $\phi$ in a static column is shown by black filled symbols in (a,e).

For rough walls, $\dot {M}$ is lower than that for smooth walls (figure 4c). However, $\dot {M}_b$ corresponds to a situation where there is a core of rapidly flowing material near the exit slot, adjacent to static or quasistatic ‘shoulders’ of material. Thus the material flows through a channel with curved, rough walls, and there is no explicit dependence on wall roughness in the correlation of Beverloo et al. (Reference Beverloo, Leniger and Van de Velde1961). For $\bar {\phi } = 0.62$, $\dot {M}/\dot {M}_b < 1$, unlike the case for smooth walls. As mentioned earlier, there is no flow for $\bar {\phi } > 0.62$, and hence the material slows down considerably as this limit of $\bar {\phi }$ is approached. The effect may be more pronounced for rough walls, leading to $\dot {M}/\dot {M}_b < 1$.

In the actual bin, $\bar {\phi }$ is not within the control of the experimenter for a fixed slot width. As there is no exit slot in the present work, the magnitude of the centreline velocity ${-u_y (x = 0)}$ is several multiples of $\sqrt {g W}$, where $g$ is the acceleration due to gravity. Bhateja & Khakhar (Reference Bhateja and Khakhar2020) have simulated slow flow in a bin with slot widths in the range 6–8 $d_p$, and with $\boldsymbol{\mathsf{D}}/(2\,W)$ in the range 0.2–0.27 (figure 1a); they obtained smaller centreline velocities (${\approx }0.13 \sqrt {gW} \text{--} 0.25 \sqrt {gW}$). In their study, the periodic b.c. in the flow direction was relaxed, but particles leaving the bin were re-introduced at the top and allowed to fall from rest at random horizontal positions. The effect of the walls was removed by applying periodic b.c.s in the other two directions, and hence the properties varied in the flow direction only. Alternatively, the emptying of a bin can be examined. At some height between the free surface and the exit slot, a quasi-steady state may be expected to prevail for a time interval before the bin is fully empty. For example, the emptying of seeds from an axisymmetric bunker was studied by Danczyk et al. (Reference Danczyk, Meaclem, Mehdizad, Clarke, Galvosas, Fullard and Holland2020) using experiments and DEM simulations. In the simulations, the velocity profile ‘reached steady state flow after about 1 s. Steady flow continued until approximately 18 s, at which point all the seeds had drained from the hopper’. Some of these strategies will be examined in future work to permit a comparison with experiments.

The slip length $l_u$ in (3.33) is evaluated in the DEM as $l_u/d_p = (u_y(x = W - 0.75\,d_p) - u_y(x = W - 2\,d_p))/(1.25\,d_p u_y(x = W - 0.75\,d_p))$, where $(x = W - 0.75\,d_p)$ and $(x = W - 2\,d_p)$ are the centres of two bins adjacent to the wall. The distance between the centres of these bins is $1.25\, d_p$. The slip length is approximately independent of $W$ and $\bar {\phi }$ for both smooth and rough walls (figure 4d). This is an encouraging result, as it suggests that $l_u$ may indeed be a material parameter that is related to the wall roughness. Further studies are needed to confirm this conjecture. The slip length is approximately 4 times higher for smooth walls than for rough walls. It may be expected that the magnitude of the velocity decreases and the magnitude of the velocity gradient increases as the wall roughness increases, leading to a lower value of $l_u$.

For the static column, the $\phi$ profile is approximately flat except near the wall, where it decreases (the black filled symbols in figure 5a,e). Because the material is sheared near the wall but not near the centre, the solids fraction $\phi$ is expected to be lower at the wall (figure 5). For smooth walls, there is a sharp decrease near the wall, and the values are less than $\phi _{min}$. Thus frictional effects are expected to be small compared with collisional effects in this region. For a fixed value of $2\,W/d_p$, as $\bar {\phi }$ decreases, the $\phi$ profiles shift downward, except near the centre. This is expected as $\bar {\phi }$ represents the average value of $\phi$. However, the value of $\phi (x = 0)$ is largely independent of the value of $\bar {\phi }$, for reasons that are not clear. For a fixed value of $\bar {\phi }$, the $\phi$ profiles are not affected much by changes in $2\, W/d_p$, except near the wall. For rough walls, the decrease is less pronounced, and the curves for smooth and rough walls cross near the wall. We do not have a simple qualitative explanation for this behaviour, even though it would be very desirable.

Consider the variation of the inertial number $I$, defined by (3.2). In the literature, the characteristic normal stress $N$ is chosen as $N = p \equiv (\sigma _{xx} +\sigma _{yy}+\sigma _{zz})/3$, where $p$ is the pressure. As mentioned in the previous section, two of the models do not provide estimates of $\sigma _{zz}$, and other two models do not predict normal stress differences. To facilitate comparison with these models, $N$ is calculated as $N = \sigma _{xx}$. It is found that there is only a slight difference in the values of $I$ calculated using $p$ and $\sigma _{xx}$ (not shown). The inertial number varies from very small values of order $10^{-5}$ (approximately zero) at the centre to large values of order 1 at the wall (figure 6). Thus the flow is quasistatic near the centre and inertial near the wall. The situation is unlike the case of steady flow down an inclined plane, where $I$ is a constant and the velocity varies in the cross-stream direction (Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006). Here, both $I$ and $u_y$ vary with $x$. For a fixed value of $2\,W/d_p$, $I$ increases as the bulk solids fraction $\bar {\phi }$ decreases (figure 6). This is in accord with intuition, as inertial effects are expected to become more important in systems that are less dense. The value of $I$ for rough walls is greater than that for smooth walls over most of the channel, probably because the shear rate is expected to be higher in the former case. However, close to the wall, the trend is reversed. As the width $2\,W/d_p$ increases, the confining effect of the walls decreases. Hence it may be expected that the flow will be faster, and the shear rate will be higher near the wall. This leads to a larger value of $I$. This is not a rigorous argument, but merely one based on intuition.

Figure 6.

Variation of the inertial number $I$ with $x/W$. The first row (a–d) and second row (e–h) correspond to $2\,W/d_p = 40$ and 80, respectively: $\bar {\phi } = 0.62$ (a,e); 0.61 (b, f); 0.60 (c,g); 0.59 (d,h). The symbols correspond to the DEM results for smooth ($\circ$) and rough ($\square$) walls, and the curves to the model predictions of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (brown solid line); Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) (red dashed line); Henann & Kamrin (Reference Henann and Kamrin2013) (green dotted line); Dsouza & Nott (Reference Dsouza and Nott2020) (magenta dash-dotted line – numerical solution, blue – approximate solution).

The $y$-component of the momentum balance shows that the shear stress $\tau$ varies linearly with $x$ if $\phi$ is a constant. This condition holds in the plug region, but is violated in the shear layer. The DEM results show an approximately linear variation, with slight deviations apparent in the shear layer (figure 7a,b).

Figure 7.

Variation of the scaled shear stress $\tilde {\tau } = \tau /(\rho _pgW)$ with $x/W$ for smooth (a) and rough (b) walls. In (c,d), the markers above and below the horizontal dashed line correspond to the scaled first $(\sigma _{xx}-\sigma _{yy})/\sigma _{xx}$ and second $(\sigma _{yy}-\sigma _{zz})/\sigma _{xx}$ normal stress differences, respectively, for smooth (c) and rough (d) walls. The data shown in (a–d) represent the DEM results for $\bar {\phi } = 0.60$. (e) Variation of the scaled normal stress $\tilde {\sigma }_{xx} = \tilde {N} = \sigma /(\rho _pgW)$ with $\bar {\phi }$. Its inset shows the DEM results for the angle of wall friction $\delta$ (in degrees) (see (3.34)). In (a–e), $\vartriangle$, $\diamond$, $\circ$ and $\square$ are for $2\,W/d_p = 30$, 40, 50 and 80, respectively. In (e), the open and filled symbols correspond to smooth and rough walls, respectively, the black symbols to the DEM results, and the curves to the models predictions of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (brown symbols solid line); Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) (red symbols dashed line); Dsouza & Nott (Reference Dsouza and Nott2020) (magenta symbols – numerical solution, blue symbols dash-dotted line – approximate solution).

For smooth walls, the scaled first normal stress difference $(\sigma _{xx} - \sigma _{yy})/\sigma _{xx}$ is positive, and varies from approximately 10 % at the centre to approximately 30 % near the wall (figure 7c). The scaled second normal stress difference $(\sigma _{yy} - \sigma _{zz})/\sigma _{xx}$ is largely negative, and is almost zero at the centre and approximately 30 % near the wall. Such differences do not appear to have been reported earlier for fully developed flow through a channel. As discussed later, they provide a stringent test of constitutive equations. Unfortunately, we do not have a simple physical explanation for their occurrence. Silbert et al. (Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001) observed the normal stress difference of approximately 15 %–20 % in inclined chute flow, but they too did not provide a simple explanation. For rough walls, the trends are similar (figure 7d), except that the values are slightly smaller in magnitude near the wall.

Values of the scaled normal stress $\tilde {N} \equiv N/(\rho _p g, W)$ are of O(1) for the DEM results (figure 7e). For a fixed value of $2\,W/d_p$, $\tilde {N}$ decreases as $\bar {\phi }$ decreases. The stress tensor is evaluated as the sum of two contributions, one arising from contacts (including collisions), and one arising from streaming (Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001). Our calculations show that the latter is small compared with the former. Small changes in $\bar {\phi }$ cause much larger changes in $\tilde {N}$. For systems with high values of $\bar {\phi }$, this may be anticipated as the frictional or contact stresses are known to be sensitive to small changes in $\phi$. Further, the material is expected to lose sustained contact with the wall as $\bar {\phi }$ decreases. Hence $\tilde {N}$ is likely to decrease. As mentioned by one of the referees, the stress depends on both the strength of the contact forces and the number of contacts. If the latter dominates, $\tilde {N}$ will decrease as $\bar {\phi }$ decreases. However, we do not have a convincing argument to prove this assertion at present. It is hoped that a more satisfactory argument may be advanced in the future. For a fixed value of $\bar {\phi }$, $\tilde {N}$ is a very weak function of $2 W/d_p$. Thus the scaling $\tilde {N} = N/(\rho _p g W)$ is reasonable but not exact. The normal stress for the smooth wall is approximately 12 %–20 % higher than that for the rough wall (figure 7e). Intuitively, as the angle of wall friction $\delta$ is expected to be higher for the rough wall, (3.35) explains the decrease of $\tilde {N}$ with the increase in wall roughness for a fixed value of $\bar {\phi }$.

The value of the angle of wall friction $\delta$, calculated from (3.34) decreases by approximately 50 % of the larger value as $\bar {\phi }$ increases, for both smooth and rough walls (see the inset of figure 7e). For a fixed value of $\bar {\phi }$, $\delta$ is approximately 15 %–20 % higher for rough walls than for smooth walls. Compared with the slip b.c. (3.33), the wall friction b.c. (3.34) appears to be more robust. This issue merits further investigation.

As the shear rate and the shear stress increase with $x/W$ and the normal stress $N$ is a constant, both $\mu \equiv \tau /N$ and $I$ increase, leading to the behaviour shown in figure 8. For very small values of $I$, $\mu$ appears to tend to a non-zero value. For intermediate values of $I$, the curves for smooth and rough walls are almost identical except for $\bar {\phi } = 0.62$. For large values of $I$, $\mu$ is larger for rough walls than for smooth walls (figure 8). As the curves for $I$ (figure 6) intersect near the wall, the value of $I$ in the wall region is lower for the rough wall . Hence, for a fixed value of $I$, the corresponding value of $x/W$ is larger for the rough wall, leading to a larger value of $\tau$. As the normal stress $N$ is also smaller for the rough wall (figure 7e), $\mu$ is higher near the wall.

Figure 8.

Variation of the friction coefficient $\mu$ with the inertial number $I$. The first row (a–d) and second row (e–h) correspond to $2\,W/d_p = 40$ and 80, respectively: (a,e), $\bar {\phi } = 0.62$; (b, f), $\bar {\phi } = 0.61$; (c,g), $\bar {\phi } = 0.60$; (d,h), $\bar {\phi } = 0.59$. The markers ($\circ$ – smooth walls, $\square$ – rough walls) correspond to the DEM results and the curves to the model predictions of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (brown solid line); Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) (red dashed line); Henann & Kamrin (Reference Henann and Kamrin2013) (green dotted line).

As $I$ increases, the shear rate increases. Hence the material is expected to dilate more leading to a lower value of $\phi$, or a higher value of $\phi _{drp} - \phi$ (figure 9). Except for $\bar {\phi } = 0.62$, the results for smooth and rough walls almost coincide, indicating that $\phi$ vs $I$ relation is independent of the wall roughness. For $\bar {\phi } = 0.62$, $\phi$ is slightly lower for rough walls. For $2\,W/d_p = 80$, there are two distinct regimes, except for $\bar {\phi } = 0.62$, (i) where $\phi _{drp} - \phi$ increases gradually with $I$, and (ii) where $\phi _{drp} - \phi$ increases more rapidly for $I > 10^{-1}$.

Figure 9.

Variation of $\phi _{drp} - \phi$ with $I$. The first row (a–d) and second row (e–h) correspond to $2\,W/d_p =$ 40 and 80, respectively: (a,e), $\bar {\phi } = 0.62$; (b, f), $\bar {\phi } = 0.61$; (c,g), $\bar {\phi } = 0.60$; (d,h), $\bar {\phi } = 0.59$. The markers ($\circ$ – smooth walls, $\square$ – rough walls) correspond to the DEM results; Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (brown solid line); Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) (red dashed line).

4.2. Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017)

For reasons explained in § 4.2.1, this model could not be solved for the entire range of parameter values used here. Using (3.23), (3.25), (3.27) and (3.31), $\phi$, $I$ and $N$ can be obtained independently of the wall roughness. This result is a consequence of matching $\bar {\phi }$ to the DEM results, and not using the alternative b.c. (3.34). Hence the velocity gradient (3.30) is the same for both smooth and rough walls, and the velocity profile corresponds to a shift along the velocity axis when the wall roughness changes. Thus the form of the velocity profile in unaffected by the wall roughness, as shown in figure 3(c,d) (the brown curves). The same conclusion holds for the model of Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019).

For $2\,W/d_p = 40$, the magnitude of the scaled velocity is qualitatively similar to the DEM results (the brown curves in figure 3c,d). The values are underestimated near the centre and overestimated near the wall for both smooth and rough walls.

For $2\,W/d_p = 80$ and $\bar {\phi } = 0.62$, the profile agrees fairly well with the DEM results for smooth walls (the brown curve in figure 3e). However, the values are overestimated near the centre and underestimated near the wall for $\bar {\phi } = 0.61$ and smooth walls, unlike the case for $2\,W/d_p = 40$ (the brown curve in figure 3f). For rough walls, the scaled velocity changes sign near the wall because of the shift, and hence leads to unrealistic profiles (not shown in figure 3e, f).

The scaled thickness of the shear layer $\varDelta /d_p$ varies with $\bar {\phi }$ and $2\,W/d_p$ like the DEM results (the brown symbols in figure 4b). It is of the same order of magnitude, but the range is larger than the latter. For reasons explained earlier, this model predicts that $\varDelta /d_p$ is independent of the wall roughness, in contrast to the DEM results.

The values of the slip length $l_u$ calculated from (3.33) vary widely (the brown symbols in figure 4d). In some cases, they exceed the half-width of the channel. Thus this b.c. is not reasonable for this model. Conversely, if $l_u$ is treated as a material parameter, the predicted mass flow rates will differ significantly from the DEM results.

The solids fraction profiles are approximately flat, as shown by the insets of figure 5(c–f) (the brown curves). It behaves like an incompressible model, at least in the present problem, and hence does not match the DEM results well. As discussed earlier, the profiles are independent of the wall roughness for a fixed value of $\bar {\phi }$.

The vertical asymptotes in the $I$ profiles correspond to the boundary between the plug and the shear layers (the brown curves in figure 6c–f), whereas the DEM results show a more gradual variation. The profiles agree reasonably well with the latter near the wall.

The scaled shear stress $\tilde {\tau }$ almost coincides with the DEM results for all the models, and hence the profiles are not shown. The first normal stress difference $\sigma _{xx} - \sigma _{yy}$ is zero, and the second normal stress difference $\sigma _{yy} - \sigma _{zz}$ cannot be evaluated as $\sigma _{zz}$ is unknown. This is at variance with the DEM results, particularly near the wall (figure 7c,d). The scaled normal stress $\tilde {N} = N/(\rho _p g W)$ underestimates the DEM results, and is independent of wall roughness as mentioned earlier (the brown symbols in figure 7e). Unlike the latter, $\tilde {N}$ varies more strongly with $2\,W/d_p$ and $\bar {\phi }$.

Hence the angle of wall friction $\delta$ cannot be treated as a material parameter (see (3.35)). Further, as $\tilde {N}$ is independent of wall roughness, the model predicts that $\delta$ is also independent of the wall roughness, contrary to expectation. This problem arises because $\bar {\phi }$ has been matched to the DEM value. On the other hand, if $\delta$ is determined from independent experiments, then the predicted value of $\bar {\phi }$ for a fixed value of $2\,W/d_p$ may not match the value used to generate the DEM results.

From (3.26), as $I \rightarrow 0$, $\mu \rightarrow \mu _s = 0.38$, which is much larger than the DEM results (the brown curves in figure 8c–f). At large values of $I$, the trends are similar to the DEM results for rough walls, but the values of $\mu$ are overestimated. For smooth walls, the latter show a different curvature. Thus a single function $\mu (I)$ cannot fit data for both smooth and rough walls. The brown curves in figure 8(c–f) are obtained by fitting (3.26) to the DEM results of Chialvo et al. (Reference Chialvo, Sun and Sundaresan2012) for homogeneous shear between parallel plates (figure 2a). As this function does not fit our DEM results well, it follows that $\mu (I)$ is not a true constitutive equation, but one that depends on the geometry. As $\phi$ does not vary much with $x/W$ (the insets in figure 5c–f), plots of $\phi$ vs $I$ are effectively horizontal lines, and hence they do not capture the variation shown by the DEM results (the brown curves in figure 9).

4.2.1. Reasons for the lack of solutions for some parameter values

Introducing the dimensionless quantities $\tilde {x}_p \equiv x_p/W$ and $\tilde {N} \equiv N/(\rho _p g W)$, (3.29) reduces to $\tilde {x}_p = (\mu _s \tilde {N})/ \phi _p$. As $\tilde {x}_p \leqslant 1$, we have

(4.1)\begin{equation} \mu_s\tilde{N} \leqslant \phi_p . \end{equation}

As $\beta (I \rightarrow 0) = 2\,\mu _s$, (3.27) implies that

(4.2)\begin{equation} \phi_p = \frac{\mu_s \tilde{\varLambda}\phi_{min} + \tilde{N}\phi_{max}}{\mu_s \tilde{\varLambda} + \tilde{N}} , \end{equation}

where $\tilde {\varLambda } \equiv \varLambda /(\rho _p gW)$. Using (4.1) and (4.2), we obtain

(4.3)\begin{equation} \mu_s \tilde{N}^{2} + (\mu_s^{2}\tilde{\varLambda} - \phi_{max})\tilde{N} - \mu_s\tilde{\varLambda}\phi_{min} \leqslant 0. \end{equation}

For (4.3) to hold, $\tilde {N}$ must be less than a critical value $\tilde {N}_c$, which is the positive root obtained by equating the left-hand side of (4.3) to zero. For a chosen value of $\tilde {\varLambda }$, the root corresponds to the point of intersection of the curve (4.2) and the line given by equality in (4.1) in the $\phi _p$–$\tilde {N}$ plane (figure 10). For a value of $2\,W/d_p$, or $\tilde {\varLambda }$, if the curve (4.2) corresponding to that $\tilde {\varLambda }$ crosses the line (4.1) at a $\phi _p = \phi _{u}$, the model cannot be solved for $\bar {\phi } > \phi _{u}$.

Figure 10.

Variation of the dimensionless normal stress $\tilde {N}$ with the solids fraction $\phi _p$ in the plug layer for the model of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017): solid curves, plots of (4.2) for specified values of $\tilde {\varLambda } = \varLambda /(\rho _p g W)$, where $\rho _p$ is the particle density, $g$ is the acceleration due to gravity and $W$ is the half-width of the channel in the $x$-direction; upper dash-dot line, equality in (4.1); lower dashed line, equality in (4.6).

An approximate solution provides an additional insight. As seen from our DEM results, the model predicts a very small variation of $\phi$ (the insets in figure 5c–f). Treating $\phi$ as a constant (${=}\phi _p$), and using (3.25) and (3.26), the second of (3.23) can be integrated to obtain

(4.4)\begin{equation} \frac{I}{I_0+I} = \frac{\phi_p}{\tilde{N} \Delta \mu} (\tilde{x}-\tilde{x}_p) , \end{equation}

where $\tilde {x} \equiv x/W$ and $\tilde {x}_p \equiv x_p/W$. Simplifying, we obtain

(4.5a,b)\begin{equation} I = \frac{I_0 (\tilde{x} - \tilde{x}_p)}{\left(\dfrac{\tilde{N} \mu_d}{\phi_p} - \tilde{x} \right)} ; \quad \mu_d \equiv \mu_s + \Delta \mu . \end{equation}

Equation (4.5a,b) shows that $I \rightarrow \infty$ as $\tilde {x}$ increases if

(4.6)\begin{equation} \frac{\tilde{N} \mu_d}{\phi_p} \leqslant 1. \end{equation}

For a chosen value of $\tilde {\varLambda }$, the limiting condition corresponds to the intersection of the curve (4.2) with the line given by the equality in (4.6) (figure 10). This gives a limiting value of $\phi _p = \phi _l$ such that there is no solution for $\bar {\phi } < \phi _l$. Thus the model cannot be solved for both small and large values of $\bar {\phi }$, but the limits depend on the value of $\tilde {\varLambda }$. The lower limit should be used with care as it is based on an approximate solution.

For typical values of a particle density $\rho _p = 2650\ \textrm {kg}\ \textrm {m}^{-3}$, $\varLambda = 555 \ \textrm {N}\ \textrm {m}^{-2}$, a particle diameter $d_p = 1$ mm and $2\, W/d_p = 30$, 40, 50 and 80, the values of $\tilde {\varLambda } = \varLambda /(\rho _p g W)$ are 1.4, 1.1, 0.8 and 0.5, respectively. For these values of $\tilde {\varLambda }$, figure 10 shows that solutions can be obtained for $\bar {\phi } = (0.6, 0.59)$ for $2\,W/d_p = 30$ and 40, (0.61, 0.6, 0.59) for $2\,W/d_p = 50$ and (0.62, 0.61) for $2\,W/d_p = 80$.

4.3. Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019)

As $\bar {\phi }$ is matched to the DEM results and $\phi _c = 0.613$ in (3.40), a solution could not be obtained for $\bar {\phi } = 0.62$. As $\phi (x = 0) \approx 0.63$ for the DEM results and $\phi (x = 0) \leqslant \phi _c$ for the model, we expect the predictions to overestimate $\phi$ near the wall for $\bar {\phi } \leqslant 0.61$ (the red curves in figure 5). The scaled velocity profiles are qualitatively similar to the DEM results, except that the values are higher near the wall and lower near the centre (the red curves in figure 3). This holds for both smooth and rough walls. The scaled thickness of the shear layer $\varDelta /d_p$ overestimates the DEM results for most of the parameter values (the red symbols in figure 4b), but is of the same order of magnitude. It varies with $2\,W/d_p$ and $\bar {\phi }$, but does not vary with the wall roughness, similar to that of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017). The range of $l_u/d_p$ values is comparable to that of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) (the red symbols in figure 4d). As in the case of the latter, $l_u$ is not a material parameter.

The $\phi$ profiles (the red curves in figure 5) are better than the flat profiles predicted by Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017). However, there is no improvement in the profiles of $I$, which are similar to those of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017), except near the wall (the red curves in figure 6). For a fixed value of $\bar {\phi }$, the $I$ profiles are independent of $2\,W/d_p$ unlike the profiles of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017). This is because the right-hand side of (3.41) involves the ratio of $x/W$ to the scaled normal stress $\tilde {N} = N/(\rho _p g W)$, and $W$ does not occur explicitly. Hence $\phi$ is independent of $2\,W/d_p$ (see (3.40)), leading to a similar behaviour for $\tilde {N}$ (the red symbols in figure 7e). Equation (3.35) shows that $\delta$ is also independent of $2\,W/d_p$, and varies only with $\bar {\phi }$. The values $\tilde {N}$ are underestimated. The first normal stress difference is zero, and the second normal stress difference is undetermined. As explained in § 4.2, the profiles of $\phi$, $I$, and $\tilde {N}$ are unaffected by the wall roughness, and hence the remarks about the use of the friction b.c. (3.35) apply here also.

The profiles of $\mu$ vs $I$ are identical to those of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017), as (3.26) has been used in both cases (the red curves in figure 8). For $I$ in the range 0.002–0.2, (3.40) has been fitted to the DEM results. Hence the agreement is good (the red curves in figure 9). For smaller values of $I$, deviations are apparent.

4.4. Henann & Kamrin (Reference Henann and Kamrin2013)

For smooth walls, the scaled velocity profiles are qualitatively similar to those for the DEM results (the green curves in figure 3), but the values are underestimated near the centre and overestimated near the wall. For rough walls, the solution could not be obtained as discussed in § 4.4.1. For smooth walls, the scaled thickness of the shear layer is comparable to the DEM results (the green symbols in figure 4b).

The value of the slip length $l_u$ is approximately independent of $\bar {\phi }$ and $2\,W/d_p$, and comparable to the DEM results (the green symbols in figure 4d). Hence it appears that the slip velocity b.c. (3.33) can be used for this model, with $l_u$ as a material parameter.

The slip length for the fluidity $l_f$ is calculated using (3.58), and is shown in figure 11(a). The DEM results of $l_f$ for smooth walls do not vary with $2\,W/d_p$ and $\bar {\phi }$, which is not true for rough walls. For smooth walls, the value of $l_f$ obtained from the model is much smaller than that for the DEM results, and shows a very weak dependence on $2\,W/d_p$ and $\bar {\phi }$. Thus (3.58) can possibly be used as a b.c., but the discrepancy in the $l_f$ values shown in figure 11(a) suggests that it may not give good predictions of various quantities.

Figure 11.

Variation of $l_f$ (a) and $\tilde {x}_s = \mu _s \tilde {N}/\phi _p$ (b) with $\bar {\phi }$, $\tilde {f} = f/f_w$ (c) and $T_1$ and $T_2$ (d) with $\tilde {x} = x/W$ for the model of Henann & Kamrin (Reference Henann and Kamrin2013). Here, $l_f$ is calculated from (3.58), $T_1$ and $T_2$ are defined by (4.9a,b). The value of $\tilde {N}$ is matched to the DEM results. In (a,b), open and filled symbols correspond to smooth and rough walls, and $\vartriangle$, $\diamond$, $\circ$, and $\square$ are for $2\,W/d_p = 30$, 40, 50 and 80, respectively. In (a), the solid lines correspond to the DEM results, and the dotted lines to the model predictions. In (c,d), $2\,W/d_p = 40$; black dash-dotted line, $\bar {\phi } = 0.62$ (smooth wall, $\tilde {x}_s > 1$); black solid line, $\bar {\phi } = 0.6$ (smooth wall, $\tilde {x}_s$ slightly $<$ 1); black dashed line, $\bar {\phi } = 0.6$ (rough wall, $\tilde {x}_s$ much $<1$). The numerical (blue) and approximate (red) solutions for $\tilde {f}$ are shown in (c). Using the numerical solution, the profiles of $T_1$ (black) and $T_2$ (magenta) are shown in (d).

As $\phi$ is assumed to be a constant, the profiles are flat, in contrast to the DEM results (the green curves in figure 5). The variation of $I$ with $x/W$ is predicted fairly well by the model (the green curves in figure 6). The normal stress differences are zero (see (3.44)), in contrast to the DEM results. The scaled normal stress $\tilde {N}$ has been matched to the DEM values, and hence there is no discrepancy in the values of $\tilde {N}$ and $\delta$. Conversely, if the value $\delta$ is taken from independent experiments, the predicted value of $\tilde {N}$ may not match the DEM results. For smooth walls, $\mu$ vs $I$ profiles are predicted much better than the other models (the green curves in figure 8). In particular, the change in curvature of the profiles near the wall for some of the parameter values is captured by the model. Solutions could not be obtained for some parameter values (e.g. for rough walls), for reasons discussed below.

4.4.1. Reasons for the lack of solutions

Introducing the dimensionless variables $\tilde {x} = x/W$, $\tilde {W} = W/d_p$, $\tilde {u}_y = u_y/\sqrt {g W}$, $\tilde {N} = N/(\rho _p g W)$, $\tilde {f} = f/f_{w}$, $\tilde {f}_{w} = f_{w}/\sqrt {g/W}$, where $f_w = f(x = W)$, and substituting $\mu$ from (3.52), (3.53) can be rewritten as

(4.7)\begin{equation} \frac{ \mbox{d}^{2} \tilde{f}}{ \mbox{d} \tilde{x}^{2}} = \frac{\Delta \mu}{A^{2}} \left( \frac{\tilde{x}_s - \tilde{x}}{\tilde{x}_d - \tilde{x}} \right) \tilde{W}^{2} \tilde{f} + \frac{\Delta \mu}{I_0 A^{2}} \, \frac{\phi_p}{\tilde{N}^{3/2}} \tilde{W} \tilde{f}_{w} \tilde{x} \tilde{f}^{2} ,\end{equation}

or

(4.8)\begin{equation} \frac{ \mbox{d}^{2} \tilde{f}}{ \mbox{d} \tilde{x}^{2}} = T_1 + T_2 , \end{equation}

where $\tilde {x}_s \equiv \mu _s \tilde {N}/\phi _p$ and $\tilde {x}_d \equiv (\mu _s + \Delta \mu ) \tilde {N} / \phi _p$, and

(4.9a,b)\begin{equation} T_1 \equiv \frac{\Delta \mu}{A^{2}} \left( \frac{\tilde{x}_s - \tilde{x}}{\tilde{x}_d - \tilde{x}} \right) \tilde{W}^{2} \tilde{f}; \quad T_2 \equiv \frac{\Delta \mu}{I_0 A^{2}} \, \frac{\phi_p}{\tilde{N}^{3/2}} \, \tilde{W} \tilde{f}_{w} \tilde{x} \tilde{f}^{2} .\end{equation}

The values of $\mu _s$ and $\Delta \mu$ are such that $\tilde {x}_d > \tilde {x}_s$ and $\tilde {x}_d > 1$ for the ranges of $\tilde {N}$ and $\phi _p (= \bar {\phi })$ used here. However, $\tilde {x}_s$ may be $> 1$ for some cases (figure 11b). If $\tilde {x}_s > 1$, the right-hand side of (4.7) is non-negative. Using the b.c. (3.56), if $\tilde {f}(\tilde {x} = 0) > 0$, then (4.7) implies that $\tilde {f}$ increases monotonically throughout the domain. Hence a solution can be obtained that matches the b.c. (3.57) for smooth walls (the blue dash-dot curve in figure 11c).

If $\tilde {x}_s$ is slightly ${<} 1$, $T_1$ is ${<}0$ near the wall, whereas $T_2$ is ${>}0$ (see (4.9a,b)). If the latter dominates, a solution can be obtained (the solid blue curve in figure 11c). If $\tilde {x}_s$ is much ${<}1$, $T_1$ is ${<}0$ for $\tilde {x}_s < \tilde {x} \leqslant 1$, and may dominate $T_2$ in magnitude. Hence the right-hand side of (4.7) may be negative for some values of $\tilde {x}$. This causes the slope $\mbox {d} \tilde {f}/\mbox {d} x$ to decrease. If the slope decreases, $\tilde {f}$ may eventually become $< 0$, as in the case of the blue dashed curve in figure 11(c). This is an unrealistic result as (3.45) implies that $\tilde {f}$ must be $\geqslant 0$. For rough walls, $\tilde {x}_s$ is significantly ${<}1$ (figure 11b), and hence a realistic solution could not be obtained.

4.4.2. An approximate solution by the method of matched asymptotic expansion

An insight into the behaviour of solutions to (4.7) may be obtained by omitting $T_2$ as discussed below, even though the numerical solution shows that $|T_1| \ll T_2$ near the wall (figure 11d). Introducing a small parameter $\varepsilon = 1/\tilde {W} = d_p/W$, $\tilde {f}$ may be expanded in powers of $\varepsilon$ as

(4.10)\begin{equation} \tilde{f} = \tilde{f}_0 (\tilde{x}) + \varepsilon \tilde{f}_1(\tilde{x}) + \varepsilon^{2} \tilde{f}_2 (\tilde{x}) + \dots . \end{equation}

Substituting (4.10) in (4.7), omitting $T_2$ and matching the terms of ${O}$($\varepsilon ^{0}$), we obtain $\tilde {f}_0 = 0$. This does not satisfy the b.c. (3.57) at $\tilde {x} = 1$. Regarding this as the outer solution $\tilde {f}_{outer} = 0$ which is not valid near the wall, we seek an inner solution $\tilde {f}_{inner}$ by introducing a new coordinate

(4.11)\begin{equation} X = (1 - \tilde{x})/\varepsilon . \end{equation}

Expanding $\tilde {f} = \tilde {f}_{inner} \equiv F$ in powers of $\varepsilon$, we obtain

(4.12)\begin{equation} F = F_0 (X) + \varepsilon F_1 (X) + \varepsilon^{2} F_2 (X) + \dots.\end{equation}

Substituting (4.11) and (4.12) in (4.7), we obtain

(4.13)\begin{equation} \frac{ \mbox{d}^{2} F}{ \mbox{d} X^{2}} = \frac{\Delta \mu}{A^{2}} \, \frac{\tilde{x}_s - (1-\varepsilon X)}{\tilde{x}_d - (1-\varepsilon X)} \, F . \end{equation}

The b.c.s required to solve (4.13) are $F (X = 0) = 1$ and $F \rightarrow \tilde {f}_{outer} = 0$ as $X \rightarrow \infty$. Substituting (4.12) in (4.13), and matching the terms of O($\varepsilon ^{0}$), we obtain

(4.14)\begin{equation} \frac{ \mbox{d}^{2} F_0}{ \mbox{d} X^{2}} = m^{2} F_0 ,\end{equation}

where $m^{2} = (\Delta \mu /A^{2}) ((\tilde {x}_s - 1)/(\tilde {x}_d - 1))$, and the leading-order terms in the b.c. are $F_0 (X = 0) = 1$ and $F_0 \rightarrow \tilde {f}_{outer}$ as $X \rightarrow \infty$.

Case 1: for $\tilde {x}_s > 1$, $m^{2}$ is $> 0$ and the inner solution is

(4.15)\begin{equation} \tilde{f}_{inner} = F = F_0 + {O}(\varepsilon) = \exp ({-}m X) + {O}(\varepsilon) = \exp \left( -\frac{m}{\varepsilon} (1-\tilde{x}) \right) + {O} (\varepsilon), \end{equation}

which ensures the solution $\tilde {f}$ is $\geqslant 0$ in the whole domain. The exponentially decaying trend of (4.15) from the wall towards the centre agrees qualitatively with the numerical results (the blue dash-dot curve and the blue solid curve in figure 11c).

Case 2: for $\tilde {x}_s < 1$, $m^{2}$ is $< 0$ and the inner solution is

(4.16)\begin{equation} \tilde{f}_{inner} = F (X) = F_0 (X) + {O} (\varepsilon) = c_1 \cos (m X) + c_2 \sin (m X) + {O}(\varepsilon) ,\end{equation}

where $c_1$ and $c_2$ are constants. Using the b.c. $F_0(X = 0) = 1$, $c_1 = 1$. The b.c. $F_0 \rightarrow \tilde {f}_{outer}$ as $X \rightarrow \infty$ cannot be satisfied. Hence the solution (4.16) does not ensure $\tilde {f} \geqslant 0$ in the whole domain for the instances $\tilde {x}_s < 1$.

The above argument very crudely suggests that there is a possibility of obtaining a solution with $\tilde {f} < 0$ in a part of the domain if $|T_1| > T_2$. The numerical solution shows that $T_2$ is ${\gg }|T_1|$ near the wall, although both the terms are comparable in a zone far from the wall (figure 11d). This leads us to seek an approximate solution by omitting $T_1$.

Introducing a small parameter $\epsilon _* = 1/(\tilde {W} \tilde {f}_w)^{1/2}$ (the DEM results show that $\tilde {f}_w$ is of the same order as $\tilde {W}$), (4.7) can be rewritten as

(4.17)\begin{equation} \frac{ \mbox{d}^{2} \tilde{f}}{ \mbox{d} \tilde{x}^{2}} = \epsilon_*^{2} \frac{\Delta \mu}{I_0 A^{2}} \, \frac{\phi_p}{\tilde{N}^{3/2}} \, \tilde{x} \tilde{f}^{2} . \end{equation}

Expanding $\tilde {f}$ in powers of $\epsilon _*$

(4.18)\begin{equation} \tilde{f} = \tilde{f}_0 + \epsilon_* \tilde{f}_1 + \epsilon_*^{2} \tilde{f}_2 + \dots . \end{equation}

Substituting (4.18) in (4.17), and matching the terms of O($\epsilon _*^{0}$), we obtain

(4.19)\begin{equation} \tilde{f}_0 = 0 . \end{equation}

The solution (4.19) satisfies the b.c. $\mbox {d} \tilde {f}/\mbox {d} \tilde {x} (\tilde {x} = 0) = 0$, but not the b.c. $\tilde {f} (\tilde {x} = 1) = 1$. Hence we regard (4.19) as the outer solution $\tilde {f}_{outer}$.

To obtain the inner solution, let

(4.20)\begin{equation} \chi = \frac{1 - \tilde{x}}{\epsilon_*} . \end{equation}

Expanding $\tilde {f} = \tilde {f}_{inner} \equiv \varphi$ in powers of $\epsilon _*$, we obtain

(4.21)\begin{equation} \varphi = \varphi_0 (\chi) + \epsilon_* \varphi_1 (\chi) + \epsilon_*^{2} \varphi_2 (\chi) + \dots . \end{equation}

Substituting (4.21) in (4.17), we obtain

(4.22)\begin{equation} \frac{ \mbox{d}^{2} \varphi}{ \mbox{d} \chi^{2}} = \frac{\Delta \mu}{I_0 A^{2}} \, \frac{\phi_p}{\tilde{N}^{3/2}} (1 - \epsilon_* \chi) \varphi^{2} . \end{equation}

The b.c.s for (4.22) are $\varphi (\chi = 0) = 1$ and $\varphi \rightarrow \tilde {f}_{outer} = 0$ as $\chi \rightarrow \infty$. Substituting (4.21) in (4.22), and matching the terms of O($\epsilon _*^{0}$), we obtain

(4.23)\begin{equation} \frac{ \mbox{d}^{2} \varphi_0}{ \mbox{d} \chi^{2}} = \alpha \varphi_0^{2} , \end{equation}

where $\alpha = (\Delta \mu \phi _p)/(I_0 A^{2} \tilde {N}^{3/2})$. Multiplying both sides of (4.23) by $\mbox {d} \varphi _0 / \mbox {d} \chi$, and using the b.c.s $\varphi _0 (\chi = 1) = 1$ and $\varphi _0 \rightarrow 0$ at $\chi \rightarrow \infty$, the inner solution is given by (Zaitsev & Polyanin Reference Zaitsev and Polyanin2002)

(4.24)\begin{equation} \tilde{f}_{inner} = \varphi = \frac{1}{\left( \sqrt{\dfrac{\alpha}{6}} \chi + 1 \right)^{2}} + {O}(\epsilon_*) . \end{equation}

The approximate solution (the red curves in figure 11c) agrees reasonably well with the numerical solution (the blue curves), except when the fluidity is negative.

4.4.3. Use of $\phi = \phi (I)$ to improve the model

To permit variation of $\phi$ with $x$, we assume that $\phi = \phi (I)$, and use the form (3.40). Equations (3.30), (3.50a,b) and (3.53) are solved using the b.c.s (3.56) and (3.57). As in the other models, $\bar {\phi }$ and $\dot {M}$ are matched to the DEM results. For $\bar {\phi } = 0.61$, the $\phi$ profiles (the blue and red solid curves in figure 12b) agree better with the DEM results than the incompressible model (the dashed line). Because $\phi (x = 0) = \phi _c = 0.613$ at the centre, which is slightly ${>}\phi _p = \bar {\phi }$ for the incompressible model, we expect the scaled velocity to be smaller near the centre (the blue and red solid curves in figure 12a). The profiles of $I$ for both the models are qualitatively similar (figure 12c). These remarks also apply to the case $\bar {\phi } = 0.60$ (figure 12d–f).

Figure 12.

Variation of the scaled vertical velocity $-u_y/\sqrt {gW}$ (a), the solids fraction $\phi$ (b) and the inertial number $I$ (c) with $x/W$ for the model of Henann & Kamrin (Reference Henann and Kamrin2013): $\circ$, the DEM results; blue and red dashed line, $\phi = \textrm {const}. = \phi _p$; blue and red solid line, $\phi = \phi (I)$; $2\,W/d_p = 40$, red curves; $2\,W/d_p = 80$, blue curves. Here, $\bar {\phi } = 0.61$ (a–c), and 0.60 (d–f).

4.5. Dsouza & Nott (Reference Dsouza and Nott2020)

The value of the parameter $\ell$ is estimated by minimizing the root mean square error $N_{rms}$ between the normal stresses $N_{dem}$ and $N_{model}$, with

(4.25)\begin{equation} N_{rms} \equiv \sqrt{ \frac{1}{n_p} \sum_{i = 1}^{n_p} \left( \frac{N_{model} - N_{dem}}{N_{dem}} \right)^{2} } , \end{equation}

where $n_p$ is the number of bins used for the DEM. For the models, $N_{model} = \sigma _{xx}$ is independent of $x$, but $N_{dem}$ varies slightly. Here, $N_{rms}$ increases with $\ell /d_p$ (figure 13a), and hence there is no obvious choice for $\ell /d_p$. Clearly, $\ell$ should be small compared with the half-width $W$ of the channel and large enough so that the averaging volume used to average the flow rule contains a reasonable number of particles. The smallest value of $\ell /d_p$ used here is 3. Choosing a value of $N_{rms}$ that is 5 % larger than the value for $\ell /d_p = 3$, and averaging the results for $2\,W/d_p = 30$, 40 and 50, we obtain $\ell /d_p = 5$. A sphere of radius $5 \, d_p$ contains ${\approx }75$ particles. The value used by Dsouza & Nott (Reference Dsouza and Nott2020) is $\ell /d_p = 10$ but the reason behind this choice is not clear. It leads to a larger discrepancy compared with the DEM results (figure 13a). As mentioned in § 3.5, the dimensionless shear rate $\tilde {S} = S' \sqrt {W/g}$ is replaced by $\tilde {S} + \tilde {\epsilon }$, where $\tilde {\epsilon }$ is chosen as $10^{-4}$. The results are unaffected for $\tilde {\epsilon }$ in the range $10^{-6}{-} 10^{-2}$, and the relative error in the results with respect to $\tilde {\epsilon } = 10^{-4}$ is in the range $10^{-4}\,\%{-} 10^{-2}$ %.

Figure 13.

Variation of the root mean square error for the normal stress $N_{rms}$ (see (4.25)) with $\ell /d_p$ (a), solids fraction $\phi$ with $x/W$ (b) and the lower limit $\bar {\phi }_l$ for $\bar {\phi }$ with $2\,W/d_p$ and $\phi ^{0} \equiv \phi (x = 0)$ (c) for the model of Dsouza & Nott (Reference Dsouza and Nott2020). Here, $\ell$ is the effective radius of averaging volume (see (3.59)). In (b), the symbols represent the DEM results for smooth ($\circ$) and rough ($\square$) walls, and the curves represent the numerical solutions based on (3.60) (black dashed line) and (3.61) (black solid line) for $2\,W/d_p = 40$ and $\bar {\phi } = 0.62$. The scaled normal stresses $\tilde {N}$ predicted by (3.60) and (3.61) are 6.26 and 4.16 $\times 10^{4}$ , respectively. In (c), $\bar {\phi }_l$ corresponds to $\phi (x = W) = \phi _{min}$ (see (3.60)).

The results are sensitive to the choice of the functional form of $\varPi (\phi )$, the local value of the mean stress at a critical state. The two forms used here (3.60) and (3.61) lead to very different profiles of $\phi$ in the shear layer and also different estimates of the normal stresses (figure 13b). The local contribution of the critical state pressure $\varPi (\phi )$ diverges faster for (3.61) than (3.60). The ratio of $\varPi (\phi )$ calculated from (3.61) and (3.60) for, say, $\phi (x = 0) = 0.63$ is ${\approx }10^{3}$, which explains the much higher prediction of the scaled normal stress $\tilde {N}$ for (3.61). There is a dire need for data on $\varPi (\phi )$ at stress levels comparable to those that occur in laboratory-scale experiments.

Equation (3.64) could not be solved except for $\bar {\phi } = 0.62$. The reason for this behaviour is discussed in § 4.5.1. As the centreline velocity is forced to match the DEM value, the prediction is good near the centre (the magenta curves in figure 3a,e). However, the DEM results are significantly underestimated except close to the wall. Consequently, the predicted thickness of the shear layer is much higher (the magenta symbols in figure 4b). The value of $l_u$ calculated from (3.33) varies from 1 to 20 as $2\,W/d_p$ increases from 30 to 80 (figure 4d). Hence (3.33) is not a suitable b.c. For rough walls, the magnitude of the scaled velocity changes sign near the wall, which is unrealistic. The reason for this behaviour is discussed in § 4.2.

For smooth walls, the $\phi$ profiles agree well with the DEM results (the magenta curves in figure 5a,e). For rough walls, discrepancies are evident for intermediate values of $x/W$. For a fixed value of $\bar {\phi }$, $\phi$ profiles depend on $\phi (x = 0)$, which is matched to the DEM results. As $\phi (0)$ is found to be insensitive to the wall roughness (figure 5a,e), the predicted $\phi$ profiles are independent of the wall roughness, in contrast to the DEM results.

The normal stress differences are zero (the first of (3.62)). The scaled normal stress $\tilde {N}$ is much higher than the DEM results (the magenta symbols in figure 7e). The local model, obtained by setting $\ell = 0$ in (3.62)–(3.64), implies that the momentum balance (3.64) cannot be satisfied, and hence it does not predict the variation of $\phi$ in the current geometry. However, the local contribution to the normal stress may be readily evaluated using the values of $\phi$. For example, if $2\,W/d_p = 40$ and $\bar {\phi } = 0.61$, $\tilde {N}_{local}/\tilde {N} \approx 0.25$ at $x = 0$. Thus large normal stresses arise mainly from the non-local term involving $\ell ^{2}$ in the first of (3.64). As $2\,W/d_p$ increases, $\tilde {N}$ decreases, in contrast to the behaviour of the DEM and the other models. Some insight can be gained from an approximate solution discussed in § 4.5.2.

4.5.2. Approximate solution

The structure of (3.64) and the b.c. permits the construction of a Taylor series solution. Expanding $\phi$ in a Taylor series about $x = 0$, and using the symmetry condition $\mbox {d}\phi /\mbox {d}x (x = 0) = 0$

(4.26)\begin{equation} \phi = \phi^{0} + \left.\frac{1}{2} \,\frac{\mbox{d}^{2}\phi}{\mbox{d}\kern0.7pt x^{2}} \right|_0 \; x^{2} + {{O}}(x^{4}) . \end{equation}

Using the first of (3.64)

(4.27)\begin{equation} \left.\frac{\mbox{d}^{2}\phi}{\mbox{d}\kern0.7pt x^{2}} \right|_0 = \left.\frac{\varPi - N}{\ell^{2} \, \dfrac{\mbox{d}\varPi}{\mbox{d}\phi}} \right|_0 . \end{equation}

Integrating (4.26) with respect to $x$ from 0 to $W$, we obtain

(4.28)\begin{equation} \bar{\phi} = \phi^{0} + \left.\frac{1}{6} \, \frac{\mbox{d}^{2}\phi}{\mbox{d}\kern0.7pt x^{2}} \right|_0 W^{2} . \end{equation}

Using (4.26) and (4.28), the approximate solution for $\phi$ is

(4.29)\begin{equation} \phi = \phi^{0} + 3 \left(\bar{\phi} - \phi^{0} \right) \left(\frac{x}{W} \right)^{2} . \end{equation}

To obtain the normal stress $N$, we use (4.27) and (4.28)

(4.30)\begin{equation} N = \varPi_{\phi = \phi^{0}} - 6 \left(\frac{\ell}{W}\right)^{2} \left(\bar{\phi} - \phi^{0} \right) \left.\frac{\mbox{d}\varPi}{\mbox{d}\phi} \right|_{\phi = \phi^{0}} . \end{equation}

Setting $q = \mbox {d}u_y/\mbox {d}x$, (3.67) reduces to

(4.31)\begin{equation} \frac{ \mu_*}{q} \left( N q - \ell^{2} \varPi \left.\frac{\mbox{d}^{2} q}{\mbox{d} x^{2}} \right) \right\vert_x = \rho_p g \int_{0}^{x} \phi \, \mbox{d}\kern0.7pt x' . \end{equation}

Expanding $q$ in a Taylor series about $x = 0$, and using the conditions $q(x = 0) = 0$ and $\mbox {d}^{2}q/\mbox {d}x^{2}(x = 0) = 0$ (see (3.66))

(4.32)\begin{equation} \frac{\mbox{d}u_y}{\mbox{d}\kern0.7pt x} = q = \left.\frac{\mbox{d}q}{\mbox{d}\kern0.7pt x}\right|_0 x + \frac{1}{3!} \left. \frac{\mbox{d}^{3} q}{\mbox{d} x^{3}}\right|_0 x^{3} + {{O}} (x^{5}) . \end{equation}

Now, applying limit $x \rightarrow 0$ to (4.31), L’Hopital's rule implies that

(4.33)\begin{equation} \frac{N \left.\dfrac{\mbox{d} q}{\mbox{d} x} \right|_0 - \ell^{2} \varPi_{\phi = \phi^{0}} \left.\dfrac{\mbox{d}^{3} q}{\mbox{d} x^{3}}\right|_0}{\left.\dfrac{\mbox{d} q}{\mbox{d} x} \right|_0} = 0 . \end{equation}

Assuming that ${\mbox {d} q}/{\mbox {d} x} |_0 \neq 0$, (4.33) reduces to

(4.34)\begin{equation} N \left. \frac{\mbox{d} q}{\mbox{d} x} \right|_0 - \ell^{2} \varPi_{\phi = \phi^{0}}\left. \frac{\mbox{d}^{3} q}{\mbox{d} x^{3}}\right|_0 = 0 . \end{equation}

Solving for ${\mbox {d}^{3} q}/{\mbox {d} x^{3}}|_0$ from (4.34), integrating (4.32) and using $u_y (x = 0) \equiv u_y^{0}$, we obtain

(4.35)\begin{equation} u_y = u_y^{0} + \left( \frac{1}{2} x^{2} + \frac{1}{24\ell^{2}} \, \frac{N}{\varPi_{\phi = \phi^{0}}} x^{4} \right) \left. \frac{\mbox{d} q}{\mbox{d} x} \right|_0 . \end{equation}

Using (3.32), (4.29) and (4.35), we obtain

(4.36)\begin{equation} \left.\frac{\mbox{d} q}{\mbox{d} x} \right|_0 = \frac{ \left(2 \rho_p B ({-}u_y^{0}) \bar{\phi}W \right) - \dot{M}}{2 \rho_p B \left( \dfrac{3}{10}\bar{\phi}W^{3} + \dfrac{1}{56\ell^{2}} \dfrac{N}{\varPi_{\phi = \phi^{0}}} \bar{\phi} W^{5} + \dfrac{1}{60} \phi^{0} W^{3} - \dfrac{1}{105\ell^{2}} \dfrac{N}{\varPi_{\phi = \phi^{0}}} \phi^{0} W^{5}\right)} . \end{equation}

The functions (4.29), (4.30), (4.35) and (4.36) provide an approximate solution. The approximate solution for the velocity is close to the numerical solution, except near the wall (figure 3a,e). Owing to the Taylor series expansion, the approximate solution can be constructed for $\bar {\phi } \leqslant 0.62$, whereas the numerical solution cannot be obtained for $\bar {\phi } < 0.62$ (see § 4.5.1). For fixed values of $\bar {\phi }$ and $x/W$, (4.29) shows that $\phi (x/W)$ is independent of $W$ (the blue curves in figure 5). The numerical solution (the magenta curves in figure 5a,e) reflects this behaviour except very close to the wall. This is in keeping with the trends shown by the DEM results.

In the absence of the non-local term involving $(\ell /W)^{2}$ in (4.30), $N$ is approximately a constant as $\phi ^{0}$ does not vary much with changes in $\bar {\phi }$ and $2\,W/d_p$ for the DEM results. When the non-local term is included, (4.30) shows that $N$ increases as $W$ decreases provided $\bar {\phi } < \phi ^{0}$, in contrast to the DEM results (the blue symbols in figure 7e). The numerical solution (magenta) is close to the approximate solution. The slope $\mbox {d} N/\mbox {d} \bar {\phi }$ of (4.30) is negative, in contrast to the DEM results. Thus higher order terms may be required to get a more accurate solution.