Numerical Modeling of Wall Pressure in Silo with and Without Insert

A numerical model was built to simulate the experimental study conducted on the metal silo, in order to study the distribution of the resulting pressure on the silo wall with and without the use of the insert. This study is based on a previous study conducted by Wojcik [15] and Wang [9–18], in which the flow was modeled only within the hopper using the FEM. The shape and dimensions of the silo and the insert used in the model were similar to those in the authors’ previous experimental study.

Structural analysis (Abaqus, V6.12-1) and FEM were used [14]. The numerical model was implemented through several steps. This started by choosing the dimensions and sections of the model and then selecting the appropriate finite elements for each element of the model. The next step involved selecting the appropriate methods for modeling the behavior of the materials, and then the type of contact between the different elements to ensure coherence and interoperability. The appropriate boundary conditions were then selected to simulate the experimental model. Analysis and extraction of results were conducted and the adaptive mesh (ALE Adaptive Mesh) defined. This will be used to analyze the large deformation resulting from the flow of stored material, which is relative to the computational mesh used in silo modeling. To adopt the numerical model, the acquired data must be compared with the results of the approved experimental model. To clarify the above, a silo with an axisymmetric geometry is considered in the present finite element simulation. The height of the silo is 1.80 m, radius 0.8 m, and the radius at the outlet is 0.025 m.

The stored material (corn) has been simulated in the numerical model, where the dilation angle is ψ = 30°, the elastic model is E = 36 MPa, Poisson's ratio is υ = 0.3, the internal angle of friction of the corn is ‘ϕ = 23°, and the volumetric weight is γ = 800 kg/m³.

The steel wall was also considered as elastic–perfectly plastic material. To model the steel wall behavior by Abaqus program, steel elasticity modulus was introduced (E = 210 GPa), and the Poisson's ratio was υ = 0.3 and the elasticity limit was Fy = 275 MPa (see Fig. 9).

Figure 9

The flow in the silo is considered as the case of stress with displacement (stress/displacement) because the movements of the material stored in the silo are displacements that result in stress.

The element CAX4R was chosen as in Fig. 10 (four-node, bilinear, reduced integration with hourglass control), which is an axisymmetric solid element for stress condition with displacement (solid stress/displacement element). It consists of four nodes, with each node having two degrees of freedom 1, 2 at the element level. It allows modeling for half of the silo without the need to perform it on the entire model [9–15–16]. This element was used to model the corn, the silo wall, and the insert elements.

Figure 10

It is necessary to form a connection between the simulated stored materials, the simulated model, and the insert walls to obtain the exact behavior between them. Therefore, the appropriate contact between the elements should been achieved. This was accomplished as follows: during the discharge of the stored material, relative transition arises in the contact surface between the wall and the stored material. This leads to a rise in pressure between the two surfaces (contact pressure).

This pressure and a transition at the same point of contact create frictional forces and shear stress. However, this stress increases by the influence of gravitational forces, self-weight, and continuous discharge and causes slipping between the two surfaces with the continuation of transfer of pressure forces, as shown in Fig. 11. In order to achieve the mechanical, contact pair surface was used from a kind of Surface to Surface as shown in Fig. 12.

Figure 11

Figure 12

To complete the simulation of the experimental model, the boundary conditions shown in Fig. 13 were adopted. Transitions were prevented at each of the silo-based points, the posts, the discharge hole, and the corners of the input element (insert).

Figure 13

The behavior of granular material during discharge was described with a simple constitutive law for granulates, which is available in the finite element code Abaqus [14], namely with an elasto–perfectly plastic constitutive law using a linear Drucke–Prager yield criterion without cap, hardening, influence of the intermediate principal stress and rate-dependent effects [15].

