Study of Static Pressures of Granular Materials in a Silo Using the Distinct Element Method t

The pressure inside a silo, especially the dynamic pressure during discharge is a very important factor in its structural design. However, numerical simulations of this problem have seldom been performed because of its complicated behaviors. In this report, the distinct element method (DEM) which is suitable for discontinuous media is used for the analysis of the static state of granular materials in silos, and the possibility of its use to analyze the silo problem is considered. Also, the characteristic distribution of pressures on walls is studied using the DEM simulation, and the influence of friction coefficients and particle arrangement is examined. The result of this study is that the pressure distributions obtained from the DEM simulation are similar to the theoretical or experimental ones. This result of static analysis is used for the initial condition of dynamic analysis, and it is expected to be used to simulate the com plicated behaviors of granular materials during discharge in silos.


Introduction
The materials stored in a silo such as coal, grain, and flour are usually treated as powders or granular materials. The behavior of such materials resemble that of fluids, but doesn't comply with the rules of fluid mechanics or those of solid mechanics. Thus, this complicated behavior in a silo was mainly studied experimentally 1 l.2), and the theory for its structural design is derived from experiments and observations or is based on the classical theory of Rankine. The design basis for a silo in most countries is founded on Janssen's formula for the static state, while for the dynamic case the incremental coefficient added to the static pressure is considered. In this way the behavior in a silo is understood experimentally but not theoretically. It is very difficult to simulate the behavior of the particles in a silo, although the inside pressure, especially the dynamic pressure during discharge, is a very important factor in its structural design.
On the other hand, the Distinct Element Method (DEM in· short) is a numerical method used to simulate the assembly of discrete rigid bodies such * 27F, Fukoku Seimei BLDG., 2-2-2 Uchisaiwai-cho, Chiyoda-ku Tokyo 100, Japan t This report was originally printed in]. Soc. Powder Technology, japan. 29, 673-86 (1992) in Japanese, before being translated into English with the permission of the editorial committee of the Soc. Powder Technology, Japan. as rock mass. The DEM is often used in the civil engineering field, for example the failure of rock slopes, the stability of rock caverns, the dynamic behavior of sands, and so on. For the problem of silos, Kiyama et al 3 ). applied it to a simple model, but the number of its particles was too small to describe the complicated behavior in a silo. Recently, the development of computers has made great strides, and massive calculations can now be performed much faster than 10 years ago. Therefore a large scale model composed of 10,000 circular elements can be calculated on the DEM. This report shows the analytical result for granular material in a silo using the DEM, and its possibility for use to simulate large numbers of elements is considered. Based on the analytical results, the relationship between the properties of granular material and the static pressures towards the walls in a silo is studied.

Distinct Element Method
The DEM is a numerical method proposed by Cundall in 1971 4 l, and is very useful for the simulation of the complicated behavior of discontinuous media. Its main advantage is that it can be used for the large dynamic deformation problem continuously from static state to dynamic failure. It is usually applied for the problem of rock mass failure and the analysis of particle motions.
The solution technique of the DEM is relatively simple. It is based on a few assumptions (see Fig. 1) and while the numerical integration of the equation of motion is performed, equilibrium conditions are considered in algorithms by which calculated unbalanced forces are reduced step by step. One of the difficulties is that the time increment for each calculation step must be defined so small that a very long execution time is required to analyze the behavior of materials in real time. In order to reduce the execution time, my program 5 l has been improved on the following two aspects. 1) The old program contains many rearrangements of variables and external file I/Os in order to reduce the necessary main memory capacity, but recently a large memory capacity can be used as a result of the extensional region of FORTRAN and the enhancement of computers themselves. Therefore, in the new program, the variables are directly set to their dimensions, and external file I/Os are concentrated. 2) To reduce the time of the searching process for the contacts between particles, adjacent particles are memorized between two boxing processes, which are also useful techniques of searching for contacts introduced by Cundall. Owing to the development of computers and the improvement of the program mentioned above, it is possible to solve a case involving 10,000 particles, while it was formerly limited to 1,000 elements. However, in order to deal with a real situation for powder technology the number of particles which must be handled would be at least 10 times higher.

Influence of friction coefficient
The dynamic pressure during discharge is usually several times higher than the static pressure, and is one of the main problems in a silo design. In this report, as a reasonable initial state of the discharge analysis, the static state in a silo is simulated using the DEM, and the static pressures acting towards the side walls and the bottom are examined from the viewpoint of the properties of the particles.
Silos are grouped into two classes as shallow and deep bins according to the ratio of their diameter to height and to the static angle of repose of the materials inside them. The pressure acting towards the side wall of shallow bins is relatively in good agreement with the pressure based on Rankine's theory, and that of deep silos can be calculated using Janssen's formula.
In this chapter, a model of randomly arranged particles (see Fig. 2) is simulated, and the friction coefficient is used as a parameter having several values. In order to compare the analytical results,

