Influence of particle size in silo discharge

Recently Janda et al. [Phys. Rev. Lett. 108, 248001 (2012)] reported an experimental study where it was measured the velocity and volume fraction fields of 1 mm diameter stainless steel beads in the exit of a two-dimensional silo. In that work, they proposed a new expression to predict the flow of granular media in silos which does not explicitly include the particle size as a parameter. Here, we study if effectively, there is not such influence of the particle size in the flux equations as well as investigate any possible effect in the velocity and volume fraction fields. To this end, we have performed high speed motion measurements of these magnitudes in a two-dimensional silo filled with 4 mm diameter beads of stainless steel, the same material than the previous works. A developed tracking program has been implemented to obtain at the same time both, the velocity and volume fraction. The final objective of this work has been to extend and generalize the theoretical framework of Janda et al. for all sizes of particles. We have found that the obtained functionalities are the same than in the 1 mm case, but the exponents and other fitting parameters are different.


Introduction
Despite it has been studied since years ago, the discharge of silos through holes is still an open problem in granular matter physics.It is well known that the flow rate of particles W is a power law of the outlet size D = 2R, so that the exponent depends on the dimensionality of the system.The most widely expression used to fit that magnitude was proposed by Beverloo et al. [1] in 1961 from a dimensional analysis, where C is a fitting constant, n is the dimensionality and ρ b = ρφ b is the bulk density of the material, being φ b the bulk packing fraction and ρ the density of the pure material.The factor √ g arises from the concept of free fall arch developed by Brown and Richards [2].This idea assumes the existence of a parabolic [3] or hemispheric [4] region over the outlet above which the particles are subject to contact stress and under which they fall freely under the action of gravity [5].Instead of D, Eq.(1) includes a reduced outlet size D − kd involving a term proportional to the particle size.This term, also supported by Brown and Richards [2] and known as empty annulus, is associated with the existence of a space of indefinite size near the outlet borders where particles are forbidden to pass through.Although the ideas above are difficult to justify from a physical point of view, they have been used in diverse systems [6] and Beverloo's expression has been proved to be useful to fit the flow rate in different circumstances [7].
e-mail: dgella@alumni.unav.esHowever, over the years the community has found some inaccuracies in the expression of Beverloo.In 1971, Hunginton and Roney [8] reported that materials dilated differently when were near the exit and defined a new flowing density ρ f to be inserted in Eq.( 1) instead ρ b .This ρ f was obtained as the ratio of the mass flow rate to the volumetric flow rate calculated at the top of the silo.In fact, packing fractions φ vary considerably from point to point within the silo [9], so there is a potential error up to 20% in the flow rate if the choice of ρ f is not the appropiate.
In 2007, Mankoc et al. [10] showed that the flow rate is not possible to be fitted to the Beverloo's expression if a wide span of outlet sizes is regarded.So, they proposed a new empirical expression, where b and C are fitting parameters.Apart from settting k to 1, the main novelty of Eq.( 2) is the inclusion of an exponential term, which was suggested to arise from a dependence of ρ f (at the exit) on the outlet size.
In 2012, Janda et al. [11] reported an experimental work where they measured vertical velocity v z and packing fraction φ profiles for several outlet lengths at the exit of a two-dimensional silo by means of high speed image adquisition.The authors found that these profiles were self-similar and could be collapsed by means of the expressions: where ν and μ are fitting exponents, and x is the horizontal position measured from the center of the outlet.v c and φ c represent, respectively, the vertical velocity and packing fraction measured at the center of the outlet and are fitted to the following expressions in function of R: (5) where γ, α 1 , α 2 and φ ∞ are fitting parameters.Despite it has been recently proved [12] that the free fall arch idea is quite different from what it was initially thought, the square root dependence of the velocity on gravity is still recovered by Eq. ( 5).Furthermore, the fact that data are well fitted by Eq.( 6) is the experimental proof of the flowing density dependence on R with the form suggested in [10].By using mass conservation arguments, the mass flow rate could be calculated assuming density and velocity profiles as continuum functions by integrating ρφ(x)v z (x)dx.Thus, the proposed expression of the mass flow rate for a two-dimensional system is where C depends on the curvature of the profiles, that is, on the exponents μ and ν.From the original Beverloo expression (Eq.( 1)), the reduced aperture D−kd p was substituted in Eq. ( 7) by the actual dependence of φ and v z on x and R. It should be noted that the particle size is not explicitly included in Eq. ( 7).This implies that the mass flow rate would be independent on the particle size unless there was some dependence through the fitting parameters of Eqs.(3)(4)(5)(6).Janda et al. found the values of them for 1 mm diameter stainless steel beads in their experimental work.In a more recent work, Zhou et al. [13] performed numerical simulations of a two-dimensional flat-bottomed silo with 2 and 6 mm diameter disks.They used the equations proposed by Janda et al. to adjust the data, encountering self-similarity in the profiles.Also, their packing fractions fitted successfully to Eq. ( 6), obtaining parameters that described some dependence on the particle size.However, they could not fit the velocity data to Eq.( 5) and proposed an alternative expression for this magnitude.The aim of our work is therefore trying to generalize Eq. ( 7), checking the validity of the equations (3-6) and analysing whether the fitting parameters change with the particle size.

