3.1 $\unicode[STIX]{x1D6F1}$ theorem

The mass flow rate out of the silo, $Q$ , depends on the density $\unicode[STIX]{x1D70C}$ , the gravity $g$ and the geometrical parameters: width $L$ , orifice height $D$ , thickness $W$ , the filling height $h_{p}$ and the grain size $d$ . Standard application of dimensional analysis or $\unicode[STIX]{x1D6F1}$ theorem gives us the relation between non-dimensional numbers (eight quantities, three units, which gives five non-dimensional numbers). The flux must scale like $\unicode[STIX]{x1D70C}\ell ^{2}\sqrt{g\ell }$ where $\ell$ is any of $D,W,h_{p},d$ (and $\ell ^{2}\sqrt{\ell }$ any combination of these lengths), and this flux must also be a function of the geometric ratios e.g. $D/W,d/D,h_{p}/W,L/D$ (note that if we use the mass of grains in the silo, this will give us an extra parameter with dimension, which can be reduced to a new dimensional number, the packing fraction $\unicode[STIX]{x1D719}_{b}$ ). Classical experiments on silos show that the flux is independent from the width $L$ (if large enough), from the filling height $h_{p}$ (if large enough) and from the grain size $d$ (if small enough). We can therefore neglect the corresponding geometric ratios. This then reduces $Q/(\unicode[STIX]{x1D70C}\ell ^{2}\sqrt{g\ell })$ (where $\ell$ is any of $D,W$ ) to a function of the aperture aspect ratio ${\mathcal{A}}=D/W$ only.

Figure 4.

Experimental mass flow rate $Q$ in the rectangular silo versus the height of orifice $D$ for two sizes of particles $d$ . (a) $W=40$  mm, the dashed (respectively full) line represents the equation $Q=c_{1}D^{3/2}$ where $c_{1}=2.79~\text{g}~\text{s}^{-1}~\text{mm}^{-3/2}$ (respectively $c_{1}=2.65~\text{g}~\text{s}^{-1}~\text{mm}^{-3/2}$ ) is obtained using a least-squares fit. (b) $W=3.5$  mm, the dashed (respectively full) line represents the equation $Q=c_{2}D$ where $c_{2}=0.58~\text{g}~\text{s}^{-1}~\text{mm}^{-1}$ (respectively $c_{2}=0.43~\text{g}~\text{s}^{-1}~\text{mm}^{-1}$ ) is obtained using a least-squares fit.

If $W$ is large, the problem becomes bi-dimensional and $D/W$ has no influence anymore. The velocity then scales like $\sqrt{gD}$ , the size of the aperture being scaled by $D$ , and the flux per transverse unit $Q/W$ is $\unicode[STIX]{x1D70C}D\sqrt{gD}$ . This gives the equivalent for this rectangular geometry of the Hagen–Beverloo relation (Beverloo, Leniger & Van de Velde Reference Beverloo, Leniger and Van de Velde1961; see a translation of the original article of Hagen in Tighe & Sperl Reference Tighe and Sperl2007)

This scaling relation was recovered experimentally for a bottom aperture, spanning the width of a rectangular silo, by severals authors (Choi, Kudrolli & Bazant Reference Choi, Kudrolli and Bazant2005; Benyamine et al. Reference Benyamine, Djermane, Dalloz-Dubrujeaud and Aussillous2014). This behaviour is also observed in our experimental set-up with a side aperture. It can be seen in figure 4(a) where the flow rate is plotted as a function of the outlet length $D$ for two particle diameters and for the larger silo thickness $W=40$ mm. For a given particle diameter the data are well described by (3.1), see the full line and the dashed line in the figure. Therefore, following the $\unicode[STIX]{x1D6F1}$ theorem, the relevant normalisation for the flow rate seems to be $\unicode[STIX]{x1D70C}\sqrt{gW^{5}}$ , giving for the former relationship

In the following section this scaling will be compared to the data obtained for the rectangular silo.

3.2 Experimental results for the rectangular silo

Figure 5.

