Role of particle size in the kinematic properties of silo flow

We experimentally analyze the effect that particle size has on the mass flow rate of a quasi two-dimensional silo discharged by gravity. In a previous work, Janda et al. [Phys. Rev. Lett. 108, 248001 (2012)] introduced a new expression for the mass flow rate based on a detailed experimental analysis of the flow for 1-mm diameter beads. Here, we aim to extend these results by using particles of larger sizes and a variable that was not explicitly included in the proposed expression. We show that the velocity and density profiles at the outlet are self-similar and scale with the outlet size with the same functionalities as in the case of 1-mm particles. Nevertheless, some discrepancies are evidenced in the values of the fitting parameters. In particular, we observe that larger particles lead to higher velocities and lower packing fractions at the orifice. Intriguingly, both magnitudes seem to compensate giving rise to very similar flow rates. In order to shed light on the origin of this behavior we have computed fields of a solid fraction, velocity, and a kinetic-stress like variable in the region above the orifice.


I. INTRODUCTION
The flow of granular matter in the discharge of silos or hoppers has been studied for many years [1,2] due to its application in industry.In fact, materials of diverse nature, such as cereals, drugs, or minerals, are stored in granular form and have to be processed accordingly.Even so, given the complexity of discrete media, the silo discharge continues to be an open problem from the point of view of its underlying physics and is, indeed, a topic under very active research [3][4][5][6][7].The most widely accepted expression to predict the mass flow rate W of granular matter in a silo was proposed by Beverloo et al. [8] in 1961, where C and k are fitting parameters, g is the gravity acceleration, ρ B is the bulk density, N is the dimensionality of the system, D = 2R is the diameter of the outlet, and d p = 2r p is the beads' diameter.Over the years, this equation has been used in diverse situations.Some examples are vibrated silos [9], a system where the material is driven by a conveyor belt [10], mixtures of granular media [11,12], water submerged grains [13], or even the flow of air bubbles through a two-dimensional (2D) slit [14].However, despite the popularity of the expression of Beverloo et al. [8], it involves some features of uncertain physical sense, such as the inclusion of the bulk density ρ B instead of the flowing density or the reduced aperture size D − kd p , which accounts for a hypothetical forbidden area of the orifice through which the beads are not allowed to pass.This idea, popularly known as the empty annulus [15], has been the traditional way to include the dependence of the mass flow rate on the particle size.
A few years ago, Mankoc et al. [16] proposed a new empirical expression in which the bulk density was substituted by an exponential form of the flowing density as a function of the outlet size, where C and b are fitting parameters.Shortly after, Janda et al. [17] reported an experimental study of a discharge of monodisperse 1-mm diameter stainless steel beads in a twodimensional silo.By means of an image processing technique, they measured vertical velocity v z and solid fraction φ profiles at the outlet for a wide range of apertures.These profiles were proved to be self-similar and could be fitted to these expressions, where ν and μ are fitting exponents whereas φ c and v c represent the value of both magnitudes at the center of the outlet.Note that, as the silo considered was two dimensional, φ(x) accounted for the projection of the three-dimensional (3D) solid fraction onto the XZ plane and R is in this case half of the outlet size.The solid fraction at the center φ c was fitted to an exponential function in R in the same way as in Ref. [16], where α 1 and α 2 are also fitting parameters.In addition, v c was fitted successfully to where γ is a fitting parameter.This √ 2gR scaling of the central velocity has been both experimentally [18] and numerically [19] demonstrated in recent times.Traditionally, it was related to the concept of the free fall arch, which suggests the existence of a hypothetical hemispheric [1] or parabolic [20] region proportional to R, just over the outlet, where the flow properties undergo a transition.Above this region, particles would experiment contact stress and descend with negligible acceleration; below it, grains would fall down freely by gravity with minimal energy dissipation.Despite the soundness of this idea, some evidence has resulted to be in contradiction with it [21,22].Recently, Rubio-Largo et al. [23] performed a micromechanical study that combined numerical simulations and the analysis of experimental data.In that way, they revised the idea of the free fall arch and proposed a modified version based on relating the spacial distribution of kinetic pressure and velocity profiles.Also, they showed that acceleration profiles describe smooth functionalities instead of the discontinuous functions expected.Finally, they justified the parameter γ included in Eq. ( 6) as the integral of the relative acceleration at the center of the orifice along the vertical (z) axis, where F (h) = a(h)/g is the normalized evolution of the acceleration along the vertical axis and h = z/R is the normalized height.
Going back to the work of Janda et al. [17], a new expression for the mass flow rate was derived from the solid fraction and velocity profiles combining Eqs. ( 3)-( 6), where ρ is the density of the material.Then, solving the integral, the expression of the mass flow rate obtained for a two-dimensional system is as follows: where C is a constant which depends on the parameters μ and ν, which determine the width of the corresponding profiles, and γ .We should emphasize that, unlike in the fit of Beverloo et al. [8] [Eq.( 1)], in Eqs. ( 3)-( 9) the particle size is not explicitly included as a variable.However, it is possible that any of the fitting parameters which appear in the expressions above depends on the particle size.Recently, Zhou et al. [24] computed velocity and solid fraction profiles in a numerical study with particles of 2-and 6-mm diameters.Therein, both profiles were fitted successfully by Eqs.(3) and (4).Nevertheless, instead of Eqs. ( 5) and ( 6), they used other functions which included the particle size as a parameter to fit the dependence on the outlet size.In particular, the velocity involved an exponential term in d p /D which is not compatible with Eq. ( 6).
Thus, the main objective of this paper is to determine whether there exists a dependence of either φ or v z (and consequently W ) on the particle size.The first part of this paper consists of reproducing the measurements performed by Janda et al. in Ref. [17] but for spheres of 4-mm diameter instead of the ones of 1-mm diameter used there.The second part of the paper deals with studying the micromechanical behavior of the system.To achieve this, we have applied a coarse-graining tool that, as will be explained in the following sections, enables us to measure velocity and solid fraction fields as well as to calculate a kinetic-stress like magnitude.

