# Stresses in silos: Comparison between theoretical models and new experiments

## Abstract

We present precise and reproducible mean pressure measurements at the bottom of a cylindrical granular column. If a constant overload is added, the pressure is linear in overload and nonmonotonic in the column height. The results are quantitatively consistent with a local, linear relation between stress components, as was recently proposed by some of us. They contradict the simplest classical (Janssen) approximation, and may pose a rather severe test of competing models.

PACS numbers: 46.10.+z, 83.70.Fn

The prediction of static stresses in dry, cohesionless granular matter has become the focus of renewed attention (see ). Surprisingly, there is no consensus on what is the basic physics involved. Some argue that the behavior is essentially elastic (ultimately justified by the slight elastic deformation of individual grains); others that it is dominated by the extremely nonlinear constraint that tensile intergranular forces are absent . Indeed, some of us  have argued that the statics of granular materials can be described, without considering elastic displacements, by assuming a local, history-dependent, relation between stress tensor components. This gives hyperbolic equations for the stress field, in contrast to the elliptic (or elliptic-hyperbolic) equations of conventional elastic (or elastoplastic) models. Our approach provides a simple continuum model of ‘force chains’ ; (physical) force chains become (mathematical) characteristics of the hyperbolic equations. In the simplest case, these form a regular array; stresses propagate through space via a wave equation . According to the model, the medium is ‘fragile’ in a precise sense : it responds linearly to a specific class of ‘compatible’ loads; all others cause plastic reorganization.

This approach accounts well  for the pressure ‘dip’ below the apex of a conical sandpile poured from a point source . (It also predicts  that the dip is absent for a pile made of successive horizontal layers, as recently confirmed by experiment .) However, it has excited strong criticism in some quarters , and certainly demands further experimental test . For example, such models predict that if a small localized overload is placed on top of a granular layer, the excess weight at the bottom is maximal, not directly beneath the weight, but on a ring . To test this directly is difficult, because of strong nonlinearity and (especially) noise effects which hinder the interpretation of data .

A more robust and practical situation, is the cylindrical granular column, or bin. Here also noise effects come into play; but ways around these (by careful ensemble averaging of experimental data) have been pioneered in . Below we report precise measurements (beyond those of ) of the effective mass $`M_e`$, supported by the bottom plate, as a function of the total mass $`M_t`$ poured into a (small) bin, with and without an added overload. With no overload, as expected, $`M_e(M_t)`$ first rises linearly, then saturates at a column height comparable to its width; for high bins, most of the mass is ‘screened’ by frictional transfer to the walls. A simple hyperbolic model (called osl for ‘oriented stress linearity’ ) gives bin results close to, but different from, the classical Janssen approximation (recalled below) . In contrast to traditional methodologies  our new ensemble-averaged experiments can distinguish these predictions; we find that osl, which has an extra fitting parameter, is discernibly better. Another classical model (ife, see below) gives wholly inadequate answers unless unphysical values of the wall and bulk friction constants are used.

There then follow, from the osl model, two important new predictions for the effect of a uniform overload of mass $`Q`$ at the top of the granular column. First, $`M_e`$ should be linear in $`Q`$; second, for large $`Q`$, $`M_e`$ should be nonmonotonic in $`M_t`$. We find that, with no further fitting, our overload experiments quantitatively confirm the osl predictions, strongly supporting the hyperbolic picture. At the end of this Letter, we comment on the challenge these new results pose to other modelling strategies.

First we recall our own approach. By stress continuity,

$$_i\sigma _{ij}=\rho g_j$$
(1)   
where $`\sigma _{ij}`$ is the (symmetric) stress tensor, $`\rho `$ is the density of the material, and $`g_j`$ is the gravitational acceleration. In general one needs extra physical assumptions to close Eq.1. For an elastic body, one assumes a (single-valued) displacement field, and a linear relation between stresses and strains (Hooke’s law). For poured cohesionless grains, the definition of a macroscopic displacement is problematic (see ). Instead we assume that the arrangement of granular contacts gives, on continuum length scales, a definite relation between components of the stress tensor . One such relation, often used in the literature, is the ife (‘incipient failure everywhere’) assumption: that the material is everywhere on the verge of Coulombic failure (see e.g. ). Then there exists a (locally varying) set of axes $`𝐧𝐦`$ such that $`\sigma _{nm}=\sigma _{nn}\mathrm{tan}\varphi `$ where $`\varphi `$ is the Coulomb angle.

Our modelling strategy instead gives a fundamental rôle to the network of force chains which, if grains are undeformable, must carry forces longitudinally . One interpretation of our equations is that the friction between parallel force chains is fully mobilized; a Coulomb-like condition, $`\sigma _{nm}=\sigma _{nn}\mathrm{tan}\psi `$, then holds (with $`\psi \varphi `$ an ‘effective’ friction angle) but the orientation $`𝐦`$, which is directed along the force chains, is now fixed by the construction history and not (as in ife) by the load . (This assumes the load is a compatible one.) For simple construction histories, like piles and bins, we assume that $`𝐦`$ is the same everywhere, up to an inversion through the central symmetry axis; $`𝐦`$ must then have a fixed angle $`\tau `$ to the vertical. In cylindrical polars ($`z,r,\theta `$) with $`z`$ downwards, we recover the osl model :

