Dynamic simulation analysis of elastic overpressure fluctuation on silo wall for the arch action

2.1. Formation process of arch

Grain is a kind of granular materials. Granular materials have characteristics like both liquid and solid, but their physical properties are not the same as solid and liquid entirely. On a mesoscopic scale, the essence of granular materials is a collection of many discrete particles, showing a liquidity in some respects. On a macroscopic or bulk scale, the storage materials can still be considered isotropic granular aggregate.

The storage materials would not arch in the static condition. The “arch action” is a process of dynamic equilibrium. The particles in the silo fall under gravity during discharge. The particles show a complex flow behavior because of the friction and collision between the particles, and between the particles and their environment. The forces and moments in the system are renewed every moment. The particles could find the best place to transfer pressure through constant adjustment of position and stress. A complete force chain connecting the two walls of silo is formed gradually. No particle vacancy is allowed to exist in the force chain, which in order to achieve the transmission of force and the impediment to the flow of upper material. In addition, it is also required that the arch feet have enough strength to support the upper material. Since the formation of the force chain affects the outflow of the particles, the decrease of the particle’s velocity causes an inertial force in the same direction with the particles velocity.

2.2. Elastic overpressure fluctuation

In the process of the arch formation, the decrease of material velocity in the silo is not an integral, continuous process. In practical engineering, the silo wall is not absolutely rigid, the wall has elasticity. The overpressure on the wall at the arch feet causes a slight lateral expansion of the silo [15], more materials will inflow to the section because of the volume increase of the expansion section, and the density of the storage in the expansion section increase because of compaction. The inflow materials also have certain velocity, they flow in the expansion section, and their velocity would decrease to zero, then the inertial force would be produced.

The wall expansion is illustrated in Fig. 1. Moreover, the velocity of particles in a micro region above the arch must be affected firstly and rapidly decreased to zero, then the inertial force appears, the silo wall pressure increases. The velocity of particles in the other micro region above the region aforementioned were further affected to decrease, the pressure on the wall section also be influenced.

The overpressure, caused by the materials continuously flowing into the expansion volume and the velocity of material decrease in order, propagates in the form of fluctuation from the arch feet.

Fig. 1Silo wall expansion

4.1. Experimental model

The experimental model was established based on the engineering example of silo accident in Shandong, Xinliang [17]. In this project example, the geometrical scale was 20 m in height, 6 m in the inner diameter and 160 mm in wall thickness, and the diameter of the outlet is 1 m.

The scale model of 20:1 was studied in the experiment, the geometrical symbols of the silo are illustrated in Table 1 and Fig. 2.

Fig. 2The geometrical symbols of the silo for experiment

The silo is made of plexiglass with a thickness of 8 mm for convenient viewing. The experimental model silo is filled with wheat with the diameter of 5 mm to 6 mm. The plate of the outlet is removed to perform the discharging experiment after charging is completed. Repeat experiments to achieve the “arch action” by adding dust, changing the moisture content, placing sandpaper on the wall, and changing the properties of the wheat. The grain arch is shown in Fig. 3.

Fig. 3Experimental model

4.2. Monitoring points

10 monitoring points are arranged along the left side silo wall, strain gauges are placed on each monitoring point to monitor the pressure. The pressure for the arch action could be obtained. The distribution of the monitoring points is shown in Fig. 4.

Fig. 4The distribution of the monitoring points on silo wall

4.3. Experimental results

The pressure at the monitoring points numbered 4, 5, 6, 7, and 8 were selected to monitor. The pressure at each monitoring points is shown in Fig. 5.

Fig. 5The pressure at the monitoring points

In the experimental results, after the unloading of the silo, there is a momentary drop in the pressure inside the 0.1 seconds, and the side pressure fluctuation reaches the maximum with the formation of the arched force chain. The formation of the “arch action” pressure could be clearly seen through the chart of the dynamic silo wall pressure, the pressure fluctuated during discharge, and the overpressure for arch action occurrences in the fluctuation peak. But limited by objective experimental conditions, it is difficult to reflect the formation of overpressure at mesoscopic scale accurately. Therefore, there would be further research by PFC simulation in the later part of the paper, and the results of experimental method would be compared with simulation to prove the accuracy of PFC simulation.

5.1. Contact model and constitutive model

5.1.2. constitutive model

1) Between particles.

