Bending back stress chains and unique behaviour of granular matter in cylindrical geometries

We analyse the general solutions for the stress field in planar annuli of isostatic media, a model often used for marginally rigid granular materials in Couette cells. We demonstrate that these solutions are much richer than in rectangular symmetries. Even for uniform media, stress chains are found to curve, broaden away from the stress source, attenuate and leak stress into a cone of influence. Most spectacularly, stress chains may bend back and transmit forces oppositely to the original direction. None of these phenomena arises in solutions for uniform media in Cartesian coordinates. We further analyse non-uniform media, which exhibit chain branching and stress leakage from the chains. These results are directly relevant to the many experiments on granular materials, carried out in Couette cells. They also shed light on, and are supported by, hitherto unexplained experimental observations of curved and back-bending chains, which we point out. In particular, we use our results to provide a new interpretation for the pattern of slip lines observed experimentally.


Introduction
The ubiquity and significance of granular matter (GM) have focused scientific and technological attention for millenia [1], but the theoretical understanding of this form of matter is far B Raphael Blumenfeld rbb11@cam.ac.uk 1 Imperial College, London, UK 2 College of Science, NUDT, Changsha, Hunan, China 3 Cavendish Laboratory, Cambridge, UK from complete. One of the most important theoretical challenges is the development of a theory of stress transmission in dry such media. Predicting stress is essential to a wide range of technological and geo-mechanical applications, as well as being a springboard for modelling the dynamic behaviour of GM. Experimental observations of nonuniform stress transmission in GM, e.g. via 'force chains', date back to the 1940s [2][3][4], with modern experiments revealing detailed features of this phenomenon [5][6][7]. The limitations of conventional theories, such as elasticity, to account for force chains led to investigations of ideal systems: isostatic media. These are statically determinate and marginally rigid aggregates, characterised by a low mean contact (or coordination) number, whose intergranular forces can be determined, in principle, from balance conditions alone. The stress equations of isostaticity theory are hyperbolic, differing markedly from the conventional elliptic equation of elasticity theory [8][9][10][11][12][13][14][15][16]. Understanding these ideal systems is an essential step to a theory of real GM, which comprise both isostatic regions and denser regions, where conventional theories are valid [16].
Many physical and numerical experiments are carried out in cylindrical setups [17][18][19][20], yet most theoretical analyses of stresses in isostatic materials are based on rectangular coordinate systems. This practice is based on an implicit, normally little discussed, assumption that the stress chains phenomenon is independent of the system symmetry. The purpose of this paper is to show that, at least in two dimensions (2D), this is not the case, namely, that isostatic stress solutions in cylindrical symmetry have features that do not arise in rectangular coordinates. This we do by solving explicitly for the stress field in annuli and then demonstrating these solutions by examples. Specifically, we show that uniform isostatic media exhibit stress chains that may curve, broaden, dissipate and even bend backwards. The characteristics of these solutions are analysed and their implications discussed. Experimental support in the literature for some of these phenomena are pointed out and interpreted in view of these solutions.

The isostaticity stress equations
In 2D Cartesian coordinates, the isostaticity stress field Eqs. are [11][12][13][14][15] where σ is the stress tensor, g includes all external and body forces and p i j are the components of a fabric tensor, P(x, y), which can be determined directly from the grain structure [15] and upscaled to the continuum [21]. It has been shown that the determinant of P is generically negative [16] and therefore that these equations are hyperbolic and yield stress chain solutions, [11][12][13][14][15]. Equations (1)- (3) have been analysed and solved for uniform fabric tensors [10][11][12][13][14]16], as well as for position dependent ones [22,23]. Yet, in spite of numerous experiments in annuli, e.g. Couette cells, exact solutions in the literature are often derived in rectangular coordinates, presumably under the assumption that the behaviour would not differ significantly in cylindrical geometries. This is not the case. Moreover, since force chains often originate in very localised regions, down to single grains, it is more appropriate to describe the field around such an origin in polar coordinates. In polar coordinates, the balance equations are The equations are closed by a stress-structure relation that, for consistency, has the same form as (3):

