This paper presents a methodology based on the finite element method to simulate the flow of granular materials. Moreover, it allows proper estimation of dynamic pressure during silo discharge since this subject is still under discussion, especially for designing silos with an insert (an input element). A 2-D simulation of the discharge process of a cylindrical silo with cone and a central discharging orifice was performed. Two cases were studied, with and without using insert in silo. Numerical analysis was carried out with the help of the uncoupled arbitrary Lagrangian–Eulerian (ALE) approach. The resulting dynamic pressure distribution on the silo wall for each of the two cases was inferred numerically. The resulting values of pressure were compared with the results of the experimental study on a cylindrical metal silo to demonstrate the accuracy of the numerical model in determining the dynamic wall pressure, especially in the case of using an insert in silo during discharge.

#### Keywords

- Numerical modeling
- silo
- silo model
- flow pressure
- insert
- wall pressure

Silos are important structures in agricultural and industrial countries. This type of structures still faces many problems causing damages or even collapse. The main reasons for silos’ failures are problems related to storage and discharge of silos [1–2]. More specifically, it can be said that the major cause of failure is the inaccuracy in determining the dynamic pressure produced during the discharge, which is mainly related to the flow pattern. It has been demonstrated that using an input element (insert) inside the silo can help to direct the flow and improve the shape of the material flow, which leads to a reduction in the resulting dynamic pressure. Several studies have been conducted on the effectiveness of the use of the element (insert) inside the silo. The flow in silos containing inserts has been investigated experimentally [3–4–5–6–7–19–20–26] and numerically using the finite element method (FEM) [21–22–28–29]. It also dealt with the effect of insert (an input element) on the distribution of the resulting dynamic pressure on the silo wall because it is of great importance to the structural state of the silo [17]. The idea of using inserts in order to improve the material flow pattern from funnel (core) flow to mass flow during emptying was originally proposed by Johanson [19] who observed that an inverted cone located at a specific position above the outlet might significantly reduce the size of stagnant zones during funnel flow. Rombach and Eibl [23] performed a dynamic finite element analysis and presented dynamic pressure profiles, dependent on space and time just at the beginning of the emptying stage. In other studies, either remesh–rezoning technique or assumed failure boundary was used to describe large deformation occurring during the discharge stage in order to avoid mesh distortions [24–25]. However, with the general lack of understanding and information concerning the discharge process with insert, the finite element FE modeling of the wall pressures during discharge is relatively rare and requires more fundamental research [8–10–11]. In addition, estimation of the dynamic horizontal pressure produced on the silo wall during discharge is still mainly based on a simple analytical model of Janssen [13]. The use of this model is limited to the fact that the stored material is at rest and cannot be applied in the case of material flow. Moreover, it cannot be applied if a proposed insert is used within the silo.

This paper presents a finite element simulation of a silo discharge to investigate the dynamic pressure during discharge with and without using an insert. The FE simulation is performed using the uncoupled arbitrary Lagrangian–Eulerian (ALE) formulation with an adaptive meshing technique using the FE software Abaqus [14]. This method allows for the simulation of almost the entire silo discharge process without involving the mesh distortion problem, which is often encountered in the modeling of granular material flow [12–15]. Thus, we can deduce the dynamic pressure on the silo wall with the case of using the insert. In addition, this method can be used to estimate the resulting pressure during discharge for different positions and dimensions of the insert.

This study is the continuation of a previous experimental research carried out by the authors to investigate the flow pressure of granular material in silos and determine the best location for the input element (insert). This would facilitate obtaining a better flow of granular material, which reduces the resultant dynamic pressure [26]. Furthermore, this research is complementary to our previous research which investigated the best shape of the insert to enhance the flow shape [6].

The main goal of this study is to build a numerical model to simulate the granular flow within the silo with a proposed insert, using FEM and uncoupled ALE approach, in order to deduce the distribution of dynamic horizontal pressure produced on the silo wall during discharge.

Based on the previous study [26], a cylindrical metal silo with conical hopper was used (Fig. 1). The dimensions of the silo were 150 cm in height and 80 cm in diameter. The hopper angle was α = 50°; the outlet hole was central and circular with a diameter of 2.5 cm (Fig. 2). The used inserted element was a top cone and a trunk cone bottom (CTC) as shown in Figs 3 and 4. The stored granular material used was corn, with diameter ranging from 3 to 4 mm; its volumetric weight was γ = 800 kg/cm^{2} (see Fig. 5).

Nine levels were identified to study the resulting dynamic pressure on the silo wall. Four levels were distributed in the hopper (H_{1}, H_{2}, H_{3}, H_{4}) and one at the transition section (T_{1}) (the section between the cylindrical part and hopper). Additionally, four levels were distributed at the cylindrical section of the silo (C_{1}, C_{2}, C_{3}, C_{4}; see Fig. 2).

