Critical numerical analysis of quasi-two-dimensional silo-hopper discharging

We present a critical comparative analysis between numerical and experimental results of quasi-two-dimensional silo and hopper flows. In our approach, the Discrete Element Method was employed to describe a single-layer mono-disperse sphere confined by two parallel walls with an orifice at the bottom. As a first step, we examined the discharge process, varying the size of the outlet and the hopper angle. Next, we set the simulation parameters fitting the experimental flow rate values obtained experimentally. Remarkably, the numerical model captured the slight non-monotonic dependence of the flow rate with the hopper angle, which was detected experimentally. Additionally, we analyzed the vertical velocity and solid fractions profiles at the outlet numerically and experimentally. Although numerical results also agreed with the experimental observations, a slight deviation appeared systematically between both approaches. Finally, we explored the impact of the system’s confinement on this process, examining the consequences of particle-particle and particle-wall friction on the system macroscopic response. We mainly found that the degree of confinement and particle-wall friction have a relevant impact on the outflow dynamics. Our analysis demonstrated that the naive 2D approximation of this 3D flow process fails to describe it accurately.


Introduction
Filling, storing, and discharging of granular media are conventional industrial procedures, which are highly automated. These processes are designed empirically and ultimately need manual intervention: any failure in these devices could eventually lead to devastating consequences [1][2][3]. Unfortunately, there are no well-founded theoretical frameworks to describe these macroscopic phenomena in terms of the particle's micromechanical properties. Moreover, particles and grains are diverse; thus, their properties may differ depending on the specific context. Hence, the mechanical response of granular materials is a topic in continuous development.
When examining granular flows, there are several practical restrictions, and sometimes, it is not possible to fully describe the 3D behavior of the grains. Quasi-two-dimensional experimental setups are often developed [4][5][6][7][8][9][10][11][12], because they capture the particle dynamics precisely. These setups are generally very flexible, and they typically consist of two transparent plates between which opaque particles are poured. Afterward, the use of non-invasive optical methods allows measuring particle motion via image processing. However, quantifying the particle-particle interaction forces is a more challenging task, which is only achieved in static conditions [13,14]. There is a real need to perform numerical simulations in that context because they provide the local particle-particle interaction and the force-chain network structure. Discrete element modeling (DEM) is widely accepted as an effective method for addressing engineering problems concerning dense granular media [15][16][17]. In DEM, the particle-particle contact laws define the macroscopic collective results with high accuracy. Moreover, it is relatively easy to implement, and it has the advantage of being parallelizable.
Particle discharging of a silo is a canonical example of granular flow. Over the centuries, experimental investigations have been carried out to explore the silo flow. Very early on, it was shown that the discharge flow rate Q of a silo is practically independent of the filling height when the orifice diameter size is large enough. Although it is well accepted that the value of Q in a silo varies as D 5∕2 (or D 3∕2 in 2D), a theoretical framework explaining this robust finding is lack. It is worth mentioning, a few groups have recently progressed in understanding the origin of this scaling by combining numerical and experimental approaches [18][19][20][21][22][23][24].
Very recently, the flow of mono-dispersed stainless-steel spheres in a quasi-two-dimensional hopper was experimentally analyzed, varying the hopper angle and the size of the outlet [25]. Here, we performed a critical numerical analysis of the particle discharge process, validating numerical algorithm guided by the experimental outcomes [18,25]. As well as this, we explored the significance of the system's confinement, examining the impact the particle-particle and particle-wall friction have on the system macroscopic response. More importantly, we analyzed the relevance of considering this "hypothetical 2D system" as a real 3D single layer of particles. The work is organized as followed: in Sect. 1.1, we summarize a theoretical framework concerning the particle flow rate in silos and hoppers, which are used further in the discussion. In Sect. 2, the simulation method is described briefly. Next, Sect. 3 summarizes the main outcomes, discussing the simulation method's feasibility to reproduce the experimentally obtained results.
Decades ago, Beverloo et al. [26] introduced a phenomenological formulation quantifying the mass flow rate Q in silos and hoppers in terms of the outlet size D = 2R, where B = M B , and M is the material density of the grains, B is the solid-fraction, while C and k are fitting constants. The factor (D − kd) is the hypothetical effective outlet size when considering grains of dimension d and the existence of an empty annulus. Moreover, Beverloo's approach is based on the assumption that the velocity of the grains at the exit scales with the outlet size as √ D . Note that Eq. (1) corresponds to the 3D case, while for the 2D case the exponent is 3/2.