Analysis
Rewriting the stress components of Eq. (6) in Cartesian coordinates and comparing to (3), a relation can be derived between the Cartesian, P, and polar, Π , fabric tensors: where, for brevity, (S, C) ≡ (sin 2θ, cos 2θ). The condition that det{P} < 0 translates to: From this condition we can determine a region in the constitutive parameter space π rr -π θθ -π r θ , in which the determinant of Π is negative for all values of θ . Note that this eliminates a wide range of fabric tensors, some of which were studied in the literature, e.g. in [11][12][13][14]22,23].
To obtain the general solution for the stress, we follow a similar procedure to the one initiated in [22,23]. Assuming first that π θθ = 0, we define q i j ≡ π i j /π θθ and substitute σ rr from the stress-structure condition (6) into Eq. (4). The resulting balance equations can be written as The characteristic variables, w 1,2 , can be expressed in terms of the stress components as where λ 1,2 = q r θ ± q 2 r θ − q rr are the eigenvalues of A and Y −1 AY = Λ is diagonal. Since q 2 r θ − q rr = −det {Π } /π 2 θθ > 0, λ 1,2 are real and distinct. In terms of w, (9) becomes where h = Y −1 g.
It is now convenient to parameterise the characteristic paths by length variables, s i : ∂ ρ s i = 1/λ i and ∂ θ s i = 1 (i = 1, 2). For spatially uniform Π , ∇π i j = 0, this reduces (11) to the linear form This equation shows that, even for spatially uniform fabric tensors, w 1 and w 2 are coupled by the off-diagonal terms of Y −1 BY . This is in contrast to the situation in rectangular coordinates, where uniform fabric tensors give rise to decoupled characteristics and to straight path solutions, on which the stress is constant whilst vanishing elsewhere. The inherent coupling in polar coordinates makes stress chains δ θ δθ r Fig. 1 The right hand side characteristic path (in blue), emanating from a source of width r 0 δθ at the internal boundary, flares out and broadens to r δθ (colour figure online) in annuli a much richer phenomenon, which we proceed to explore.
Consider the response to a localised stress source at the inner boundary, r 0 . Clearly, the linearity of the equations guarantees that the response to any stress source distribution can be found by superposition. Along the characteristic paths, we have ∂ ρ θ = 1/λ i . Therefore, a point source at (r, θ) = (r 0 , θ 0 ) gives rise to a characteristic path, whose trajectory satisfies It follows that a stress source of width r 0 δθ at the internal boundary generates two stress chains whose trajectories 'flare out' and broaden.This is illustrated in Fig. 1 for the right hand side characteristic path.
For completeness, we present an alternative analysis that allows π θθ = 0. In this case, π rr = 0 (π θθ and π rr cannot both vanish lest the equations are no longer hyperbolic). Following the same procedure as above, we now scale the fabric tensor,q i j ≡ π i j /π rr , substitute from relation (6) into (4)-(5) and rewrite these as whereÃ In terms of v, the characteristic variables, rr > 0, these eigenvalues are also real and distinct. In terms ofw, (14) becomes whereh ≡Ỹ −1 g. Unsurprisingly, the forms of Eq. (16) and (11) are very similar. Again, we see that the characteristics are coupled byỸ −1BỸ even for spatially uniform fabric tensors. The length parameters along the characteristic paths, s i , are now: ∂ ρ s i = 1 and ∂ θ s i = 1/λ i (i = 1, 2). Along the paths ∂ ρ θ i =λ i and the two paths flare out following the relation For spatially uniform fabric tensors, Eq. (16) reduces to

Example solutions
The coupling between the characteristics makes it difficult to derive general analytic solutions and we resort below to insight derived from numerical solutions. However, it is instructive to analyse first the special case when Π is diagonal. In this case,q θθ < 0,λ 2 = −λ 1 = −q θθ and a straightforward calculation yields from which the stress is found exactly An example of such a solution for σ rr is shown in Fig. 2, when a narrow Gaussian σ rr stress source is applied at (r, θ) = (r 0 , 0) (where s = 0 for both characteristics). This solution exhibits several interesting features. One is the aforementioned curving of the stress paths. Another is flaring out of paths with r . A third is a broadening of each path with r . A fourth is 'leaking' of stress to the region between the paths, known as the cone of influence [22,23]. This 'leak' is also the cause of attenuation of the stress along the paths. Again, in solutions in rectangular coordinates, these phenomena cannot occur for spatially uniform fabric tensors, which would exhibit only straight path trajectories carrying constant stress. A fifth, and a spectacular feature, is that the stress chains can curve backwards! This underlines the difference between isostaticity and strain-based theories, where such a phenomenon cannot occur. In the isostatic medium, the 'back-bending forces' are balanced by the stress that leaks to the cone of influence. We have checked that the back bending of the stress chains does not lead to sign change of the stress anywhere in the entire system. This means that tensile forces do not develop anywhere in the system, which would have destabilised the structure of dry granular materials. Therefore, such solutions are physically viable. Indeed, back-bending force chains have been observed experimentally [24], as can be seen in Fig. 2. We can also predict the conditions for this phenomenon to be observed. At the point where bending back first occurs, the tangent to the path trajectory makes an angle of π/2 with the original orientation, which we designate as the x-axis. Writing this condition as dx/dy = 0 and converting it to polar coordinates we get that this point corresponds to θ > π/4 on the left branch and θ < −π/4 on the right, or θ >| π/4 |. Using then (17), we find that the ith characteristic path starts bending back at a critical radius, r ci , satisfying In other words, to observe stress chains bending back, the stresses along the paths should not attenuate to invisibility before the critical radius is reached.
For more involved fabrics, whenq r θ = 0, the characteristic paths may no longer be symmetric. We plot an extreme example of such a solution, obtained numerically for σ rr and σ θθ , in Fig. 3. For this solution we usedq r θ = 2 and q θθ = −1. Note the bending back of one of the paths. The stress in this solution also remains compressive throughout the system.
The full richness of the solutions in cylindrical geometries emerges for fabric tensors that vary spatially across the medium. Analytic solutions for such media are difficult to obtain and we resort to numerical solutions for insight. In the following, we keep to narrow Gaussian stress sources at the inner boundary. The example illustrated in Fig. 4 is of a localised perturbation to the fabric at r = 3r 0 : q θθ = −1 − ex p −100 (r − 3r 0 ) 2 /r 2 0 . The perturbation gives rise to a clearly observed path branching. This is reminiscent of the branching observed in rectangular systems [22,23].
The effects of spatial non-uniformity are nicely isolated and illustrated in the symmetric case and, for clear visualisation of the branching effect, we solve for the fabric tensor (q rr ,q r θ ,q θθ ) = (1, 0, −0.2 cos(5r/r 0 ) − 0.3). The solutions for σ rr and σ θθ , plotted in Fig. 5, show clearly the periodic branch-like behaviour induced by the periodicity in the medium. The stress outside the cone of influence is identically zero whilst stress leaks from the main paths into the cone of influence at periodic intervals, where the gradient of the fabric tensor are largest. The leakage is along secondary characteristic paths.