Experiments in the rectangular silo for a given particle diameter, $d=190~\unicode[STIX]{x03BC}\text{m}$ , and various silo widths, $W$ . (a) Dimensionless mass flow rate $Q/(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}\sqrt{gW^{5}})$ as a function of $D/W$ . The dashed line represents equation (3.2) with $c_{D}=0.51$ , the dashed-dotted line represents equation (3.3) with $c_{W}=0.68$ and the solid line represents (3.4) with the same parameters. (b) Corresponding dimensionless discharge velocity $u_{D}/\sqrt{gd}$ as a function of $D/d$ . The dashed lines represent equation (3.5) with the same parameters.

Figure 5(a) represents the dimensionless flow rate $Q_{{\mathcal{A}}}=Q/(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}\sqrt{gW^{5}})$ versus the aperture aspect ratio ${\mathcal{A}}=D/W$ , for various thicknesses $W$ and lengths $D$ of the orifice and for a given diameter of particle $d=190~\unicode[STIX]{x03BC}\text{m}$ . In this representation the data collapse as suggested by the $\unicode[STIX]{x1D6F1}$ theorem.

As expected, for large thicknesses ( $D/W\ll 1$ ), $Q_{{\mathcal{A}}}$ follows a power law with an exponent $3/2$ corresponding to the Hagen–Beverloo relation (3.2), see the dashed line in the figure. However, for large values of $D/W$ (i.e. thin enough silos), an important experimental result is that the dimensionless flow rate of particles depends linearly on the aperture aspect ratio (see the dashed-dotted line in figure 5 a):

This can also be seen in figure 4(b) where the flow rate is plotted versus the aperture length $D$ for the smallest thickness explored $W=3.5$ mm and for two particles diameters $d$ . For a given particle diameter the flow rate indeed exhibits a linear trend $Q\propto D$ . The transition between these two regimes occurs around a specific aperture ratio ${\mathcal{A}}_{c}\approx 2$ . According to these two asymptotic regimes, we can thus propose a fitting formula for the flow rate:

As illustrated in figure 5(a) (solid lines), this formula describes well the flow rate dependency on aperture dimensions in any regime, with the fitting parameters $c_{D}=0.51$ and $c_{W}=0.68$ .

These two regimes can be interpreted in term of the (horizontal) discharge velocity $u_{D}=Q/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}WD$ , which can be expressed using (3.4) as

This relation again describes the data very well whatever the silo width as can be seen in figure 5(b) (dashed lines). The first regime corresponds to the classical Hagen–Beverloo relation with a velocity which scales with the aperture length:

Whereas in the second regime the discharge velocity scales with the aperture width:

This gives a velocity scaling with the smallest of the two lengths $W$ and $D$ .

Historically the Hagen–Beverloo relation was explained with the concept of a ‘free-fall arch’ located at the outlet, from which the particles fall freely. Several recent experiments and numerical simulations question this hypothesis. From an experimental point of view, Janda et al. (Reference Janda, Zuriguel and Maza2012) and Rubio-Largo et al. (Reference Rubio-Largo, Janda, Maza, Zuriguel and Hidalgo2015) showed that the granular medium remains dense, with a small dilation, close to the outlet and that the particles do not undergo a free fall. These observations were validated using discrete simulation (Rubio-Largo et al. Reference Rubio-Largo, Janda, Maza, Zuriguel and Hidalgo2015). Moreover Navier–Stokes simulations, assuming a continuum frictional rheology of the granular media, have been successfully used to recover the Hagen–Beverloo scaling as in Staron et al. (Reference Staron, Lagrée and Popinet2012, Reference Staron, Lagrée and Popinet2014), Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015) and Davier & Bertails-Descoubes (Reference Davier and Bertails-Descoubes2016), with very different numerical methods.

The fact that we recover the Hagen–Beverloo relation with a vertical aperture also tends to refute this concept, as pointed out by Sheldon & Durian (Reference Sheldon and Durian2010). Based on these observations, we performed a continuum simulation, using our Navier–Stokes solver with the granular rheology, to test whether we can reproduce at least qualitatively the experimental behaviour.

3.3 Comparison with continuum numerical simulations

Figure 6.