Flow rate in hoppers
More recently, a new formulation of Q(R) was introduced for flat bottom silos, in terms of the outlet radius [18,19]. This formulation assumes that the vertical velocity v z (r) and solid-fraction (r) are self-similar functions of the horizontal radial coordinate r/R. Importantly, both magnitudes are self-similar functions at the silo exit, that can be scaled by their values at the orifice center. Such an idea was recently generalized for 2D-hoppers by Darias et al. [25]. Accordingly, solid-fraction and vertical velocity profiles at the hopper exit can be expressed as: where v c and c are the velocity and solid-fraction at the center of the orifice, respectively. The exponents a and b control the curvature of the dome shape profiles, whereas the two R dependent coefficients read as, Equation 4 accounts for the dependency of the material dilatancy on R, assuming an exponential saturation set by two fitting parameters, 1 and 2 , that controls its magnitude in the region where the system can clog. More importantly, the expression captures that the discharged material passing through the exit is loose for very large orifices. Accordingly, ∞ is the solid-fraction in the limit of very large orifices R → ∞ [18] which is typically lower than the bulk solidfraction. The effective acceleration along the vertical axis, accounted by , determines the magnitude of v c , which is obviously influenced by the hopper angle. Finally, as shown in [25] the hopper also has a relevant impact on the shape of the solid-fraction and velocity profiles accounted by the parameters. Altogether, this set of parameters fix the mass flow rate Q in quasi-2D hoppers, resulting in: where accounts for the front-back wall separation. The numerical factor is a consequence of considering a single layer of spheres as a flowing material. -function results from the shape of the velocity and density profiles, while accounts for the influence of the internal stresses. We emphasize that Eq. (6) was derived only assuming the selfsimilar properties of velocity and density profiles, and the corresponding scaling functions. Consequently, any numerical approach must reproduce these distinctive features as we discuss in the next sections.

Numerical model
In our simulations, we used discrete element modeling (DEM), describing the evolution of an ensemble of contacting particles. The model implementation is a 3D-version of a previous work reported in [27]. For the particle-particle interaction, it adopts a linear-viscoelastic force model, which is based on the Cundall-Strack approximation [28,29]. The equations of motion are integrated via a velocity-Verlet algorithm with a predictor-corrector modification [27,30]. A complete description of this theoretical framework can be found in [15,27,30].
Here, our numerical procedure resembles the experimental setup reported in Ref. [25]. Thus, the system consisted of a single layer of spherical beads of diameter d = 0.1 cm , confined by two walls at a distance , which we will from here on refer to as the confinement width. It had a width of W = 40 cm and silo-hopper outlet size D = 2R was configured. Specifically, three different tilt hopper angles, = 0 • , 30 • and 60 • were implemented. In all cases, the silo was high enough to contain N = 1.2 × 10 5 particles.
To mimic steel spheres we set density = 7.52 g/cm 3 , particle stiffness n = 1.94 ⋅ 10 6 g/s 2 , restitution coefficient e n = 0.9 , and friction coefficient = 0.25 , which is a very good estimation for a steel-steel surface [31]. The ratios between normal and tangential visco-elastic constants were taken as s n = 2 7 and n s = 3 , with n = √ 2mk n 1+( ∕ln(e n ) 2 ) = 4.15 g/s. The walls were rigid and obeyed the same rules of interaction as the spheres. Frontal and rear planes had a density p = 2.5 cm 3 (mimic glass walls) and a friction coefficient w , to be adjusted.
The protocol initialization places the beads within the silo in a hexagonal arrange with an inter-particle distance slightly larger than the particle diameter. Then, a randomly oriented initial velocity of 10 cm/s is assigned to each bead. Then the system evolves through its initial stationary configuration. Finally, the silo exit outlet is opened, and the particles start to fall by gravity towards the outlet. After crossing the outlet, they travel a distance of 20d and are removed from the simulation. The results were visualized using the open-source tool OVITO [32]. As example, Fig. 1 illustrates a silo-hopper in discharge after 2.5 s with approximately 1.2 × 10 5 spheres, = 30 • and D = 2 cm.
In all cases explored here, the systems quickly evolve to a steady-state characterized by a stable macroscopic flow rate. The particle flow rate through the outlet was computed by taking the elapsed time Δt i , in which ΔN = 80 spheres pass through the outlet, namely, Q i = ΔN Δt i . Then, the mean flow rate Q was computed, averaging 64 Q i consecutive samples that imply less than one second simulation time. The standard deviations found in the complete range of parameters explored were always under 5% of its mean value.
The consistency of the numerical findings was also examined, computing the spatial profiles of the vertical-velocity v z (x) and particle solid-fraction (x) at the orifice location.
The v z (x) profiles were sampled by taking the particle vertical-velocity just after it crosses the exit line, at a given horizontal position x. The solid-fraction profiles (x) were obtained as the time average of occupation histograms at z = 0 . The exit line was divided into 2048 bins, where each bin is 1 if there is mass and 0 if it is an empty space. It is worth mentioning that the protocol implemented to determine these magnitudes is the same as the one used in the corresponding experimental setup.