The unit vector model and the contact deformation of two particles are obtained by Newton’s Second Law of motion. Two particles $i$ and $j$ are in contact with each other, as shown in the Fig. 7. The central coordinates of two particles are $(x_i,y_i)$, $(x_j,y_j)$ respectively. There is little difference in the particle diameter distribution of wheat, therefore, it can be considered that the radius of both particles is $r_c$ normal vector at the point of contact is $\overrightarrow{n}=(n_x,n_y)$ the tangential vector is $\overrightarrow{t}=(t_x,t_y)$.

Which can be obtained:

3

${\overrightarrow{n}}_x=\overrightarrow{n}\mathrm{c}\mathrm{o}\mathrm{s}\alpha ,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\overrightarrow{n}}_y=\overrightarrow{n}\mathrm{s}\mathrm{i}\mathrm{n}\alpha ,$

4

${\overrightarrow{t}}_x=\overrightarrow{t}\mathrm{c}\mathrm{o}\mathrm{s}\alpha ,\mathrm{}\mathrm{}{\overrightarrow{t}}_x=\overrightarrow{t}\mathrm{s}\mathrm{i}\mathrm{n}\alpha ,$

where ${\overrightarrow{v}}_i$, ${\overrightarrow{v}}_j$ are used to represent the velocity at the center point of particles called $i$, $j$ the relative velocity of the two particles at the contact point can be expressed as:

5

${\overrightarrow{dv}}_c={\overrightarrow{v}}_{ci}-{\overrightarrow{v}}_{cj}.$

Fig. 6Linear contact model: a) tangential contact model, b) normal contact model

a)

b)

Fig. 7Contact between two particles

The velocity of particles $I$ and $j$ contact points can be expressed by Eq. (6), (7):

6

${\overrightarrow{v}}_{ci}={\overrightarrow{v}}_i+{\overrightarrow{w}}_i\times \overrightarrow{r},$

7

${\overrightarrow{v}}_{cj}={\overrightarrow{v}}_j+{\overrightarrow{w}}_j\times \overrightarrow{r},$

which can be obtained:

8

${\overrightarrow{dv}}_c=\left({\overrightarrow{v}}_i-{\overrightarrow{v}}_j\right)+\left({\overrightarrow{w}}_i\times \overrightarrow{r}-{\overrightarrow{w}}_j\times \overrightarrow{r}\right).$

The normal and tangential components of the relative velocity is:

9

${\overrightarrow{dv}}_n={\overrightarrow{dv}}_c\overrightarrow{n}=\left({\overrightarrow{v}}_i-{\overrightarrow{v}}_j\right)\overrightarrow{n},$

10

${\overrightarrow{dv}}_t={\overrightarrow{dv}}_c\overrightarrow{t}=\left({\overrightarrow{v}}_i-{\overrightarrow{v}}_j\right)\overrightarrow{t}+\left({\overrightarrow{w}}_i\times \overrightarrow{r}-{\overrightarrow{w}}_j\overrightarrow{r}\right).$

For the particles $i$, $j$ the normal and tangential velocity components of velocity at the contact point both cause displacements in the corresponding directions. According to the velocity component calculated in the above manner, the variation of displacement component within a moment could be obtained as:

11

$\mathrm{\Delta }{\overrightarrow{r}}_n={\overrightarrow{dv}}_n\mathrm{\Delta }t,$

12

$\mathrm{\Delta }{\overrightarrow{r}}_t={\overrightarrow{dv}}_t\mathrm{\Delta }t.$

According to the law of force displacement:

13

$\begin{array}{l}\mathrm{\Delta }{\overrightarrow{F}}_n=k_n\mathrm{\Delta }{\overrightarrow{r}}_n,\\ \mathrm{\Delta }{\overrightarrow{F}}_s=k_s\mathrm{\Delta }{\overrightarrow{r}}_s.\end{array}$

Could be obtained, the motion equation of the particle $i$ could be obtained by Newton’s Second Law of motion:

14

$\sum \overrightarrow{F}=m\frac{\overrightarrow{dv}}{dt}=m\frac{\stackrel{}{\overrightarrow{d^{2}r}}}{dt},$

15

$\sum \overrightarrow{M}=I\frac{\overrightarrow{dw}}{dt}=I\frac{\overrightarrow{d^2\theta }}{dt},$

where $I$ is the particles’ moment of inertia. In the linear contact model, it is generally considered that the moment at particle contact point is always zero).