II. EXPERIMENTAL SYSTEM A. Setup
A photograph of the experimental setup is shown in Fig. 1.It consists of a quasi two-dimensional silo of 160 × 61 cm 2 built with two glass sheets.Between them there are two aluminum blanks of 4-mm thickness stuck at both sides of the glasses to serve as the silo lateral walls.These blanks are supplemented by two thin strips of paper in such a way that there is only space for a single layer of particles between the glass plates.The particles are monodisperse 4-mm diameter spheres made of stainless steel AISI 420.
At the bottom of the silo there are two movable bladeshaped pieces of stainless steel of 4-mm thickness.In this way, it is possible to change conveniently the size of the outlet through which the beads leave the silo and fall into a bucket.The process of filling the silo is carried out manually through a hopper placed at the top.All of the structure rests on a framework built mainly with aluminum which also provides rigidity to the glass sheets, avoiding that they bend and ensuring that the material remains in a single layer.
The way to proceed in this experiment has been basically to record the material flowing out of the silo for different sizes of the outlet.The recording was performed by means of a "Photron FastCam-1024 PCI 100K" high-speed camera placed in front of the silo, focusing on the outlet and the surrounding area.Behind the silo a light emitting diode light panel illuminates the outlet region producing a great contrast between the beads and the background.This facilitates the subsequent image processing and data analysis.

