## Abstract

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 *k*_{f} and rotational springs of stiffness *k*_{b} at the hinge positions as shown. This model has just *N* degrees of freedom, *q*_{i}*L*, one for each joint displacement, while the granular model supplements this with an additional *N* rotational 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 *k*_{b}. If the vertical shortening of the (*i*+1)th link is
the total potential energy can be written as
2.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.

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 *q*_{i}*R* 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 *k*^{r} and *k*^{t}, respectively. Similarly, an elastic spring of stiffness *k*^{s} provides lateral support to each disc.

From the general position of figure 2*a*, the centre-lines can be returned to alignment by first rotating the top disc anticlockwise through the *rolling angle*,so the centre-lines become parallel as seen in figure 2*b*, and then sliding the discs relative to one another for them to line up. With reference to figure 2*a*, *α* and *β* can be written asand the misalignment of centre-lines seen in figure 2*b*, as measured by the difference in arc-lengths around the two circumferences, is given byThe end shortening Δ_{i} between cells *i* and (*i*+1) isand so the total nonlinear potential function for the complete system is
2.2where *F* is the load. Summations from 0 to *N* enable different options for boundary conditions. In the following, zero displacement (*q*_{i}=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 *q*_{i} 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 *k*_{b}. 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 *k*_{b} 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,
3.1Following Thompson & Hunt (1973, p. 100), we next introduce the diagonalizing transformation,
3.2producing in effect a set of *N* orthogonal buckling modes with wavenumbers *m* and amplitudes *u*_{m}(*m*=1,…,*N*). Substitution into *V*_{lin} 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 *V*_{lin}-function is positive definite and stability is ensured. A critical load *P*^{C} 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 *u*_{i}, 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 *k*_{b} and *k*_{f}. Critical load plots for two typical parameter sets are given in figure 3. As the foundation stiffness *k*_{f} increases in comparison with *k*_{b}, 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.

### (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 §2*b* becomes
3.3involving in general a total of 2(*N*+2) unknowns, *q*_{i} and *ω*_{i} (0≤*i*≤*N*+1).

Following the earlier approach of equation (3.2), we next expand *q*_{i} and *ω*_{i} as Fourier modes,
3.4and substitute into the linearized form of *V* to give a new energy function *W* of the 2(*N*+2) modal amplitudes *u*_{m} and *ϕ*_{m}. *W* is in a near-diagonal quadratic form; terms in *u*_{i}*u*_{j}, *ϕ*_{i}*ϕ*_{j} and *u*_{i}*ϕ*_{j} are zero for *i*≠*j*, but those in , and *u*_{i}*ϕ*_{i} in general exist.

The above forms are useful for boundary conditions where *q*_{i}=0 and *ω*_{i} is symmetric (such that *ω*_{i−1}=*ω*_{i+1}) at the ends. The alternative of zero *ω*_{i} and symmetry in *q*_{i} is obtained by switching sines and cosines in these expressions (Tordesillas *et al.* in preparation). Also, the Fourier expansions for *q*_{0} and *q*_{N+1} ensure that *u*_{0} and *u*_{N+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 2*N* unknowns *u*_{i} and *ϕ*_{i}: 1≤*i*≤*N*. Next, the *ϕ*_{i} can be eliminated via the *N* first-order equilibrium equations ∂*W*/∂*ϕ*_{i}=0, giving (Thompson & Hunt 1973)
3.5where superscript 0 indicates evaluation in the undeflected state. As no *u*_{i}*ϕ*_{i} cross-terms appear in the work done by load, these equations are independent of *F*. The resulting expressions for *ϕ*_{i} in terms of *u*_{i} are then readily substituted into the *W*-function to produce an *N* degree of freedom diagonalized energy function, written as *A*(*u*_{i},*F*) to distinguish it from *W*(*u*_{i},*ϕ*_{i},*F*).

Double differentiation of *A* with respect to each of the *u*_{i} in turn leads to a set of stability coefficients which can be set to zero to obtain *N* critical loads *F*^{C}, corresponding to the harmonically varying critical modes or eigenfunctions defined by the integer values of 1≤*m*≤*N* 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 *k*^{s}. The plots have much in common with the link model response of figure 3. Again, increasing the support stiffness (*k*^{s} or *k*_{f}) raises the minimum critical load and lowers its wavelength.

## 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 *u*_{i}*u*_{j}*u*_{k}*u*_{l}, 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):
4.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.

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 *k*_{b} in comparison with the foundation *k*_{f} has a generally destabilizing effect. When *k*_{b}=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 *u*_{i}*u*_{j}, *ϕ*_{i}*ϕ*_{j} and *u*_{i}*ϕ*_{j} (*i*≠*j*) are absent, the equivalent form of equation (4.1) is (Thompson & Hunt 1973)Figure 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.

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 §4*b*(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 §4*a* 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.

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 (–) has localization occurring at one end. This is not particularly surprising as, in the spirit of simply supported boundaries, rotational springs *k*_{b} 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 in figure 7*c*, which localize internally and show much in common with a shear band.

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 8*a* shows load/end shortening responses emerging from the three lowest bifurcations of figure 3*b*. Solution with *m*=51 develops into a symmetric localized form as shown in figure 8*c* and ends up in an apparent tangle (not shown). The two with next lowest critical loads, with *m*=50 and with *m*=52, grow into the antisymmetric forms shown in figure 8*b*. 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 in figure 9.

The figure depicts joint displacement and link rotation plotted along the length, at four different positions on the snaking path 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 *k*^{t} 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 *k*_{b}, 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.

## Acknowledgements

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

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society