Analytical Model and Material Properties
The model used in this chapter is visualized as a 2-dimensional test model measuring 40 em in width and 50 em in height. Within this area 5000 particles of 6mm-diameter are arranged at random. Table 1 shows the material properties of the particles which are defined as constant except for the friction coefficient. The friction coefficient between two particles is defined as being the same as that between a particle and wall in this analysis, and 4 values, 0.3, 0.5, 0. 7, and 0.9, are used. The sequence of this analysis is as follows. First, the particles are arbitrarily arranged in the analytical area with small spaces between them. Second, the gravity force is applied to each particle, and the Table 1  calculation is continued until the motion of the particles is sufficiently converged. However, it is very difficult to control the motion of the particles perfectly by the algorithm of the DEM, therefore in practice the total calculation is performed up to 10000 steps with a time increment of 5.0 x I0-5 sec.

Results of Analysis and Discussion
The final state of each friction coefficient is shown in Fig. 2. In this figure circles show the particles themselves, and short lines indicate the contact forces between the particles and the wall. The filled circles are merely the marks at each 10 em interval.
When the friction coefficient is high the highest position of the particles is slightly higher. The difference is smaller than a 6 mm diameter particle between the friction coefficients 0.3 and 0.9. The direction of the contact force for the larger coefficient is not normal to the wall because of the greater tangential force than that for the smaller coefficient. The maximum contact force given on the top left of the figure is not significantly different, and it acts on the bottom. Therefore, it is assumed that the contact force acting towards the bottom is not seriously affected by the variation of the friction coefficient in this condition.
The contact forces acting towards the wall tend to decrease at higher locations of the particles but in the case of the adjacent forces the relationship is sometimes different. Therefore, in order to examine the distribution of the pressures towards the walls, the contact forces are integrated at each 4 to 5 em-width interval and the distributions are shown in Figs. 3 and 4. Fig. 3 shows the pressure dis- which shows the calculated pressure exerted by the weight of all the particles. The pressures towards the bottom are lower near the side walls, and this tendency is similar to that in practical or experimental conditions. As the pressure distribution towards the upper part of the side walls has a nearly triangular form, it is proportional to the depth and the dispersion is small. Contrariwise, the lower part shows a large variation of the pressures and the distribution has a form other than triangular. Each case is shown with its corresponding result and the difference between the cases of each friction coefficient is not large. The straight lines in Fig. 3 represent the Rankine's active pressure for the friction coefficient 0.3, and if the coefficient is larger, the Rankine's pressure must be smaller. However, the results obtained from the DEM analysis are larger than the Rankine's pressure, therefore the pressure calculated using Janssen's formula must be smaller than indicated by these straight lines. The forms of the pressure distributions resemble those of a deep bin, therefore the analytical results must be compared with the Janssen's. However, a greater difference between theory and analysis would be shown even if the result of the DEM is compared with the Janssen's.
The map of the contact forces is shown as short lines in Fig. 5. The forces describe a random spidery figure because of the concentration of the forces by the local arching phenomenon, and the stress concentration on the walls shown in Fig. 3 is also related to the arching. In the case of particles arranged arbitrarily, the stress distribution is strongly influenced by their arrangement. It is considered that when the size of the particles is relatively large with respect to the bin, the distribution of the forces in the bin is not homogeneous and the arrangement of the particles is closely related to the concentration of the stress. Also, the stress reduction at the comer of the model shows the influence of the friction between particles and walls. It can be said that the arrangement of the particles and the friction force at the walls are important factors which affect the stress distribution.

Influence of the arrangement of particles
The result of the last chapter indicates that the influence of the arrangement of the particles on the pressure is relatively large, and the partial concentration of the contact forces leads to an unbalanced stress distribution. In this chapter a slightly different arrangement of the particles from the previous analysis is examined, when the other conditions are kept constant. Regularly arranged models are also introduced to clearly demonstrate the influence of the arrangement.

Random Arrangement
The material properties and the conditions of the analysis except the arrangement of the particles are the same as in the case of the friction coefficient 0. 7 in Chap. 3. The arrangements before the calculation are changed and the vertical distance of the particles is set 0.01 mm and 0.1 mm higher in two cases, so that the appearances of the models are almost identical.
The lateral pressure distributions on the side walls are shown in Fig. 6 and the vertical pressure distributions on the bottom are shown in Fig. 7. The total distributions are not significantly different, but the difference between each case is as large as those of Figs. 3 and 4. The map of the contact forces does not appear so different (see Fig. 8), but there are some differences in the connectivities of the forces, and they are the cause of the variation of the pressure distribution. The maximum values of the contact forces are not equal, and that of Case 1 which is the result of the last chapter is close to the maximum value in Fig. 5 for the friction coefficient 0.5, and the form of map also resembles that of Fig. 5, but not the other two cases. From these results it can be observed that the pressures towards the walls are very sensitive to the arrangement of the particles, and when the arrangements of the particles are similar, the result of the analysis does not vary significantly even if the friction coefficients differ.

