The role of the particle aspect ratio in the discharge of a narrow silo

The time evolution of silo discharge is investigated for different granular materials made of spherical or elongated grains in laboratory experiments and with discrete element model (DEM) calculations. For spherical grains, we confirm the widely known typical behavior with constant discharge rate (except for initial and final transients). For elongated particles with aspect ratios between 2 ⩽ L/d ⩽ 6.1, we find a peculiar flow rate increase for larger orifices before the end of the discharge process. While the flow field is practically homogeneous for spherical grains, it has strong gradients for elongated particles with a fast-flowing region in the middle of the silo surrounded by a stagnant zone. For large enough orifice sizes, the flow rate increase is connected with a suppression of the stagnant zone, resulting in an increase in both the packing fraction and flow velocity near the silo outlet within a certain parameter range.


Introduction
The discharge rate of inelastic granular material out of a container is known to be independent of the filling height, thus constant in time, which is opposed to the case of a liquid where the flow rate decreases with filling height. This feature of granular flow is useful for industrial applications. It enabled humankind to construct the hourglass, a reliable device to measure time, known already in the middle ages. Various aspects of granular discharge have been investigated including the flow rate as a function of the orifice size (Beverloo law [1][2][3][4][5]), clogging phenomena for small orifices [6][7][8][9], the geometry of the flow field and possible stagnant zones [10][11][12][13][14][15], fluctuations in the flow velocity and density [12,16,17], resonance phenomena (silo music) [18][19][20], influences of the particle properties, such as surface friction [21][22][23][24], stiffness [24][25][26] or shape [9,15,[27][28][29][30][31][32][33]. Recently, the discharge of soft particles with very low surface friction was shown to be rather similar to the case of a liquid with decreasing flow rate during the process [24].
Small-scale laboratory experiments and numerical modeling are often used to quantify the above-described effects. Exploring the silo discharge, the studies typically use simplified systems, e.g., with spherical particles, two-dimensional or quasi-two-dimensional geometries, and limited system size. Despite these simplifications, many of the basic features of silo discharge can be captured and quantified, e.g., the dependence of the flow rate on the orifice size [3,4], the typical stress distribution within the whole silo and especially in the vicinity of the orifice [5,34,35], etc. The flow of elongated particles leads to orientational ordering and realignment of the grains in the sheared regions [36][37][38]. Various studies investigated the effects of grain shape on the clogging properties [9,39], the geometry of the flow field [10,12,15,30,32,40], or discharge rates [9,27,28,31]. Such studies are even more limited for the system size because modeling of elongated grains in numerical investigations is more time-consuming (either for multi-sphere grains, spherocylinders or ellipsoids), and the determination of particle orientations in laboratory experiments requires more effort than only detecting particle positions.
In the present work, we show that in a typical laboratory-scale experiment with a narrow silo, the discharge process of elongated grains has a peculiar time-dependent behavior. Namely, the flow rate notably increases near the end of the process. Our findings were reproduced numerically by a 3D model with the same dimensions containing about 100 000-400 000 particles. We emphasize that previous work using similar system sizes and grain shapes did not explore this behavior [9,[27][28][29][30][31][32]40]. Here, we thoroughly investigate the cause of the increasing flow rate by exploring the spatial structure of the packing fraction and flow velocity inside the silo.

Experimental setup
The granular material is filled in an acrylic cylindrical silo with an inner diameter of D silo = 144 mm, a height of H = 800 mm, and a circular hole at the center of the bottom (see figure 1). The diameter D of the hole can be adjusted using replaceable inserts.
The basal plate of the silo is separated from the silo wall with a gap of approximately 1 mm. The basal plate is supported by two load cells in order to monitor the vertical basal force. The discharged material is collected in a container, while its mass is measured with another load cell.
In the experiments presented in this work, five granular samples were used: glass beads (Sigmund-Lindner GMBH, Germany) and plastic rods (cylinders with length L and diameter d) with elongations L/d = 2, 3, 4.5, and 6. The plastic rods were obtained by cutting d = 2.4 mm nylon trimmer lines (Oregon) with a cable cutting machine (Model GlobalCut 100, Navia International GMBH, Germany). Table 1 shows photographs, of the particles and also indicates their length L and diameter d. The equivalent diameter d * is defined as the diameter of a sphere with the same volume as the rod. Table 1 also indicates the value of d * for each particle.
After filling the silo manually with a filling rate of about 500 cm 3 s −1 directly from the top of the silo, the material was discharged, and the evolution of the discharged mass and the silo basal force were continuously recorded during the whole discharge process with a sampling rate of 1 kHz. A typical discharge process lasted for about 50-100 s. The data series were smoothed by convolution with a Bartlett Table 1. Pictures of the granular samples (right): glass beads and plastic rods with elongations of L/d = 2, 3, 4.5, and 6, and sketches of the particles used in the DEM simulations (left). The grains are characterized by their length L, diameter d, the equivalent diameter d * and their surface friction coefficient μ. The equivalent diameter is defined as the diameter of a sphere with the volume as the elongated particle. For the numerical system, the number of grains used in the simulations is also indicated.
(triangle) window of a width of 1 s to reduce noise and fluctuations. The flow rate was obtained by differentiating the curves of the discharged mass and smoothing the result again with the same smoothing procedure. We mention that due to the impact of grains, there is a momentum transfer from the particles to the container. If the flow rate is constant, this effect only causes an offset in the measured mass, thus has no influence on the time derivative, except for an initial and final transient. When the flow rate is changing, the impact momentum is also changing, but in our case its contribution is minor compared to the main effect of the increasing mass in the container. This effect was thus neglected when calculating the flow rate from the load cell data. The value of the flow rate is presented in grains/s, using the average unit mass of the individual grains, which was determined by measuring the mass of 200 particles.

Numerical model
Numerical simulations were carried out in order to gain further insight into the dynamics of silo discharge with elongated grains. For this purpose, a custom hybrid GPU-CPU discrete element modeling (DEM) code was used, taking advantage of the computing power of several NVIDIA graphical processing units [41][42][43]. The simulations resemble the experimental setup, modeling a cylindrical container of diameter D silo = 144 mm with a flat basal plate, in which there is an orifice centered at the cylinder z-axis. The particles are defined as spherocylinders with aspect ratio L/d, where L is the end-to-end length and d is the diameter. Initially, N particles were loosely inserted into the simulation domain at random locations and with random orientations. N was chosen such that the height of the dense packing is near the experimental value of 800 mm, in the range of [750-800] mm (see table 1). Note that this quantity depends on the particle aspect ratio, since the packing fractions and particle volumes differ. Next, the particles fell freely under gravity, forming a dense granular system as the initial configuration. After the packing became static, the orifice was opened so that the grains were allowed to discharge until the process ended.
The DEM algorithm numerically computes the translational and rotational motions of each particle i (i = 1, . . . , N). As input, it utilizes all pairwise interactions between contacting particles i and j. The particle-particle interaction was computed using an algorithm for interacting spheropolyhedra [44,45], implemented on GPU architecture [42]. The particular case of a spherocylinder is defined by a sphero-radius r = d/2, and the location of its two vertices. Its surface is delineated by all points at distance r from the edge delimited by the two vertices. The contact detection between two spherocylinders consists in finding the closest point between two edges. It results in an overlap distance δ n , defined by the overlap of two spheres of radius r, located at each of the edges.
In the model, the contact force F ij has a tangential and a normal component, We use a linear spring-dashpot model, where the normal force is F n ij = −k n δ n − γ n v n rel . Here, k n is the spring constant in the normal direction, δ n denotes the overlap between particles, γ n is the normal damping coefficient, and v n rel is the normal relative velocity between particles i and j. The tangential force F t ij models the friction between the grains. It takes into account Coulomb's friction constraint, namely, where μ is the friction coefficient. At each contact, δ t represents the elastic deformation of the tangential spring with elastic constant k t , which increases as d δ t /dt = v t rel while the particles are overlapped. Moreover, γ t is the damping coefficient in the tangential direction, and v t rel is the tangential component of the relative velocity of the overlapping pair.
The translational and rotational equations of the particle i read as where m i and I i are the mass and the moment of inertia tensor of the particle i, respectively. In the implementation, the sums only account for the N c (i) contacts of each particle. τ ij are the torques corresponding to each contacting force F ij . The acceleration g represents the acting gravitation. Moreover, v i and ω i are the translational and rotational velocities. A velocity Verlet numerical algorithm [46] is used to integrate the translational equations of motion of each particle, while the rotational degrees of freedom were resolved using a modified leap-frog algorithm [47].
The main model parameters of the numerical analysis correspond to plastic particles with ρ p = 1120 kg m −3 , k t = 2/7k n , γ n = γ t , e n = 0.9, g = 9.8 m s −2 and k n = 2 × 10 5 (m p g/d). The particle friction coefficients were set to μ = 0.3 (plastic particles) and μ = 0.5 (glass particles). Similarly to the experimental procedure, we executed a systematic study of the system response, varying the particle aspect ratio L/d = 2, 3, 4.5, 6 and the size of the orifice D = 2R in terms of the equivalent diameter d * , which is defined as the diameter of a sphere with the same volume as the spherocylinder.

Coarse-graining (CG) micro-macro technique
The continuous description of the granular flows often brings useful insight into their complex behavior. DEM algorithms provide the position and velocities of every particle, as well as the interaction forces. To access the relevant macroscopic fields, we use a well-accepted coarse-graining technique [48][49][50][51]. This methodology requires the introduction of a non-negative integrable function ϕ i ( r ), which maps the microscopic details into macroscopic fields. Here, we use a truncated Gaussian function , with A −1 ω chosen such that the integral of ϕ i ( r ) over all the space results in 1, the cutoff distance is r c = 4ω and ω = L/4. Thus, the macroscopic fields of density ρ( r, t), momentum P( r, t), and velocity V( r, t) = P( r, t)/ρ( r, t) are obtained from the particle positions r i and velocities v i , respectively. Complementarily, the density gradients and the components of the shear rate tensor γ ij are computed via central difference.
The flow-induced orientation of the non-spherical particles is very relevant in our context. The particle orientation field is characterized by the orientation tensor O( r, t), defined in terms of the particle orientation l i = [l x , l y , l z ] i , which is a unit vector parallel to the particle long axis. It reads The stress field σ( r, t) is computed using the formulation introduced in references [48,49]. Following this approach, the stress σ( r, t) = σ k ( r, t) + σ c ( r, t) has two contributions, the kinetic stress field σ k ( r, t) and a contact stress field σ c ( r, t).
The contact stress field σ c ( r, t) is computed from the interaction forces and branch vectors r ij by the line integral of ϕ( r ) along r ij as follows, In the implementation, the sum over j only accounts for the N c (i) contacts of each particle. On the other hand, the kinetic stress σ k ( r, t) accounts for the velocity fluctuations relative to the mean velocity field V( r, t). In our case, however, σ k ( r, t) resulted in contributions several orders of magnitude lower than σ c ( r, t).
Based on the previous formulation, a post-processing data analysis tool was implemented. It allows the computation of all the relevant macroscopic fields, using DEM outcomes as input data.

Macroscopic response experimental approach
In the experiments, the flow rate and the basal force F b were measured during the whole discharge process. Figure 2 shows the evolution of these quantities as a function of the mass in the silo for glass beads and plastic rods with elongations L/d = 2, 3, 4.5, and 6. As we see, the nature of the discharge is gradually changing with increasing elongation of the grains: for spherical grains, the flow rate is basically constant during the discharge process (except for initial and final transients), while for elongated grains, this is not always true, as for large orifice sizes the flow rate is clearly increasing in time. For a smaller orifice, the flow rate remains nearly constant also for elongated grains, similarly to the case of beads. Note that we chose D/d * ≈ 11 as the lower limit of orifices due to the clogging of L/d = 6 particles below this orifice size.
The basal force F b was normalized by the weight W D silo of the material, corresponding to a filling height of the diameter of the silo. As we see in figure 2, the value of F b is much smaller than the weight of the granular column for all cases. For glass beads, the evolution of F b shows a plateau (right after the initial transient) and then a gradual decrease in the second half of the experiment. This behavior is in accordance with the Janssen effect, i.e., part of the weight of the granular column is held by the silo walls due to arches in the contact network. Interestingly, for elongated grains, the plateau is much less pronounced, and the value of F b is larger. This is similar to our previous observations with lubricated (reduced friction) glass beads (see figure 5(d) of reference [24]). This might be an indication that the friction between the rods (which are mostly aligned parallel with boundaries) and the silo wall is smaller than that in the case of glass beads. Note that for each material, the curves measured for different orifice sizes overlap, which shows that the basal force does not depend on the flow rate (at least in the examined range of orifice sizes.) The measurements presented in figure 2 were obtained in the orifice range of 42 mm D 64 mm corresponding to 11.14 D/d * 16.18. Increasing the orifice size to very large values would lead to an accelerated flow regime, basically the same as without constriction at the bottom of the cylinder. Are the curves with increasing flow rate in figure 2 already the precursors of such a transition? Figure 3 confirms that this is not the case. Here, we plot the flow rate curves for glass beads and rods with L/d = 4.5 up to the orifice size of D = 102 mm, i.e., D/d * = 22.5. The sequence of curves for elongated particles shows that the increasing flow rate scenario is present in a limited range of orifice sizes, and at around D/d * = 22, we recover the constant flow rate scenario. The flow rate as a function of time is shown in figures 3(c) and (d), presenting the typical timescales of the process. Another important observation is that the flow rate increase is not particularly sensitive to the initial filling height. In figure 3(b), three measurements are presented for the orifice of D/d * = 18.5, with the initial fill heights of H/D silo = 5.5, 4.2, and 2.8. We see that the three curves nicely overlap, i.e., the flow rate increase is the same for all three measurements.
We can characterize the flow rate increase by two quantities: the peak flow rate and the initial flow rate. The initial flow rate is defined as the flow rate right after the transient when the material in the silo is between 5m D silo > m > 4m D silo . Figures 4(a) and (b) show the initial and final flow rates as a function of the reduced orifice size D/d * for all five samples. We see that while the final flow rate increases similarly for all materials (at least up to about D/d * = 15), the increasing trend of the initial flow rate is getting significantly weaker with increasing aspect ratio L/d. This leads to an increasing value of the ratio of the final and initial flow rates with increasing aspect ratio of the particles, as we see in figure 4(c).

Macroscopic response numerical approach
DEM simulations were performed for a direct comparison with the experimental system. In figure 5, we plot the flow rate and the basal force F b as a function of the mass in the silo. For clarity, the data corresponding to the same cases discussed in figure 5 are shown.
In general, the numerical and experimental results are in quantitative agreement with a difference of up to 20%-30%, which can be attributed to the slightly different particle shapes (spherocylinders vs rods with slightly irregular ends) or surface friction of the grains. As we can see, DEM simulations reproduce the most  relevant features of the discharge process, with constant flow rate for spherical particles and a complex behavior for elongated grains. Thus, the flow rate for elongated grains shows an increasing trend for large orifices. However, the behavior changes as the orifice size reduces, resulting in a nearly constant value. For very small orifices, and after half of the process, the flow rate slightly decreases towards the end of the process, denoting the increasing risk of getting clogged.
The basal force F b in the numerical simulations also shows a similar trend as the experimental data, with overlapping curves for different orifice sizes (no flow rate dependence) and after a nearly constant value in the first half of the discharge process, there is a gradual decrease in the second half.

Coarse-graining analysis of the granular flows
In order to calculate the relevant macroscopic fields from the DEM data, we execute a coarse-graining analysis, which was introduced in section 3.1. As a first step, our analysis focuses on the behavior of the mean packing fraction φ , grain velocity in the vertical direction v z , and vertical component of the orientation tensor O zz , averaged in the orifice region. We define this region as a cylindrical volume centered around the silo axis with its base at the orifice level, with the same diameter as the orifice size and a height of L/2. The evolution of these quantities is presented in figure 6. We observe an overall non-monotonic trend of φ by increasing the particle aspect ratio L/d with the largest value of φ for particles with L/d = 2. This is coherent with earlier observations on the packing fraction of disordered jammed granular materials, where a similar non-monotonic dependence was found as a function of the particle elongation L/d, with the largest value of φ at around L/d = 1.5 [52,53]. For particles with L/d > 2, the value of φ monotonically decreases, since particles with larger elongation tend to form more dilute packings with a lower degree of orientational ordering [53].
In detail, examining spheres (top row, first panel), we see that φ is practically constant during the discharge, and it is slightly growing with increasing orifice size D/d * , in agreement with previous numerical findings [54].
Focusing the attention on moderately elongated particles, i.e., for L/d < 4.5, φ varies only slightly during the discharge process, regardless of the orifice size D/d * . Typically, we obtain a slightly increasing value of φ throughout the discharge process, with the largest φ values are always detected near the end. This trend also holds for very long rods L/d > 4.5 and small orifices, but a distinct behavior is depicted for large orifices. In those cases, the flow is getting significantly diluted in the beginning, and a growing value of φ is observed only in the second half of the process.
When inspecting the evolution of v z (figure 6, middle row) for spheres, we obtain v z practically constant during the discharge. As expected, the value of v z typically rises with increasing the relative orifice size D/d * . For elongated particles, however, near the end of the process, we find increasing velocity for large orifices but decreasing velocity for small orifices. Consistently, this behavior is in line with the obtained trends for the particle flow rate when varying D/d * , both experimentally and numerically. Here, we can argue that changes in the flow patterns lie behind these different trends.
The behaviors of φ and v z suggest that there are significant changes in the morphology of the flow during the discharge. In figure 6, bottom row, we present the mean value of the vertical component of the orientation tensor O zz , which shows that the particle orientation also changes during discharge. In general, its behavior correlates with the trend of φ in all cases. This indicates that the level of alignment and order of the flowing structures determines the macroscopic volume fraction at the orifice. Accordingly, these microstructural changes also correlate with the macroscopic particle flow rate obtained numerically and experimentally.
In order to better understand the complex nature of the data presented above, in the following, we extend our analysis by investigating the spatial profiles in the whole silo. Figure 7 illustrates packing fraction φ (r, z), velocity v(r, z), contact pressure p c (r, z) = 1 3 Tr(σ c (r, z)), and shear rateγ(r, z) by presenting 4 The results shown in figure 7 indicate that the nature of the discharge process is very different for beads and rods. For beads (column I), the macroscopic fields are relatively homogeneous, while for rods (II-IV) strong gradients are observed. For rods, we see a fast flowing central region surrounded by a slowly moving (or stagnant) zone, with a strongly sheared region in between. This is coherent with our earlier observations in a quasi-two-dimensional hopper [12]. Looking for the cause of the flow rate increase near the end of the discharge process in case III, we now compare cases II and III, corresponding to the same particle elongation L/d = 6.1 but different orifice sizes D/d * = 10.95 and D/d * = 12.97, respectively. Focusing the attention on the φ colormaps at the orifice region, it is noticeable that in the last stage, the values are comparable in both cases. However, in case III, φ is significantly lower in the first three snapshots than in case II, while in case II such flow-dilation is only observed in the first snapshot, after which the packing fraction at the orifice quickly increases. Examining the velocity fields, they also indicate important differences in the nature of the flow field. For case II, a well developed stagnant zone is formed, and it also expands in time. It seems to increase the systems' effective friction, leading to a decrease in the flow velocity in the orifice region (see the blue velocity curve for L/d = 6.1 in figure 6). This compensates for the slight increase in the packing fraction and altogether leads to an approximately constant flow rate for case II. In contrast, it does not occur in case III where the material around the surging inner region is also flowing, and it is hard to identify a stable stagnant zone. As a result, there is a significant increase of the flow velocity near the end of the discharge process (see the green velocity curve for L/d = 6.1 in figure 6). The above described two ingredients: (i) the increase in the packing fraction in the orifice region and the suppressed stagnant zone together lead to an increased flow rate for case III near the end of the discharge process.
As we have seen, the scenario of case III is only observed for a relatively large orifice (in relation to the inner silo diameter), when the stagnant zone is suppressed. Can we get back to the constant flow rate scenario by increasing the silo width and thereby recovering the stagnant zone? This configuration is shown in our last test (case IV in figure 7) in a silo with a 40% larger diameter. As we see, the velocity field of case IV resembles that of case II with a clearly developed stagnant zone. With the reappearance of the stagnant zone, the dilation in the region above the orifice also disappeared. Altogether, this results in a constant flow rate throughout the discharge process in case IV, similar to case II. We also note that when the orifice size is further increased we have observed a high constant flow rate as shown in figure 3 for D/d * = 22.5. The mesoscopic fields for this case (not shown here) are very similar to the ones observed for spherical grains (case I in figure 6) without strong gradients.
In figure 8 we visualize the experimentally determined discharge flow regimes in the explored range of the parameters. The black circles and red crosses indicate where the system exhibits constant or increasing flow rates, respectively. The increasing flow rate scenario is only observed in a finite range of the parameters. As expected, in the limit of a very large silo diameter compared to the orifice size, the classical constant flow rate response is achieved, regardless of particle shape and size. In the opposite limit, when the orifice size is comparable with the silo dimension, the constant flow rate is recovered.

Summary
We have shown that the discharge of elongated particles from a narrow silo can lead to strongly time-dependent outflow rates with a significant increase before the end of the discharge process. This is not observed for spherical particles of similar size in the same silo. The peculiar scenario for elongated grains is observed both in a typical size lab experiment and in corresponding numerical (DEM) simulation using 100 000-400 000 particles. The analysis of the numerical data revealed that the difference between the case of beads (constant flow rate) and elongated particles (time-dependent flow rate) originates from the difference in the flow field. Namely, for spherical grains, we observe smooth, homogeneous fields for the packing fraction, velocity, or contact pressure, representing a mass flow type discharge. For elongated grains, however, these quantities show large spatial variations inside the silo, with a strong gradient representing funnel flow behavior. In particular, a fast-flowing region is observed in the middle of the silo, which is surrounded by a slowly moving (or stagnant) zone. This is in accordance with our earlier observations in a quasi-two-dimensional silo [12], and it leads to slower discharge than for beads. For a large enough orifice size, however, the stagnant zone is suppressed, leading to an increase in the packing fraction and in the flow velocity before the end of the process. Altogether, the increasing flow rate scenario is observed in a finite range of the parameters.