$$\sigma _{rr}=\eta _1\sigma _{zz}+\eta _2\sigma _{rz}$$
(2)   
with $`\eta _1=\mathrm{tan}\tau \mathrm{cot}(\tau \psi )`$ and $`\eta _2=\mathrm{tan}\tau \mathrm{cot}(\tau \psi )`$. Eq.2 closes the problem in two dimensions ($`d=2`$): inserting it into Eq.1, gives an anisotropic wave equation, with one characteristic along $`𝐦`$, and another along a direction $`𝐦^{}`$ at angle $`\tau \psi \pi /2`$. (These can be interchanged without affecting Eq.2; so $`𝐦^{}`$ describes a second family of force chains .) For $`d=3`$, a further closure equation is needed. Our choice here is $`\sigma _{rr}=\sigma _{\theta \theta }`$; but from work on conical piles, we expect insensitivity to this choice . In the bin geometry, the osl model can then be solved exactly ($`d=2`$)  or numerically ($`d=3`$). Note that ife, like osl, gives propagative (hyperbolic) equations; but these are nonlinear, unlike our wave equation.

For nonzero $`\eta _2`$, the force chain network distinguishes between inward and outward radial directions. This does not contradict the axial symmetry present . But if as well the medium is locally symmetric, then $`\eta _2=0`$; in Eq.2, this recovers the model of Ref.. The latter can be viewed as a local version of the classical Janssen hypothesis . Janssen proposed a constant ratio between horizontal and vertical stresses, $`\sigma _{rr}=K\sigma _{zz}`$, but neglected altogether their dependence on $`r`$. Assuming also that friction at the wall is fully mobilized, with a friction coefficient $`\mathrm{tan}\varphi _w`$, he found the equation:

$$M_e=M_{\mathrm{}}\left(1\mathrm{exp}\left[M_t/M_{\mathrm{}}\right]\right)$$
(3)   
with $`M_{\mathrm{}}=\rho D^2/2K\mathrm{tan}\varphi _w`$ for $`d=2`$, and $`M_{\mathrm{}}=\rho \pi D^3/16K\mathrm{tan}\varphi _w`$ for $`d=3`$; $`D`$ is the bin diameter.

We turn now to the experimental procedure, described in detail in . The bin is a tube of diameter $`D=3.8`$ cm, filled with beads of glass (density $`\rho _b=2.6`$ g/cm<sup>3</sup>, diameter $`2`$ mm). The bottom comprises a very stiff scale plate ($`2\times 10^4`$ N/m). Initially, the tube is filled with a low packing density; this is increased by giving it small taps. The bottom plate is then lowered (by a few tens of microns) and the effective mass decreases monotonically to an asymptotic value; $`M_e`$ and the mean density $`\rho `$ are measured. The density is again increased by tapping, the plate lowered and further measurements taken. This entire procedure is done about $`30`$ times – each run giving results for the whole range of densities. The measured results for $`M_e`$ show a certain ($`M_t`$-dependent) ‘error bar’: not a measurement error of the mass, but arising from intrinsic fluctuations in the packing. This protocol is a major advance because (a) an ensemble average value for $`M_e`$ is found, improving accuracy; (b) due to the downward motion of the base, that wall friction is fully mobilized, which might not otherwise be the case . The wall friction angle is measured separately as $`\varphi _w=22^o\pm 2^o`$ , thus eliminating one fit parameter.

The experimental results, for a packing density $`\rho =1.53`$ g/cm<sup>3</sup>, are compared in Fig.1 with three models: ife (which has no adjustable parameter once the internal friction angle $`\varphi =25^o\pm 2^o`$ is known); Janssen’s equation (one adjustable parameter); and the osl model (two adjustable parameters). Each plotted datapoint is itself a mean value, with an error bar $`\mathrm{\Delta }`$ shown in inset (a). (This is small at small $`M_e`$ but then grows rapidly.) To find the best fits, we have minimized the following:

$$E^2=N^1\underset{i}{}(\delta M_e^i/\mathrm{\Delta }^i)^2$$
(4)   
where $`\delta M_e^i`$ is the difference between the $`i`$th experimental datapoint and the theoretical $`M_e`$ value, $`\mathrm{\Delta }^i`$ the observed error bar, and $`N`$ the number of datapoints.

For our data, the (active) ife approach, using the measured friction values $`\varphi `$ and $`\varphi _w`$ is plainly inadequate. Better agreement with ife is found by taking $`\varphi `$ and/or $`\varphi _w`$ as fit parameters. Even then, the fit remains poor (e.g. $`E=4.43`$ for $`\rho =1.53`$ g/cm<sup>3</sup>) ; and the fitted values, $`\varphi =\varphi _w=30^o`$ are incompatible with those found by direct experiment. For given $`\varphi _w`$, ife systematically overpredicts the asymptotic stress; so the fitted $`\varphi _w`$ exceeds the real one. In systems where the wall friction is not fully mobilized, the error is harmlessly absorbed by the fit. In our system, the fitted value is higher than the fully mobilized $`\varphi _w`$ measured separately, which is unphysical.