Continuum simulations results in two dimensions and in pseudo three dimensions for various $W$ with $H=360d$ , $L=90d$ and $D=45d$ . (a) Flow rate $Q$ normalised by $\sqrt{gW^{5}}$ as a function of $D/W$ . The dashed line represents equation (3.2) with $c_{D}=0.76$ , the dashed-dotted line represents equation (3.3), with $c_{W}=1.49$ , and the solid line represents equation (3.4), with the same parameters. (b) Mean horizontal velocity at the orifice $\overline{u}/\sqrt{gd}$ , as a function of $D/d$ . The dashed lines represent equation (3.5), with the same parameters.

We first performed continuum numerical simulations in the 2-D case. To compare with the experimental results we have plotted in figure 6(a) the dimensionless flow rate $Q_{{\mathcal{A}}}=Q/(\sqrt{gW^{5}})$ versus the aperture aspect ratio ${\mathcal{A}}=D/W$ , using the same series of data but rescaled several times with each width $W$ used for the pseudo-3-D simulations. In a pure bi-dimensional flow (red circles), the Navier–Stokes simulations recover the Hagen–Beverloo scaling, equation (3.2) (dashed line). This suggests that the $\unicode[STIX]{x1D707}(I)$ fluid rheology provides consistent results for the discharge of a silo with a side-located aperture.

Then, to mimic the effects of the lateral walls from a Hele-Shaw point of view, we added the friction term proportional to the pressure and aligned with the velocity, equation (2.4), and we varied the pseudo-width $W$ (see the different symbols in figure 6 a). We first notice that when the aperture ratio is small, similar to the experimental case, the data of the pseudo-3-D simulations are superimposed onto the 2-D simulations, and we recover the Hagen–Beverloo scaling, equation (3.2). When increasing the aperture ratio, we observe a departure from this scaling towards a regime where $Q_{{\mathcal{A}}}\propto {\mathcal{A}}$ (see the dashed-dotted line). We can see that we do not completely reach this regime due to the numerical upper limit of achievable values for $\boldsymbol{f}_{w}$ , however the data are well described by (3.4) with $c_{D}=0.76$ and $c_{W}=1.49$ (solid line).

Following the experimental section, the dimensionless mean horizontal velocity at the outlet, $\bar{u}/\sqrt{gd}$ , is plotted versus the dimensionless outlet length $D/d$ in figure 6(b). We observe the same trends as for the experimental data: in the regime controlled by the outlet length $D$ , the velocity tends toward the Hagen–Beverloo scaling ((3.6) and 2-D data) whereas when the ratio ${\mathcal{A}}$ increases, the velocity tends toward the regime controlled by the silo width given by (3.7). Again each series of data for a given $W$ are fairly well described by (3.5) (dashed lines).

In the continuum simulation, the width of the silo appears only in the term describing the side wall friction, which suggests that the second regime is controlled by this term. It is interesting to note that in the regime dominated by the lateral friction the flow rate per unit length at a given $D$ is lower than in the first regime. Even if we cannot fully reach the second regime ( ${\mathcal{A}}\gg 1$ ), we carry on with the analysis and present some comparisons of the details of the internal fields at the limit of the numerical model.

Figure 7.

Velocity field (a,b) and streamlines (c,d) for the experiment (a,c) with $D=20$  mm,  $d=1129~\unicode[STIX]{x03BC}\text{m}$ and $W=[3.5,10,20,40]$  mm, and for the continuum simulation (b,d) with $D=56.25d$ and $W=[13.5d,45d,180d]$ .

Figure 8.

Limit of the stagnant zone for various $D$ : (a) experiments with $W=5$  mm and (b) continuum simulations with $W=17.4d$ .

Figure 7 shows the velocity field and the streamlines of both the experimental and the numerical runs for increasing silo width. The numerical velocity field (figure 7 b) is qualitatively very similar to the experimental velocity field (figure 7 a). In both configurations, the streamlines tend to vertical lines upward from the orifice at a distance which decreases when $W$ increases. On the same figures, we observe that the flowing zone is found to be thinner when the lateral friction increases (i.e. when $W$ decreases). The limit of the stagnant zone for various outlet sizes $D$ is plotted in figure 8 both for the experiment (a), and the continuum simulation (b), in the regime dominated by the lateral friction. Interestingly this position does not depend on $D$ in either cases, but depends strongly on $W$ .