The applied constitutive law includes only four material parameters for cohesionless solids: modulus of elasticity E, Poisson's ratio n, internal friction angle b, and dilatancy angle y. The material exhibits unlimited flow (perfectly plastic material) when the stress reaches yield depending on pressure. The yield surface in the deviatoric plane is circular. A modified linear Drucker–Pruger yield criterion is described by the following equation as suggested by Abaqus [14]: F=q-p.tan b{\rm{F}} = {\rm{q}} {\text{-}} {\rm{p}}.\tan \,{\rm{b}} In turn, the plastic flow potential is chosen as G=q-p. tany{\rm{G}} = {\rm{q}} {\text{-}} {\rm{p}}.\,{\rm{tan\,y}} wherein p is the mean stress and s_ij is the stress tensor p=−13σij{\rm{p}} = - {1 \over 3}{{\rm{\sigma }}_{{\rm{ij}}}} and q denotes the Mises equivalent stress (s_ij is the deviator stress tensor and d_ij is the Kronecker delta): q=32(sijsji) with sij=σij+p⋅δij{\rm{q}} = \sqrt {{3 \over 2}\left( {{{\rm{s}}_{{\rm{ij}}}}{{\rm{s}}_{{\rm{ji}}}}} \right)} \;\;\;{\rm{with}}\;\;\;{{\rm{s}}_{{\rm{ij}}}} = {{\rm{\sigma }}_{{\rm{ij}}}} + {\rm{p}} \cdot {\delta _{{\rm{ij}}}} The slope of the linear yield surface is described by the parameter tan b (b is the internal friction angle of the Drucker–Pruger condition) and y is the dilatancy (dilation) angle (tan y – slope of the potential surface). The parameter b is connected in plane strain with the internal friction angle f of the material described by a Mohr–Coulomb failure surface (Abaqus [14]) (being more realistic for granular bodies) as follows: tan β=9sin ϕ3(9−tan2ψ) +sinϕtanψ\tan \,{\rm{\beta }} = {{9\sin \,\phi } \over {\sqrt {3\left( {9 - {{\tan }^2}{\rm{\psi }}} \right)\,} + \sin \phi \tan {\rm{\psi }}}} For the parameters f = 23° and y = 30°, a plane strain matching of Drucker–Pruger and Mohr–Coulomb models results in b = 33.4°. The model used is very simple for describing very complex flow behavior of cohesionless granular materials.

Algorithms of continuum mechanics usually use two descriptions of motion: the Lagrangian description (where the state variables are a function of the material coordinate) and the Eulerian description (where the state variables are a function of the current coordinates). However, each description has some disadvantages. In the pure Lagrangian formulation, the numerical analysis usually looses accuaracy, so the size of the time increment has to be significantly reduced due to convergence problems. The formulation of the pure Eulerian leads to difficulties when free surface conditions, prescribed boundary conditions, or deformation history–dependent material properties are considered, as the element mesh is not connected to the material. The ALE formulation has been developed in order to combine the advantages of both formulations (Lagrangian and Eulerian), so the state variables are functions of the referential coordinates (not connected to the material points). In the ALE formulation, the mesh is neither connected to the material nor fixed to the spatial coordinate system. However, it can be prescribed in an arbitrary manner; thus, the ALE formulation keeps distortions of the spatial mesh under control to ensure the quality of the numerical simulation [15] (Fig. 14). The Abaqus/Explicit software takes advantage of the so-called uncoupled ALE method based on an operator split.

Figure 14

The ALE field was defined based on the previous study [2–3], on the entire model. The contact surfaces between the corn and both the wall of the model and the surface of the insert are defined as the Lagrangian boundary. The free top surface of the corn is named as the Lagrangian boundary, where the material will remain bound to the mathematical mesh and deform with it, undergoing deformation of the elements. The formed surface of the nodes at the discharge outlet (where the stored material leaves and separates the computational mesh) is considered as the Eulerian boundary (see Fig. 15).

Figure 15

As a result, the mesh velocity has to be determined in order to compute the mesh. Grid points on the surface move with the mesh velocity, but these points must remain on the free surface. Since the mesh is not connected to the material, a remap of state variables must be performed.