Results
As mentioned above, we examined the discharging of silos and hoppers, resembling the experimental conditions of Ref. [18,25]. After a careful analysis of the numerical and experimental flow rate results [18,25], the confinement width = 1.10 d and friction coefficient w = 0.23 were chosen as the best set of parameters to reproduce experimental conditions. In a later section, we demonstrate that these parameters have the most significant impact on numerical outcomes. Figure 2a illustrates the correlation between time-averaged flow rate Q obtained numerically and the experimental observations (data taken from [25]). The graph includes outflow data corresponding to the entire range of explored outlet sizes R and two hopper angles ( = 0 • and = 30 • ). As can be noticed, in Fig. 2, the agreement between both (a) A system view notes that crystallization domains coexist with grain boundaries defects due to the particle monodispersity, as occurs in the experiment. (b) Zoom view around the outlet during the discharge. (c) A lateral view; note the spheres can slightly move between the frontal and rear planes, in y-direction due to the wall separation approaches is excellent, resulting in Pearson correlation coefficients larger than 0.99 in both cases.
The analysis of Fig. 2b supports, even more, the validity of the numerical approach. It shows the numerical mass flow rate data Q(R) vs. the outlet sizes R, for hoppers with = 0 • , 30 • , and 60 • . The agreement with the mass-flow expression, Eq. (6), introduced recently by Darias et al. [25] is excellent. Note that, in each case, the used parameters a, b, and are the ones reported in Ref. [25], which correspond to experimental data. Significantly, the numerical outcomes also recover a slightly non-monotonic dependence of Q when varying the tilt angle [33,34]. It is worth mentioning that this trend has been very recently found in 3D hoppers, both experimentally and numerically [35].
Given the model's feasibility to reproduce the macroscopic system response, it is crucial to verify that the kinematic fields that determine the resulting flow are also compatible with those reported in the experiments. Figure 3 shows the velocity (vertical and horizontal) and solid fraction values, obtained for the simulations and the experiments (taken from [25]) for a flat bottom silo and R = 2 cm. Insets in Fig. 3a and b illustrate the data cloud of the velocity components and their corresponding coarse-grained profiles v z (x) , v x (x) , which are comparable with the experimental ones introduced in Ref. [25]. The main panels show a comparison between numerical and experimental coarse-grained profiles and the outcomes of Eqs. (2) and (3). In general, we found a reasonable agreement between the numerical results and the results introduced in [18]. However, a systematic deviation seems to exist between numerical results and the solution given by the set of parameters introduced to fit the experimental results.
To explore more carefully such a difference, we also analyzed the outcomes corresponding to the full range of explored outlet sizes in the cases of = 30 • (Fig. 4a ) and (b)). The insets shows a sample of the numerical magnitudes. Red circles are obtained by the the same smoothing protocol introduced to present the experimental results (blue squares) in [18]. The agreement between experimental and numerical mean field is remarkable. (c) 2D-solid fraction profiles just at the outlet. Again, the agreement between numerical (red dots) and experimental results (blue squares) is excellent. Continuous lines correspond to the fitting functions discussed in [18,25], see text for details and b) and = 60 • (Fig. 4c and d). The vertical velocity values are shown in Fig. 4a and c, and solid fraction values Fig. 4b and d. For comparison, the graphs at the left illustrates the numerical results, while the right ones their corresponding experimental values (taken from [25]). It is noticeable that minor discrepancies systematically appeared. Whereas numerical velocity profiles slightly overestimate the experimental data, the solid-fraction displays an inverse trend. Moreover, the smaller the size of the outlet, the more significant the differences between numerical and experimental solid fractions. Given the excellent agreement when comparing the flow rate values in all cases, we decided to explore the origin of the mean value profile discrepancies carefully.
With this aim, our analysis focuses on the maximum velocity v c (R) , solid fraction c (R) , and the characteristic particle flux J c (R) = c (R) ⋅ v c (R) detected at the center of the orifice, see Fig. 5a-c respectively. For comparison, the experimental findings reported in Ref. [25] are also included in the figures. In each case, the data values were fitted using Eqs. (4) and (5), respectively. In performing this analysis, we used the = 30 • case, where the differences are noticeable. For v c (R) , we obtained that the functional dependence Eq. (5) is in excellent agreement with the numerical data, resulting in num = 1.18 , whereas exp = 0.98 was the reported experimental value [25]. More importantly, the numerical solid-fraction reproduced the significant dilatancy of the granular flow, even at R → ∞ (numerical num ∞ = 0.77 , experimental exp ∞ = 0.84 ). The latter indicates conclusively that the phenomenological formulations, Eqs. (4) and (5), are good quantitative descriptions of the micro-mechanical interactions, which is behind the observed macroscopic dynamic. However, we obtained that the numerical values of v c and c at the orifice center differed systematically from the experimental result. In fact, the experimental results of v c were systematically overestimated by the numerical tool (see Fig.5a), while the experimental c values appeared to be underestimated (see Fig.5b). The reasons behind the discrepancies between both approaches are not obvious. Firstly, we could argue against the accuracy of the used contact model as the major discrepancy source. Secondly, we could also mention possible experimental bias due to the experimental resolutions. Moreover, in the experiments, the gap between walls at the silo exit was not accurately reported.
Importantly, when examining the characteristic current density J c (R) vs. R curves (see Fig.5c) we recovered the consistency between both approaches, in line with the macroscopic results regarding the flow rate. The two data sets J c (R) show a significant correlation (see inset of Fig.5c). These outcomes are very consistent with the behavior obtained for the macroscopic particle flow rate, which was presented earlier (see Fig. 2).
Developing the previous analysis, we have chosen the = 30 • case to emphasize the differences between experimental observations and numerical results. The data analysis of v c (R) and c (R) corresponding to = 60 • showed very similar trends. Nevertheless, despite the similarities, it is worth mentioning that the deviation of the numerical results seems to be more pronounced for = 30 • as occurs in another alternative approach [35]. This feature could be related to the interplay between collisional and sliding regimes against the hopper wall, but this will be discussed elsewhere.

