The Friction Stress Distribution of Loose Materials to Silo Walls in Deep Pyramidal Silos †

Referring to the classification of vertical silos, the classification rules for deep silos and shallow silos on pyramidal silos and conical silos have been established. Based on the basic assumptions of Janssen Formula, the friction stress distribution of loose materials to silo walls was discussed in detail on deep pyramidal silos. An analytic equation was established used the slice method. Finite element model was developed to simulate the friction stress distribution. The analytic result was almost the same with the simulation result. However, the equation yielded a minor error mainly because there were some differences between the assumptions presented in the equation and the real state of materials. Based on the simulation results, the analytic equation was revised. The results showed that the friction stress distribution of loose materials to silo walls in deep pyramidal silos was related to the material characteristics and design parameters of silo walls, as well as the material powder attribute. The revised equation precisely indicated the contact friction stress behaviors of loose materials and provided a shortcut method of calculation for related research.


Introduction
As a special industry thin-shell structure, silos have many advantages, ranging from having a simple structure that is low cost, space efficient, and stable, having an easy loading and unloading process that can be conveniently managed.Silos are widely used in many fields to stockpile loose materials, such as mineral, grain, foodstuff, coal, chemical fertilizer, and concrete.In the state of stockpiling loose materials, the mechanics behavior of the silos and the loose materials are all very complex (Jansseune A. et al., 2016a,b;Park H.W. et al., 2016;Jagtap P. et al., 2015;Colonnello C. et al., 2014).However, with the mature development of the finite element method and discrete element method, the numerical simulation method managed to contribute greatly to the silo research (Ding S. et al., 2015;Wang Y. et al., 2014;Gallego E. et al., 2010) and the loose materials or the solid particle behavior, including the powder friction angles, forces, stresses and so on (Hua X. et al., 2013;2015).Basing on relevant literature, many scholars are interested in study-ing the lateral pressure of stored loose material in silos.Although the famous static deep-silo lateral pressure formula has already been put forward by JANSSEN over 100 years ago (Janssen, H.A. 1895), the JANSSEN formula is still being widely cited, used and discussed (Matchett A.J. et al., 2015;LI Z.F. et al., 2014).In the meantime, many new research results on lateral pressure have been put forward (Wang Y. et al., 2014;Couto A. et al., 2013;Nascimento J.W.B.d. et al., 2013).
In general, these research works are mainly concerned with cylindrical silos (Iwicki P. et al., 2016;Jansseune A. et al., 2016a,b;Sadowski A.J. et al., 2014), particularly on the pressure (Wang Y. et al., 2014;Couto A. et al., 2013;Nascimento J.W.B.d. et al., 2013), stability (Sondej M. et al., 2016;Iwicki P. et al., 2015;Lozano C. et al., 2015), and fluidity effect or discharge flow (Ding S.et al., 2015;Gallego E. et al., 2015;Benyamine M. et al., 2014) of the loose material on the silos.The square platform silos have already been studied by some scholars (Goodey R.J. et al., 2006;Guines D. et al., 2001;Brown C.J. et al., 2000).However, there are only a few studies so far which analyze the mechanical behavior of loose materials in pyramidal silos.Comparing with the general cylindrical silos, the structure of the pyramidal silos is special.It has been applied in some engineering fields because of its simple structure and easy manufacture.The friction stress distribution of loose materials in these silos is helpful to study the force obtained in storing and unloading, optimize the design structure and cut down the manufacturing cost.
In this paper, an analytic equation of loose material friction stress distribution was established used the slice method in deep pyramidal silos with reasonable assumptions and simplifications.The friction stress was solved by a finite element contact analytic method, and the contactfriction state was discussed in deep pyramidal silos.Lastly, a more reasonable equation was published after verification and modifying to provide foundation for the deep pyramidal silos' structure design and storage efficiency.