Regular Arrangement
It has been shown that the arrangement of the particles is a very important factor affecting the pressure distribution on the wall. Therefore, a regular arrangement model was analyzed in order to compare it with the random arrangement model. Here, regular means that the particles are horizontally and equidistantly arranged with no empty space. therefore in the case of shallow bins the result of the DEM analysis is similar to that of the traditional theory, but in the case of deep bins the D EM results in larger pressures that derived from experience. As mentioned above, it is assumed that the interaction between the particles is mainly slip when the friction coefficient is small, and when it is larger, the interaction gradually changes to include toppling. As in the traditional theory in which the interaction is assumed to be only slip, the result of the DEM is different from the theoretical one. On the other hand, the pressure towards the bottom is relatively constant with the exception of the reduction near the side walls, and the extent of this reduction depends on the friction coefficient. These phenomena indicate that the friction forces between the particles and the side walls are closely related to the reduction of the pressure near the lower corners. Otherwise, the forces between the particles may have arbitrary directions, but the direction of the forces acting on the wall is almost normal to the wall. Therefore the pressures near the lower corners are lower than at other parts of the wall. These two reasons are believed to be the explanation for the reduction of the pressures.
The distributions of the contact forces in Fig. 11 are relatively homogeneous and the maximum value is smaller than in the case of the random arrangement. This fact shows that the local concentration of the stress is caused by the localization of the arrangement or a vacancy in the rows. However, in real situations the diameter of the particles is much smaller, and the influence of the arrangement seems to be smaller.

Influence of the diameter of particles
The pressures in a silo calculated using the DEM analysis are fairly varied depending on the particle arrangement. However, the size of the particles in this report compared to the size of the silo is much larger than the particles in real silos, and this deviation must be of strong influence on the distribution of the pressures. Therefore, in this chapter a model of small particles having almost half the area of the previous models is used in order to evaluate the effect of the size of the particles. The diameter of the particles is taken as D = 4 mm, and the number of the particles as N = 10000. Therefore, the total mass of the model is not so different from that of the previous model in which D = 6 mm and N = 5000. The other analytical conditions are the same as for the previous models with the exception of the time increment which is defined as 2.0 x 10·5 sec (see Table 2).
The map of the contact forces in Fig. 12 shows the local concentration as in the cases of D = 6 mm,  in Fig. 13 show a well-balanced and smooth form compared to those in Fig. 6. The pressure itself is larger than that of Janssen's formula, but the distribution is similar to the other results. The form of the curve in Fig. 14 is uneven like that in Fig. 7, but the standard deviation for D = 4 mm in Table 3 is smaller than that for D = 6 mm. As mentioned above, when the diameter of the particles is smaller, the result becomes more homogeneous and well-balanced. However, in practical situations the diameter of the particles is much smaller and their number is huge. From such a point of view in order to analyze a practical problem using the DEM, a massive and high-speed program is required, and the diameter of each particle which is related to the size of the whole model must be reasonably determined or a regular arrangement must be properly used with the appropriate properties.

Conclusions
The behavior of granular materials in a silo, especially the pressure distribution on the walls, is studied using the DEM simulation. From this study it is concluded that: 1) The form of the distribution curve of the pressure towards the side walls is identical to that of the traditional theory. However, the value obtained by simulation is larger than the theoretical one, and one reason for this deviation is believed to be the fact that the rotation motion is not considered in the theory. 2) When the particles are arranged at random, the simulation is significantly influenced by the variations of the arrangement, and especially the size of the particles is so large with respect to the bin that a vacancy in the rows causes stress concentration. This is also one reason for the irregular form of the pressure distribution on the walls. 3) When the size of the particles is smaller, the pressure distribution becomes more homogeneous. However, to deal with the real size of particles is very difficult due to the limitations of the capability of computers.

222
4) The pressure on the walls decreases near the lower comer of the bin, and it is considered that the friction force between the particles and the wall causes this phenomenon. However, the effect of this phenomenon on the total behavior must be examined. As mentioned above, the static state in a silo can be simulated using the DEM, even with some remaining problems. Taking this static state as a reference, the simulation of the discharge can be analyzed, and the complicated behavior during the discharge could be explained using the DEM analysis.