Particle velocity fluctuations
Recently, Rubio-Largo et al. [19] correlated the discharge dynamics with the stresses inside a flat bottom silo. These authors showed that kinetic stress displayed comparable qualitative features regarding the free fall arch idea [36]. They suggested that kinetic pressure plays a central role in controlling the silo discharge rate. The magnitude of kinetic pressure mainly depends on the particle velocity fluctuations relative to the coarse-grained velocity fields. Consequently, not only mean-field outcomes should be calibrated to guarantee the predictive strength of the numerical approach. When doing predictive analysis using DEM numerical tools, we should also ensure it reproduces the particle velocity fluctuations with sufficient accuracy. Hence, the quantification of the velocity fluctuations is needed to extrapolate the numerical predictions to the actual experimental situation. Figure 6 compares the magnitude of the particle velocity fluctuations with respect to the mean-field given by Eq. 3 for the experimental and numerical data corresponding to R = 1.95 (see continuous lines in Fig. 4a). The figure illustrates the components of the fluctuations of the vertical, v z , and horizontal, v x , components. Firstly, it should be highlighted that the magnitude of the fluctuations effectively covered the same ranges. Besides this, fluctuations in the orthogonal directions are uncorrelated, and more importantly, numerical and experimental qualitative features are indistinguishable (see Fig. 6a ). Both experimental and numerical results fluctuated against the fitted profiles. They passed the Shapiro-Wilk normality test at a p=0.05 level of significance; nevertheless, although both cases have nearly equal standard deviations, like the data set shown in Fig. 6b, a few of the samples studied do not pass equally of variances tests. More importantly, the numerical protocol captures the spatial distribution of the fluctuation at the outlet region. Figure 6c shows that the velocity fluctuations become more prominently at borders of the outlet, and the numerical protocol also describes such tendencies with accuracy. Altogether, these results stress the importance of choosing the adequate averaging protocol to construct kinematic fields like the discussed in Ref. [19].

