# Effect of Friction Coupling on Discharge Velocity Profiles and Force Chain Distribution of Maize Particles in Silos

Authors: Yong Feng; Caihua Yu; Fan Pan
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

The evolution mechanism of discharge velocity profiles and force chain distribution of maize particles in silos was studied based on the interaction between internal and external rolling friction of particles. Through EDEM, the silo and maize grain models were established for unloading simulation, whose flow pattern was compared with the silo unloading test to verify the rationality of the simulation. By slice observation, we compared and analyzed the time evolution rules of particle mesoscopic parameters under different friction conditions. The results show that the larger the interparticle friction coefficient is, the longer the total discharge time is and the smaller the coefficient of rolling friction between particles, the earlier the particle flow from mass flow to funnel flow. For silos with the funnel, the reduction of interparticle friction will change the limit between the mass flow and the funnel flow, thus increasing the area of the funnel flow. When the coefficient of rolling friction increases, the vertical velocity and angular velocity of the particle near the silo middle increase. However, the effects of internal and external friction coupling on the vertical velocity of the side particle, the horizontal velocity of the whole particle, and the spatial distribution and probability distribution of the force chain are more significant.


---

## Body

## 1. Introduction

The frictional properties of coarse grains are much more complex than those of engineering materials or solid dry friction. In particular, the irregularity of grains has a great influence on the friction characteristics during the dynamic unloading process in silos. Moreover, the friction characteristics of particles directly or indirectly affect the flow pattern, thus affecting the storage and unloading safety of granary.Filling and discharge of grains from silos are very frequent in industrial and agricultural applications. The dynamic characteristics of granular flow [1–3] and velocity are of great significance to the safety of silos. Weinhart [4] used the discrete element method to capture particle velocity and stress in DEM simulation. Diego [5] modified the unloading rate formula through 2D simulation and explored the effect of particle size on unloading rate. Through quantitative analysis of discrete element simulation, Tian [6] concluded that the unloading flow state of silo is related to the loading height. Wan [7, 8] studied the effect of different curvature of the funnel on the unloading velocity and efficiency. Besides the external factors such as silo shape and diversion body, the importance of friction characteristics is increasingly found during studying the macroscopic changes of particle flow.R. Kobyłka et al. [9, 10] explored the evolution of the pressure fluctuation and friction characteristics by simulating the velocity and stress distribution. Yuan et al. [11–13] used particle image velocimetry to study the displacements of the soils around the pile. Yung and Xu [14] found that the rolling friction of a flexible ball moving on a rigid plane has a nonlinear relationship with the rolling speed. Different tests to detect particle rolling friction and conducted comparative verification of repose angle through the discrete element method were designed by Liu et al. [15–17], which provided basic parameter standards for numerical simulation. When conducting numerical simulation, Weizman [18] proposed to convert the rotation velocity into the translational velocity and determined the variation coefficient of the rolling friction of the rigid ball on the flexible surface. Balevičius [19–21] studied the influence of the shape of the hopper on the discharge flow of particles based on circular particle simulation and pea experiment and studied the influence of coarse granulation parameters on the particle flow state and obtained the best CG parameters. Silbert [22] found that the spatial distribution of particles in the particle heap depends on the particle friction coefficient. Goniva [23] found that the introduction of rolling friction into a single nozzle fluidized bed simulation could greatly improve the simulation results. Golshan [24] studied the velocity profiles and residence time distribution of nonspherical particles in silos by experiments and DEM. Zeng [25] combined the fluctuation of the unloading velocity with the disappearance in cycles of the internal force chain in silos. These studies correspond the macroscopic changes to the microscopic properties of particles. Yet, even in the monoball particle case, the variation of particle velocity and force chain distribution due to the frictional coupling effect in a silo is not well studied or understood.Based on the above analysis, it is of great significance to study the friction characteristics of irregular maize particles in silos during the dynamic unloading process. We use discrete element analysis software (EDEM) to study the dynamic process of mesoscopic variables, building the bridge macro-mesoscopic. The main work of this paper is to compare and analyze the temporal and spatial variation of the velocity field and force chain under different frictional coupling effects.


## 2. Discrete Element

### 2.1. Modeling

#### 2.1.1. Particle Model

In this paper, the simulated grain prototype is a variety of Demeiya No. 1 maize with a wide planting area in China. The three-dimensional geometric dimensions of the selected corn seeds are shown in Figure1, with the upper base width l, the lower base width L, the upper base thickness d, the lower base thickness D, and the height h.Figure 1
Contour size of corn granules.In the simulation, we adopted the multiball method [26] to approximate the shape of corn particles, and a total of six ball models were combined with each other, whose size and position are shown in Figure 2.Figure 2
Shape and size of the corn model.


#### 2.1.2. Silo Model

The small silo in the laboratory is formed by modeling and pouring with acrylic materials. Because of its good transparency, we can directly observe the movement of particles in the silo during filling and discharging. The top and bottom sides of a silo are square, with side lengthL = 120 mm and height H = 250 mm. The hopper of the silo is in the form of a truncated pyramid, the walls of which subtended an angle of 60° to the horizontal, and the square side of the funnel hole is D = 50 mm. The size of the test silo is shown in Figure 3(a). According to the size of the experimental silo, 3Dmax modeling software is used to draw the silo model diagram of the same size, and its model is shown in Figure 3(b).Figure 3
Dimensions of silos.
(a)(b)


### 2.2. Particle Contact Model