Figure 9.

Continuum simulation results: angle of inclination $\unicode[STIX]{x1D703}$ of the central streamline at the orifice relative to vertical for various $W$ (a) versus $D/L$ , (b) versus $D/W$ . The dashed line represents equation (3.11) with $\unicode[STIX]{x1D6FE}_{1}=0.56$ and $\unicode[STIX]{x1D6FE}_{2}=0.25$ .

We may explain these behaviours at small $W$ by looking at the Navier–Stokes equations. Let us consider a steady flow (implying that we neglect inertia). As the friction at the walls increases when $W$ decreases, the friction term (of magnitude $\unicode[STIX]{x1D707}_{w}p/W$ ) will be at most as large as gravity ( $\unicode[STIX]{x1D70C}g$ ). The gradients of the stress tensor are of order $\unicode[STIX]{x1D707}_{s}p/L$ . This order of magnitude may be rewritten as $(\unicode[STIX]{x1D707}_{s}/\unicode[STIX]{x1D707}_{w})(W/L)\unicode[STIX]{x1D707}_{w}p/W$ . It is clearly smaller than $\unicode[STIX]{x1D707}_{w}p/W$ as $(\unicode[STIX]{x1D707}_{s}/\unicode[STIX]{x1D707}_{w})=O(1)$ and $(W/L)\ll 1$ . Hence, as the wall friction increases to balance gravity, the momentum equation can be approximated as:

where $\boldsymbol{e}$ is a unit vector in the main direction of the force resulting from internal friction. The velocity $\boldsymbol{u}$ is aligned with gravity at first order, the more so as the thickness decreases. This is noticeable in figure 7(c,d) where for small $W$ the streamlines are clearly more vertical than for larger $W$ . To quantify this effect, we have plotted the angle of inclination $\unicode[STIX]{x1D703}$ of the central streamline at the orifice relative to the vertical in figure 9. If we note $u_{0}$ the horizontal velocity at the centre of the outlet and $U_{0}$ the norm of the velocity at this position, we can write $\sin \unicode[STIX]{x1D703}=u_{0}/U_{0}$ (note that $u_{0}$ is proportional to the previously defined discharge velocity $u_{D}$ , as we will see in § A.2). We observe in figure 9(a) that this angle decreases slightly with $D$ and strongly with $W$ as expected. More interestingly, if plotted versus $D/W$ (figure 9 b), the data superimpose. In order to study the scaling relation of $\sin \unicode[STIX]{x1D703}$ , we then turn to the dependency of the norm of the velocity $U_{0}$ and horizontal velocity $u_{0}$ on $W$ and $D$ . In figure 10 we have represented the norm of the velocity $U_{0}/\sqrt{gd}$ on the central streamline at the orifice versus $D/d$ for various $W$ . Surprisingly we can see that the norm of the velocity does not significantly depend on $W$ . The data can be described by

with $c_{U}=1.2$ . This suggests that the kinetic energy always scales like $\unicode[STIX]{x1D70C}gD$ whatever the lateral friction. For the horizontal velocity $u_{0}$ , since the flow rate $Q\propto u_{0}WD$ , from (3.4), we can write

We can thus obtain

Moreover we see that the data in figure 9(b) are well described by this formulation with $\unicode[STIX]{x1D6FE}_{1}=0.56$ and $\unicode[STIX]{x1D6FE}_{2}=0.25$ , with the fitting parameter $\unicode[STIX]{x1D6FE}_{2}=0.25$ corresponding to $(c_{D}/c_{W})^{2}=0.26$ .

Figure 10.

Continuum simulation results: norm of the velocity $U_{0}/\sqrt{gd}$ on the central streamline at the orifice versus $D/d$ for various $W$ . The dashed line represents equation (3.9) with $c_{U}=1.2$ .

Figure 11.