B. Image recording and processing
The experiments have been carried out in two stages.First, we focused on measuring magnitudes just at the hole in order to reproduce the profiles obtained in Ref. [17].Hence, the recorded region is just the outlet and a little surrounding area of around three particle sizes.We have taken videos for 11 outlet sizes in a range beyond the region where clogging starts to appear.The films, recorded at a frame rate of 1000 fps (frames per second), span durations from around 30 s for the smallest outlet sizes to around 90 s for the largest ones.Then, the images are processed by means of a program developed using MATLAB.Basically, the program core binarizes each frame and detects the centroid of all beads.After that, the contribution of each bead to the solid fraction at the outlet has been calculated in each frame by using a simple trigonometric relationship.In this way, given the centroid position (x c ,z c ) of a bead with a radius r p , the length of the chord that overlaps the outlet line , where r p is the particle radius.The outlet is divided in as many cells as needed, and the program builds an array whose elements are these cells assigning a value of 1 to the cells occupied by a bead and 0 to the ones unoccupied.Thus, in each frame, 1's are placed in the positions of the chords, which are segments with the center in x c and length c.The final solid fraction is achieved by time averaging the array for all frames.
Particle velocities have been calculated from the difference of centroid positions in two consecutive frames which is divided by the time interval between them.To this end, we had previously ensured that the beads' displacement in the time interval between two frames is less than a particle radius.Since we wanted to focus on the velocity profiles just at the orifice, we only have taken account of the particles that were above the outlet line z = 0 at one frame and under it at the next one.Furthermore, the flow rate has been obtained simply by counting the number of beads passing through the outlet and dividing by the total time.
In the second part of the paper we calculated space-time averaged magnitudes [25] in order to investigate the dynamics of the grains in the surrounding area above the orifice.To this end, we have taken complementary films with a larger visual field (see Fig. 1).This is performed for only five orifices but spanning the same size range as before.The frame's height covers nearly 25 bead diameters in such a way that the outlet is placed about two or three diameters above the bottom of the frame.In this case, for each outlet size, we have recorded one single video of around 4 s of duration with a frame rate of 1000 fps.After that, the images are processed obtaining the centroids and velocities of all beads within the field of view.From these data, we implemented a coarse-graining mathematical protocol which was introduced by Goldhirsch [26] and recently studied by Weinhart et al. [27] for its application in silos.This technique, which is explained in detail in the Appendix, makes it possible to calculate continuous fields of density, velocity, and stress from discrete systems.

A. Profiles at the exit
The experimental velocity profiles at the exit line (z = 0) for 11 outlet sizes are presented in Fig. 2(a).These profiles are self-similar in agreement with Ref. [17] as they have collapsed into a single master curve governed by Eq. ( 4) [Fig.2(b)].Remarkably, we have found an exponent (μ = 2.7) slightly greater than observed for 1-mm particles (μ = 2, see Table I).This indicates that our curves are somehow wider than the previously reported circular profiles.Interestingly, μ = 2.7 is similar to the figure given by Zhou et al. [24].Figure 2(c) shows the experimental data of solid fraction profiles at z = 0. First, the rise in resolution with respect to previous works has allowed evidencing a smooth ripple in the shape of the curves, which is more marked as we get close to the edges.This feature takes account of the discreteness of the granular material.Passing over this second order effect, all curves can be collapsed into a single one as shown in Fig. 2(d).The exponent obtained when fitting the data to Eq. ( 3) has been ν = 5, greater than the one reported in Ref. [17] (ν = 4.55) and slightly lower that the one found in Ref. [24] (ν = 5.26).These results evidence the validity of Eqs. ( 3) and (4) to describe the profiles studied, yet the exponents reveal some small dissimilarities.Now, we will check if the scaling parameters (φ c and v c ) also depend on R with the proposed functionalities in Eqs. ( 5) and (6).
The experimental results of v c vs R are represented by black squares in Fig. 3. To establish a comparison, the data corresponding to 1-mm diameter beads [17] have also been included in the chart.Both outcomes have been fitted successfully by Eq. ( 6), nevertheless, the experimental velocities for FIG. 3. Comparison between the experimental data of the velocity at the center of the outlet for 4-mm (this paper) and 1-mm [17] diameter beads.The solid lines are the corresponding fitting functions with Eq. ( 6) obtaining γ = 1.45 and γ = 1.07 for d p = 4 and d p = 1 mm, respectively.Expt.