Experimental support
Several experiments support these results. Flaring of force chains has been observed, e.g. in figure 7 of [25] and in our Fig. 2 [24]. Back-bending forces have also been observed, as Fig. 2 shows. Another intriguing and potentially related experimental observation has been reported in [26]. The experiments consisted of shearing sand-filled Couette cells by a very slowly rotating inner boundary and observing formation of patterns of Mandala-like slip lines in the medium. The slip lines were narrow, well defined and appeared in pairs, flaring out almost exactly symmetrically from points along the inner boundary (their Fig. 4). Moreover, the slip lines pass through one another with very little interaction, if The σ r θ -chain solutions (left) for the fabric tensor measured from Fig. 4 of [26]. For a clear comparison with the experimentally observed Mandala-like slip lines (centre, courtesy of Bobryakov and Revuzhenko [26]), only ten chain pairs are shown. The rotation in the experiment breaks the left-right symmetry and gives slight preference to one family of characteristics over the other (right). The chains geometries compare well to the experiment despite the rough estimate from their figure at all. The pattern and shapes of these slip lines are identical to the stress chain solutions derived here and it is tempting to relate the two. Our interpretation of their observations is that, on the verge of slipping, their medium is almost perfectly isostatic and the low shear rate makes it possible for the medium to remain very close to this state, making possible the appearance of our solutions. The slow shear generates distinct stress sources along the inner cylinder, which give rise to characteristic paths in the material, following the analysis presented here. The slip lines then form at the zones of highest shear stress. The slow rotation of the inner axis breaks somewhat the left-right symmetry, which is why one characteristic family is more evident than the other.
We can further use their Fig. 4 and data to obtain information about the structural constitutive properties of the granular medium in their experiments. Assuming that the structure can be approximated by a spatially uniform fabric tensor, the fact that there is very little interaction between the slip lines suggests that π r θ = 0. This also accounts for the near-symmetry of the slip line pairs. We can also estimate, from the figure and the reported data, an approximate relation between the flaring out angle and radius along every slip line. Using this estimate in our Eq. (17) we deduce that their effective fabric tensor satisfiesλ 1,2 ≈ ± π 6 ln 2.5 and hence thatq θθ = −λ 2 i = −0.33 ± 0.06, where the error stems from the uncertainty of measurements from their Fig. 4. In Fig. 6, we plot the shear stress chains, σ r θ , for this fabric tensor. For clarity, we show only ten pairs, as the experimentally produced twenty plus pairs would clutter the figure. We also plot the left hand side characteristic paths family on its own for clear comparison. The strong similarity between the experimental and theoretically derived chain geometries supports both our analysis, as well as our interpretation of their observations.
It should be noted, however, that this Mandala-like pattern can only arise when the fabric tensor is very close to uniform across the system and π r θ = 0. It is very likely that this is the case in the experiment because of the careful initial preparation of the material and the very slow shear rate. More general fabric tensors are expected to be more disordered, having have local gradients, which would give rise to secondary disordered stress chains.

