Effects of Wall Inclinations and Wall Imperfections on Pressures during Granular Flow in Silos †

The effects of wall inclinations and wall imperfections at the onset of granular flow in a silo are numerically studied. The calculations were carried out for the onset of quasi-static mass flow in a silo with a controlled outlet velocity along the entire bottom. In the analysis, a finite element method and a polar hypoplastic constitutive law were used. A polar hypoplastic law describes the salient properties of granular bodies. FE calculations were performed with a plane strain silo with parallel walls, slightly convergent walls, slightly divergent walls, and parallel and convergent walls. The influence of a wall imperfection directed inwards and a wall imperfection directed outwards in a silo with parallel walls was also analysed. The FE results showed a large sensitivity of stresses in bulk solids due to the change of the direction of shear deformation along the silo wall.


Introduction
Experiments with bulk solids during silo flow show a large effect of the wall inclination and wall imperfections on wall stresses [1]- [7]. In general, wall stresses increase when silo walls become convergent, and they decrease when silo walls are divergent. A significant growth of wall stresses is usually obser ved at the transition between bin and hopper (the so-called "switch"). A wall imperfection directed to the inside of the bulk solid causes an increase of wall stresses. In turn, a wall imperfection directed to the outside of the bulk solid diminishes wall stresses.
The wall inclinations and wall imperfections inf luence wall stresses due to the fact that bulk solids are very sensitive to every change in the direction of shear deformation [8]- [11]. Thus, a realistic calculation of wall pressures due to the change of the wall inclination and the occurrence of imperfections in a silo is possible using a constitutive law taking into account this important granular property.
The intention of the paper is to analyse numerically the effect of the wall inclination and wall imperfections on the evolution of wall stresses in granular materials at the beginning of confined granular f low in a plane strain model silo with very rough walls. The calculations were carried out for quasi-static mass f low in a silo with a constant outlet velocity along the entire bottom (inertial forces were neglected). In a theoretical study, a finite element method on the basis of a polar hypoplastic constitutive law was used. The constitutive law can describe the salient properties of granular bodies during deformation with allowances made for shear localisation [12]- [15]. The effect of density, pressure, deformation direction, mean particle diameter and particle roughness is also captured. The FE calculations were carried out within a silo research study aimed at determining an optimal shape and roughness of silo walls to reduce an increase of wall stresses after outlet opening. In the paper, compressive stresses were considered negative (as in soil mechanics).