Expt.
FIG. 4. Comparison between the experimental data of the packing fraction at the center of the outlet for 4-mm (this paper) and 1-mm [17] diameter beads.Both sets of data are fitted to Eq. ( 5) obtaining in each case different values of the parameters α 1 and α 2 but the same φ ∞ as shown in Table I. d p = 4 mm are systematically above the ones for d p = 1 mm.This is captured by the γ parameter which is greater in the case of 4-mm diameter beads (γ = 1.45) as is shown in Table I.
In Fig. 4 the data of the solid fraction at the center of the orifice for 4-and 1-mm diameter particles are compared and fitted according to Eq. ( 5).First, the agreement with this expression has been satisfactory in both cases, and the asymptotical value obtained for large R is the same (φ ∞ = 0.83).However, for 4-mm particles, the growth of φ c with R is rather slow compared to the case of 1-mm beads.Therefore, the parameters obtained when fitting the data to Eq. ( 5) are considerably different as reported in Table I.Unfortunately, it was not possible to establish a straightforward relationship among these parameters and the particle size.

B. Coarse-grained fields
Aiming to further investigate the origin of the different behaviors obtained when changing the size of the particles we have calculated the maps of solid fraction φ, vertical velocity V z , and a kinetic-stress like magnitude σ k (see the Appendix for more details about the coarse-graining technique used to calculate these fields).The set of the solid fraction maps is shown in Fig. 5.At first sight, it can be noted that, mostly for small outlet sizes, there are two triangular regions at both sides of the outlet where the density is considerably higher than in the central part of the silo.This suggests the development of a highly crystallized natural hopper that behaves like a solid.Moreover, for the smallest hole it is possible to discern some yellowish lines crossing these regions which correspond to shear bands.
Regarding the central region of the map, we observe that the material fluidizes as approaching downwards to the outlet line.Also, for big orifices, there are two roughly triangular areas of low solid fractions (blue in the color map) near the edges at the inner part of them.In fact, their sizes seem to be independent of the outlet size so that the material falls more compacted at the center as the outlet gets larger.This could be the cause of the asymptotical value that φ c reaches in the limit of large R [as captured by Eq. ( 5)].For small apertures, however, both areas overlap, and φ c is reduced.
In Fig. 6 we display the set of velocity fields for the same orifices as in Fig. 5. Interestingly, the shape of the pattern obtained above the outlet is not as round as its counterpart for 1-mm diameter beads (reported in Ref. [23]).Also, the fields appear to be self-similar, given that the pattern size grows with R but their shape barely changes.
Finally, Fig. 7 depicts the fields of a magnitude analogous to the kinetic stress σ k .In all cases, it is possible to observe that on their way out, the particles collide with the boundary of the motionless hopper giving rise to a sharp change in σ k from a finite value (in white) to nearly zero (in blue).Apart from that, it can be appreciated that the lobes of maximum σ k close to the edges point upwards and have stretched shapes.
In general, the fields of the different magnitudes reveal that a modification of the beads' sizes induces changes in the whole configuration of the system.In the following we try to establish a connection between these changes in the macroscopic fields and the alterations observed in the velocity and volume fraction profiles reported in Sec.III A.

