Royal Society Publishing

Force-chain buckling in granular media: a structural mechanics perspective

Giles W. Hunt, Antoinette Tordesillas, Steven C. Green, Jingyu Shi


Parallels are drawn between the response of a discrete strut on a linear elastic foundation and force-chain buckling in a constrained granular medium. Both systems buckle initially into periodic shapes, with wavelengths that depend on relative resistances to lateral displacement, and curvature in the buckled shape. Under increasing end shortening, the classical structural model evolves to a localized form extending over a finite number of contributing links. By analogy, it is conjectured that the granular model of force-chain buckling might follow much the same evolutionary route into a shear band.

1. Introduction

When viewed from the structural mechanics perspective, the response of a constrained granular medium under load has a number of familiar characteristics. In particular, loss of initial stability is often at a point of bifurcation. This means that such systems are open to study from the point of view of nonlinear bifurcation theory, for example, as applied to highly nonlinear plate and shell problems and stemming notably from the work of Koiter (1945). Such theories focus in the first instance on a perfect structural system, often carrying a high level of symmetry. Interest is then drawn to the disproportionate effects of small imperfections, which can have a dramatic influence on load-carrying capacity and led Bernard Budiansky in the 1960s to coin the phrase imperfection sensitivity. A perfect system for a granular medium would comprise identical spheres in an ordered array. Under uniaxial loading, resistance is then provided by parallel force chains. From the modelling perspective, it is useful in the first instance to focus attention on the response of a single chain, supported laterally by a system of confining forces that are due to other force-chains and/or particles in the complementary weak network (Tordesillas & Muthuswamy 2009). This is the concept behind the two-dimensional models which follow.

We shall consider a two-dimensional rectangular array of identical flat circular discs, and model the buckling of a single internal column of cells arranged orthogonally to the direction of the load (Tordesillas & Muthuswamy 2009). This is compared and contrasted with the model of figure 1, comprising N+1 rigid links of length L linked by N hinges and supported by linear springs of stiffness kf and rotational springs of stiffness kb at the hinge positions as shown. This model has just N degrees of freedom, qiL, one for each joint displacement, while the granular model supplements this with an additional N rotational degrees of freedom.

Figure 1.

Spring and link model. (a) Undeflected state with unstressed springs marking the perfect system. (b) Deflected state showing degrees of freedom.

In keeping with the cross-disciplinary spirit of this Theme Issue, each system is expressed in terms of a conservative potential energy function. Although friction has a clear role to play in granular media, this approach follows a tradition of applying such methods to non-conservative problems where load reversal is either absent, or is covered by appropriate switching between potentials (Hunt & Bolt 1986; Budd et al. 2003). In the granular context, here it is assumed that deflections and rotations are always increasing and accordingly in the small displacement range are given linear elastic descriptions. It is anticipated that nonlinear effects, such as those associated with friction, for example, can be brought in at a later stage.

2. Total potential energy functions

(a) Link model

The model of figure 1 is effectively that of Thompson & Hunt (1973, pp. 99–102), with the addition of bending stiffness in the springs kb. If the vertical shortening of the (i+1)th link is Embedded Imagethe total potential energy can be written as Embedded Image2.1

(b) Granular model

A two-dimensional discrete-element granular model of a force chain (Tordesillas & Muthuswamy 2009) would have three degrees of freedom for each cell, two translational and one rotational. Here this is reduced to two: one for sideways or lateral displacement as in figure 1 and one for cell rotation. Elimination of the third freedom is justified by the following reasoning.

The buckling responses of two-dimensional plate and shell structures are marked by strong interactions between out-of-plane (bending) and in-plane (stretching) contributions, but in one-dimensional struts such interactions are absent. Thus, a model like that of figure 1 but with (high) axial stiffness in the links would behave almost identically to its rigid counterpart, but would seemingly need an extra set of degrees of freedom for a full description. The same reduction in the requisite number of degrees of freedom can be obtained in a discrete granular model by assuming the discs to be rigid. Further, recent discrete-element simulations and photoelastic disc experiments show that the failure of contacts, i.e. irreversible relative motion at contacts, along the force-chain column is predominantly due to the large relative rotations developed between force-chain particles (Tordesillas 2007; Tordesillas et al. 2009). Tordesillas & Muthuswamy (2009) recently explored the complex interplay between the various resistances to force-chain buckling. Findings from this analysis and past experimental studies discussed in Tordesillas & Muthuswamy (2009) suggest that the two major resistances to force-chain buckling are: (i) resistance to relative rotations between particles at contacts along the chain and (ii) lateral support to the force-chain column.