Significance of the degree of confinement and wall friction
For completeness, we also thoroughly analyzed the significance of the system confinement and the particle-wall friction w . To this end, we performed a systematic numerical analysis, varying both the separation between the walls and w . The data has been analyzed using the reported experimental data of Ref. [18] as a baseline reference.
For the sake of simplicity, we examine a flat silo case ( = 0 • and R = 2.05 cm ), systematically changing = [1.05, 1.10, 1.15, 1.20 and 1.30] in terms of the particle diameter d, and w from 0.0 to 0.8 in steps of 0.2. Figure 7 shows the particle flow rates Q obtained in all cases. First, it is evidenced that variations of the confinement are not very relevant in comparison to changes induced by the wall friction w , which affect the outflow significantly. However, interesting features also arise when varying the , keeping w constant. As expected, for friction-less particles, the flow rate slightly increases when increasing . This is simply because the effective size of the orifice increases with increasing . However, the trend is inverted for systems with moderate wall friction w , i.e where w is comparable with the particle-particle friction = 0.2 . We remark that this system is the one that better reproduces the experimental results [25]. In that case, we obtain that Q slightly decreases when increasing the confinement . Interestingly, the trend is inverted again when exploring significant wall friction w ≥ 0.6 . This fact indicates that the degree of confinement width and the wall friction play a non-trivial role in defining the particle outflow.
A well-founded explanation of this distinct behavior can be extracted from a deeper analysis of the numerical data in the region of the orifice. Figure 8 illustrates the solid fractions (x∕R) and vertical velocity v z (x∕R) profiles, obtained changing and w . Firstly, the outcomes show that when using friction-less walls (see Fig. 8a and d), the confinement width has no significant impact on the outflow dynamics. The solid fraction at the orifice practically does not change, and the particles move slightly faster when rising . As a result, the particle flow rate has a slightly increasing trend with . Including the wall-friction force leads to a significant reduction of the solid fraction and particle velocity in the orifice region. When w is comparable to the particle friction, both (x∕R) and v z (x∕R) decreases with increasing (see Fig. 8b and e). As expected, it correlates with the decreasing trend of the particle flow rate (see Fig.7). However, when exploring the system response for large w values (see Fig. 8c and f), we detected that the particle-wall collisions effectively control the outflow dynamic, and the profiles show very different features. Firstly, the particles are notably slowed down, and a uniform velocity profile develops. Moreover, significantly lower solid fraction values are obtained for the lowest confinement width, while the system is better packed with increasing . As a result, the particle flow rate slightly increases with increasing .