A cell pressure gauge was used to obtain the pressure values on the silo wall at the specified levels. A special device was used to connect the cell pressure gauge to a computer with a special program to draw pressure curves resulting from the flow during the discharge (Fig. 6).

The resulting pressure on the model wall during the discharge was investigated for two cases:

the flow pressure without using insert and

the flow pressure while using insert located at the transition section.

The dimensions of the insert were chosen according to a previous study [6]. This study was conducted by the authors, in which it was found that the insert whose diameter is equal to quarter of the silo diameter and the angle of inclination of the insert wall with the horizon 60° give the best flow. The resulting dynamic pressure change was plotted versus the height of the model in a chart (see Fig. 7).

The chart presents the distribution of the resulting dynamic horizontal pressure on the silo model wall during discharge at the nine defined levels without using of insert. Figure 8 shows the dynamic pressure distribution that resulted on the silo wall model during discharge while using the insert inside the silo model. It is clear from the two charts that the value of dynamic pressure recorded at the transition section for the case of the silo without the insert element is greater than the pressure recorded at the same level for the silo with the insert, considering that the transition section is the most critical section in silo.

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^{3}.

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).

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.

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.

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).

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]:
_{ij} is the stress tensor
_{ij} is the deviator stress tensor and d_{ij} is the Kronecker delta):

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.

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).

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]:
^{*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:
^{*g} = v^{g} is the velocity of a grid point coinciding with x at time t. Solving the spatial rate of change
^{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).

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 verify the present numerical model and demonstrate its ability to model discharge without insert, we consider comparisons with the results of the experimental study and adopting the values of the variables of this analysis to perform the required modeling when using the insert. Therefore, the numerical model verification process is divided into two stages: the first stage is the model verification stage without using the insert and the second stage is the model verification stage using the proposed insert.

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.

The numerical study variables for the case of discharge without insert were adopted and applied to the model with the use of insert located at the transition section. As mentioned earlier, transfers were prevented at silo-based points on the columns, at the discharge hole, and at the corners of the insert to simulate the arrangements of the experimental model used in this research (see Figs 19 and 20).

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.

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.

This research has shown that the finite element model using the uncoupled (ALE) approach is an effective technique for the simulation of silo discharge process that involves large plastic deformations. Moreover, by defining an adaptive mesh for the granular solid and suitable boundary conditions for the mesh region, it is possible to model in principle almost the entire process of the discharge in a satisfactory manner.

The modeling of granular material flow was carried out inside the silo. The stored material and the silo wall were modeled. The appropriate element in the software program was chosen in order to represent elements that were the subject of the current study and to ensure the appropriate contact between the elements. Additionally, the necessary boundary conditions were adopted. To improve the mesh quality during discharge and reduce distortion, adaptive contact was used between the mathematical mesh and the physical mesh. This enabled faster and more accurate solutions to be reached through the ALE adaptive mesh technology. Numerical model was analyzed using Abaqus/Explicit with nonlinear analysis. Adaptive mesh deformation was observed over time, and the resulting dynamic pressure values over time were extracted on the model wall, which were compared with the results of the experimental study. The values of the model variables were adopted to perform modeling of granular material flow within the silo with insert. Finally, the following conclusions can be made:

The results revealed that the model can capture qualitatively and quantitatively several aspects of the resultant dynamic pressure.

The model provides the potential of being used as a predictive tool to find accurate values for the resulting pressure on the silo wall in the case of using the insert, which is an essential issue in the design process.

The mathematical model provides accurate information for the resultant pressure, especially at the transition section of the silo, which is the most critical region in which a peak pressure is always recorded.

The model can be used to obtain information for the internal stresses on the inserts, which is difficult to obtain experimentally and can be relevant for design purposes.

The model can provide the resulting pressure during discharge for different positions and dimensions of the insert in order to have optimal flow of granular materials, which leads to less pressure. This is due to the significant relationship between the flow pattern and the resultant wall pressure.

However, several aspects are yet to be explored and further experimental and numerical work is required to better define the physics that underlies some of the parameters of the model.

Safety analysis of the Żelazny Most tailings pond: qualitative evaluation of the preventive measures effectiveness Impact of subgrade and backfill stiffness on values and distribution of bending moments in integral box bridge Laboratory tests and analysis of CIPP epoxy-resin internal liners used in pipelines – part I: comparison of tests and engineering calculations Behaviour of eccentrically inclined loaded rectangular foundation on reinforced sand Application of the thermoporoelasticity model in numerical modelling of underground coal gasification influence on the surrounding medium Load capacity of the mixed bench and slab foundation. Numerical simulations and analytical calculation model Termite Mound Soils for Sustainable Production of Bricks Laboratory tests and analysis of CIPP epoxy resin internal liners used in pipelines – part II: comparative analysis with the use of the FEM and engineering algorithms Geostatistical analysis of spatial variability of the liquefaction potential – Case study of a site located in Algiers (Algeria) The effect of the pre-wetting of expanded clay aggregate on the freeze-thaw resistance of the expanded clay aggregate concrete