C. Connecting coarse-grained fields with profiles at the outlet
We will start by analyzing the results of σ k at x = 0 as a function of the height z scaled by the outlet size R (Fig. 8).We observe a certain self-similarity, which indicates that R is the parameter that determines the σ k values.Although the curves are not identical, they have the same shape and seem to collapse in a single one.For the three intermediate outlet sizes the correspondence is noticeable.It is just for the largest and smallest exit sizes for which there are some little deviations.Let us remark that all plots exhibit a maximum at the same values of z/R, except for the smallest outlet where the maximum appears for a higher value.
In Fig. 9 we represent the height z c within the silo at which σ k is maximum in each horizontal position (both distances are rescaled with the outlet size).First, it is observed that most of the values collapse on a single curve when normalizing by the outlet size.The only discrepancy has been observed for the outcomes of the smallest outlet size (R 4r p ) for which the curve is higher than the others and could not be scaled successfully.According to the modified concept of the free fall arch it is possible to connect the region of maximum σ k with the spatial distribution of the velocity profiles shown in Fig. 2(a).This is performed by means of z c ∼ [1 − (x/R) 2 ] 2/μ , which would be the shape of the region from which the particles would initiate a free fall to produce a velocity profile governed by Eq. ( 4).Note that, for 1-mm particles we obtained μ = 2, a scenario compatible with z c displaying a parabolic arch giving rise to a circular velocity profile.For 4-mm particles, the arch shape emerging when using μ = 2.7 [as obtained when fitting the profiles of Fig. 2(b)] agrees quite well with the experimental data as shown in Fig. 9.
Finally, the acceleration of the material just above the center of the orifice has been calculated from the velocity at x = 0 (V zc ) by means of the expression, The results of this variable with respect to the height are presented in Fig. 10.Again, the acceleration curves have been FIG.6. Vertical velocity V z maps in cm/s for different outlet sizes as indicated above each chart.Note that the color bar is the same for all plots.FIG. 7. Maps of σ k (in arbitrary units) for different outlet sizes as indicated above each chart.Note that the color bar is the same for all plots.
resulted to be self-similar when normalizing by R. Incidentally, it seems that the material experiences an acceleration slightly larger than gravity for some hole sizes.This fact was not seen in the 2D experimental work with 1-mm particles, yet similar behavior was reported for 3D simulations [23].More importantly, the acceleration profiles provide an alternative way of linking the velocities at the outlet with the behavior above it.As explained in Sec.I this is performed by calculating the γ parameter as the area under the curves of Fig. 10 [Eq.( 7)].We should note that the integral goes from 0 to ∞ and our range of measurement is limited by the visual field of the camera, therefore, the value of γ obtained by this method will always be below the real one.This problem becomes specially relevant for big outlets where the integral is far from convergence.For this reason, the results of R 16r p and R 23r p have not been taken into account.For the other smaller outlets we have obtained values of γ = 1.34, γ = 1.35, and γ = 1.08 for R 12r p , R 8r p , and R 4r p , respectively.Indeed, the two first values of γ are remarkably greater than 1 and compatible with the one achieved in the fitting of Fig. 3 (γ = 1.45).However, for the case of R = 4r p the system seems to behave differently as revealed by the trends of all the magnitudes analyzed.In fact, if we look in detail at the two black squares of lower R in Fig. 3, we note that they are slightly under the proposed fit, in agreement with the smaller value of γ obtained (γ = 1.08).This feature could be attributed to the fact that, for such low values of R, the flux is no longer stationary and intermittencies start to appear [28][29][30].