Discussion and outlook
Our critical comparative analysis of the numerical and experimental data led to several well-founded findings. Firstly, the used DEM implementation predicted the particle flow rate quantitatively in quasi-2D hoppers with three different angles and several aperture sizes. We remark that this was done by identifying realistic values for the model parameters, namely, the stiffness, restitution coefficient, and friction coefficient of the steel beads. Such quantitative agreement has not been reported previously, using pure 2D approaches (simulation using disks). Our outcomes indicate that it is necessary to implement the rear and front planes as rigid and frictional walls to reproduce the experimental findings accurately. The role of confinement and the friction coefficient between the grains and the walls were the most significant parameters determining the particle outflow dynamics. We remark again that when implementing a hypothetically exact 2D approach, namely, the w = 0 and = d the flow rate values were significantly overestimated. It is worth mentioning that similar disagreement was reported by other numerical studies of quasi-2D silos, where the system confinement was not considered [7], and also that the packing properties of a confined particle mono-layer behave neither as a 2D nor as a fully 3D limit [11]. Interestingly, it has been found that the transit from the quasi-2D to a 3D configuration (decreasing the confinement width ) induces the presence of defects, which play a significant role due to the rearrangement of the grains spatial configuration [11].
Focusing on the kinematic profiles at the orifice region, we obtained that all of the simulated regimes are in reasonable agreement with the phenomenological expressions introduced very recently, predicting the hopper angle's role on the flow rate [25]. Thus, the velocities and volume fraction profiles at the exit line agreed with the experimental data. The agreement is also valid for the magnitude of the correlations and even for the spatial distribution of the fluctuations against their corresponding coarse-grained profiles. Nevertheless, we also encountered some systematic discrepancies between experimental and numerical results in Sec. 3. These differences are particularly noticeable for = 30 • hoppers. It that case, we found that the wall angle has a counter-intuitive impact on particle dynamics, reducing the mass flow rate. Remark that this has been reported experimentally in quasi-2D [25] and 3D hoppers [35].
The numerical approach allowed us to perform a systematic study exploring the consequences of the dynamics of the particles combined with the particle-wall influence. We found that the wall friction affects the outflow dynamic, reducing the particle flow rate significantly. Moreover, the impact of the degree of confinement on the outflow is also dependent on the wall friction. Hence, the outflow dynamic seems to be strongly influenced by the particle-wall collisions in the limit of considerable wall friction. This issue requires a more careful investigation, exploring, among other things, the evolution of the force chains, the contact network, and the developments of stable clogs.
Summarizing, the exposed results provide clear pieces of evidence about the complex interplay existing between friction, confinement, and particle collisions in quasi-2D silos and hoppers. All these factors influence the local values of particle velocities and solid-fraction, even when the resulting flow rates fit the experimental results closely. Thus, our analysis demonstrates that the naive 2D approximation of this 3D flow process fails to describe it accurately. Moreover, our outcomes strongly suggest that speed and solid-fraction can not be controlled independently by fine-tuning model parameters used to reproduce the macroscopic observations. Thus, new variables or more sophisticated particle-particle and particle-wall interaction rules must be analyzed to reproduce the microscopic In each case, outcomes corresponding to three different confinement width are shown dynamics more precisely. Nevertheless, we found that the disagreement between the numerical and experimental velocity fluctuations was not relevant. It indicates that the mean-field approximation introduced in [19] calculating the material kinetic stress seems to be justified, being very useful in analyzing the hopper angle's role on the discharging flow rate.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was partially funded by Universidad de Navarra (Grant number: PIUNA) and Ministerio de Economía y Competitividad (Spanish Government) through Project No. FIS2017-84631 MINECO/AEI/FEDER,UE. R B-R thanks the CONACYT fellowship program. G P-A and R B-R thank SEP-CON-ACYT project No. A1-S-46572 from the Mexican Government.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.