Figure 2.

Discs i and i+1. (a) General deflected state. (b) Rolling without sliding so that the original centre-lines are again parallel.

With reference to figure 2, we therefore take a model comprising N+2 identical rigid discs of radius R, with lateral displacements of the centres given by qiR and clockwise (positive) rotations from vertical by ωi. In keeping with the concept of a perfect structural system, the discs are assumed initially to be perfectly aligned. Resistances to relative rotation as well as relative tangential motion or slip operate at contacts between the discs: these are governed by linear elastic spring forces with stiffnesses kr and kt, respectively. Similarly, an elastic spring of stiffness ks provides lateral support to each disc.

From the general position of figure 2a, the centre-lines can be returned to alignment by first rotating the top disc anticlockwise through the rolling angle,Embedded Imageso the centre-lines become parallel as seen in figure 2b, and then sliding the discs relative to one another for them to line up. With reference to figure 2a, α and β can be written asEmbedded Imageand the misalignment of centre-lines seen in figure 2b, as measured by the difference in arc-lengths around the two circumferences, is given byEmbedded ImageThe end shortening Δi between cells i and (i+1) isEmbedded Imageand so the total nonlinear potential function for the complete system is Embedded Image2.2where F is the load. Summations from 0 to N enable different options for boundary conditions. In the following, zero displacement (qi=0) and reflective symmetry (ωi−1=ωi+1) are taken at the ends (i=0 and N+1). An alternative might be to assume reflective symmetry in qi and zero ωi. The difference is explored in more detail in a companion paper (Tordesillas et al. in preparation).

3. Linear eigenvalue solutions

(a) Link model

The linear eigenvalue analysis of Thompson & Hunt (1973) is readily extended to include bending energy resulting from the springs kb. That analysis showed that with simply supported boundary conditions, each of the resulting N critical buckling modes can be inscribed within a sine wave. In the absence of bending stiffness, the lowest critical load is associated with the highest wavenumber such that the system buckles initially with alternating right/left displacements at the hinges. With the appearance of the springs kb this is no longer the case. The sinusoidal shapes remain, but the lowest critical load is now not necessarily linked to the highest wavenumber.

First-order expressions for the curvature and end shortening of the V -function (2.1) give the linearized form for the model of figure 1, Embedded Image3.1Following Thompson & Hunt (1973, p. 100), we next introduce the diagonalizing transformation, Embedded Image3.2producing in effect a set of N orthogonal buckling modes with wavenumbers m and amplitudes um(m=1,…,N). Substitution into Vlin leads to a diagonalized quadratic form with P-dependent coefficients. The stability of the undeflected equilibrium state now rests with the signs of these stability coefficients (Thompson & Hunt 1973): if all are positive, then the Vlin-function is positive definite and stability is ensured. A critical load PC for each of the modes in turn can be obtained by setting the appropriate coefficient, effectively the second derivative of V with respect to its corresponding amplitude ui, to zero. The lowest critical load and its related mode—decided by the criterion of minimum strain energy stored per unit end-shortening—now depends crucially on the relative values of kb and kf. Critical load plots for two typical parameter sets are given in figure 3. As the foundation stiffness kf increases in comparison with kb, the lowest critical load gets larger while its characteristic wavelength drops. This indicates a growing preference to store energy in bending rather than in the supporting medium.

Figure 3.

Critical loads for the link model with N=100 and kb=1: (a) Embedded Image at m=22; (b) Embedded Image at m=51.

(b) Two-dimensional granular model

Performing the necessary substitutions and taking first-order expressions for the arcsines and square roots of equation (2.2), the linearized potential function for the model of §2b becomes Embedded Image3.3involving in general a total of 2(N+2) unknowns, qi and ωi (0≤iN+1).

Following the earlier approach of equation (3.2), we next expand qi and ωi as Fourier modes, Embedded Image3.4and substitute into the linearized form of V to give a new energy function W of the 2(N+2) modal amplitudes um and ϕm. W is in a near-diagonal quadratic form; terms in uiuj, ϕiϕj and uiϕj are zero for ij, but those in Embedded Image, Embedded Image and uiϕi in general exist.

