Stress pulsation mechanism during filling and discharging granular materials form silos

Abstract: In this paper, we investigate numerically the dynamic pulsation nature of stresses acting on hopper walls during filling and discharging granular materials from silos. From the results, at the end of filling process, it can be concluded that the larger the frictional wall coefficient is, the less the ultimate stress values acting on the hopper wall will be. Furthermore, for smaller values of internal friction, a consistent stress fluctuation appears for a slightly long period of time before the stress trace turns to a straight line. For the discharging process, the required time for the discharging of the silo is almost the same for different values of wall frictions. Moreover, for large values of internal friction coefficients, the discharging process takes more time.


Introduction
The silo system has been considered as one of the emerging issues in the granular materials research over many years. Past and current research efforts focus on investigating many aspects relevant to the silo system. Problems in silo, such as flow blockage, jamming, segregation, dead zones, and silo collapse attract the interest of researchers and raise the importance to study the storage vessels. Great efforts have been made to closely study aspects such as flow pattern of granular material through silos, stress field propagation along silo walls, force fluctuation at the boundary of granular flow in silo, and the effect of particle shape on silo discharge.

PUBLIC INTEREST STATEMENT
The silo system has been considered as one of the emerging issues in the granular materials research over many years. It is known that granular materials cannot be characterized as gas and fluid or solid. Discrete element methods (DEM's) are numerical techniques capable of simulating the entire behavior of systems of discrete interacting elements. In the present paper, by using this numerical simulation technique, we study the effects of the friction between the particles themselves and between the particles and the boundary walls in a silo system during filling and discharging processes. It is found that, for large values of friction between the particles' system, the discharging process takes more time. Moreover, the required time for the discharging of the silo is almost the same for different values of particle/ wall frictions.

Implementation of DEM model
The present numerical simulations are based on the soft-particle discrete method which originally was proposed by Cundall and Strack (1979). In our simulation model, the granular particles are subjected, during simulation, to two types of forces including contact forces and gravitational force. Unlike the gravitational body force, the particle contact forces act only when the particle is in contact with other particles and/or with a boundary wall. The normal contact force is determined by a damped linear spring model, while the tangential contact force is determined by a linear spring in series with a frictional sliding element model. The interaction forces developed between two particles or between particle and boundary wall are calculated based on the physical properties and the relative velocities. The dynamic equations for the motion for each of the particles under these forces are (1) where m i , I i , ⃗ r i , and ⃗ i (t) are, respectively, the mass, rotational moment of inertia, position, and rotational vectors of the centre of particle i, ⃗ F gi , ⃗ F c,ij and ⃗ F c,iw k are, respectively, gravitational body force and contact forces acting on particle i due to particle j and boundary wall w k , N p and N w are, respectively, the number of particles and walls in the simulation. For simplicity, we denote Therefore, Equations (1 and 2) become Unlike the gravitational body force, the particle contact forces act only when the particle is in contact with other particles and/or with a boundary wall. In the analysis of the system, any contact force is decomposed into normal and tangential components. For any typical pair of particles i and j, if they are in contact, then the contact force is determined by the following formula where n and ŝ are unit vectors in the normal and shear directions of the contact plane, ⃗ F̂n ,ij (t) and ⃗ F̂s ,ij (t) are, respectively, the magnitudes of the normal contact force and shear contact force, namely, where k̂n and k̂s are, respectively, the particle-particle normal and tangential spring coefficients, ⃗ F e ij,n (t) and are, respectively, the elastic contribution of the contact force between the particles i and j in the normal direction (n direction) and the friction coefficient of the granular particles, respectively, the normal compression and the tangential displacement between the particles and over the time step Δt = t -t 0 , R i and R j are the radii of the particles i and j. In order to reduce the number of contact checks for each particle, a neighboring cell technique is implemented in the present simulations. Once all forces and moments acting on particle due to contacts are identified, new velocities and positions of all particles are computed by numerical integration of Newton's equation of motion. Setting Therefore, ∀i = 1, 2, 3, …, N p , we have a system of first-order ordinary differential equations as follows Thus, setting yields a system of 4N p first-order ordinary differential equations that represents the dynamic equations for the motion of the N p particles and can be expressed by By using the central difference scheme and considering that ⃗ d⃗ The new translational and rotational velocities ⃗ v i and ⃗ i , at time t = N + 1, can be determined by, In order to reduce the number of contact checks for each particle, neighboring cell technique is implemented in the present simulations. Once the particle contact forces are calculated, the new velocity and position can be obtained by integrating the differential equations for particle motion over small simulation time step. The input parameters in our simulations include: physical properties of the granular particles under consideration, initial conditions, and boundary conditions. The initial conditions include: the initial positions and velocities of all particles and the geometry of the boundary conditions. All parameters are non-dimensionalized using the density of particle, gravitational acceleration, and particle diameter. A list of used mechanical and environmental dimensionless simulation parameters are given in Table 1.

Model setup and simulation schemes
The silo system under consideration consists of two connected silos one above the other as shown in Figure 1. To simulate the real silo filling process, the particles are modeled as a conglomeration of discs with a uniform distribution of diameters, 0.81 ≤ d ≤ 0.99 mm. The discs dispersed inside the upper silo, with given random initial positions and velocities, and are allowed to fall under gravity down into the upper silo. During the filling process, the outlet of the upper silo is kept closed. Hence, due to damping forces, the particles ultimately come to rest at the bottom of the silo. Once the upper silo is almost filled with particles and the system has achieved a nearly steady-state condition, the outlet of the upper silo is opened and the particles start to flow out of the silo under gravity toward the lower silo. To study the nature of normal stresses exerted on the hopper wall of the lower silo, we divide the hopper wall into five equidistance segments as shown in Figure 1. The segment length is approximately five times the particle diameter.
During typical simulation, each segment element of the hopper's wall bears stresses that are exerted by the granular material. To analyze these stresses, we suppose that segment j (where 1 ≤ j ≤ 5) has m contact points c 1 , c 2 , ..., c m as shown in Figure 2.  Particle density, ρ/ρ 0 (where ρ 0 = 1.0 × 10 3 kg / m 3 ) 2.5 Particle-particle normal spring stiffness, k n,pp ∕ 0 ⃗ gd 2 0 1.266 × 10 5 Particle-particle normal dashpot coefficient, C n,pp ∕ 0 ⃗ g 0.5 d 2.5 0 2.524 × 10 2 Particle-particle tangential spring stiffness, k s,pp ∕ 0 ⃗ gd 2 0 1.266 × 10 5 Then the corresponding contact forces acting on segment j are  We define the average contact force, ⃗ F ave j (t), acting on segment j(with length l seg j = 4.66 mm) due to the m contact points by such that and are, respectively, the normal and tangential components of ⃗ F ave j (t) and n i ,ŝ i represents normal and tangential unit vectors corresponding to the segment j at the contact points. Therefore, the total stresses acting on segment j, namely, ⃗ S T j (t) can be obtained by where ⃗ F j,g x (t) and ⃗ F j,g y (t) are, respectively, the gravity force acting on segment j in the x and y directions.

Flow pattern and velocity field
Once the outlet of the upper hopper is opened, the particles start to flow out of the silo under gravity toward the lower silo. The flow pattern and the velocity field of the granular material at different simulation time steps are shown in Figure 3. Mass flow can be observed to occur in the entire upper silo. At the point of transition where the silo wall meet the hopper wall, the particle moving direction changes sharply from straight down to a direction parallel to the hopper inclined wall. Consequently, during the discharging process, the stress on the hopper wall reaches the maximum at these two points and decreases toward the outlet. Hence, we expect that the value of the acting stress on segment five of the hopper wall to be the maximum. However, it has been observed that the more particles reach the point close to the outlet the more they gain room to move. Thus, the stress profile of segment one may more likely show an oscillating nature than other segments of the hopper's wall.
Even though mass flow can be seen to occur in the entire silo as can be seen in Figure 4(a), it is possible to recognize different zones with different flow patterns. These flow zones can be divided into: pipe, pipe feed, and plug flow zones. Figure 4(b) shows the three different flow zones.

Simulation results and discussion
In this section, we investigate the stress fluctuating nature developed along the hopper wall during filling and discharging of a silo. A simulated filling and discharging processes were carried out to study the variation in the stress that acts on the hopper wall, with time. Moreover, we investigate the dynamic nature of these stresses at various values of the wall and the internal friction coefficients.

Stress pulsation mechanism
We, partially, filled both the upper and the lower silos by allowing particles, having random initial positions and velocities, to fall down under gravity. The upper silo was filled with 2,500 particles while 4,500 particles were placed in the lower silo. Then the outlet of the upper silo was opened and the particles start to flow down from the upper silo under gravity toward the lower silo. Snapshots of the process at different time steps are shown in Figure 5.
At t = 0 s, each of the five segments of the lower silo wall, bears an initial stress value due to the initial existing particles in the lower silo. The particle packing plays an important role on the initial stress profile on these segments. Once the particles that flow out from the upper silo strike the static surface of the granular material in the lower silo, the stress value acting on each segment starts to increase. Figure 6 shows the stress profile acting on each segment of the hopper's wall of the lower silo throughout the simulation time.
At t ≅ 12.2 s, particles from the upper silo started to hit the static surface of granular material in the lower silo. Consequently, stress trace of each segment starts to fluctuate as the height of bulk material in the lower silo starts to increase. For each stress's profile, two fluctuation traces can be identified, namely, inconsistent and consistent stress fluctuation traces. The inconsistent stress trace dominates as long as there are particles falling down from the upper silo and hitting the granular surface in the lower silo. As shown in Figure 6, it extends from t ≅ 12.2 -64.6 s. Once all the particles of the upper silo reach the lower silo and settle down, the stress trace on each segment does not turn to a straight line directly. A temporary consistent fluctuation which has a harmonic motionlike trace starts to dominate before the stress profile turns to a straight line. Therefore, the system possesses a temporary dynamic pulsation stress-like nature on the hopper wall after it appears to reach a near static state. The stress profile on segment five which is the closest to the surface of bulk material appears to have larger trace fluctuation amplitude than that on other segments of the hopper wall during the process. Figure 7 shows the dimensionless mean stress acting on the hopper wall of the lower silo at two different simulation stages, namely, at t = 0 s and t = 64.6 s. It is obvious that the stress value on each segment increases and almost on all segments the stress increases monotonically throughout the simulation time. In both cases, the stress acting on segment one appears to be the largest and the wall stress decreases gradually along the hopper wall toward segment five. This result is, generally, in agreement with existing experiment results (Cleary & Sawley, 2002;Goda & Ebert, 2005;Masson & Martinez, 2000).

Influence of wall friction
In this section, we investigate the effect of the wall friction on the dynamic nature of the wall stress developed during filling and discharging of the lower silo. A sequence of test simulations was carried out. In each simulation, we vary the value of the wall frictional coefficient µ pw and maintain other simulation parameters the same. The experiments were conducted using the parameter values shown in Table 1, except for the value of µ pw which varies from 0.1 to 0.9. The methodology was to, completely, fill the upper silo with 7,000 particles and let them to settle down till they come to rest. Then the outlet was opened and the particles flow down under the gravitational force into the lower silo. The stresses acting on each segment of the hopper's wall, during filling the lower silo for various values of µ pw , are shown in Figure 8.
Our results reveal that, after the filling process, granular systems with small wall fictional coefficients, in general, achieve static state faster than those with large values. For the cases of µ pw = 0.1 and 0.3, there are no stress fluctuation after t = 100 s which means that the two systems are stable and there is no stress variation after this time. However, the stress fluctuation continues after t = 100 s for the cases of and µ pw = 0.5 and µ pw = 0.9, which means that there is still a stress variation occurring beyond this time. Nevertheless, the amplitude of the stress fluctuation tends to be smaller as the wall frictional coefficient increases. This is because of the rapid dissipation of the kinetic energy of particles in direct contact with the hopper wall with high friction. At the end of the filling process, in general, we observe that the ultimate stress values acting on each segment decrease as the value of µ pw increases. This can be explained by the mechanism that the larger is µ pw , the more percentage of material weight the silo walls will carry. Figure 9 shows the mean stress acting on each segment of the hopper wall in the lower silo at the end of the filling process. Except the case when μ pw = 0.5, stresses exerted on segment one were the maximum while those exerted on segment five were the minimum. However, no considerable changes were observed on the stresses that acting on other segments. To investigate the influence of µ pw on the discharging process, another series of simulations was performed for discharging the particles from the lower silo this time. In these simulations, we also use the mechanical and environmental parameter values listed in Table 1. Different values of µ pw were considered in each simulation, namely, 0.1, 0.3, 0.5, and 0.9. During the process, the stresses acting on the hopper's segments were measured. At t = 0 s, the lower silo outlet opened and the particles started to flow out of the silo under gravity. As the particles flow out form the silo, the values of the stresses acting on each segment start to fluctuate and decrease gradually till the end of the discharging process. Figure 10 shows the stress traces for each segment along the hopper's wall of the lower silo during the discharging process. Each graph corresponds to a different value of µ pw while other simulation parameters are kept the same. The results show that, for different values of µ pw , the required time for the discharging of the silo is almost the same. Moreover, the fluctuation amplitude of the stress that acting on segment five is larger than that acting on the other segments.
In general, the graph shows that the stress fluctuation amplitude decreases as the value of µ pw increases.
According to our simulation outputs, the stress fluctuation profile, during the discharging process, is not considerably affected by the value of µ pw . However, in general, the fluctuation amplitude for larger values of µ pw tends to be smaller than that for small values of µ pw . On segment five the stress fluctuation trace has the largest amplitude compared to those on all the other segments. Therefore, the points of transition from vertical silo wall to the hopper wall bear the maximum value of the acting stresses during the discharging process.

Influence of internal friction
To investigate the effect of the internal friction coefficient µ pp on the stresses profile, we follow the same procedure as that in 5.2. We run a series of simulations in which all simulation parameters are fixed while µ pp is given different values each time. We consider four different values of µ pp , namely, 0, 0.1, 0.5, and 1. For the filling process, the outlet of the filled upper silo was opened and the particles started to fall down into the lower silo. Figure 11 shows the resulting stresses profiles on the segments of the hopper's wall for the four different values of µ pp . Results show that granular systems with small values of µ pp comes to rest faster than those with large values of µ pp . Moreover, at the end of the filling process and when the assembly eventually comes to rest at the base of the lower silo, stress profiles appear to have different traces. For the case of µ pp , a consistent stress fluctuation that has a simple harmonic motion-like trace appears for a slightly long period of time before the stress trace turns to a straight line. On the other hand, for other cases i.e. µ pp = 0.1, 0.5, and 1, no long consistent stress fluctuation has been observed and the stress trace turns to a straight line rapidly. This phenomenon can be explained as flows, for µ pp = 0, there is no internal friction between particles and the system has a fluid-like properties. In this case particles take more time to settle down and achieve a real steady state. Another observation, that can be obtained from Figure 11, is that the stress fluctuation amplitude for the systems with smaller value of µ pp is greater than that for the systems with larger values of µ pp . This is reasonable for the same reason that has been mentioned earlier. At the end of the filling process, the mean stress acting on each segment is plotted in Figure 12. For µ pp = 0, the ultimate stress profile is reasonably close to a straight line as the bulk material has a fluid-like property. It is also found that at the end of filling process, for pp = 0.5 and 1, stresses exerted on segment five were the maximum while those exerted on segment one were the minimum. In these cases, the transition's points of the silo system carry most of the load of the bulk material due to the large internal friction coefficients. Conversely, segment 1 bears the maximum stress for pp = 0 and 0.1 due to the low internal friction between particles. For the effect of the internal friction coefficient µ pp on the stresses during the discharging process, we performed a series of simulations using the same four different values of µ pp . Other parameters were preserved the same. Figure 13 shows the stresses profile developed on each segment of the hopper's wall during discharging of the lower silo for the four different values of µ pp .
According to our results, the discharging process takes more time for large values of µ pp as would be expected. Moreover, due to dilation of the particles at the outlet zone of the silo system during flow, segment one almost maintains the same fluctuation amplitude throughout the discharging process. For other segments, the fluctuation stress trace gradually decreases till it turns to a straight line at the end of the process. On the transition between the hopper wall and the vertical bin wall (segment 5), the stress appears to have the maximum fluctuation.

Concluding remarks
Two-dimensional soft particle discrete element simulations have been conducted to simulate the filling and discharging of granular particles from a silo system. The simulated filling and discharging processes have been carried out to study the variation in stress, which acts on the hopper wall, with time. The effects of wall and internal friction coefficients on the dynamic nature of stresses acting on hopper wall have also been investigated. From the results, it can be concluded that the stresses that act on the hopper wall have a temporary dynamic pulsation-like nature at the end of the filling process when the granular particles appear to reach a near static state. Accordingly, under the condition that the particle-wall contact force, in the normal direction, is modeled as visco-elastic (damped harmonic oscillator), the proposed soft particle discrete element model is capable of simulating the silo system. The influences of wall and internal friction coefficients on the oscillating nature have been investigated during filling and discharging processes. At the end of filling process, it can be concluded that the larger the µ pw is, the less the ultimate stress values acting on the hopper wall will be. Moreover, for smaller values of µ pp , the system has fluid-like properties and hence, the particles take more time to settle down and achieve a real steady state. Therefore, a consistent stress fluctuation appears for a slightly long period of time before the stress trace turns to a straight line. For the discharging process, the results show that for different values of µ pw , the required time for the discharging of the silo is almost the same. Furthermore, the fluctuation amplitude of stress acting on the upper part of the hopper wall is larger than that on the other parts. Also, it is found that the discharging process takes more time for large values of µ pp . Due to the particles' dilation at the zone of the silo's outlet, the lower part of the hopper wall maintains almost the same fluctuation amplitude throughout the discharging process. On the transition between vertical bin and hopper's wall, the stress appears to have the maximum fluctuation during the discharging process for all different values of μ pp .