IV. DISCUSSION
In this paper we have presented a detailed analysis of the flow of grains through orifices driven by gravity.In particular, we have tested the validity of previously proposed expressions by Janda et al. [17] when modifying the particle size.Despite the fact that the obtained data can be represented by the same functionalities, we observe a discrepancy in the dependence of both v z and φ on x and R, captured by the different fitting parameters encountered for Eqs. ( 3)- (6).In order to explain the origin of these differences, the fields of the packing fraction, vertical velocity, and σ k have been calculated above the outlet.The most characteristic feature is the formation of a natural hopper of static grains at both sides of the orifice, which is suggested to produce an alteration in the profiles.In fact, in a recent a study [4], similar density variations in the stagnant zone have already been found when changing the shape of the particles.Furthermore, some of the changes in the profiles' parameters of 4-mm particles are explained from the macroscopic fields above the outlet.In particular, the higher value of μ encountered is compatible with a slight modification of the region of maximum σ k above the outlet.Similarly, the higher γ obtained has been related to an alteration of the acceleration evolution along the vertical axis at x = 0. Surprisingly, despite the differences in the bulk behavior and the profiles of density and velocity at the outlet, the experimental mass flow rates (normalized with the thickness of the silo [31] to account for the 2D projection of the problem) resulted to be almost identical for beads of different sizes.Both sets of flow rate data (for d p = 1 and d p = 4 mm) are depicted as functions of R in Fig. 11.Also, we represent the calculated flow rate by means of the product of the velocity and packing FIG.11.Experimental mass flow rates normalized to the bead size with respect to the outlet size for spheres of 4 mm (this paper) and 1 mm [17].The continuous lines correspond to the expressions obtained with Eq. ( 9) using the parameters obtained by Eqs. ( 5) and ( 6), indicated in Table I. fraction fittings according to Eq. ( 8), which, as expected, agree very well with the experimental values.
As stated, within the experimental measurement range, the flow rates of both cases are approximately similar.The reason is that, whereas the velocity for the 4-mm diameter particles is larger, the solid fraction is smaller, so they seem to compensate.Nevertheless, as the outlet size becomes larger, the extrapolation of both curves results in a small divergence above the experimental data range.This is due to the fact that both systems tend to have similar solid fractions for large orifice sizes.Hence, the deviation in the limit velocity and, consequently, the difference in the values of γ would lead to the separation of flow rates for large outlet sizes.Although more investigations are needed to confirm this result, we believe that the natural hopper developed for the 4-mm particles might be behind the higher velocities observed which, for sufficiently large R, will lead to higher flow rates.Even though this hypothesis should be confirmed, it is sustained by findings using two-dimensional hoppers where the flow rate has been shown to be an increasing function of the bottom angle with respect to the horizontal axis [2,32]; hence, the appearance of the natural hopper could generate similar effects.
Finally, let us note that it is difficult to definitively state what is the specific property of the particles responsible for the macroscopic changes observed in the system.By changing the particles' sizes we also have altered their masses, so it is hoped that any of these variables is behind the appearance of the natural hopper and the differences in the kinematic properties.Nevertheless, we should not discard other features which are related, directly or not, with the beads' masses and dimensions.Some examples could be the nature of the contacts between the beads, the friction, or the difference in kinetic energy per unit area.Anyway, the relationship among all these magnitudes is not trivial, and further research is necessary to clarify these questions.

APPENDIX: COARSE-GRAINING CALCULATIONS
Macroscopic continuum fields are obtained by means of the application of an integrable coarse-graining function .
In particular, we have chosen a two-dimensional Gaussian function with the form where ω, the Gaussian rms width, has been selected to be equal to the particle radius r p in all cases.This election is a compromise between an exceedingly small value (which would not have the desired smoothing effect) and a too large one (which will require major corrections [33]).Thus, the two-dimensional solid fraction field φ of the system can be

FIG. 1 .
FIG. 1.A photograph of the experimental setup.The red parallelogram indicates approximately the area where films are recorded to calculate the fields.An image from the film corresponding to an outlet size of R 16r p is shown on the right.The origin of the Cartesian coordinate system has been set at the center of the exit orifice.

FIG. 2 .
FIG.2.Experimental profiles of different magnitudes at the exit line (z = 0) for several outlet sizes as indicated in the legend.(a) Vertical velocity, (b) normalized vertical velocity where the solid line is the best fit according to Eq. (4) with μ = 2.7, and (c) solid fraction.(d) A normalized solid fraction where the continuous line is the best fit according to Eq. (3) with ν = 5.

FIG. 5 .
FIG.5.Solid fraction φ maps for different outlet sizes as indicated above each chart.Note that the color bar is the same for all plots.

FIG. 8 .
FIG. 8. Values of σ k above the center of the orifice (x = 0) vs the height (z) in units of outlet size (R).The results for different outlet sizes are reported as indicated in the legend.

FIG. 9 .
FIG.9.Heights where σ k is maximum for each horizontal position x (both rescaled by the outlet size R).Results are obtained from the data of Fig.7for different outlet sizes as indicated in the legend.The solid line is a function compatible with the modified concept of the free fall arch (using μ = 2.7) as explained in the text.

TABLE I .
[17]arison between the fitting parameters obtained in this paper for particles of d p = 4 mm and in the work of Janda et al.[17]for d p = 1 mm.