Conclusion and discussion
To conclude, we analysed the isostaticity stress field equations in polar coordinates. We derived the equivalence relation between the constitutive fabric tensors in cartesian and polar coordinates and the constraint that the components of the latter must satisfy for the equations to be hyperbolic and yield stress chain solutions. The stress equations were analysed and an explicit formula has been obtained for the flaring out of the stress chains and their broadening away from the inner boundary. This relation was then used to show that stress chains can bend back and exert force components in a direction opposite to the original loading! This striking phenomenon, which is impossible in strain-based theories, such as elasticity, is a fingerprint of the arching effect. This phenomenon is not dissimilar to chains of dominos, which can be made to fall, and thus have momentum, in opposite direction to the initial domino.
We emphasise that this analysis is significant beyond the relevance to cylindrical geometries and Couette cells. It is often the case that an external load is localised almost at the particle level, say due to grains being pressed by the boundaries more than their neighbours. Since the size of the system is normally much larger than the localised length scale, the stress field around such a source has locally unavoidably a cylindrical symmetry. Consequently, the stress field near a source is better described by Eqs. (4)-(6) then by (1)-(3). This suggests that some of the effects we derived here may also be observed in experiments other than in Couette cells.
It is interesting that the flaring out phenomenon does not arise in most of the solutions in the literature, where rectan-gular coordinates were used and straight stress chains were found for uniform fabric tensors. We believe that the reason for the apparent discrepancy between those and the solutions studied here is that we constrain our fabric tensors to have a negative determinant for all azimuthal angles θ , a constraint that is missing from those initial studies. Moreover, in view of the above discussion, not imposing this constraint leads to regions in the plane where the fabric tensors, used in the literature, may have positive determinants, as shown by our Eq. (8). In reality, this means that stress chains may disperse when incident on a region with a positive determinant of the fabric tensor, where the equations become locally elliptic.
Numerical solutions were then derived for several uniform fabric tensors in annuli, supporting the analytical results and demonstrating that stress chains indeed flare out, broaden and bend back. These solutions also demonstrated the further rich behaviour in cylindrical systems: chains leak stress into the region between the chain pairs -the cone of influence. That these phenomena occur even for perfectly uniform media, when ∇π i j = 0, is in stark contrast with solutions in rectangular geometries, where such fabric tensors can only give rise to constant stress along straight stress paths and zero stress elsewhere.
We then studied several numerical solutions for spatially varying media and showed that large gradients in the local fabric tensor lead to stress chain branching, a phenomenon seen also in rectangular geometries [22,23]. This branching is in effect a strong 'leak' from a localised region on a specific path, whose trajectory is along a conjugate secondary characteristic path into the cone of influence. Effects of general spatially varying structures were also illustrated by solving for a fabric tensor with a component periodic in the azimuthal angle. The periodicity in the fabric induces periodic leaking stresses from the main paths, again via secondary paths into the cone of influence. A similar phenomenon would be observed for periodicity in the r -direction.
Finally, we pointed out experimental observations that not only support our results but also can be explained afresh in view of them. To the best of our knowledge, back-bending forces, although observed in experiments, have been neither discussed nor studied in the literature. It would be interesting to study this phenomenon in more detail in light of the predictions made here. Our results also suggest a new quantitative interpretation of the shearing experiment in [26].
This work can be extended in a number of directions. Theoretically, the next natural step is to obtain solutions for nonuniform fabric tensors and test how the statistics of their local gradients affect the statistics of the main stress chains, as well as those of the secondary and tertiary branching ones. Following the discussion above, concerning the differences between the solutions in rectangular and polar coordinates, another direction to explore is a more detailed understanding of the behaviour of stress chains for the fabric tensors, conjec-tured in the literature for rectangular coordinates, once used as input for our equations here. In particular, it would be interesting to observe how straight stress chains [11][12][13][14]22,23] disperse upon entering an elliptic region of the fabric tensor. Another natural extension is to three dimensional systems, but this extension has to wait until a proper such theory is constructed even for rectangular coordinates. These theoretical directions are being taken currently in our group. We are also looking forward to real and numerical experiments to test our theory, in particular of very slow shear in Couette cells. For example, related simulations [27] found an inclination of the principal stress direction relative to the radial direction, which initially increases with radius and then saturates to a constant of about π/4. Our theory predicts that this inclination, which is the direction of the stress chain, would increase continuously. Those experiments, while relevant, are inconclusive as a test. Firstly, the increase found there is over a 2-particles thick shear band, which cannot be considered a continuum. Secondly, their observed saturation to a constant angle is probably due to the increasing connectivity away from the shear band, which takes the medium away from the isostatic state, where the theory is valid. Thirdly, our theory predicts two such principal stresses emanating from the inner cylinder, while they observe only one. This could also be a result of the broken symmetry by the shear. Nevertheless, in view of our results, it would be useful to modify such experiments accordingly to provide a rigorous test of our solutions, e.g. by generating wider isostatic-like shear bands.

Compliance with ethical standards
Conflict of interest The authors declares that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.