The moisture content of maize was 14.23% by the drying method. It is assumed that the stress, displacement, and velocity of the particles are determined by the difference of the tiny overlap between the particles and the contact body. Based on this hypothesis, the Hertz–Mindlin (no slip) contact mechanics model related to relative velocity was established in EDEM software to define particle contact. In this model, the collision force of particles at the contact point is decomposed into normal force and tangential force, as shown in Figure4, and both tangential and normal contact forces are simplified to parallel connection of the spring and damper.Figure 4
Particle contact model.According to Newton’s second law, the translational equations of particles are as follows:(1)midvidt=mig+∑j=1niFn+Fnd+Ft+Ftd,where mi is the mass of particle i,kg; vi is the velocity of particle i,m/s;t is the time, s; g is the acceleration of gravity,m/s2; Fn is the normal stress and Fnd is the particle normal damping force; and Ft is the shear stress and Ftd is the particle tangential damping force:(2)Fn=43E∗R∗1/2δn3/2,where E∗ is the equivalent modulus of elasticity, Pa; R∗ is the equivalent radius, m; and δn is the normal overlap, m. Also,(3)Fnd=−256lneln2e+π2Snm∗vnrel,Ft=−Snδt,where e is the recovery coefficient; Sn is the normal stiffness, N/m; m∗ represents the equivalent mass, kg; vnrel is the relative normal velocity, m/s; and δt is the tangent overlap, m.The shear stress follows Coulomb’s law of friction, and the maximum shear stressFt,max is as follows:(4)Ft,max=δFn+c.σ is the coefficient of static friction. In EDEM simulation, if the absolute value of Ft is greater than Ft,max, take Ft=Ft,max.(5)Ftd=−256lneln2e+π2Stm∗vnrel,where St is the tangential stiffness, N/m and vtrel is relative tangential velocity, m/s. Particles are also affected by tangential torque Tt and rolling friction torque Tr, and the motion velocity is as follows:(6)Iidwidt=∑j=1niTt+Tr,Tt=RiFt+Ftd,where Ii is the moment of inertia of particle i, kg⋅m2 and Ri is the distance from the center of mass of particle i to the contact point, m.In the above particle collision process, the antirolling friction torque [27–29] is(7)Tr=−μrRrFnwrelwrel,where wi is the angular velocity, rad/s; Tt is the tangential torque, N⋅m; Tr is the moment of rolling friction, N⋅m; Rr is the effective rolling contact radius, m; wrel is the relative angular velocity between two particles, rad/s; and μr is the coefficient of rolling friction.


### 2.3. Parameters of Material Performance

The performance indicators of the materials are listed in Table1. Poisson’s ratio of maize particles was used as the reference value according to the results described by Tao et al. [30]. According to the standard ASAE S 368.4 (2006) [31], the shear modulus (G) of maize was determined by shear test analysis on a single grain. The density of corn was calculated according to the mass and volume of corn, and the silo density was determined according to the acrylic material standard.Table 1
Performance index values of materials.
ParametersValuePoisson’s ratio of the maize0.4Poisson’s ratio of the silo0.28Shear modulus of the maize (Pa)1.37 × 108Shear modulus of the silo (Pa)8.10 × 1010Density of the maize (kg/m3)1154Density of the silo (kg/m3)7850The interaction property parameters of materials are listed in Table2. The static friction coefficient means that the maximum static friction force is proportional to the positive pressure between contact surfaces, independent of contact surface size, and its direction is opposite to the direction of the movement trend. The recovery factor is usually defined as the absolute ratio of the velocity after the collision to the velocity before the collision. Gonzalez [26] designed a device for measuring the static friction coefficient between corn particles and between particles and contact materials and a device for measuring the recovery coefficient. The correlation coefficient was measured, and the measured data were simulated and verified. On this basis, we use the method to measure the static friction coefficient and recovery coefficient between corn particles and acrylic materials and introduce them into our simulation.Table 2
Interaction property parameters.
ParametersInteractionParticleSiloCoefficient of restitution0.3620.613Coefficient of static friction0.3490.461Coefficient of rolling friction0.06070.0311


## 2.1. Modeling

### 2.1.1. Particle Model

In this paper, the simulated grain prototype is a variety of Demeiya No. 1 maize with a wide planting area in China. The three-dimensional geometric dimensions of the selected corn seeds are shown in Figure1, with the upper base width l, the lower base width L, the upper base thickness d, the lower base thickness D, and the height h.Figure 1
Contour size of corn granules.In the simulation, we adopted the multiball method [26] to approximate the shape of corn particles, and a total of six ball models were combined with each other, whose size and position are shown in Figure 2.Figure 2
Shape and size of the corn model.


### 2.1.2. Silo Model

The small silo in the laboratory is formed by modeling and pouring with acrylic materials. Because of its good transparency, we can directly observe the movement of particles in the silo during filling and discharging. The top and bottom sides of a silo are square, with side lengthL = 120 mm and height H = 250 mm. The hopper of the silo is in the form of a truncated pyramid, the walls of which subtended an angle of 60° to the horizontal, and the square side of the funnel hole is D = 50 mm. The size of the test silo is shown in Figure 3(a). According to the size of the experimental silo, 3Dmax modeling software is used to draw the silo model diagram of the same size, and its model is shown in Figure 3(b).Figure 3
Dimensions of silos.
(a)(b)


## 2.1.1. Particle Model

In this paper, the simulated grain prototype is a variety of Demeiya No. 1 maize with a wide planting area in China. The three-dimensional geometric dimensions of the selected corn seeds are shown in Figure1, with the upper base width l, the lower base width L, the upper base thickness d, the lower base thickness D, and the height h.Figure 1
Contour size of corn granules.In the simulation, we adopted the multiball method [26] to approximate the shape of corn particles, and a total of six ball models were combined with each other, whose size and position are shown in Figure 2.Figure 2
Shape and size of the corn model.


## 2.1.2. Silo Model

The small silo in the laboratory is formed by modeling and pouring with acrylic materials. Because of its good transparency, we can directly observe the movement of particles in the silo during filling and discharging. The top and bottom sides of a silo are square, with side lengthL = 120 mm and height H = 250 mm. The hopper of the silo is in the form of a truncated pyramid, the walls of which subtended an angle of 60° to the horizontal, and the square side of the funnel hole is D = 50 mm. The size of the test silo is shown in Figure 3(a). According to the size of the experimental silo, 3Dmax modeling software is used to draw the silo model diagram of the same size, and its model is shown in Figure 3(b).Figure 3
Dimensions of silos.
(a)(b)


