A differential equation for the flow rate during silo discharge: Beyond the Beverloo rule

We present a differential equation for the flow rate of granular materials during the discharge of a silo. This is based in the energy balance of the variable mass system in contrast with the traditional derivations based on heuristic postulates such as the free fall arch. We show that this new equation is consistent with the well known Beverloo rule, providing an independent estimate for the universal Beverloo prefactor. We also find an analytic expression for the pressure under discharging conditions.


Introduction
The discharge of grains though an opening at the base of a silo has been considered in a number of studied since the 19 th century (see for example [1][2][3][4][5][6][7][8][9] and references therein). The most salient feature is the fact that the flow rate does not depend on the height h of the material in the container, in clear contrast with the behavior of viscous fluids. This has been sometimes attributed to the pressure saturation observed in static silos, however, discharging silos have a continuously evolving pressure [10].
If the discharge orifice is circular and large enough to avoid clogging [11], the mass flow rate Q is described by the so called Beverloo rule [1,12] where D o is the diameter of the opening, ρ b the bulk density of the granular sample, g the acceleration of gravity and d the diameter of the grains. Here, k and C are two fitting dimensionless constants. Interestingly, while k may vary up to a factor of 2, depending on the grains used, C ≈ 0.58 for virtually any material tested [12]. However, some deviations are observable for very low friction materials [13]. The Beverloo rule is generally explained based on heuristic models such as the "free fall arch" model and the concept of "empty annulus" [12]. However, these two models have been recently challenged [14,15]. Staron et al have also shown that the Beverloo rule can be obtained if the Navier-Stokes equations are solved for a plastic fluid including a constitutive equation for the effective friction based on the µ(I)-rheology [16].
In this work, we use the global energy balance for the granular material inside a discharging silo and a constitutive relation to derive a simple differential equation for the mass M(t) in the silo. We show that for a simple discharge the equation is consistent with the Beverloo rule, providing an independent estimate for the prefactor without fitting data. In the process of validating the differential equation, we also find a functional form for the pressure in the silo in the dynamic regime which is rather different to the well known Janssen equation.

Energy balance
Consider a cylindrical silo of radius R s , discharging through an opening of radius R o (see Fig. 1). The energy balance requires that, at any time, where W g is the power injected by the force of gravity acting on the grains,K in is the rate of change of the kinetic energy of the grains inside the silo, W out is the power loss due to the grains that leave the silo at a velocity v out ,Ė is the rate of change of the elastic energy of the grains, and W D is the dissipated power due to the non-conservative interactions between grains and between grains and walls. During the discharge, the mass M(t) in the silo can be written as where m is the mass of one grain, N(t) is the number of grains in the silo, ρ b is the apparent density in the bulk, A s = πR 2 s is the cross section of the silo, z(t) is the head of material, and z cm (t) = z(t)/2 is the center of mass of the granular column. We have assumed that the density is homogeneous throughout the column. Therefore, the flow rate q(t) in particles per unit time is where v cm (t) is the velocity of the center of mass of the granular column.

Internal kinetic energy (K in )
The kinetic energy of the grains inside the silo is where the sum includes all particles in the silo at time t and v i (t) is the velocity of particle i. In terms of the center of mass K in (t) can be expressed as where the term [u i (t) − u cm (t)] 2 can be neglected according to results obtained from DEM simulations (data not shown), hence Then, the rate of change of K in iṡ which can be written, using Eqs. (3) and (4), aṡ

Gravitational energy (W g )
The gravitational potential energy of the particles inside the silo using Eq. (3) is The power injected is therefore

Discharge energy (W out )
The kinetic energy that is removed from the system due to the particles that exit through the opening is where v out is the velocity of the grains that exit the system and N 0 is the initial number of grains in the silo. Hence, the power removed from the system is If we consider that the mass flow rate is mq(t) = ρ o A o v out , with ρ o the apparent density at the opening and A o the cross section area of the opening, then

Dissipated energy (W D )
It has been shown that for a shear cell of thickness L the tangential stress τ necessary to develop a flow at velocity v of a granular sample can be put in terms of the inertial number I as [17] τ = µ(I)P, where P is the confining pressure, µ(I) is the effective friction coefficient, is the inertial number that characterizes the flow if the gains are stiff and L d, and ρ is the density of the material the grains are made of.
If we assume that the flow in the silo can be represented as a shear flow, where v is the velocity of the free surface, L is the silo radius and P is the mean hydrostatic pressure in the column (i.e., 1/3 of the trace of the stress tensor), then, the power dissipated by friction is where A(t) is the area of frictional contact between the grains and the silo. This area includes the lateral walls, which is connected to the height of the column at time t, plus a fixed area that should be proportional to the base of the silo, i.e., A(t) = 2πR s z(t) + απR 2 s = 2πR s M(t) ρ b A s + αA s . Then, using Eq. (4) and v(t) = 2v cm (t), we obtain