Constitutive law for granular bodies
A polar hypoplastic law [12]- [15] is an alternative to a polar elasto-plastic formulation [16]- [17] for continuum modelling of granular materials. It was formulated within a polar (Cosserat) continuum [18]. A polar (Cosserat) continuum differs from a non-polar one by additional rotations. For plane strain or axial symmetry, each material point has three degrees of freedom: two translations u 1 and u 2 , and one rotation ω c . The polar rotation ω c is related with the microrotation and is not determined by displacements as in a non-polar continuum The gradients of the polar rotation are connected with curvatures which are associated with couple stresses. The deformation and stress tensor become non-symmetric. As a consequence of the presence of rotations and couple stresses, the constitutive equation is endowed with a characteristic length corresponding to the mean particle diameter. Thus, numerical results are independent of the spatial discretisation, and boundary value problems remain mathematically well posed. Due to the presence of a characteristic length, a polar approach can model the thickness of shear zones (which are inherent characteristics of granular flow in silos) and related particle size effects. The polar hypoplastic law can reproduce essential features of granular bodies. It takes into account such important properties as: dependence on pressure level, on density, and on the direction of deformation rate, dilatancy and contractancy during shearing with constant pressure, and increase and release of pressure during shearing with constant volume. The constitutive law is characterised by simplicity and a wide range of applications. The material constants can be found by means of standard element tests and simple index tests. They are simply correlated with particle properties. Thus, they can be estimated from granulometric properties (encompassing particle size distribution curve, shape, angularity and hardness of particles) [11], [19]- [21]. The capability of this law has been already demonstrated in solving boundary value problems involving localisation such as biaxial tests, monotonous and cyclic shearing of a narrow granular layer, trap door, silo filling, silo f low, furnace f low, footings, retaining walls and sand anchors. A close agreement between calculations and experiments was achieved. The FE calculations showed also that the thickness of shear zones did not depend upon the mesh discretisation if the size of finite elements in the shear zone was not more than five times the mean particle diameter when using triangular finite elements with linear shape functions for displacements and a Cosserat rotation [22].
Stress changes and couple stress changes during the plane strain deformation of granular bodies are expressed by wherein σ ij and m i are the Cauchy stress tensor and the Cauchy couple stress vector, respectively, e denotes the void ratio, the Jaumann stress rate tensor, and the Jaumann couple stress rate vector, d kl c the polar rate of deformation, k k the rate-of-curvatures vector and d 50 the mean particle diameter. The following representations of the general constitutive equations (Eqs. 2 and 3) are used: wherein the normalised stress tensor ŝ ij and the normalised couple stress vector m i are defined by The scalar factors f s ҃f s (e, σ kk ) and f d ҃f d (e, σ kk ) take into account the influence of the density and pressure level on the stress and the couple stress rates. The stiffness factor f s is proportional to the granulate hardness h s and depends on the mean stress and void ratio: with The constant c 1 is defined by Eq. 16. The granulate hardness h s represents a density-independent reference pressure and is related to the entire skeleton. It should not be confused with the particle hardness which refers to single particles. The density factor f d , a kind of a pressure-dependent relative density index, is represented by Here e is the current void ratio, e c is the critical void ratio, e d denotes the void ratio at maximum densification due to cyclic shearing, e i is the maximum void ratio, and α and n are constants. The void ratio e is thus bounded by e i and e d . The values of e i , e d and e c are assumed to decrease with pressure Ҁσ kk according to the equations [21]: wherein e i0 , e d0 and e c0 are the values of e i , e d and e c for σ kk ҃0, respectively. For the tensor and vector s ij s kk functions L ij , L i c , N ij and N i c , the following representations are used where φ c is the critical (residual) angle of internal friction, θ denotes the Lode angle, i.e. the angle on the deviatoric plane σ 1 ѿσ 2 ѿσ 3 ҃0 between the stress vector and the axis σ 3 (σ i is a principal stress component).
The coefficient a 1 Ҁ1 determining the shape of the stationary stress surface lies empirically in the range 3 to 4.5. The dimensionless polar constant a c controls the inf luence of the Cosserat quantities on the material behaviour. It lies in the range of 1.0-5.0 and is correlated with the particle roughness; the higher the constant a c , the smaller are the polar effects. σ ij denotes the deviatoric part of σ ij . For an isotropic stress state with σ ij ҃0, the expressions in Eq. 15 become cos(3θ)҃0 and a 1 Ҁ1 ҃c 1 . The first terms are linear and the second terms in Eqs. 4 and 5 are nonlinear in d kl c and kd 50 , respectively. The second terms in Eqs. 4 and 5 describe the modulus of the deformation. Thus, a polar hypoplastic law depends on the direction of deformation.
The FE analyses of silo flow for the case of plane strain were carried out with the following nine material constants (for so-called cohesionless Karlsruhe sand [21]): e i0 ҃1.3, e d0 ҃0.51, e c0 ҃0.82, φ c ҃30°, h s ҃190 MPa, α҃0.3, n҃0.5, d 50 ҃0.5 mm and a c ҃a 1 Ҁ1 . The parameters h s and n are determined from a single oedometric compression test with an initially loose specimen, h s ref lects the slope of the curve in a semi-logarithmic representation, and n its curvature. The constant α is found from a triaxial test with a dense specimen. It ref lects the height and position of the peak value of the stress-strain curve. The angle φ c is estimated from the angle of repose or measured in a triaxial test with a loose specimen. The values of e i0 , e d0 , e c0 and d 50 can be obtained with index tests (e c0 Ȃe max , e d0 Ȃe min , e i0 Ȃ(1.1-1.5)e max ).