## 2.2. Particle Contact Model

The moisture content of maize was 14.23% by the drying method. It is assumed that the stress, displacement, and velocity of the particles are determined by the difference of the tiny overlap between the particles and the contact body. Based on this hypothesis, the Hertz–Mindlin (no slip) contact mechanics model related to relative velocity was established in EDEM software to define particle contact. In this model, the collision force of particles at the contact point is decomposed into normal force and tangential force, as shown in Figure4, and both tangential and normal contact forces are simplified to parallel connection of the spring and damper.Figure 4
Particle contact model.According to Newton’s second law, the translational equations of particles are as follows:(1)midvidt=mig+∑j=1niFn+Fnd+Ft+Ftd,where mi is the mass of particle i,kg; vi is the velocity of particle i,m/s;t is the time, s; g is the acceleration of gravity,m/s2; Fn is the normal stress and Fnd is the particle normal damping force; and Ft is the shear stress and Ftd is the particle tangential damping force:(2)Fn=43E∗R∗1/2δn3/2,where E∗ is the equivalent modulus of elasticity, Pa; R∗ is the equivalent radius, m; and δn is the normal overlap, m. Also,(3)Fnd=−256lneln2e+π2Snm∗vnrel,Ft=−Snδt,where e is the recovery coefficient; Sn is the normal stiffness, N/m; m∗ represents the equivalent mass, kg; vnrel is the relative normal velocity, m/s; and δt is the tangent overlap, m.The shear stress follows Coulomb’s law of friction, and the maximum shear stressFt,max is as follows:(4)Ft,max=δFn+c.σ is the coefficient of static friction. In EDEM simulation, if the absolute value of Ft is greater than Ft,max, take Ft=Ft,max.(5)Ftd=−256lneln2e+π2Stm∗vnrel,where St is the tangential stiffness, N/m and vtrel is relative tangential velocity, m/s. Particles are also affected by tangential torque Tt and rolling friction torque Tr, and the motion velocity is as follows:(6)Iidwidt=∑j=1niTt+Tr,Tt=RiFt+Ftd,where Ii is the moment of inertia of particle i, kg⋅m2 and Ri is the distance from the center of mass of particle i to the contact point, m.In the above particle collision process, the antirolling friction torque [27–29] is(7)Tr=−μrRrFnwrelwrel,where wi is the angular velocity, rad/s; Tt is the tangential torque, N⋅m; Tr is the moment of rolling friction, N⋅m; Rr is the effective rolling contact radius, m; wrel is the relative angular velocity between two particles, rad/s; and μr is the coefficient of rolling friction.


## 2.3. Parameters of Material Performance

The performance indicators of the materials are listed in Table1. Poisson’s ratio of maize particles was used as the reference value according to the results described by Tao et al. [30]. According to the standard ASAE S 368.4 (2006) [31], the shear modulus (G) of maize was determined by shear test analysis on a single grain. The density of corn was calculated according to the mass and volume of corn, and the silo density was determined according to the acrylic material standard.Table 1
Performance index values of materials.
ParametersValuePoisson’s ratio of the maize0.4Poisson’s ratio of the silo0.28Shear modulus of the maize (Pa)1.37 × 108Shear modulus of the silo (Pa)8.10 × 1010Density of the maize (kg/m3)1154Density of the silo (kg/m3)7850The interaction property parameters of materials are listed in Table2. The static friction coefficient means that the maximum static friction force is proportional to the positive pressure between contact surfaces, independent of contact surface size, and its direction is opposite to the direction of the movement trend. The recovery factor is usually defined as the absolute ratio of the velocity after the collision to the velocity before the collision. Gonzalez [26] designed a device for measuring the static friction coefficient between corn particles and between particles and contact materials and a device for measuring the recovery coefficient. The correlation coefficient was measured, and the measured data were simulated and verified. On this basis, we use the method to measure the static friction coefficient and recovery coefficient between corn particles and acrylic materials and introduce them into our simulation.Table 2
Interaction property parameters.
ParametersInteractionParticleSiloCoefficient of restitution0.3620.613Coefficient of static friction0.3490.461Coefficient of rolling friction0.06070.0311


## 3. Physical Test

### 3.1. Test Process

The experimental materials were mainly yellow corn kernels, and a small number of black corn kernels were added as the spatial compartment so as to clearly observe the flow state of the particles in the process of unloading. The particles are loaded with the central charge method. During the charging process, due to the falling speed of particles, violent collisions will occur between particles of different colors and the particles of different colors will mix slightly. Therefore, the color boundary is not completely horizontal.After loading, let the full load model sit for two hours. Then, open the digital camera and draw out the baffle at the funnel hole to make the particles unload naturally, recording the whole process of unloading. The test flow is shown in Figure5.Figure 5
Flow chart of the test.In order to reduce the accidental error of the experiment, we repeated the experiment many times (5 times).


### 3.2. Authenticity Verification of Simulation

The material performance parameters listed in 2.3 were input into EDEM 2018 software as basic parameters to conduct loading and unloading simulation of corn granules. Five parallel repeated experiments were performed, and the average results were compared with the average results as follows.In terms of loading, the test loading height in the silo is 245.6 mm, and the simulated loading height is 250 mm. From the perspective of unloading time, the average time from the bottom plate extraction to the full discharge of particles from the funnel is 6.3 s in the test, while it of simulation is 6.6 s, and the simulation error is 4.7%.By comparing the flow behavior of particles in the experiment and simulation (Figure6), it can be observed that the height of the free surface decreases with the progress of unloading. In the picture at 0.25 T, the layered staining particles showed an obvious V-shape at the lower part of the silo, and the particle aggregation movement was obvious, while the upper layer was relatively flat and even. When the height of the stock bed decreases with time, the flow pipe that diffuses to the surface of the stock bed can be observed in the upper part of the silo. At 0.75 T, particle retention can be observed in both simulated and experimental silos and hopper joints, i.e., dead material area. The veracity of the simulation was verified by comparing the flow behavior of dye tracer particles in the experiment and simulation.Figure 6
Comparison of experimental and simulated flow patterns: (a) 0T, (b) 0.25T, (c) 0.5T, and (d) 0.75T.


