Simulation study of the discharge characteristics of silos with cohesive particles

In many industrial applications the silo for bulk materials is an important part of an overall process. Silos are used for instance to buffer intermediate products to ensure a continuous supply for the next process step. This study deals with the discharging behaviour of silos containing cohesive bulk solids with particle sizes in the range of 100-500 μm. In this contribution the TOMAS [1,2] model developed for stationary and non-stationary discharging of a convergent hopper is verified with experiments and simulations using the Discrete Element Method. Moreover the influence of the cohesion of the bulk solids on the discharge behaviour is analysed by the simulation. The simulation results showed a qualitative agreement with the analytical model of TOMAS.


Introduction
In industry silos are mainly used to store intermediate products in an overall process.In order to ensure a continuous supply for the next process step it is important to know the discharging behaviour of a silo.The discharging becomes more complex, when the silo contains cohesive bulk solids with fine particle sizes in the range about 100 µm or smaller.Since there are a large number of interfering factors, like the influence of the surounidng air on particles with d < 500 µm, affecting the discharging behavior of a silo an analytical calculation is difficult.However, there are many empirical approximations and models based on dimensional analysis.In this contribution the model of TOMAS [1,2] for stationary and non-stationary discharging velocity of a convergent hopper (1) is verified with simulations using the Discrete Element Method (DEM).
In the model in Eq. (1) t 76 is a time parameter, which describes the time needed to reach 76 % of the stationary outflow velocity and b min is a geometrical parameter, which corresponds to the minimum width of the silo outlet to ensure that no bridging occurs.Both parameters can be estimated with the granular flow properties of the bulk solids as given in [2].Moreover particle parameters, like particle diameter d p and density U s , the bulk porosity H as well as the fluid viscosityK f are included.By using a contact model which includes also attractive forces in DEM the influence of the cohesion of the bulk solids on the discharging behaviour is analysed in the simulations.Thereto the hopper was completly filled at the beginning and discharged afterwards.

Simulation methods
The DEM, introduced firstly in 1979 by Cundall and Strack [3], can be used to simulate the dynamics of particulate systems.The method describes each particle in the system and its interactions with other particles or rigid walls via contact forces F c .The contact forces can be devided in a normal component F n and a tangential component F t .In this study the particles are assumed to be spherical and indestructible.Based on the force balance for every single particle the equations of motion (2) and (3) are numerically solved.
In the equations m p , and J p are the mass of a particle and its moment of inertia while v p and w p are the velocity and the angular velocity.The tangential force F t,i due to the interactions with contact partner i is acting at the contact point with the radius r p,i and causes the torque M p,i .The contact forces need to be modeled with an appropriate contact model.For particles with a size less than 500 µm cohesive forces can become important.For that reason a model for the cohesion force F coh,i on the particles was included in the force balance.Also the drag force ‫ܨ‬ ௗ caused by the interaction of the particles with the surrounding air and the gravitational force ‫ܨ‬ is considered.The contact forces F c,i were modelled in DEM according to the soft sphere approach, which is based on the calculation of the forces depending on the overlap of the contact partners.For the calculation of the repulsive force in normal direction a non-linear elastic model, based on the Hertz theory [4], is used.To account for the energy dissipation during a collision viscous damping terms for the normal and tangential contact were considered, taken from the work of Tsuji et al. [6].
In the equation ( 4) for the normal force, E* is the equivalent Young´s modulus, R* the equivalent radius of the particles and δ n the normal overlap.The damping term which considers the energy dissipation depends on the equivalent mass m*, the relative velocity in normal direction δ ̇n and a damping factor α, which depends on the restitution coefficient e [6].The tangential force is modelled based on the non-slip solution of the theory provided by Mindlin and Deresiewicz [5] and is limited by Coulomb friction μ•F n : Equation ( 5) contains the tangential overlap δ t and the relative velocity in tangential direction δ ̇௧.Also the shear modulus G* is included.The attractive forces are modelled using a linear term as follows: where k is the energy density, which contais all cohesive effects, and A is the contact area.In this approach the particle cohesion in a contact increases proportionally to the applied force and forming contact area.A behaviour like this was qualitatively confirmed for silica particles by atomic force instruments (AFM) in [7].This effect will be dominant if the yield point of the particles is reached and the plastic contact area increases [8].The influence of the cohesion on the forcedisplacement curve is shown in Figure 1.

Fig. 1. Force-displacement curves with different cohesion
To account for the influence of the surrounding air on the particle motion, the drag force was calculated according to the work of Lord Rayleigh [9]: where U F and v F are the fluid density and fluid velocity, respectively.A p is the projection area of the particle and v p its velocity.Based on the assumption that the particles do not affect the air, a constant fluid velocity of v F = 0 is set.The drag coefficient c D depends on the particle Reynolds number Re P and can be calculated based on the equation of Kaskas [10]: ) As fluid properties the dynamic viscosity ߤ ி and the fluid density influence the particle Reynolds number.