Differential equation
If the grains are stiff, the variation in the elastic energy are expected to be small and we can neglectĖ. Notice also thatK in [Eq. (9)] falls quadratically with A s in comparison with the remaining terms of Eq. (2) that fall linearly with A s or are independent of A s . Therefore, for a wide siloK in vanishes. Plugging in Eqs. (11), (14), and (17) into Eq.
(2) and considering Eqs. (3) and (4), we obtain Solving forṀ and using A o = πD 2 o /4, being D o = 2R o , the diameter of the openinġ We note that Eq. (19) is similar to the Beverloo rule Eq. (1). The exponent 5/2 is not apparent since D o is squared in Eq. (19). We will come back to this point below. Before, it is important to mention that Eq. (19) includes in the prefactor ρ o instead of ρ b . However, if we approximate roughly ρ o ≈ ρ b /2, then the proportionality constant becomes C = π √ 2 8 ≈ 0.56. This has to be compared with the value obtained by fitting the Beverloo experiments that led to C = 0.58, which has been shown to be remarkably insensitive to the material properties [12].
Equation (19) is a first order differential equation for M(t) that can be closed with an initial condition such as M(t = 0) = M 0 . To solve this equation it is necessary to know µ(I), P(t) = P(M(t)) and α. In the rest of this paper, instead of directly solving Eq. (19), we estimate µ(I) and then extract an expression for P(M(t)) based on the condition that the flow rate must be constant throughout the discharge, as observed experimentally. This allows to fit the value of α using pressure data from DEM simulations of the discharge.
To estimate µ(I), consider that if R s → ∞, then I → 0. Da Cruz et al. have shown that in the quasistatic limit (i.e., I → 0) µ(I = 0) ≈ 0.26 for all material properties if the grain-grain friction coefficient is above 0.4 [17]. We have carried out DEM simulations using YADE [18], for a silo where R s = 15 mm and D o = 15 mm, and a granular sample with d = 1 mm, ρ = 2000 kg/m 3 , ρ b = 1160 kg/m 3 and the friction coefficient of the grain-grain and grainwall interaction is 0.55. The inertial number estimated as varies during the discharge in the range 9.8 × 10 −3 < I < 1.6 × 10 −2 , which is in the quasistatic limit consistent with an effective friction coefficient 0.26 [17].
As we mentioned, Eq. (19) does not show the classical 5/2 exponent for D o , nor the correction −kd for the "empty annulus". However, the expression under the square root in Eq. (19) must compensate the difference between D 2 o and the traditional Beverloo factor (D o − kd) 5/2 . During the discharge it is known that Q is constant. Hence, the radicand in Eq. (19) must be also a constant γ. Therefore, the pressure can be written as .   Figure 2 shows a fitting of Eq. (20) to the aforementioned DEM data of the silo discharge. The pressure is calculated as the average of Tr(σ)/3 over the entire volume of the granular column a any given time during the discharge. Here, σ is the stress tensor. We took, following the discussion above, µ(I) = 0.26 and γ = 9.2 × 10 −3 . As we can see, the only fitting parameter is α whose best fit value is α = 3.6 ± 0.2. The fit is only applied to a section of the process since the assumption of a constant flow rate is not justified in the final stage of the discharge. The fit is qualitatively fair. In particular, Eq. (19) has a saturation of the pressure that is much slower than the traditional Janssen law for static silos. This is consistent with the results found here for the simulations and also results reported in experiments of discharging silos [19].
We speculate that the value of α found may be "universal" in the sense that it does not depend significantly on D o , R s , and material properties. However, this remains a subject of investigation.

Conclusions
We have shown that an energy balance, coupled with a constitutive equation for the rheology of a granular sample leads to a simple differential equation for the flow rate of a discharging silo. The equation resembles the Beverloo rule and yields a proportionality constant C = π √ 2 8 ≈ 0.56 which is very close to the experimental fitted value of 0.58.
The assumption of a constant flow rate leads to an expression for the pressure in the silo that is more suitable for these dynamic conditions than the Janssen law derived for static silos. We have to bear in mind, however, that the pressure in the silo is not homogeneous and a more elaborated analysis should include the local pressure in the de-scription rather than the global pressure as we have done in this work.
The derived differential equation allows to explain the known phenomenology of silo discharge without the need of heuristic postulates such as the "free fall" arch and the "empty annulus". Moreover, this equation also opens the way to consider more complex situations in silo discharge, such as forced flow or vibrated discharges, that have not been explored from a theoretical perspective.