The above forms are useful for boundary conditions where qi=0 and ωi is symmetric (such that ωi−1=ωi+1) at the ends. The alternative of zero ωi and symmetry in qi is obtained by switching sines and cosines in these expressions (Tordesillas et al. in preparation). Also, the Fourier expansions for q0 and qN+1 ensure that u0 and uN+1 make no appearance analytically, and equilibrium equations ∂W/∂ϕ0=∂W/∂ϕN+1=0 then show also that ϕ0=ϕN+1=0. This reduces the problem to the 2N unknowns ui and ϕi: 1≤iN. Next, the ϕi can be eliminated via the N first-order equilibrium equations ∂W/∂ϕi=0, giving (Thompson & Hunt 1973) Embedded Image3.5where superscript 0 indicates evaluation in the undeflected state. As no uiϕi cross-terms appear in the work done by load, these equations are independent of F. The resulting expressions for ϕi in terms of ui are then readily substituted into the W-function to produce an N degree of freedom diagonalized energy function, written as A(ui,F) to distinguish it from W(ui,ϕi,F).

Double differentiation of A with respect to each of the ui in turn leads to a set of stability coefficients which can be set to zero to obtain N critical loads FC, corresponding to the harmonically varying critical modes or eigenfunctions defined by the integer values of 1≤mN in equation (3.4). Here m=1 denotes a single half wave over the full length, and m=N is the alternating left/right displacement mode of the Thompson & Hunt (1973) model. The lowest critical load again depends on the relative values of the supporting, sliding and rolling stiffnesses. Figure 4 shows the spectrum of critical loads for a realistic set of parameters (Tordesillas & Muthuswamy 2009) with two different values of ks. The plots have much in common with the link model response of figure 3. Again, increasing the support stiffness (ks or kf) raises the minimum critical load and lowers its wavelength.

Figure 4.

Critical loads for a 100 cell granular model with R=0.00114 m, kt=0.3kn and kr=0.5kn, where kn=105 N m−1. (a) ks=0.016kn, Embedded Image N at m=20. (b) ks=0.28kn, Embedded Image N at m=50.

4. Post-buckling behaviour

(a) Symmetry tests

In the post-buckling range, behaviour will be governed by the full nonlinear potential functions (2.1) and (2.2), rather than their linear counterparts (3.1) and (3.3). Each can be seen as a Taylor series expansion about the undeflected equilibrium state, with the leading terms being the quadratics of the previous section (Thompson & Hunt 1973). The following symmetry tests apply equally to link and granular models and highlight the underlying similarities between the two systems.

All cubic terms of energy can be shown to be absent by left/right reflective symmetry. For each model, if the amplitudes of all modes are reversed simultaneously, a perfect mirror-image reflection about the vertical centre-line is obtained. This must have the same potential energy as the original and implies that no cubic terms can exist; if one does, it must bring with it a change in V (Hunt 1986).

Certain quartic terms can also be seen to be absent by reversing the amplitudes of all even-numbered waves in the presence of their odd-numbered counterparts. This again leads to an identical shape, reflected in this case about the horizontal axis, and hence identical potential energy. By the same reasoning, it implies the absence of all terms in uiujukul, where i+j+k+l is odd. In turn, this also suggests that secondary bifurcations from initially symmetric or antisymmetric post-buckled shapes into more general forms are likely, as discussed later.

(b) Stability of initial bifurcation

(i) Link model

Stability of initial post-buckling rests in the first instance on the existence or otherwise of cubic terms of energy (Koiter 1945; Thompson & Hunt 1973). Here all cubics are absent, so the bifurcation has an initial slope, on a plot of load against amplitude of buckling, of zero. It is next classified as either stable symmetric or unstable symmetric, depending on whether the load increases or decreases with buckling amplitude. The difference is found in the quartic terms of energy.

The diagonalization means the initial post-buckling curvature for mode m of the link model can be written in terms of derivatives evaluated at its critical point C as follows (Thompson & Hunt 1973): Embedded Image4.1Figure 5 plots this result for the two situations of figure 3. In a normal loading sequence, the denominator would be negative, so stability rests with the numerator, representing the fourth derivative of V with respect to the amplitude of buckling, evaluated at the critical point in question. Its sign then determines directly whether the post-buckling is stable-symmetric or unstable-symmetric.