Finite element data
The plane strain calculations of granular flow were performed with a bin with a height of h҃0.50 m and a width of b҃0.20 m (h/b҃2.5, b/d 50 ҃400) [23]. The FE mesh with quadrilateral finite elements composed of four diagonally crossed triangles was applied to avoid volumetric locking and spurious element behaviour. A total of 3000 triangular elements with linear shape functions for the displacements and the Cosserat rotation were used. Symmetry with respect to the centre line was taken into account. In order to realistically describe the interface behaviour along the wall, the FE mesh was strongly refined at the wall. The width of three quadrilateral elements close to the wall was equal to 0.5 mm, 1.5 mm, 3 mm and the width of the next five elements was 5 mm, respectively. The height of all quadrilateral elements was 10 mm. Quasi-static flow (without inertial forces) was initiated through constant vertical displacement increments prescribed to all nodes along a smooth bottom.
Thus, mass f low was obtained in the silo.
The FE calculations were carried out mainly with an initially dense solid (low initial void ratio) between very rough walls (r w ͧd 50 ). In this case, the largest stress increments occur during silo flow after bottom displacement [22]- [24]. The wall roughness r w is depicted as a relative height between the highest and the lowest point along the surface at a length of (2-3)҂d 50 [22]. For modelling of very rough walls in silos, full shearing of the material along a wall was assumed [22]: As a result of the assumed polar boundary conditions along the wall (Eq. 17), the wall friction angle can be derived and no special interface elements are needed [22]. The remaining boundary conditions in the silo fill were along the top traction and momentfree, and in the symmetry axis: u 1 ҃0, ω c ҃0 and σ 21 ҃0 (σ 21 Ҁ vertical shear stress).
To investigate the influence of the wall inclination and wall imperfection, the analyses were carried out with the same initial stress state, namely with the socalled K o state (frequently used in soil mechanics): σ 22 ҃γ x 2 , σ 11 ҃σ 33 ҃K 0 γ x 2 , σ 12 ҃σ 21 ҃m 1 ҃m 2 ҃0 (σ 11 Ҁ horizontal normal stress, σ 22 Ҁ vertical normal stress, σ 33 Ҁ normal stress perpendicular to the plane of deformation, σ 12 Ҁ horizontal shear stress, σ 21 Ҁ vertical shear stress, m 1 Ҁ horizontal couple stress, m 2 Ҁ vertical couple stress, γ Ҁ initial unit weight of sand, x 2 Ҁ vertical coordinate measured from the top of the solid, K 0 ҃0.40 Ҁ pressure coefficient at rest). To study the effect of the initial stress state on wall stresses, a comparative calculation was also performed with initial stresses after filling according to a slice method by Janssen [6]: where N K is the mean pressure coefficient, ϕ w denotes the mean wall friction angle, x 1 is the horizontal coordinate measured from the wall, b denotes the silo width, and σ 33 is the normal stress perpendicular to the plane of deformation. Both model experiments [22] and numerical calculations [12] show that stresses in silos during filling can be approximated with a slice method provided that N K and ϕ w are empirically estimated. The choice was made on the basis of model tests [22]: N K҃0.29 and ϕ w ҃43°(initially dense sand between parallel and very rough walls).
Four different wall inclinations in a silo were assumed: parallel walls (Fig. 1a), slightly convergent walls with a wall inclination from the bottom equal to Ҁ2N Ktanj w x 2 b gb 2N Ktanj w α҃88° (Fig. 1b), slightly divergent walls with a wall inclination from the bottom equal to α҃88° (Fig. 1c), and parallel walls (h҃0.25 m) and slightly convergent walls (h҃0.25 m) with a wall inclination from the bottom equal to α҃88° (Fig. 1d).
Two different wall imperfections at the mid-height of the wall (h҃0.25 m) were taken into account: a small imperfection directed to the inside of the bulk solid (Fig. 1e) and a small imperfection directed to the outside of the bulk solid (Fig. 1f ). The height of both imperfections was h 1 ҃40 mm (h 1 /b҃0.2) and the width b 1 ҃2.5 mm (b 1 /b҃0.0125).
The calculations were carried out with large deformations and curvatures. In this case, an Updated Lagrangian formulation was applied using both the Jaumann stress rate and the Jaumann couple stress rate [22]. At the same time, the changes of the element configuration and the element volume were taken into account.
The integration was performed in 3 sample points placed in the middle of each triangular element side. The density of the silo fill was kept constant. For the solution of the non-linear equation system, a modified Newton-Raphson scheme with line search was applied using an initial global stiffness matrix calculated with only two first terms of the constitutive equations (Eqs. 4 and 5). To accelerate the calculations in the softening regime, the initial increments of displacements and the Cosserat rotation in each calculation step were assumed to be equal to the converged increments from the previous step. In addition, to prevent possible inadmissible stress states (in particular in the corner of the wall at the bottom), a sub-stepping algorithm was used (deformation and curvature increments were divided into smaller parts within each step). The iteration steps were performed using translational and rotational convergence criteria (found by means of preliminary FE calculations). For the time integration of stresses and couple stresses in finite elements, a one-step Euler forward scheme was applied.