Profiles of the top surface of granular materials for various $W$ at constant time steps: (a) experiments with $D=25$  mm,  $d=190~\unicode[STIX]{x03BC}\text{m}$ , time interval $\unicode[STIX]{x0394}t_{e}=0.5$  s, (b) continuum simulation with $D=56.25d$ and time interval $\unicode[STIX]{x0394}t_{s}/\sqrt{d/g}=9.5$ .

Finally figure 11 shows the time evolution of the top surface of granular material for various silo widths for (a) the experiment and (b) the continuum simulation. From experiments we can see that the surface of particles begins to tilt in the early stage for a small width $W$ , whereas it remains constant till the last stage for the largest $W$ . The tilted interface exhibits a slope starting from the wall opposite to the outlet and reaches a flattened surface, or sometimes even a small bump, on the wall containing the outlet. The slope of the interface is higher than the angle of repose for $W=3.5$  mm, and smaller than the angle of repose for $W=10$  mm. In continuum simulations, we explored a smaller range of $W$ , the smallest $W$ case $W=13.5d$ , is comparable to the case $W=20$  mm in the experiments, but we clearly recover the same trend: the larger the thickness $W$ , the longer the interface remains symmetrical. Note also that the no-slip conditions imposed at the wall are not exactly the same as in the experiments, where a sliding velocity is observed. Nevertheless, the same qualitative profiles are obtained. This behaviour is consistent with the previous observation: for large $W$ , the flow far from the outlet is symmetrical. For small $W$ , the lateral friction breaks this symmetry and localises the flow on the side of the outlet, which inclines the surface in this direction.

3.4 Experimental results in the cylindrical silo

The previous results were given for a simplified geometry, with a rectangular silo and an outlet which spans the width of the silo. Experimentally we also performed measurements of the flow rate for a cylindrical silo, with an outlet located at its side as schematised in figure 1(b). This situation generates a flow in a fully 3-D geometry, and is of practical interest. To characterise the different roles played by the length $D$ and the width $W$ of the outlet, we have plotted in figure 12 the dimensionless mass flow rate $Q/(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}\sqrt{gW^{5}})$ as a function of the aperture aspect ratio $D/W$ , as done in the rectangular configuration.

Figure 12.

Experimental results in the cylindrical silo: dimensionless mass flow rate $Q/(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}\sqrt{gW^{5}})$ as a function of the aperture aspect ratio $D/W$ . The dashed line represents equation (3.2) with $c_{D}=0.49$ , and the dashed-dotted line represents equation (3.3) with $c_{W}=0.46$ .

As suggested by the $\unicode[STIX]{x1D6F1}$ theorem in § 3.1 the data superimpose in this representation. Moreover, we observe the same behaviours as in rectangular silos. For large $W$ , the flow rate follows the Hagen–Beverloo relation ((3.2), dashed line in the figure) and for small $W$ , the flow rate follows (3.3) (dashed-dotted line in the figure) suggesting that the flow is dominated by the lateral friction. Over the whole range, the data are well correlated by (3.4) (solid line in the figure) with $c_{D}=0.49$ and $c_{W}=0.46$ . The fitting parameter $c_{w}$ is found to be slightly smaller than that in the rectangular silo. This suggests that the dissipation is higher in the cylinder, and more specifically in the regime controlled by the lateral friction where we can assume that the friction with a smooth lateral wall is smaller than that with an erodible granular media. Consequently, the transition between the two regimes occurs for a smaller aperture aspect ratio ${\mathcal{A}}_{c}\approx 1$ .

3.5 Dependency of the flow rate on the particle diameter