Unlike the ife model, Janssen’s model gives a fair approximation ($`E2`$; Table 1) but, as shown in inset (b), there is a clear systematic deviation: screening by the walls is in turn over- and underestimated for small and large $`M_t`$ values. (Note also that our $`K`$ parameters are higher than those usually reported : but as with ife, low fitted values might compensate for incompletely mobilized of wall friction.) This has led two of us  to propose elsewhere an empirical model (not shown) where an excess contribution from grains at the bottom of the pile is added to the Janssen result. As shown in Table 1, the best-fit osl model does as well as this empirical model, with an error $`E1`$ : the systematic deviations are reduced, in particular in the first part of the curve. This can be understood by noting that within the osl model, the grains contained within a ‘light-cone’, resting on the bottom plate, cannot interact with the walls ; the mass of these grains is completely unscreened.

Note the values found for $`\eta _2`$. The minimum of $`E(\eta _2)`$ is not sharp, but positive $`\eta _2`$ is always preferred (as for other types of grains ). A positive $`\eta _2`$ means that most of the weight follows the ‘inward’ characteristic thus reducing the screening effect of the walls. Conversely, in sandpiles (created from a point source) $`\eta _2`$ is negative ; this ‘outward’ transfer of weight is responsible for the pressure dip underneath the apex. Positive $`\eta _2`$ could be caused by slight inward avalanches of material as the base is lowered. Its decrease at higher densities might indicate a diminished susceptibility to this effect; alternatively the tapping procedure could progressively erase a local assymmetry induced by the initial fill.

We now turn to the key results of this paper, for the response to an overload $`Q`$ placed on top of the granular column. (This is a solid piston, just narrower than the cylinder.) This is taken into account within the osl model by modifying the boundary conditions to include a uniform downward stress at the top surface. Such a load is found to be compatible. From the linearity of the osl model (also true of Janssen’s model) we then have:

$$M_e=M_{\mathrm{}}f_0\left(\frac{M_t}{M_{\mathrm{}}}\right)+Qf_Q\left(\frac{M_t}{M_{\mathrm{}}}\right)$$
(5)   
In Janssen’s description, $`f_0(x)=1e^x`$ and $`f_Q(x)=e^x`$, so $`M_e`$ is monotonic in the poured mass $`M_t`$ (and constant when $`Q=M_{\mathrm{}}`$). The result of the osl model are more surprising: $`f_0(x)`$ and $`f_Q(x)`$ have different $`x`$ dependences. Hence $`M_e`$ is not monotonic in $`M_t`$; at intermediate $`Q`$ there is an ‘overshoot’ (Fig.2). In addition, both functions have a (slight) oscillatory character, caused by ‘resonances’: these are standing-wave modes of the wave equation, damped by ‘absorption’ arising from wall friction (see  and ). In Fig.2, we show the experimental results obtained for various overloads $`Q`$. As shown in the inset, these results do indeed obey the linear relation, Eq.5, to good accuracy. A clear overshoot effect is also seen, although any further ‘resonant’ oscillations are small (even theoretically). Note that the osl predictions in Fig.2 use the same parameters as determined previously for $`Q=0`$. Thus osl, with no further fitting, gives a good quantitative account of the data for all $`Q`$.

We have shown that simple hyperbolic models , encoding the presence of linear force chains , can be used to reproduce quantitatively the observed stress response of cohesionless granular media, not only in piles , but in bins. The same is not true of the traditional Janssen analysis. Nor is it true of ife; this does predict resonant behaviour (at least in local stresses ), but our results, even without overload, rule it out entirely as a physical model. Any expectation of nonmonotonicity in Fig.2 based on ife would thus have been misplaced.

What of other continuum modelling strategies? Much recent work on bins and silos has studied elastoplastic constitutive models (also widespread in soil mechanics), often by a finite-element method. There are many such models, and a recent comparative study found little consensus among them . But we wonder whether these approaches can, with reasonably few fit parameters, reproduce the results of Figs.1 and 2. For example, the observed linearity in $`Q`$ (seen even for $`Q/M_{\mathrm{}}1`$) may set a challenge, although one finds numerically that, after summing stresses over the base, the (non-linear) ife model obeys to a good precision the linear relation (5). Linearity is, of course, also recovered if the material is entirely Hookean. The challenge is then to explain within a purely elastic theory the nonmonotonic (if not oscillatory) curves of Fig.2. The investigation of these important questions is underway .

Finally it is important to map out more clearly the domain of validity of the hyperbolic approach (see e.g. ). In particular, our granular columns are tiny: only twenty grains or so across. These data clearly do not rule out a crossover to more conventional elastic or elastoplastic behavior at larger scales (e.g. where the grains start to deform) , although the hyperbolic approach also works well in conical piles up to 1 metre wide . Careful overload experiments on larger bins could be very valuable, as well as local stress measurements, which would reveal more clearly the oscillatory nature of the response.

We thank P.G. de Gennes, J. N. Roux and G. Combe for very useful discussions. E.C. and L.V. thank J. Lanuza for technical assistance.