The analytic formula of friction stress distribution
2.1 The silo structure EN 1991-4.Eurocode 1: actions on structures, Part 4: silos and tanks" (CEN, 2004) prescribed that a deep cylindrical silo was the one with aspect ratios (the ratio of h to D) greater than or equal to 2.0; a shallow cylindrical silo was the ratio of h to D greater than or equal to 0.4 and less than 1.0; otherwise, it was a medium deep silo [see Fig. 1(a)].However, there was no rule for conical silos and pyramidal silos.Here, such a rule was prescribed for them with reference to "EN 1991-4.Eurocode 1: actions on structures, Part 4: silos and tanks" (CEN, 2004), as shown in Fig. 1.
A pyramidal silo has some features as follows: (1) it has a simple structure (see Fig. 2) and is easy to manufacture; (2) the friction stress of the materials varies with the obliquity of silo walls; (3) the contact friction behaviors of loose materials in deep silos as presented in this paper is different from those of the shallow silos.
In this paper, the friction stress distribution of loose materials to silo walls was investigated for the case of the deep pyramidal silos.

The friction stress distribution
In 1895, according to two basic assumptions, the German civil engineer JANSSEN gave the famous calculation formula of the static pressure in a cylindrical silo, namely, the JANSSEN Formula (Janssen H.A., 1895).The two as- sumptions were expressed in the following: (1) the vertical pressures were equally distributed in the same depth; (2) the ratio of horizontal pressure to vertical pressure was a constant k, which was the lateral pressure coefficient.
Below are some of the assumptions in this paper based on above: (1) Loose materials are in an ultimate stress state.Before or after the state, the external force on the loose materials approached a limiting value and the loose materials will suddenly collapse.At this moment, the lateral pressure coefficient k reflects the force exerted from another direction when the force is acted upon by the loose materials either horizontally or vertically.
(2) Particles are minute (the diameter is less than 1 mm), indicating that the loose materials are homogeneous.
(3) The loose materials must have unified physical properties and be in a uniformly filling state to assure a constant internal friction factor.
(4) The vertical pressure is equal on the same plane.
As shown in Fig. 2, we supposed that the horizontal cross section was rectangular, the height was h, the obliq-uities of the four walls to the axis vertical curve were ϕ 1 , ϕ 2 , ϕ 3 and ϕ 4 , respectively, and the bulk density of loose materials was ρ.Starting at the lowermost peak, a thin layer unit with vertical thickness dy was the research object.The thin layer unit was approximated as a cuboid.
Let the areas of the four walls of the thin layer unit be dA 1 , dA 2 , dA 3 , and dA 4 , and the friction forces of the loose materials in walls be f dA1 , f dA2 , f dA3 and f dA4 .Thus, along the vertical direction, the force balance of the loose materials in the thin layer unit could be expressed in the following: Upon combining Eq. ( 1), (2), and (3), a relationship relating dp to dy could be expressed in the form of Eq. ( 4): Note that p = p 0 if y = h from Fig. 2. By solving the differential Eq. ( 4), the vertical pressure p of the loose materials in deep pyramidal silos could be determined in Eq. ( 5): In most cases, λ ≠ 1.
If there were no loose materials or other pressures at the top of the silos, p 0 = 0.
In the engineering practice, most pyramidal silos were symmetrically designed, namely ϕ 1 = ϕ 3 , and ϕ 2 = ϕ 4 .Thus, the λ might be simplified in Eq. ( 6) in the following: In some cases, the design parameters of the pyramidal silos were the length A and B of silo top opening as Fig. 2 Fig. 2 The diagram of the vertical pressure distribution of loose materials in deep pyramidal silo.
showed rather than the dip angle of the silo walls.If the pyramidal silos were symmetrically designed, namely ϕ 1 = ϕ 3 and ϕ 2 = ϕ 4 , the λ might be expressed as Eq. ( 7) in the following: The vertical pressure shown in Eq. ( 5) must be disassembled to study the friction stress between loose materials and silo walls.In the current study, the lateral pressure coefficient k was an important constant when the vertical force was changed into the horizontal force.Furthermore, the lateral pressure coefficient k could be sampled with different methods prescribed by the code for the design of silos in different countries (SAG, 2001;ACI, 1997;Jaky J., 1948;CEN, 1995;DIN 1055DIN , 1987)).In this paper, the lateral pressure coefficient k was prescribed in Eq. ( 8) according to Rankine's earth pressure theory (SAG, 2001).
In Eq. ( 8), ϕ f is the internal friction angle of loose materials.
Taking a piece of triangular prism as an object for investigation (Fig. 3).We supposed that plane A was an internal wall, ϕ i (i = 1, 2, 3, 4) was the dip angle of plane A to the vertical plane, p was the vertical pressure, pA i was the pressure of p to A i , and pA i sin ϕ i was the disassembled pressure perpendicular to plane A. Similarly, the lateral pressure kp might be disassembled.The disassembled pressure of the perpendicular pressure to surface A was kpA j cos ϕ i .Therefore, the friction force on plane A was Eq. ( 9): Thus, the friction stress distribution on plane A could be indicated as Eq. ( 10):