## 3.1. Test Process

The experimental materials were mainly yellow corn kernels, and a small number of black corn kernels were added as the spatial compartment so as to clearly observe the flow state of the particles in the process of unloading. The particles are loaded with the central charge method. During the charging process, due to the falling speed of particles, violent collisions will occur between particles of different colors and the particles of different colors will mix slightly. Therefore, the color boundary is not completely horizontal.After loading, let the full load model sit for two hours. Then, open the digital camera and draw out the baffle at the funnel hole to make the particles unload naturally, recording the whole process of unloading. The test flow is shown in Figure5.Figure 5
Flow chart of the test.In order to reduce the accidental error of the experiment, we repeated the experiment many times (5 times).


## 3.2. Authenticity Verification of Simulation

The material performance parameters listed in 2.3 were input into EDEM 2018 software as basic parameters to conduct loading and unloading simulation of corn granules. Five parallel repeated experiments were performed, and the average results were compared with the average results as follows.In terms of loading, the test loading height in the silo is 245.6 mm, and the simulated loading height is 250 mm. From the perspective of unloading time, the average time from the bottom plate extraction to the full discharge of particles from the funnel is 6.3 s in the test, while it of simulation is 6.6 s, and the simulation error is 4.7%.By comparing the flow behavior of particles in the experiment and simulation (Figure6), it can be observed that the height of the free surface decreases with the progress of unloading. In the picture at 0.25 T, the layered staining particles showed an obvious V-shape at the lower part of the silo, and the particle aggregation movement was obvious, while the upper layer was relatively flat and even. When the height of the stock bed decreases with time, the flow pipe that diffuses to the surface of the stock bed can be observed in the upper part of the silo. At 0.75 T, particle retention can be observed in both simulated and experimental silos and hopper joints, i.e., dead material area. The veracity of the simulation was verified by comparing the flow behavior of dye tracer particles in the experiment and simulation.Figure 6
Comparison of experimental and simulated flow patterns: (a) 0T, (b) 0.25T, (c) 0.5T, and (d) 0.75T.


## 4. Results

### 4.1. Macroscopic Phenomenon Analysis

With the friction coefficient between the silo wall and the particles unchanged (a nonzero constant value of 0.0311), we changed the rolling friction coefficient between the particles. We found that the value ofμ not only affected the flow pattern of the particles in the silo but also affected the time step length of unloading.Increasing interparticle friction (μ) affects the final time step (T) in two ways. First, it will produce a small volume density value, which indirectly leads to the increase in T; in addition, Balevicius [20] et al. found that the rolling friction between particles would affect the constant C of the Beverloo [32] equation, which affects the unloading rate. That is, with the increase in friction, the unloading rate decreases and the unloading time step increases.In addition, by comparing the flow state of particles under the three conditions (Figure7), it can be seen that the smaller the friction coefficient is, the freer the movement of particles is, i.e., the more violent the mixing between different color particles is. The smaller the coefficient of friction, the earlier the black particle belt bends, i.e., the earlier the particle flow shows the characteristics of funnel flow. Especially, at 2 s and 3.5 s listed in Figure 7, the first layer of black particles with low friction coefficient has begun to bend, while the larger the friction coefficient is, the less obvious the bending is. After entering 5 s, the particles have basically entered the tubular flow state. Therefore, for silos with the funnel, the reduction of interparticle friction will change the limit between mass flow and funnel flow, thus increasing the area producing funnel flow. In addition, this explanation would be more valid for real silos in which the relationship between the dimensions of the silo and the particles stored was much bigger.Figure 7
Comparison of particle flow patterns under three rolling friction coefficients.


### 4.2. The Change of Velocity Field in Space

#### 4.2.1. Slicing Parameters

A series of 7-1-10 detection grids inx-y-z direction was established at the center of the silo body to detect the flow state of particles and the variation of their mesoscopic parameters in each spatial orientation inside the silo. The specific grid division and number between Figure 8.Figure 8
Specific division of the grids.


#### 4.2.2. Vertical Velocity

During the whole unloading process, the particle velocity profile at the height ofz in the silo (corresponding to Z = 1, 3, 5, 7, and 9 in the figure) graphically represents the average vertical velocity of all particles. Figure 9 shows an example of radial distribution under three rolling friction coefficients and provides the relevant velocity distribution, which is intended to highlight the difference caused by the modification of the friction coefficient.Figure 9
Contrast diagram of the vertical velocity of particles.In the lower middle of the silo, the vertical velocity presents a relatively obvious change rule: with the increase in the coefficient of rolling friction, the vertical velocity of the central particle increases, and the closer it is to the bottom, the more obvious this phenomenon is. WhenZ = 12 mm, as the friction coefficient increases, the corresponding central vertical velocity is 0.103 m/s, 0.114 m/s, and 0.123 m/s, which increases by 10.6% and 7.9%, respectively. When Z = 62 mm, the center vertical velocity under three friction coefficients is 0.0532 m/s, 0.0591 m/s, and 0.0637 m/s, which increases by 11.1% and 7.8%, respectively. When Z = 112 mm, the corresponding central vertical velocity is 0.0423 m/s, 0.00434 m/s, and 0.0437 m/s, which increases by 2.6% and 0.7%, respectively, corresponding to the characteristics of the fast flow of the central particle in Figure 7. However, the correlation between vertical velocity and friction coefficient of particles in the upper part of silo is not obvious. Due to the way of central unloading, the extrusion pressure on the particles within the unloading port and its upper vertical range is small. In this case, the more fully the particles rub, the faster the particles roll down. The vertical velocity of particles on the wall is obviously lower than that near the central at all heights. For the particles in the upper part of the silo, the friction coefficient effect is not obvious due to the faster loss rate.