Experimental Setup
The experimental setup has been designed to be a rescaled version of that used by Janda et al. in [11].Instead of the 1 mm diameter beads, we have used 4 mm diameter spheres of the same type of stainless steel.Therefore, the only variable that has been changed is the beads size, implying a projected area 16 times larger than in the previous case, and a weight that increases 64 times.We have built a two-dimensional flat-bottomed silo from two safety glass sheets of 1600 × 700 mm 2 .Between them, we have stuck two 4 mm-width aluminium flat bars at its lateral edges.They have been supplemented with a small piece of paper in such a way that just a single layer of 4 mm diameter beads fits inside.At the bottom of the silo two 4 mm width movable blade-shaped pieces of stainless steel have been placed to adjust the size of the outlet conveniently.Finally, there is a hopper at the top of the silo that facilitates manual loading.The experiment is protected by means of an aluminium framework to make it stable and to avoid glass bending, ensuring that the grains arrange in a monolayer.
The magnitudes related to the particles flow rate have been measured by analysing a set of films taken by means of a 'Photron FastCam-1024 PCI 100K' high speed camera, which is able to record up to 10 6 frames per second at low resolution.The camera has been placed in front of the silo, and a 30 W led panel light has been placed behind it to maximize the contrast between the beads and the voids.As an example, in Fig. 1 we show a frame of one of the several videos we have taken.These videos have been processed by means of a developed program in matlab, whose code consists in the following: first, the images are binarized and the centroid positions of every particle have been determined at each frame.In this work we are interested only in measuring magnitudes at the position of the exit, so for calculating the velocity of the beads we have focused on the particles that cross the outlet.In practice we find particles whose centroid is over the line of pixels shown in Fig. 1 at one frame and below it at the following.Thus, the velocity of each ball is the result of computing the displacement of the particles divided by the time interval between two frames.Furthermore, the flow rate has been calculated simply by counting the centroids that have passed through that pixel line and dividing it by the total acquisition time.
In order to obtain the packing fraction we also take account of every particle that overlaps the pixel line at each frame.Since we know the size of the beads and the pixel length, the program calculates, from the centroid positions, how much each ball contributes to the outlet pixel line by a simple trigonometric relationship; that is, it obtains the length of the chord of each ball that is onto the outlet line.Then, the outlet is divided in as many cells as is required, and a value of 1 is assigned to each cell if it is occupied and a value of 0 if not.In Fig. 1 it is possi-  3) with ν = 2.7.c) Dependence of the velocity at the center of the exit with the outlet size.The solid line corresponds to the best fit of Eq. ( 5) with γ = 1.45.4) with μ = 5. c) Dependence of the packing fraction at the exit center with the outlet size.The solid line corresponds to Eq.( 6) with the fitting parameters showed inside the graph.ble to observe in yellow the segments which overlap the exit line.The packing fraction profile is finally calculated by time-averaging these amounts over all frames in each position.
In comparison to the work of Janda et al. [11], this way of processing the films has supposed an important enhancement in spatial resolution, as it will be appreciated in the results.