The previous results are shown for a given small particle diameter $d=190~\unicode[STIX]{x03BC}\text{m}$ . In this section we focus on the influence of the particle diameter. We have shown that the flow rate of discharge of a silo with a lateral outlet depends on the aperture aspect ratio ${\mathcal{A}}=D/W$ , and exhibits two regimes of flow described by (3.2) for small ${\mathcal{A}}$ and by (3.3) for large ${\mathcal{A}}$ . When the particle diameter is varied, these scaling relations seem to remain valid, as can been seen in figure 4(a,b) where the flow rate is plotted versus $D$ for two particle diameters, $d=190~\unicode[STIX]{x03BC}\text{m}$ and $d=1129~\unicode[STIX]{x03BC}\text{m}$ , in the two regimes. In each regime the same scaling is found for the two grain sizes, however we observe a shift of these relations towards smaller flow rates when the particle diameter increases. This behaviour can be seen in figure 13(a,c) where the dimensionless flow rate $Q_{{\mathcal{A}}}=Q/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}\sqrt{gW^{5}}$ is plotted as a function of $D/W$ for both the rectangular silos and the cylindrical silo for all the batch of particles we used. The two regimes are recovered whatever the particle diameter, see the dashed line and the dashed-dotted line. However there is a significant scattering of the data above the relation fitted for the smallest particle size $d=190~\unicode[STIX]{x03BC}\text{m}$ .

Figure 13.

Mass flow rate $Q$ as a function of $D/W$ for various $d$ for (a,c) rectangular silos and (b,d) cylindrical silos, normalised in (a,b) by $(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}\sqrt{gW^{5}})$ and in (c,d) by $\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}\sqrt{gW^{5}}G_{D}(D/d)G_{W}(W/d)$ (see (3.12)) with the fitting parameters $\unicode[STIX]{x1D6FD}_{D}=\unicode[STIX]{x1D6FD}_{W}=0.1$ , $\unicode[STIX]{x1D6FC}_{D}=1$ and $\unicode[STIX]{x1D6FC}_{W}=0.46$ (see appendix A). The dashed and dashed-dotted lines are the same as in figures 5 and 12.

To model the dependency of the flow rate on the particle diameter, we follow the work of Janda et al. (Reference Janda, Zuriguel and Maza2012). In this framework we suppose that the particle size has to be accounted for through the number of beads in the aperture via a geometrical function, significant when this number is small. We take the two dimensions of the aperture into account and thus we suppose that the flow rate depends on two geometrical functions based on the number of beads in the aperture length $D/d$ and in the aperture width $W/d$ . The details of this approach and its validation using a 2-D discrete numerical simulation are given in the appendix. This yields the following expression for the flow rate:

with $G_{D}(D/d)=[1-\unicode[STIX]{x1D6FC}_{D}\text{e}^{-\unicode[STIX]{x1D6FD}_{D}(D/d)}]$ , $G_{W}(W/d)=[1-\unicode[STIX]{x1D6FC}_{W}\text{e}^{-\unicode[STIX]{x1D6FD}_{W}(W/d)}]$ and $F(D/W)=\sqrt{1/[1+\unicode[STIX]{x1D6FE}_{2}D/W]}$ .

Figure 13 shows the normalised flow rate using the geometrical functions, $Q/[\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}_{b}G_{D}(D/d)G_{W}(W/d)\sqrt{gW^{5}}]$ , as a function of the aperture aspect ratio $D/W$ for the whole range of particle diameters for (c) rectangular silos and (d) cylindrical silos. In this representation the data collapse and can be well adjusted by the same asymptotic relations in two regimes as in figures 5 and 12.

We propose additional analysis in the appendix. The 2-D discrete numerical simulation shows that the velocity at the outlet may follow the same geometrical relation than the volume fraction. As the continuum simulation described correctly most of the observations on the discharge of a silo with a lateral aperture for a given particle diameter we wonder if the continuum simulation is able to recover partially the dependency on the particle size of the flow rate. Indeed, the $\unicode[STIX]{x1D707}(I)$ rheology contains the information on the particle diameter in the definition of the inertial number $I=d\sqrt{2}D_{2}/(\sqrt{p/\unicode[STIX]{x1D70C}})$ . Even though a dependency on the particle diameter is observed in the 2-D case, the profiles present a larger asymmetry than the one predicted by the 2-D discrete simulations. This can be due to the fact that the volume fraction has been assumed to be uniform, whereas the complete $\unicode[STIX]{x1D707}(I)$ -rheology, that includes a $\unicode[STIX]{x1D719}(I)$ relation, would predict a dilation of the granular medium when it is sheared.