Contact Analysis and Verification
Eq. ( 10) could be verified numerically using finite element contact analysis.

Contact conditions and iterative processes
Contact conditions established when two separated surfaces touched at normal direction while maintaining tangential relative movement.In general physical meaning, the surfaces in contact status had the following characteristics: (1) No penetration between each other.
(2) The normal pressure and tangent friction force might be transferred during contact.(3) Generally the normal pulling force could not be transferred when surface separation occurred.All the above three characteristics introduced significant physical variations of silo-loose material system in terms of normal and tangential stiffness.In fact, this highly nonlinear contact mechanics between silo and loose materials formed a major concern of silo engineers or civil engineers.This paper studied the contact friction stress.
Nonlinear contact problems could be solved through a series of repetitive iterative processes to accomplish accurate contact analysis.The procedure could be divided into four steps: (1) a particular contact status type was initially presumed; (2) finding out the contact stresses and displacement of contact point based on the presumed contact status; (3) checking if the contact condition could be met for all contact points; (4) if the contact condition was not satisfied, then the contact status was modified and another iteration started until all contact conditions could be met at all contacting points.

Contact analysis
Due to the symmetric property of silo structures, simplified two-dimensional contacting simulations were carried out in this paper, nonlinear finite element software ANSYS was used and the contacting surfaces between loose materials and silos were defined with rigid-to-flexible surface-to-surface contact pair.Target surfaces of silos were modeled with TARGE169 element and the contact surfaces of loose materials were modeled with CONTA171 element.
The ANSYS TARGE169 element could be used to rep-Fig. 3 The diagram of the vertical pressure disassembling.
resent various 2-D "target" surfaces together with the associated contact elements (CONTA171).The contact elements, CONTA171, were overlaid on top of the relevant solid elements in order to describe deformable loose material boundary and possible contact silo surface with target elements, TARGE169.This target surface was discretised by a set of target segment elements (TARGE169) and was paired with its associated contact surface via a shared real constant set.Any translational or rotational displacement or even temperature could be imposed on the target segment element.Moments and Forces could also be applied on target elements.CONTA171 was used to represent contact and sliding between 2-D "target" surfaces (TARGE169) and a deformable surface, defined by this element.The element was applicable to 2-D structural and coupled field contact analyses.This element was located on the surfaces of 2-D solid, shell, or beam elements without midside nodes (such as PLANE42).It had the same geometric characteristics as the attached solid, shell, or beam element.In ANSYS numerical model, contact occurred when surface elements penetrated into target segment elements (TARGE169) at allowable normal distance.Coulomb and shear stress friction could be used to simulate the tangential contact behavior.
Although no penetration occurred in real silos, numerical technique only allowed a stiffness-penetration law to simulate the contacting physics (see Fig. 4).Hence, penetration Δ must be small enough so that high accuracy could be obtained, which meant larger contact stiffness must be used.But large contact stiffness added simulation cost and convergence difficulties.To overcome this, contact stiffness in this paper was defined as the silo element stiffness multiplied by a normal contact stiffness factor.It was found that accurate results could be obtained when the normal contact stiffness factor was set to 1.0.It should be noted that this factor also indicated the contact stiffness equaled to the silo element stiffness.The maximum penetration was only δ = 0.001 m (see Table 1) when the normal contact stiffness factor was set to 1.0.Since only a brief summary of analysis method was given here, more detailed discussions and method could be found elsewhere (Wang X.W. et al., 2015).
During FE analysis, silos and loose materials were meshed with 2D solid element, PLANE42.Such element could be used either in 2D plane model (plane stress or plane strain) or axi-symmetric model.The element was defined by four corner nodes, each had two degrees of freedom: translations in the nodal x or y directions.The element could simulate plasticity, creep, swelling, stress stiffening behavior with large deflection and large strain capabilities.
The full Newton-Raphson solver was used in this study, due to its improved convergence efficiency for highly nonlinear problems, especially to the contact between granulate materials and silo surfaces.Moreover, ANSYS Newton-Raphson solver could automatically improve.Its convergence efficiency by tracing back system residual errors.
In FE contact analysis, correlation parameters on the loose materials and the silo were indicated in Table 2. Table 1 The contact parameters.