#### 4.2.3. Horizontal Velocity

The longitudinal variations (with different silo heights) and the transverse fluctuations (under different friction coefficients) of the horizontal velocity of particles are compared in Figure10. For the same rolling friction coefficient, the fluctuation of horizontal velocity decreases with the increase in the height of the silo. The horizontal velocity of particles at the bottom fluctuates the most and presents the obvious characteristics of funnel flow, while the particles in the middle and upper part (7 and 9 points) flow as a whole and hardly move horizontally.Figure 10
Contrast diagram of the horizontal velocity of particles.By comparing the horizontal velocity of particles under different friction conditions, it is found that the effect is more obvious for particles near the bottom of the warehouse. We take the horizontal velocity whenZ = 1 as the specific research object and calculate the standard value of the horizontal velocity variation with time in the process of grain unloading. When μp=0.001 and μp/μw=0.03,σ=0.0073; when μp=0.0607 and μp/μw=2, σ=0.0207; when μp=0.102 and μp/μw=3.3, σ=0.0122. The violent degree of horizontal fluctuation of silo bottom particles is no longer a linear change with the rolling friction coefficient between particles but is related to the coupling effect of internal and external friction. When the coefficient of internal and external friction is relatively high or low, the fracture zone appears near the wall surface, leading to the formation of the boundary layer [33]. Among them, the movement of particles near the wall is slower than that in the central region, which cannot promote the rolling of particles. Only within the appropriate range (the ratio is 1.8∼2.1) can the rolling friction effect of particles be best achieved.


#### 4.2.4. Angular Velocity

The data of 1-1-1 (3, 5, 7, 9) and 4-1-1 (3, 5, 7, 9) in the section were detected in the grid, and two sets of angular velocity comparisons were made under different friction coefficients. The comparison of angular velocity of particles at the wall of the silo is shown in Figure11, and the comparison at the center is shown in Figure 12. Under the same friction condition, obviously, the angular velocity at different heights fluctuates around its mean value (C value), and the lower the heights, the greater the C value. By comparing the angular velocity of particles in the central area and the wall, it can be clearly obtained that the floating average of the angular velocity in the central area at the same height is larger, especially when Z = 1, the effect is most significant. When μp=0.001, the C value near the wall is 4.0, and the C value at the center is 6.8, which is 70% higher than that near the wall. When μp=0.06007, the C value near the wall is 2.3, and the C value at the center is 5.8, which is 152% higher than that near the wall. When μp=0.102, the C value near the wall is 2.0, and the C value at the center is 4.7, 135% higher than that. It can be seen from the comparison that although the angular velocity of particles in the central region under various friction coefficients is increased somewhat than that near the silo wall, there is still a difference in the amplitude of growth. When the ratio of internal and external friction coefficients is 2, the amplitude of growth is the largest. The above effect confirms the conclusion that the particles flow in the upper part of the silo as mass flow turns into funnel flow in the lower part of the silo [4], which is also consistent with the effect of internal and external friction coupling on horizontal velocity in 4.2.3 above.Figure 11
Comparison of angular velocity near the wall of the silo.Figure 12
Comparison of angular velocity at the center of the silo.Under different friction conditions, both the particles near the wall and in the central region show the phenomenon that the angular velocity decreases gradually with the increase in the rolling friction. The data ofZ = 1 were compared. When μp=0.001, C1 = 4.0 and C01 = 6.8; when μp=0.06007, C1 = 2.3 and C01 = 5.8; when μp=0.102, C = 2.0 and C01 = 4.7. With the increase in μp, the angular velocity of particles near the wall was 4.0, 2.3, and 2.0, respectively, which decreased by 42.5% and 13%, respectively, and the fluctuation range that its decrease is bigger. The angular velocity of particles at the center was 6.8, 5.8, and 4.7, respectively, which decreased by 14.7% and 19%, respectively, and the reduced fluctuation range is small. This phenomenon is consistent with the characteristics of the mixed flow of particles [34].


## 4.1. Macroscopic Phenomenon Analysis

With the friction coefficient between the silo wall and the particles unchanged (a nonzero constant value of 0.0311), we changed the rolling friction coefficient between the particles. We found that the value ofμ not only affected the flow pattern of the particles in the silo but also affected the time step length of unloading.Increasing interparticle friction (μ) affects the final time step (T) in two ways. First, it will produce a small volume density value, which indirectly leads to the increase in T; in addition, Balevicius [20] et al. found that the rolling friction between particles would affect the constant C of the Beverloo [32] equation, which affects the unloading rate. That is, with the increase in friction, the unloading rate decreases and the unloading time step increases.In addition, by comparing the flow state of particles under the three conditions (Figure7), it can be seen that the smaller the friction coefficient is, the freer the movement of particles is, i.e., the more violent the mixing between different color particles is. The smaller the coefficient of friction, the earlier the black particle belt bends, i.e., the earlier the particle flow shows the characteristics of funnel flow. Especially, at 2 s and 3.5 s listed in Figure 7, the first layer of black particles with low friction coefficient has begun to bend, while the larger the friction coefficient is, the less obvious the bending is. After entering 5 s, the particles have basically entered the tubular flow state. Therefore, for silos with the funnel, the reduction of interparticle friction will change the limit between mass flow and funnel flow, thus increasing the area producing funnel flow. In addition, this explanation would be more valid for real silos in which the relationship between the dimensions of the silo and the particles stored was much bigger.Figure 7
Comparison of particle flow patterns under three rolling friction coefficients.


## 4.2. The Change of Velocity Field in Space

### 4.2.1. Slicing Parameters

A series of 7-1-10 detection grids inx-y-z direction was established at the center of the silo body to detect the flow state of particles and the variation of their mesoscopic parameters in each spatial orientation inside the silo. The specific grid division and number between Figure 8.Figure 8
Specific division of the grids.


### 4.2.2. Vertical Velocity