Mesoscopic parameters are changing at every moment in the operation. Taking particle $i$ as an example, the force is updated in a time stepping like:

16

${\overrightarrow{F}}_n(t+\mathrm{\Delta }t)={\overrightarrow{F}}_n\left(t\right)+\mathrm{\Delta }{\overrightarrow{F}}_n\left(t\right),$

17

${\overrightarrow{F}}_t(t+\mathrm{\Delta }t)={\overrightarrow{F}}_t\left(t\right)+\mathrm{\Delta }{\overrightarrow{F}}_t\left(t\right).$

The velocity is updated in a time stepping like:

18

${\overrightarrow{v}}_i(t+\mathrm{\Delta }t)={\overrightarrow{v}}_i\left(t\right)+\frac{\sum \overrightarrow{F}}{m}\mathrm{\Delta }t,$

19

${\overrightarrow{w}}_i(t+\mathrm{\Delta }t)={\overrightarrow{w}}_i\left(t\right)+\frac{\sum \overrightarrow{M}}{I}\mathrm{\Delta }t.$

The position of particle is updated in a time stepping like:

20

${\overrightarrow{r}}_i(t+\mathrm{\Delta }t)={\overrightarrow{r}}_i\left(t\right)+\overrightarrow{v_{\dot{i}}\left(t+\frac{\mathrm{\Delta }t}{2}\right)}\mathrm{\Delta }t,$

21

${\overrightarrow{\theta }}_i(t+\mathrm{\Delta }t)={\overrightarrow{\theta }}_i\left(t\right)+\overrightarrow{v_{\dot{i}}\left(t+\frac{\mathrm{\Delta }t}{2}\right)}\mathrm{\Delta }t.$

2) Particles and walls.

The wall $j$ revolve around $(x_j,y_j)$, and the central coordinates of the particle $i$ is $(x_i,y_i)$. The coordinates of the contact point between the wall and the particle is $(x_c,y_c)$. The normal vector at the point of contact is $\overrightarrow{n}=(n_x,n_y)$, while the tangential vector is $\overrightarrow{t}(t_x,t_y)$${\overrightarrow{v}}_i$, ${\overrightarrow{v}}_j$ is used to express the velocity at the center point of particle $i$ and the rotation center of wall $j$, respectively the relative velocity at the contact point could be expressed by:

22

${\overrightarrow{dv}}_c={\overrightarrow{v}}_{ci}-{\overrightarrow{v}}_{cj},$

23

${\overrightarrow{v}}_{ci}={\overrightarrow{v}}_i+{\overrightarrow{w}}_i\times \overrightarrow{r},$

24

${\overrightarrow{v}}_{cj}={\overrightarrow{v}}_j+{\overrightarrow{w}}_j\times {\overrightarrow{r}}_w,$

where $w_j$ is the angular velocity of the wall rotation. The wall would not translate, so $v_j=$ 0:

25

$\therefore {\overrightarrow{dv}}_c={\overrightarrow{v}}_i+\left({\overrightarrow{w}}_i\overrightarrow{r}-{\overrightarrow{w}}_j{\overrightarrow{r}}_w\right).$

The component of relative velocity in the normal and tangential direction is:

26

${\overrightarrow{dv}}_{cn}={\overrightarrow{dv}}_c\overrightarrow{n}={\overrightarrow{v}}_i\overrightarrow{n}+{\overrightarrow{w}}_j\times {\overrightarrow{r}}_w\overrightarrow{n},$

27

${\overrightarrow{dv}}_{ct}={\overrightarrow{dv}}_c\overrightarrow{t}={\overrightarrow{v}}_i\overrightarrow{t}+{\overrightarrow{w}}_i\times \overrightarrow{r}\overrightarrow{t}.$

The updating methods of contact force normal and tangential increment, position and velocity are all mentioned above.

Fig. 8Contact between particle and wall

5.2. Simulation model and parameter

According to the existing experimental model size [11], the simulation model uses the same scale model to study the dynamic simulation. A silo model with an internal diameter of 300 mm, a height of 1000 mm, a outlet diameter of 50 mm, and the cone half angle of 45° is built to simulate particle flow with PFC2D.The particles gradually settle on to the silo, the physical parameters of the wheat and walls, as shown in Table 2, were selected on the basis of the grain properties provided in document. The charged particles are shown in Fig. 9. Then delete the wall of the outlet, the materials will discharge, there are arched stress chains in the vicinity of the outlet under the “arch action”.