The item The parameters
The mesh element Plane42 The contact elements Targe169 and conta171 The contact form Rigid-to-flexible surfaceto-surface contact The normal contact rigid factor 1 The maximum penetration depth 0.001 m Table 2 Applied parameters of a pyramidal silo.

The results
Taking an engineering silo as an example (Fig. 5), the contact analysis were carried out.The analysis object is the narrow-wall of the engineering silo (Fig. 5).
Correlation parameters on the narrow-wall and the loose materials in the engineering silo were indicated in Table 2.The k and λ were determined by these parameters, where k = 0.49, λ = 1.49.The results were shown in Fig. 6, Fig. 7, and Table 3.In Fig. 6, different width rep-resents different friction zone (sliding and sticking).The curves of Fig. 7 shown the contact friction stress value and its location.

Modification
Fig. 6 indicated that the contact friction zone could be divided into two parts: sliding contact friction zone and sticking contact friction zone.The interval dividing point between the sliding zone and sticking zone was approximately in the middle.The location varied with different parameters.In Fig. 6, the location is located in the middle of node 103 and node 104.The contact friction stress reached the minimum (near 0) at the top of the pyramidal silo, and the stress gradually increased from the sliding contact friction zone to the interval dividing point (Fig. 7).The contact friction stress was at the max in the interval dividing point, and it gradually decreased from the interval dividing point to the bottom.
Comparing the FEA results value with the analytic value of the 11 nodes in Table 3 (node location shown in Fig. 6), it showed that above the interval dividing point, the FEA results was a little greater than the analytic solution for most of the nodes.The FEA results were about 1.1 times that of the analytic solution.However, the FEA results below the interval dividing point were about 0.95   times or less that of the analytic solution.Near the discharge opening (node No. 116), the FEA results were a little greater because of the dynamic discharge pressure.For dynamic discharge, more complex phenomena, such as material sticking or blocking in the form of a rat hole, arching and bridging were likely to happen.This paper did not discuss the complex phenomena near the discharge opening.As a whole, the analytic results roughly conformed to the FE contact calculation results.It had been pointed out that the analytic method was still reasonable.The error was in 10 % or less, because there were some differences between the assumptions and the real state.
Moreover, different states in the sliding contact friction zone and sticking contact friction zone had not been fully taken into account in the assumptions.In the sliding contact friction zone, the friction was introduced by external forces.The energy dissipation by friction was mainly due to the relative movement between the loose materials and the silo at the contact surfaces.In the sticking contact friction zone, the friction was internal friction.The friction energy dissipation was mainly due to the elastic deformation because of the contact.In general, the impact of external forces was a little greater than internal friction.
Taking the wide-wall of the engineering silo as analysis object, and changing attribute parameters of the loose materials, a similar contact friction stress curve and similar conclusion were obtained (Fig. 9).Thus, Eq. ( 10) might be modified.Based on Table 3 and Fig. 6, taking the discharge opening as a starting point, supposing the height of pyramidal silo was h, the modified formula of the combination of Eq. ( 5) and Eq. ( 10) were as follows: Here, y = h max when p is at the max, and ϕ i is the obliquity (see Fig. 2 and Fig. 3).
In addition, the force of loose materials in the discharge opening could not be regarded as the static force.In here, dynamic pressure would be acted upon during the unloading process of the loose materials.Eq. ( 11) would not be applied in the discharge opening where the friction stress distribution was more complex.Hence, further research would be needed, or the values near the discharge opening would increase by 1.5-2 times that of the results obtained by Eq. ( 11) in the practical engineering calculations.
The modified distribution curve (shown in Fig. 8: after modifying) was drawn according to Eq. ( 11).In fact, there was a mutation near h max .The curve in Fig. 8 was made a smooth processing.