Numerical results
Influence of initial void ratio and initial stress state Figure 2 shows the evolution of normalised hori-zontal normal stresses σ 11 /(γb) at the mid-height of the wall (h҃0.20-0.30 m) during the normalised vertical bottom displacement u 2 /h in a model silo with parallel walls of Figure 1a at the onset of quasi-static granular f low. The onset of silo f low is crucial with regard to maximum wall stresses [22]- [24]. Later, wall stresses have a tendency to oscillate with wandering small peaks due to the appearance of shear zones along a wall and inside of the flowing material [22]- [24]. The presented range of u 2 /h for dense sand (e o ҃0.60) and loose sand (e o ҃0.90) in Figure 2 is different.
The shear stresses were not taken into account during filling (K o Ҁ initial stress state), thus wall stresses decreased at the beginning of flow (Fig. 2). The wall element '21' is located at a height of 0.20 m and the wall element '30' at a height of 0.30 m above the bottom (Fig. 2). Other wall elements are located between '21' and '30' at a distance of 10 mm. The wall stresses in dense sand decrease first, then they increase and reach an asymptote or they reach an asymptote. The wall stresses in loose sand decrease and very quickly reach a residual state (u 2 /h҃0.0005). The maximum normalised horizontal stresses σ11/(γb) in the middle of the wall are about 0.70 (dense sand) and 0.55 (loose sand). The stresses in loose sand increase with increasing depth. However, the stresses in dense sand are not distributed proportionally with a silo depth due to interior shear zones inside the flowing solid [22]- [24]. The maximum horizontal wall stresses for dense sand are higher by about 30% than those calculated by the German Silo Code (plane strain, dense sand, very rough walls): σ11/(γb)҃0.55 kPa (with γ҃17.0 kN/m 3 , N K҃0.70, tanϕ w ҃0.60, b҃0.20 m, h҃0.25 m). In large silos, the normalised wall stresses become smaller due to a decrease of d 50 /b and an increase of pressure (reduction of the internal friction angle, dilatancy angle and stiffness) [23], [24].
The effect of the initial stress state on the evolution of stresses σ 11 /(γb) in a silo of Figure 1a is demonstrated in Figure 3. The calculations were performed with dense sand and initial stresses according to a slice method by Janssen (Eqs. [18][19][20][21]. The results show that the initial stress state inf luences the horizontal wall stresses only up to u 2 /h҃0.0005. Later, the evolution of stresses is the same, independent of the initial stress state. The maximum normalised stress increment after outlet opening is about ∆σ11/(γb) ҃0.40 (dense sand). The stresses during f low increase more than twice when compared to silo filling.

Influence of wall inclination
The effect of the wall inclination on horizontal wall stresses during f low (at different scales of u 2 /h) is shown in Figures 4 and 5 (K o Ҁ initial stress state). Even a small change of the wall inclination influences the stresses. The maximum normalised horizontal normal stresses on the wall σ11/(γb) (Fig. 4) in the middle region are 8% higher in a silo with slightly convergent walls of Figure 1b, and by 2% lower in a silo with slightly divergent walls of Figure 1c compared to the silo with parallel walls (Fig. 2a). The increase of horizontal normal stresses (Fig. 5) below the transition between a parallel and convergent part of the silo of Figure 1d is similar to that in the silo with convergent walls at mid-height (Fig. 4a).

Influence of wall imperfections
The evolution of normalised horizontal normal stresses σ 11 /(γb) during vertical bottom displacement in the middle region of the wall in a model silo with parallel imperfect walls (dense sand with e o ҃0.60, K o Ҁ initial stress state) is depicted in Figure 6.
The FE results demonstrate that a small imperfection in the wall increases (by about 40%) or decreases (by about 30-50%) the maximum horizontal wall stresses compared to a perfect wall (Fig. 2a). The stresses increase significantly on the convergent part of the imperfection (α҃75°) and decrease significantly on the divergent part of the imperfection (α҃75°). They also increase directly above the divergent part. Other stresses at imperfections are approximately in the range of stresses for a silo with perfect walls.

Conclusions
The following conclusions can be drawn from FE calculations at the onset of granular flow in a silo with very rough walls using a polar hypoplastic law: Bulk solids during silo f low are sensitive to every change in the direction of shear deformation.
The effect of the initial void ratio has a large inf luence on horizontal normal stresses on the wall. The wall stresses in a dense solid are larger during f low than in a loose one.
The initial stress state has an insignificant effect on maximum wall stresses.
Horizontal wall stresses decrease along divergent walls and increase along convergent ones.
Even small imperfections in a silo wall can cause a significant increase of local horizontal stresses.
The wall stresses during f low of dense solids in silos with very rough walls can be larger during f low than standard values.
The FE calculations will be continued. The inf luence of the wall roughness, the initial void ratio in the solid, and the magnitude and shape of wall imperfections will be analysed. To realistically describe the localisation of shear zones inside a flowing solid, the calculations will be performed for the entire silo. The initial void ratio will be stochastically distributed in the solid with a random generator.