During the whole unloading process, the particle velocity profile at the height ofz in the silo (corresponding to Z = 1, 3, 5, 7, and 9 in the figure) graphically represents the average vertical velocity of all particles. Figure 9 shows an example of radial distribution under three rolling friction coefficients and provides the relevant velocity distribution, which is intended to highlight the difference caused by the modification of the friction coefficient.Figure 9
Contrast diagram of the vertical velocity of particles.In the lower middle of the silo, the vertical velocity presents a relatively obvious change rule: with the increase in the coefficient of rolling friction, the vertical velocity of the central particle increases, and the closer it is to the bottom, the more obvious this phenomenon is. WhenZ = 12 mm, as the friction coefficient increases, the corresponding central vertical velocity is 0.103 m/s, 0.114 m/s, and 0.123 m/s, which increases by 10.6% and 7.9%, respectively. When Z = 62 mm, the center vertical velocity under three friction coefficients is 0.0532 m/s, 0.0591 m/s, and 0.0637 m/s, which increases by 11.1% and 7.8%, respectively. When Z = 112 mm, the corresponding central vertical velocity is 0.0423 m/s, 0.00434 m/s, and 0.0437 m/s, which increases by 2.6% and 0.7%, respectively, corresponding to the characteristics of the fast flow of the central particle in Figure 7. However, the correlation between vertical velocity and friction coefficient of particles in the upper part of silo is not obvious. Due to the way of central unloading, the extrusion pressure on the particles within the unloading port and its upper vertical range is small. In this case, the more fully the particles rub, the faster the particles roll down. The vertical velocity of particles on the wall is obviously lower than that near the central at all heights. For the particles in the upper part of the silo, the friction coefficient effect is not obvious due to the faster loss rate.


### 4.2.3. Horizontal Velocity

The longitudinal variations (with different silo heights) and the transverse fluctuations (under different friction coefficients) of the horizontal velocity of particles are compared in Figure10. For the same rolling friction coefficient, the fluctuation of horizontal velocity decreases with the increase in the height of the silo. The horizontal velocity of particles at the bottom fluctuates the most and presents the obvious characteristics of funnel flow, while the particles in the middle and upper part (7 and 9 points) flow as a whole and hardly move horizontally.Figure 10
Contrast diagram of the horizontal velocity of particles.By comparing the horizontal velocity of particles under different friction conditions, it is found that the effect is more obvious for particles near the bottom of the warehouse. We take the horizontal velocity whenZ = 1 as the specific research object and calculate the standard value of the horizontal velocity variation with time in the process of grain unloading. When μp=0.001 and μp/μw=0.03,σ=0.0073; when μp=0.0607 and μp/μw=2, σ=0.0207; when μp=0.102 and μp/μw=3.3, σ=0.0122. The violent degree of horizontal fluctuation of silo bottom particles is no longer a linear change with the rolling friction coefficient between particles but is related to the coupling effect of internal and external friction. When the coefficient of internal and external friction is relatively high or low, the fracture zone appears near the wall surface, leading to the formation of the boundary layer [33]. Among them, the movement of particles near the wall is slower than that in the central region, which cannot promote the rolling of particles. Only within the appropriate range (the ratio is 1.8∼2.1) can the rolling friction effect of particles be best achieved.


### 4.2.4. Angular Velocity

The data of 1-1-1 (3, 5, 7, 9) and 4-1-1 (3, 5, 7, 9) in the section were detected in the grid, and two sets of angular velocity comparisons were made under different friction coefficients. The comparison of angular velocity of particles at the wall of the silo is shown in Figure11, and the comparison at the center is shown in Figure 12. Under the same friction condition, obviously, the angular velocity at different heights fluctuates around its mean value (C value), and the lower the heights, the greater the C value. By comparing the angular velocity of particles in the central area and the wall, it can be clearly obtained that the floating average of the angular velocity in the central area at the same height is larger, especially when Z = 1, the effect is most significant. When μp=0.001, the C value near the wall is 4.0, and the C value at the center is 6.8, which is 70% higher than that near the wall. When μp=0.06007, the C value near the wall is 2.3, and the C value at the center is 5.8, which is 152% higher than that near the wall. When μp=0.102, the C value near the wall is 2.0, and the C value at the center is 4.7, 135% higher than that. It can be seen from the comparison that although the angular velocity of particles in the central region under various friction coefficients is increased somewhat than that near the silo wall, there is still a difference in the amplitude of growth. When the ratio of internal and external friction coefficients is 2, the amplitude of growth is the largest. The above effect confirms the conclusion that the particles flow in the upper part of the silo as mass flow turns into funnel flow in the lower part of the silo [4], which is also consistent with the effect of internal and external friction coupling on horizontal velocity in 4.2.3 above.Figure 11
Comparison of angular velocity near the wall of the silo.Figure 12
Comparison of angular velocity at the center of the silo.Under different friction conditions, both the particles near the wall and in the central region show the phenomenon that the angular velocity decreases gradually with the increase in the rolling friction. The data ofZ = 1 were compared. When μp=0.001, C1 = 4.0 and C01 = 6.8; when μp=0.06007, C1 = 2.3 and C01 = 5.8; when μp=0.102, C = 2.0 and C01 = 4.7. With the increase in μp, the angular velocity of particles near the wall was 4.0, 2.3, and 2.0, respectively, which decreased by 42.5% and 13%, respectively, and the fluctuation range that its decrease is bigger. The angular velocity of particles at the center was 6.8, 5.8, and 4.7, respectively, which decreased by 14.7% and 19%, respectively, and the reduced fluctuation range is small. This phenomenon is consistent with the characteristics of the mixed flow of particles [34].


## 4.2.1. Slicing Parameters

A series of 7-1-10 detection grids inx-y-z direction was established at the center of the silo body to detect the flow state of particles and the variation of their mesoscopic parameters in each spatial orientation inside the silo. The specific grid division and number between Figure 8.Figure 8
Specific division of the grids.


## 4.2.2. Vertical Velocity