Figure 5.

Initial post-buckling curvatures for the link model with N=100 and kb=1: (a) kf=0.2, Embedded Image at m=22; (b) kf=4, Embedded Image at m=51.

Because unstable post-buckling is a natural precursor of localization (Hunt 2006), it is interesting to observe that in each case of figure 3 the mode with the minimum critical load is initially unstable. It has become clear from these and other examples that reducing the bending stiffness kb in comparison with the foundation kf has a generally destabilizing effect. When kb=0, the outcome at zero load is thorough localization to a pair of links, as seen in fig. 12 of Hunt et al. (1997).

(ii) Granular model

Although the granular model has the added complication of cell rotations, the initial post-buckling can be found with much the same approach. For the particular form of energy function W in which all cross-terms uiuj, ϕiϕj and uiϕj (ij) are absent, the equivalent form of equation (4.1) is (Thompson & Hunt 1973)Embedded ImageFigure 6 plots this result for the two situations of figure 4. The outcome is remarkably similar to that of figure 5 for the link model. In each case, the mode with the minimum critical load is initially unstable and suggests that the system has the potential to localize.

Figure 6.

Initial post-buckling curvatures for the 100 cell granular model with R=0.00114 m, kt=0.3kn and kr=0.5kn, where kn=105 N m−1: (a) ks=0.016kn, Embedded Image at m=20; (b) ks=0.28kn, Embedded Image at m=50.

To summarize, figures 5 and 6 demonstrate the common characteristics of the two systems that post-buckling modes with large m are stable and those with small m are unstable. For all cases considered, the minimum critical buckling modes are unstable.

(c) Routes to localization

The initial post-buckling studies of §4b(i),(ii) describe the development of buckling into fixed-wavelength, periodic, modes. As an initial buckled shape develops into the post-buckling range, it is likely to be contaminated by other modes, leading eventually to localized shapes. The symmetry tests of §4a suggest certain cubic energy terms are absent, and this has the effect of splitting the system up into symmetric and antisymmetric forms. Initial buckling would be into either one or the other, and as post-buckling progresses, we would expect progressive contamination of the initial buckling mode from others with the same symmetry property, leading ultimately to localization. Secondary bifurcations can also occur into the form with the opposing symmetry property, as seen in Hunt et al. (1997, fig. 12). However, ultimately the choice of post-buckling path would depend less on symmetry properties than on least-work/minimum-energy considerations.

To explore these issues, we have performed numerical runs using the path-following continuation code AUTO (Doedel et al. 1997). This software enables equilibrium paths for boundary-value problems to be traced under parametric variation, while identifying all possible bifurcations. Figure 7 shows the load P for one of our link model examples, plotted against its corresponding end shortening Δ. As the stored energy equates to the area under this graph, such plots allow comparisons between alternative equilibrium states. The expected ‘least energy’ solution would be the one lying closest to the load axis on this plot.

Figure 7.

(a) Numerically obtained load/end shortening plots for the system of figures 3a and 5a. N=100;kb=1;kf=0.2;m=22. (b,c) Developing localizations.

The figure shows only some of the available post-buckled solutions, specifically those developing from the lowest critical bifurcation load of figure 3. An initial periodic buckling form goes through a sequence of secondary bifurcations, each leading to a different localized state. The figure shows that the least energy solution (Embedded ImageEmbedded Image) has localization occurring at one end. This is not particularly surprising as, in the spirit of simply supported boundaries, rotational springs kb have been omitted at the pin positions, i=0 and i=N+1. Such solutions are clearly sensitive to boundary effects, and a realistic experiment in a granular material (Tordesillas 2007) might well have more constraint at the ends than in the middle. We add no further comment here, but to avoid unnecessary distractions will henceforth limit discussion to solutions such as Embedded Image in figure 7c, which localize internally and show much in common with a shear band.

Figure 8.

(a) Load/end shortening plots for the lowest three critical points of figure 3b. N=100; kb=1; kf=4: Embedded Imagem=51; Embedded Imagem=50; Embedded Imagem=52. (b) Localized solutions at cross positions. (c) Route to localization for solution Embedded Image.