Results
The experimental data of vertical velocity profiles at the outlet line for a range of apertures from R = 0.86 cm to 4.46 cm are shown in Fig. 2.a.It is clearly noted that the curves shape is the same for each hole size.In fact, they have resulted to be self-similar and have been successfully collapsed in the curve governed by Eq. ( 3), just as in Ref. [11] as it is depicted in Fig. 2.b.However, our curve is slightly wider, since we have found an exponent ν = 2.7, a little greater of that of Janda et al. [11] and similar to that obtained in [13] (see Table 1).The velocities at the center as a function of R are represented in Fig. 2.c.The correspondence between the data and the solid line confirms the validity of Eq. ( 5), although in our case the parameter γ has resulted to be 1.45, a bit greater than those reported in [11] or [12].On the other hand, the experimental packing fraction profiles for the same range of apertures are plotted in Fig 3 .a.The increase in resolution enables us to appreciate a small ripple in each curve at a position that coincides with a distance to the border of approximately one particle radius.In the curves corresponding to small outlets, some secondary ripples are observed with a wavelength in the order of the particles size, evidencing that this property is consequence of the discrete nature of our system.Moreover, all the curves exhibit the same shape and they have been collapsed by Eq. ( 4), as it is shown in Fig 3 .b.The exponent obtained in this case has been μ = 5, not much different from the exponents of the Refs.[11] and [13], as it is shown in Table 1.The dependence of the packing fraction at the exit center is represented in Fig. 3.c.Although the experimental data seem to be compatible with the form of Eq.( 6) there are some features that should be considered.First, the asymptotic limit (φ ∞ ) is similar to the ones previously reported, hence appearing to be universal.Second, the parameters obtained (α 1 = 0.38 and α 2 = 3.03 cm) are quite different than in [11], in agreement with a slower rise of the value of φ with R for larger particles.Thus, our value of α 2 does not reproduce the linear behavior with the particle size that is suggested in [13] (see Table 1).In general terms, regarding the comparison among the parameters obtained in different works, it is not easy to establish a direct relationship of the parameters to the particle size.7) with the parameters calculated in the analysis of the velocity and packing fraction profiles as shown in Figs. 2 and 3.
Finally, the experimental results of mass flow rate in function of R are represented by the dots in Fig 4 .The solid line corresponds to Eq.( 7) using the parameters included in the first column of Table 1 and fits quite well to the experimental data.That proves the validity of Eq.( 7) for 4 mm stainless steel spheres.Curiously, although the parameters obtained here and in [11] are rather different, the flow rate values are similar, at least in the range of measurement.

Conclusions
The study of the flow rate in a two dimensional silo with 4 mm stainless steel beads has served to extend the validity of the expressions proposed by Janda et al.However, the parameters found seem to show up a dependence on the particle size which origin is not clear yet.Although for big apertures the density tends to the same asymptotic value in both cases, the velocity scaling seems to depend on the beads size.This fact suggests that for sufficiently large orifices the dependence of the of the mass flow rate on the particle size will become evident.Thus, more work is needed to study this question and to unveil the origin of the issue.

Figure 1 .
Figure 1.Example of a frame of a video taken for an outlet size R = 1.72 cm.The parts of the spheres occupying the exit line are highlighted in yellow.

Figure 2 .
Figure 2. a) Vertical velocity profiles for several outlet sizes.b) Normalized vertical velocities.The solid line is the best fit according to Eq. (3) with ν = 2.7.c) Dependence of the velocity at the center of the exit with the outlet size.The solid line corresponds to the best fit of Eq. (5) with γ = 1.45.

Figure 3 .
Figure 3. Packing fraction profiles for different outlet sizes.b) Normalized packing fraction profiles.The solid line is the best fit according to Eq. (4) with μ = 5. c) Dependence of the packing fraction at the exit center with the outlet size.The solid line corresponds to Eq.(6) with the fitting parameters showed inside the graph.

Figure 4 .
Figure 4. Experimental results of the flow rate.The solid line is the function of Eq.(7) with the parameters calculated in the analysis of the velocity and packing fraction profiles as shown in Figs.2 and 3.