During the whole unloading process, the particle velocity profile at the height ofz in the silo (corresponding to Z = 1, 3, 5, 7, and 9 in the figure) graphically represents the average vertical velocity of all particles. Figure 9 shows an example of radial distribution under three rolling friction coefficients and provides the relevant velocity distribution, which is intended to highlight the difference caused by the modification of the friction coefficient.Figure 9
Contrast diagram of the vertical velocity of particles.In the lower middle of the silo, the vertical velocity presents a relatively obvious change rule: with the increase in the coefficient of rolling friction, the vertical velocity of the central particle increases, and the closer it is to the bottom, the more obvious this phenomenon is. WhenZ = 12 mm, as the friction coefficient increases, the corresponding central vertical velocity is 0.103 m/s, 0.114 m/s, and 0.123 m/s, which increases by 10.6% and 7.9%, respectively. When Z = 62 mm, the center vertical velocity under three friction coefficients is 0.0532 m/s, 0.0591 m/s, and 0.0637 m/s, which increases by 11.1% and 7.8%, respectively. When Z = 112 mm, the corresponding central vertical velocity is 0.0423 m/s, 0.00434 m/s, and 0.0437 m/s, which increases by 2.6% and 0.7%, respectively, corresponding to the characteristics of the fast flow of the central particle in Figure 7. However, the correlation between vertical velocity and friction coefficient of particles in the upper part of silo is not obvious. Due to the way of central unloading, the extrusion pressure on the particles within the unloading port and its upper vertical range is small. In this case, the more fully the particles rub, the faster the particles roll down. The vertical velocity of particles on the wall is obviously lower than that near the central at all heights. For the particles in the upper part of the silo, the friction coefficient effect is not obvious due to the faster loss rate.


## 4.2.3. Horizontal Velocity

The longitudinal variations (with different silo heights) and the transverse fluctuations (under different friction coefficients) of the horizontal velocity of particles are compared in Figure10. For the same rolling friction coefficient, the fluctuation of horizontal velocity decreases with the increase in the height of the silo. The horizontal velocity of particles at the bottom fluctuates the most and presents the obvious characteristics of funnel flow, while the particles in the middle and upper part (7 and 9 points) flow as a whole and hardly move horizontally.Figure 10
Contrast diagram of the horizontal velocity of particles.By comparing the horizontal velocity of particles under different friction conditions, it is found that the effect is more obvious for particles near the bottom of the warehouse. We take the horizontal velocity whenZ = 1 as the specific research object and calculate the standard value of the horizontal velocity variation with time in the process of grain unloading. When μp=0.001 and μp/μw=0.03,σ=0.0073; when μp=0.0607 and μp/μw=2, σ=0.0207; when μp=0.102 and μp/μw=3.3, σ=0.0122. The violent degree of horizontal fluctuation of silo bottom particles is no longer a linear change with the rolling friction coefficient between particles but is related to the coupling effect of internal and external friction. When the coefficient of internal and external friction is relatively high or low, the fracture zone appears near the wall surface, leading to the formation of the boundary layer [33]. Among them, the movement of particles near the wall is slower than that in the central region, which cannot promote the rolling of particles. Only within the appropriate range (the ratio is 1.8∼2.1) can the rolling friction effect of particles be best achieved.


## 4.2.4. Angular Velocity

The data of 1-1-1 (3, 5, 7, 9) and 4-1-1 (3, 5, 7, 9) in the section were detected in the grid, and two sets of angular velocity comparisons were made under different friction coefficients. The comparison of angular velocity of particles at the wall of the silo is shown in Figure11, and the comparison at the center is shown in Figure 12. Under the same friction condition, obviously, the angular velocity at different heights fluctuates around its mean value (C value), and the lower the heights, the greater the C value. By comparing the angular velocity of particles in the central area and the wall, it can be clearly obtained that the floating average of the angular velocity in the central area at the same height is larger, especially when Z = 1, the effect is most significant. When μp=0.001, the C value near the wall is 4.0, and the C value at the center is 6.8, which is 70% higher than that near the wall. When μp=0.06007, the C value near the wall is 2.3, and the C value at the center is 5.8, which is 152% higher than that near the wall. When μp=0.102, the C value near the wall is 2.0, and the C value at the center is 4.7, 135% higher than that. It can be seen from the comparison that although the angular velocity of particles in the central region under various friction coefficients is increased somewhat than that near the silo wall, there is still a difference in the amplitude of growth. When the ratio of internal and external friction coefficients is 2, the amplitude of growth is the largest. The above effect confirms the conclusion that the particles flow in the upper part of the silo as mass flow turns into funnel flow in the lower part of the silo [4], which is also consistent with the effect of internal and external friction coupling on horizontal velocity in 4.2.3 above.Figure 11
Comparison of angular velocity near the wall of the silo.Figure 12
Comparison of angular velocity at the center of the silo.Under different friction conditions, both the particles near the wall and in the central region show the phenomenon that the angular velocity decreases gradually with the increase in the rolling friction. The data ofZ = 1 were compared. When μp=0.001, C1 = 4.0 and C01 = 6.8; when μp=0.06007, C1 = 2.3 and C01 = 5.8; when μp=0.102, C = 2.0 and C01 = 4.7. With the increase in μp, the angular velocity of particles near the wall was 4.0, 2.3, and 2.0, respectively, which decreased by 42.5% and 13%, respectively, and the fluctuation range that its decrease is bigger. The angular velocity of particles at the center was 6.8, 5.8, and 4.7, respectively, which decreased by 14.7% and 19%, respectively, and the reduced fluctuation range is small. This phenomenon is consistent with the characteristics of the mixed flow of particles [34].


## 5. Spatial Distribution and Probability Analysis of Force Chains