In addition, the simulation parameters about the walls and particles are listed in Table 3. the wall and particles in the simulation are represented as ‘Wall’ and ‘Ball’ respectively, where the kn and the ks are their normal and tangential stiffness.

20 virtual monitoring points are arranged along the wall to monitor the silo wall pressure. The monitoring results would be compared with the experimental method result, it is found that the fluctuation tendency of pressure in the experiment is consistent with the simulation result, which can prove the accuracy of the PFC simulation.

Fig. 9The charged particles

5.3. Monitoring points

The particles present a complex, asymmetrical, non-uniform flow structure. Therefore, according to the above experimental model, 20 monitoring points are uniformly distributed on the wall at intervals of 100 mm, to monitor the pressure at each points. The distribution of monitoring points is shown in Fig. 10.

Fig. 10The distribution of monitoring points in simulation

6.1. Stress chains

The velocity, position and stress of particles change constantly and a irregular force chain in arch on both wall sides is gradually formed during discharge. If the particles can transfer stresses well and the arch feet have enough strength, the arch, which hinders the outflow of the granular materials, would form. The formation of force chain in arch is shown in Fig. 12.

Fig. 11The model of the arch in silo

Fig. 12Force chain formation process

a) 105000 timestep

b) 106000 timestep

c) 107000 timestep

6.2. Simulation results

In this simulation, the arch feet form in the lower 1/3 part of the silo wall. Therefore, the pressure at monitoring points in the upper 2/3 part in silo is affected by the arch action, while the pressure of the monitoring points 1, 2, 3 was not affected. And when the arch forms, a part of particles in silo have been discharged, a zero pressure area on the wall have appeared in the upper 1/5 part of the wall. In summary, monitoring points numbered 4, 5, 6, 7, 8 were selected for observation and analysis.

The pfc software could directly demonstrate the dynamic pressure, but it was difficult to observe the pressure at each time step accurately. Therefore, in this article, pfc was used to output the Excel, which illustrates the pressure at each step, form a linear graph of dynamic pressure on silo wall. The dynamic pressure of the five points on the left wall is shown in Fig. 13. For the arch action, the pressure distribution and overpressure coefficient are shown in Fig. 14 and Table 4.

The dynamic pressure of the five points on the right wall is shown in Fig. 15. For the arch action, the pressure distribution and overpressure coefficient are shown in Fig. 16 and Table 5.

Fig. 13The dynamic pressure on the left wall

Fig. 14Pressure for each on the left wall

Fig. 15The dynamic pressure at each monitoring point on the right wall

Fig. 16Pressure at each monitoring point on the right wall

It could be clearly seen that the pressure of silo wall and the overpressure coefficient for the arch action increase with depth, in addition, the overpressure of different monitoring points for the arch action would not occur at the same time. It perhaps related to the particle velocity decrease in order and related to the elastic deformation of wall propagating from the arch feet to the upper part.

6.3. Analysis about elastic fluctuation

The simulation results demonstrated the order of different points overpressure occurrence for the arch action on the left wall, as shown in Fig. 17.

The simulation results demonstrated the order of different points overpressure occurrence for the arch action on the right wall, as shown in Fig. 18.

In view of the symmetry of the silo structure, most scholars deem that the distribution of mesoscopic parameters in silos should also be symmetrical. However, in the actual project, the flow of the particulates in the silo is complex, showing non-uniformity and asymmetry, the time step of different points overpressure occurrence for the arch action and propagation of elastic overpressure fluctuations on different wall sides are all slightly different. In case of the left wall, as could be demonstrated from Fig. 17, overpressure for arch action at the No. 4 monitoring point occurs earliest, the overpressure for arch occurs when 105850 steps are simulated. As the simulation continues to calculate, the pressure for the arch action at 4, 5, 6, 7, 8 monitoring points produced in order. The overpressure for the arch action at the No. 8 monitoring point occurs when 105920 steps are simulated. Although the elastic fluctuation on the right wall is slightly different from the result on the left wall, it could be demonstrated from Fig. 18. Fig. 9 that the propagation of overpressure fluctuation is still almost uniform.

Fig. 17The time step of overpressure occurrence for the arch action on the left wall

Fig. 18The time step of overpressure occurrence for the arch action on the right wall