Simulation results
For the simulations of the discharge behaviour of a convergent hopper, a monomodal particle system with a particle diameter of d = 500 µm was used.The energy density in the cohesion force acting on the particles was varied in the contact model.1.To compare the results obtained by the simulations of the discharging behaviour with the empirical approximation in Eq. 1 the minimum outflow width b min has to be known.b min can be calculated according to Jenike model [11] with paramters which can be measured by a shear test.However, in this study simulations were performed to obtain b min .The outlet diameter b of the geometry used in DEM was gradually reduced at a constant angle of the cone until no outflow of particles occured.This was determined by considering the average axial velocity of all particles.If the average velocity drops below a critical value for at least 10 seconds, the simulation was terminated.The results of b min determined in this way are shown in Table 2.The horizontal cross-section, displayed at the right side in Figure 3 at the bottom, shows that the particle velocity is minimal in the area of the wall (green) and increases towards the center.

Table 1. Simulation setup
In Figure 4 the averaged particle velocities at the outflow of the hopper are shown over the time.The data points show the results of the simulation while the lines are the solution of the empirical equation (1).The colours represent the different energy densities.The simulation results show a decreasing stationary outflow velocity with increasing cohesive forces.Comparing the course of the simulated velocities with the empirical results, it can be seen that the time needed to reach the stationary velocity is in a good agreement.Furthermore, it can be seen that the velocities obtained from the simulations are generally higher than the velocities calculated with Eq. ( 1).This deviation is in the same order of magnitude as the deviation showed by Tomas in [1] where the analytical results of his model for stationary velocity was compared with experiments.

Fig. 4. Averaged particle velocity at the outflow over time for different cohesion forces
The high fluctuation of the velocity at k 3 (Fig. 4) is caused by agglomeration of particles to clumps in the silo outlet.This effect is illustrated in Figure 5 which shows the particle velocities (left) and the porosity (right) in rzcross sections at t = 0.03 s.Due to bridging, agglomerates are breaking out of the bulk, what leads to an increase in the porosity between the agglomerate and the bulk.Figure 4 also shows, that for k 3 the maximal outflow velocity at t = 0.21 s is higher than the stationary outflow velocity.This effect is also a result of the breaking out agglomerates based on bridging at the beginning of the outflow.the porosity increases in the lower regions of the hopper (Figure 6, left).At a certain height of the hopper the width is getting big enough that briding not longer occurs.Because of the higher porosity in the lower regions of the hopper the particles achieve high velocities.This in turn leads to a decrease of the porosity at t = 0.22 s (Figure 6, middle) due to the compaction at the outlet.Finally, as can be seen in (Figure 6 right), at t = 0.27 s the porosity slightly increases again until the outflow is getting stationary.The simulations also allow the determination of the velocity at the outflow as a function of the radial position of the particles (Fig. 7).The maximum of the velocity is reached in the centre of the hopper.With increasing radial position the velocity drops to smaller values.The velocity gradient is getting higher with decrasing distance to the wall at r = 10 mm.The curves for the different energy densities are qualitatively similar.However, with increasing energy density the values are shifted towards smaller velocities.At the high energy density k 3 in which b min is in the range of b the averaged velocities are significantly lower then in the other cases.

Conclusions
In this work the discharge characteristics of a silo hopper which contains a cohesive bulk solids was investigated by DEM simulations.By varying the cohesive forces between the particles, the effect of the cohesion on the outflow was analysed.To model the contact force the Hertz-Mindlin model combined with a linear cohesion term was used.A comparison of the outflow velocities obtained in the simulations with velocities calculated with an empirical equation of TOMAS shows a qualitative agreement.Howeever, the simulated velocities are higher than the velocities of the empirical model.At a very high energy density in which b min is in the range of b, the simulation showed that the outflow velocity does not increase steadily up to the stationary velocity.The sudden changes are caused by the breaking out of agglomerates from the bulk.Such phenomena can not be depicted with the equation of TOMAS.Furthermore, the simulation results showed, that the influence of the cohesion on the stationary massflow is bigger than on the stationary velocity because of the changing porosity.In the next step the discharging of silos containing a highly cohesive bulk material will also be investigated in experiment.The velocities and mass flow rates which can be found in these experiments will also be compared.For this purpose, the flow properties of the bulk solid and particle contact properties were determined in shear tests using the methods described in [12].

Fig. 2 .
Fig. 2. Geometry of the simulated hopper In Figure 2 the geometric dimensions of the convergent hopper are shown.The hopper has an outflow width of b = 20 mm and an angle Θ = 10°.The simulation paramters are given in Table1.

Table 2 . 730 Figure 3
Figure 3 exemplarily shows a snapshot of the particle velocities after a simulation time of t = 0.02 s for k 1 .At the right side a rzcross section is displayed at the top.The maximum velocities (red) are reached directly above the outlet, while particles in the upper regions close to the wall have still a velocity of v = 0 m/s (blue).

Fig. 3 .
Fig. 3. Snapshot of the particle velocities after a simulation time of t = 0.02 s for k 1

Fig. 7 .Fig. 8 .
Fig. 7. Averaged particle velocity at the outflow over radial position for different energy densities in stanionary case

Figure 8
Figure8shows a comparison of the mass flow rate at the outlet over time for different energy densities.It can be seen that the mass flow rate decreases with increasing values for the energy density.On the one hand this is a result of the decreasing outflow velocities, but on the other hand it is also caused by an increase in porosity.