In view of the rapid loss of particles in the upper part of the silo and the relatively large particle velocity in the lower part near the funnel, we took the particle contact force in the 3, 4, and 5 areas for analysis and comparison.The left of Figure13 shows the spatial distribution of the particle normal force/average normal force ratio f at Z = 3 at time 0T, and the right shows the probability distribution of f.Figure 13
Distribution diagram of the force chain (0T).In the spatial distribution of the particle force chain, the middle force chain is stronger and the surrounding force chain is weaker.f>1 is the strong chain; f<1 is the weak force chain. Particles forming a strong chain constitute the main skeleton of the particle body, which can bear a large load, and particles forming a weak chain as a necessary filling, ensuring the overall stability of the particle system.In the probability distribution off on the right, the probability density near f=1 is relatively high, which is mainly due to the arching effect inside the particles. Under the weight of the upper particle, an anisotropic force chain network appears inside the particle, which can be divided into strong and weak as the load transfer path. Corresponding to the variation of the force chain in the figure, the closer the normal force is to the average normal force, i.e., the closer f is to 1, the wider its distribution is and the greater its distribution probability is. When the normal force is equal to the average normal force, its distribution probability reaches its peak value. However, the distribution of weak force chains and strong force chains on both sides of f=1 is less in the particle, making the distribution probability of force chains far away from f=1 smaller.In addition, the strong chains in the three regions are complementary to each other, and the proportion of the strong chain and the weak chain varies with the height of the particles. WhenZ = 3, 4, and 5, the strong chain accounts for 64%, 55%, and 37%, respectively, i.e., the closer the particles are to the bottom of the silo (Z = 3), the higher the pressure between particles is, the higher the proportion of the strong chain is, and the higher the probability of arch formation is.Changing the friction coefficient, the spatial distribution and probability distribution of force chains of particles at different moments atZ=3 are shown in Figure 14.Figure 14
Distribution diagram of the force chain (Z = 3).When the friction coefficient is changed at 0.2 T, the spatial distribution of the force chain is significantly different. It can be clearly seen that when the friction coefficient is 0.06, the strong chain of particles (f≥1) is widely distributed which accounts for 92%, while the weak chain only appears in the middle gap. When the friction coefficient is 0.001 and 0.102, there is a strong similarity in the force chain between them; in terms of probability distribution, the ratio of the weak strong chain is 51%: 49%; in the spatial distribution, and the strong chains are concentrated near the silo wall. When the friction coefficient is changed at 0.4 T, the distribution law of the force chain of maize particles is similar to that at 0.2 T, and the strong chain accounted for 87% when the friction coefficient is 0.06. Although it is 5.4% lower than that at 0.2 T, it is 33% higher than that at 58% when the friction coefficient is 0.001 and 64% higher than that at 31% when the friction coefficient is 0.102.The distribution of the force chain in this section is in accordance with the above characteristics of horizontal velocity and angular velocity variation, i.e., whenμp=0.0607, μp/μw=2, the proportion of the strong chain is higher, and the contact friction between particles is more sufficient.


## 6. Conclusions

In this paper, a silo model and the maize particle model are established to analyze the contact relationship and simulation of the particles, and the validity and rationality of the model and simulation are verified by laboratory experiments. On this basis, we carried out slice observation on the simulated silo and analyzed the data by using the method of reasonable grid division to explore the velocity field evolution mechanism and force chain distribution law of particles under the action of different internal and external friction coefficients. The following conclusions are drawn:(1)
Changing the coefficient of rolling friction between particles will affect the flow pattern of particles in the silo. The smaller the coefficient of rolling friction between particles is, the earlier the particles flow changes from mass flow to funnel flow. For a silo with the funnel, the reduction of interparticle friction will change the limit between mass flow and funnel flow, thus increasing the area of funnel flow.(2)
Near the center and lower part of the silo, the vertical velocity presents a V-shaped distribution, i.e., the vertical velocity in the middle of the silo is significantly higher than that in the shear zone near the silo wall. After changing the friction coefficient, the larger the rolling friction coefficient is, the larger the vertical velocity of the center particle is. When the ratio of internal and external friction coefficients is within the appropriate range, the flow of particles can be better promoted.(3)
The effect of particle friction characteristics on horizontal velocity is significantly reflected in the particles near the bottom of the silo. The violent degree of horizontal fluctuation of silo bottom particles is no longer a linear change with the rolling friction coefficients but is related to the coupling effect of internal and external friction. When the ratio of internal and external friction coefficients is relatively high or low, the fracture zone appears near the wall surface, leading to the formation of the boundary layer and preventing the transverse movement of particles. When the ratio of internal and external friction coefficients isμp/μw=2, the horizontal movement reached a greater peak value, and the standard value also reached a maximum of 0.0207.(4)
Under various friction coefficients, the angular velocity of particles in the central region is higher than that near the wall, but there is still a difference in the amplitude of increase. When the ratio of internal and external friction coefficient is 2, the maximum increase is 152%. Under different friction conditions, both the particles near the wall and in the central area show the phenomenon that the angular velocity decreases gradually with the increase in the rolling friction. The difference is the degree of impact. With the increase inμp, the angular velocity near the silo wall decreases by 42.5% and 13%, respectively, and the effect of friction is significant. However, the angular velocity near the center decreased by 14.7% and 19%, respectively, and the friction has little influence relatively on it.(5)
At 0T, the strong chains in the lower three areas of the silo are complementary. The proportion of the strong chain and the weak chain varies with the height of the area where the particles are located, i.e., the closer the particles are to the bottom of the silo (Z = 3), the higher the proportion of the strong chain is and the higher the probability of arch formation is.(6)
WhenZ=3, changing the friction coefficient, the spatial distribution and probability distribution of force chains of particles at different moments have great difference. When the friction coefficient is 0.06, μp/μw=2, the contact friction between particles is more sufficient, and the strong chain of particles is widely distributed, while the weak chain only appears in the middle gap. When the friction coefficient is 0.001 and 0.102, there is a strong similarity in the force chain between them; the proportion of the strength chain decreases, and the strength chain is concentrated near the silo wall. The probability distribution of force chains at Z = 3 is related to the characteristics of velocity variation.



---

*Source: 6655054-2020-11-24.xml*