The freedom of the mesh movement helps to handle greater distortions compared to the distortion allowed by a Lagrangian method, with more resolution than is afforded by an Eulerian approach. The material rate of any quantity f is [27]: (1)f*m=(dfdt)m= ∂f∂t+x*m∇f{{\rm{f}}^{{\rm{*m}}}} = {( {{{{\rm{df}}} \over {{\rm{dt}}}}} )^{\rm{m}}} = {{\partial {\rm{f}}} \over {\partial {\rm{t}}}} + {{\rm{x}}^{{\rm{*m}}}}\nabla {\rm{f}} where x^*m = v^m is the velocity of a material particle at spatial point x at time t. The rate of change in grid point x^g can be expressed as: (2)f*g=(dfdt)g= ∂f∂t+x*g∇f{{\rm{f}}^{{\rm{*g}}}} = {( {{{{\rm{df}}} \over {{\rm{dt}}}}} )^{\rm{g}}} = {{\partial {\rm{f}}} \over {\partial {\rm{t}}}} + {{\rm{x}}^{{\rm{*g}}}}\nabla {\rm{f}} where x^*g = v^g is the velocity of a grid point coinciding with x at time t. Solving the spatial rate of change ∂f∂t{{\partial {\rm{f}}} \over {\partial {\rm{t}}}} from Eq. (1) and substituting into Eq. (2) results in: (3)f*g=f*m+(x*g−x*m)∇f=f*m+c∇f{{\rm{f}}^{{\rm{*g}}}} = {{\rm{f}}^{{\rm{*m}}}} + \left( {{{\rm{x}}^{{\rm{*g}}}} - {{\rm{x}}^{{\rm{*m}}}}} \right)\nabla {\rm{f}} = {{\rm{f}}^{{\rm{*m}}}} + {\rm{c}}\nabla {\rm{f}} wherein c is the convective velocity. Thus, the variation of a quantity f for a given particle is the local variation plus a convective term, taking into account the relative motion between the spatial system and the material. However, determining the convective term numerically is challenging due to its discontinuous character [15]. In the ALE approach, the state variables are calculated at the points x^g. This operation consists of two steps. First, an updated Lagrangian step is performed, which results in calculating the material displacements. The material body deforms from its material configuration to its spatial one. Secondly, the mesh velocity v^g is calculated during a smoothing phase performed at each time step in order to reduce mesh distortions. The boundary after the remeshing approximately coincides with boundary obtained with the Lagrangian calculation, during the movement of the inner nodes. The mesh displacement and mesh velocity are taken as constants during a time step. A transfinite mapping method can be frequently used for the mesh management. Finally, a remap of state variables of the updated Lagrangian phase into new mesh, including the calculation of the convective term is carried out. The linear distribution of the convective term in two directions (which is required when using rectangular elements) is found in two steps: a) quadratic interpolation constructed at the integration points of the middle element and b) a trial linear distribution found by differentiating the quadratic function in the middle element. This distribution is limited by reducing its slope until it reaches the minimum and maximum values within the range of the original constant values in the adjacent elements.

The numerical model was analyzed using Abaqus/Explicit, with nonlinear analysis. Adaptive mesh deformation was observed over time, and the resulting dynamic pressure values were plotted on the model, for the case “without insert” (input element) (see Figs 16 and 17).

Figure 16

Figure 17

The explicit analysis is characterized by the fact that it is computationally effective for analyzing large models with relatively short dynamic response times where a large number of short-time increments can be accomplished with high efficiency. In addition to the analysis of discontinuous phases or processes, it also allows for the definition of general communication conditions, and it is beneficial when using models with wide distortions, such as the flow state we are studying in this paper.

To validate the results of the numerical model, they were compared with the results of the experimental model of the present study for the silo without using insert. Figure 18 shows the comparison of the pressure distribution on the wall of the numerical model with that of the experimental study. At the ninth level, the pressure values are generally converged between the experimental and numerical studies. The maximum relative errors of 30%, 12%, and 29% are obtained in the hopper, transition section and cylindrical part, respectively. This difference in pressure ratios along the height of the silo can be attributed to the complex behavior of the stored material exhibited during the discharge, which is difficult to completely enumerate. In addition to the effect of the insert on the discharge process and the distribution of the resulting pressure on the wall, and considering that the wall has not deformed. Within the limits of acceptance, the results of the numerical model can be adopted and considered appropriate for the completion of the parametric study later.

Figure 18

In the second stage of modeling, and to simulate the experimental silo model using the insert, the results of the numerical model were compared with the ones of the experimental model in the current study to verify their validity.

Figure 21 shows comparison of the resulting pressure distribution on the silo wall of the numerical model during discharge with the pressure distribution resulting from the experimental study. It is observed that the pressure values are generally converged between the current experimental and numerical study cases, where the maximum relative errors of 32%, 13%, and 22% are obtained in the hopper, transition section, and cylindrical part, respectively. Within the limits of acceptance of this difference, the results of the numerical model can be adopted.

Figure 21

One can easily observe the extent of convergence between the pressure distribution resulting from the experimental study and the resulting pressure in the numerical model, especially at the transition section (the section between the cylindrical part and the hopper). The same thing can be noticed in the conical hopper. On comparing the pressure values in the current study (experimental and numerical), it was found that the difference reached 36% within the hopper, 18% in the transition section, and 32% within the cylindrical part.