Figure 9.

Joint displacement and link rotation at different positions along the snaking path Embedded Image of figure 8.

It is interesting to observe that the wavelength of the critical buckling mode at the point of initial instability (between 9 and 10 links in the case of figure 7) seems to have a strong bearing on the length of the localization. Experience with a related continuous problem (Budd et al. 2001) tells us that the corresponding wavelength may grow as the load drops and localization develops. Here there is perhaps some evidence of this effect, but it does not appear to be of major significance.

Figure 8a shows load/end shortening responses emerging from the three lowest bifurcations of figure 3b. Solution Embedded Image with m=51 develops into a symmetric localized form as shown in figure 8c and ends up in an apparent tangle (not shown). The two with next lowest critical loads, Embedded Image with m=50 and Embedded Image with m=52, grow into the antisymmetric forms shown in figure 8b. One closely resembles the localized shape of figure 7, while the other has three separate humps. Each develops a form of the snaking response (Burke & Knobloch 2007; Dawes 2008) seen in a number of general situations when cells buckle locally and restabilize in turn (Hunt et al. 2000). This sequence is illustrated in more detail for solution Embedded Image in figure 9.

The figure depicts joint displacement and link rotation plotted along the length, at four different positions on the snaking path Embedded Image of figure 8. As end shortening increases, a localized buckle first forms and then grows in width. Meanwhile, the load oscillates between two values, one positive and the other negative. The localization develops into two loops, which pull apart leaving the in-between links in the undeflected position but inverted. The link rotations seen at the right mark this effect with relatively sharp transitions from zero to −π and back. With rotation levels set to something rather less than π, and oscillation about a positive rather than a zero load value, there would again be much in common with a shear band.

5. Concluding remarks

This paper has been written in something of a speculative vein, yet the outcome does look promising. The similarities between the two models seen in both the critical load and the initial post-bucking plots are striking. This is perhaps not so surprising in the case of the critical loads of figures 3 and 4, since after reduction to eliminate rotations ωi from the granular model the quadratic energy terms controlling the linearized responses are essentially the same. It is more of a surprise for the initial post-buckling of figures 5 and 6; these are governed by higher order energy terms which, even with kt set to infinity, fail to match in the two systems. The link model response develops a localized form that is characterized by rotation of a finite number of links over a clearly defined length. The speculation is of course that similar forms might be found in an extended granular model comprising a number of parallel force-chains, leading to something similar to a shear band on a length scale of a finite number of grains.

Given the linear elastic assumptions for the stiffnesses, it would indeed be surprising if realistic shear band shapes arose spontaneously from the model. But the images of figures 7 and 9 are evocative, and as a consequence the path-following techniques used to find these solutions are presently being modified to suit the granular model (Tordesillas et al. in preparation). The key characteristic that has emerged for the link model is that the wavelength of initial buckling is setting a length scale that is neither that of the grains nor that of the sample but somewhere in between, and this translates into the length of the localization. Thus, although in any realistic experiment small perturbations and imperfections, together with interactions with other parallel columns of grains, would decide the evolution in any particular instance, the initial elastic response could have more bearing on the general outcome than might be thought at first.

At the higher of the two values of stiffness kb, the link model also demonstrates what has come to be known as snaking (Burke & Knobloch 2007; Dawes 2008). This is a common phenomenon seen in a number of different circumstances, from the development of Luders lines in samples of plastic materials in tension to cellular buckling in axially compressed cylinders. In this case, destabilization followed by restabilization at a local level is allowing a deflection pattern to emerge that has a finite-length, large-rotation, internal region with sharply defined transitions, much like a shear, or indeed a kink (Wadee et al. 2004), band. And the sequence of destabilization followed by restabilization does resonate strongly with experimental responses for granular media (Tordesillas 2007). For a number of reasons, however, caution must naturally be exercised. First, the link rotations are not directly analogous to those of the granular model, although they may in some way be related. Second, the constitutive laws would be anything but linear. Third, these effects are bound to be influenced by transfer of load to other chains. Nevertheless, we believe that there is sufficient encouragement here to continue down this route, and that the way forward is clear.


A.T. acknowledges the support of the Australian Research Council (DP0772409) and the US Army Research Office (W911NF-07-1-0370).


View Abstract