Discussions and prospects
(1) For the purpose of engineering application, this study focused on the pyramidal silos which other researchers paid less attention.The lateral pressure coefficient k and the assumptions based on Janssen formula were also applied during the calculation to provide a convenient way for structure calculation of pyramidal silos and improve the design efficiency of field application engineers.Obviously, the stress conditions of pyramidal silos could have been some different from the cylindrical silo that Janssen formula applied on, and the mechanical  properties of loose materials were complex at the edge and corner part.Therefore, this study only did the approximate research on silo walls.More precise behavior and mechanism research of loose materials, and the friction stress distribution around the edges of the pyramid silos were beyond the scope of this study.
(2) Pyramidal silos are mainly used on the field of storing industrial bulk material (such as coal, gravel, concrete).According to the data obtained from field engineers, the main problems of this kind of silo were the poor flow during the unloading process and the wear of walls.To solve these problems, this study mainly researched on the friction stress distribution and the contact state between the loose materials and silo walls.For the normal pressure (or stress) distribution and the behavior study, the related papers were recommended (Goodey R.J. et al., 2017;Li Z.F., et al., 2014;Rotter J.M. et al., 2002;Brown C.J. et al., 2000).
(3) The slice method, based on the limit equilibrium theory, took the loose materials as the rigid body.When establishing equilibrium equations, some acceptable assumptions were used to increase the known conditions and transformed the statically indeterminate problem into the static problem.The complex boundary conditions of non-continuous, heterogeneity, anisotropy and external load of loose materials could be conveniently processed and finally derived the quantitative calculation model for field applications.
However, the slice method ignored the deformation and strength of the material itself and the complex behavior was simplified which would inevitably cause some error.In addition, the slice method could not provide accurate friction stress distribution graphics and images.The finite element analysis and software technology could make up for this defect and validate and correct the quantitative model.The combination of them could correct the calculation error to a certain extent.
(4) The finite element method was used in this study to simulate the contact and friction between the bulk material and the pyramidal silo.The discrete element method could be applied combining test researches to study the discrete characteristic of material, particle behavior and the influence mechanism of material on the silo in the future.

Conclusions
(1) The analytic equation of loose material friction stress distribution was established in deep pyramidal silos with reference to the basic assumptions of the JANSSEN Formula.
(2) The friction stress distribution of loose materials in deep pyramidal silos was studied with the finite element contact analysis.We concluded that the friction behaviors included sliding contact friction and sticking contact friction, and that the friction stress was maximum at the interval dividing point between the sliding zone and sticking zone.The friction stress decreased gradually to the roof and the bottom separately.
(3) The FEA results and the analytic result were compared and analyzed.The error was 10 % or less.The main reason was that some differences existed between the assumptions and the real state.Moreover, the different states in the sliding contact friction zone and sticking contact friction zone had not been fully taken into account in the assumptions.
(4) The analytic equation had been modified by comparing the results.It showed that the friction stress distribution of loose materials in deep pyramidal silos was related to the material and design parameters of silos, as well as the powder attribute of the loose material.The modified equation reflected precisely the contact-friction behaviors of loose materials in deep pyramidal silos and provided a shortcut method for calculation in related research.[see Eq. ( 11)]

Fig. 1
Fig.1The diagram of deep silo and shallow silo.

Fig. 4
Fig. 4 The penetration tolerance δ and the penetration Δ (a) In non-contact state; (b) In contact state.

Table 3
Results comparison.