## Abstract

Similarities and differences between the phenomena of kink banding in compressed layered structures and shear banding in compressed granular media are explored. Simple models are introduced for both, and the focus is directed onto how they can nucleate from the perfectly flat state. A convincing scenario is found for each in which a mode develops from an initial bifurcation into a periodic state, followed by rapid localization under falling load, while retaining decaying but wavy tails. At a certain lower critical load, the tails lose their waviness, and the expected form of the kink or shear band appears. In each case, good numerical evidence is provided for the existence of this form of behaviour. A second potential instability for the layered case is also explored, linked to the appearance of a critical force dipole that overcomes bending stiffness locally at some point along the length. This mode, which should appear with non-wavy decaying tails at the lower of the two critical loads mentioned earlier, proves somewhat elusive. Evidence is found for its existence in the linearized approximation to the layered model, but the search for numerical solutions to the underlying nonlinear equation is hindered by a shortage of suitable boundary conditions.

## 1. Introduction

The distinctive phenomena of kink and shear banding appear in many different structural and material systems, on scales ranging from micro to macro and beyond [1]. Figure 1 shows an experiment on half-sheets of A4 paper that have been constrained vertically, and then squashed horizontally under conditions of controlled end-shortening [2]. A black-edged sheet is introduced every 25 layers or so to aid visualization. Starting at the loaded end, kink bands appear at points B, D, F and G, in turn, shown in figure 1*b* as the output of axial load against end-displacement, each accompanied by a drop in load implying an underlying instability. Similar phenomena, known as shear bands, also appear in granular media in tri-axial tests, and softening elasto-plastic materials [3].

Kink and shear bands have been modelled from a number of perspectives, including finite-element analysis for the former [4] and discrete elements for the latter [5,6]. Comparisons with experiments have shown varying degrees of success [7,6]. The two mechanisms are clearly closely related, but there is little in the literature drawing them together. This is perhaps because of the diverse range of systems in which they appear, and possibly also the many different approaches taken to their modelling. Everyone has his or her own favourite view.

The view expressed here emphasizes the strong role that nonlinear geometry can play in such systems. We do this by limiting ourselves as far as possible to linear elastic behaviour. In so doing, we are not trying to deny the importance of material nonlinearity—merely drawing the reader's attention to what geometry can do on its own. We start by reviewing classical multi-layered folding in structural geology, as originally envisaged by Ramsay [8] and later described by Price & Cosgrove [9]. This study highlights the crucial interplay between geometry and mechanics in the selection of allowable configurations for the folding of two-dimensional multilayered structures.

If thickness is to remain constant as layers fit together without voiding, then there must be an ever-tightening curvature on the inside of a fold, as seen in figure 2*a*, a process that eventually results in a singularity [10]. If, on the other hand, layers remain the same shape but thickness variations are allowed, then the situation can be like that of figure 2*b*; we could perhaps think of (*a*, parallel folding) as being solid-like in terms of behaviour and (*b*, similar folding) as fluid-like, although, in the geological context, the actual rheology would more likely be a complex nonlinear mix [11]. These two circumstances have considerable bearing on our modelling processes, as seen below.

Figure 3 shows the nonlinear geometry associated with rotation of constrained layers, using the kink bands of figure 1 as an example. If the orientation *β* is constant and the layers remain of constant thickness outside the band, they are obliged to accommodate a change in thickness within it as they rotate through angle *α*, as shown. For *α*<2*β*, there is dilation, reaching a maximum at *α*=*β*. At *α*=2*β*, the layers have the same thickness inside and outside the band, and for *α*>2*β*, there is extra contraction in addition to the initial squash provided by the transverse constraint. It is thus proposed in Wadee *et al*. [2] that the mechanics of kink banding for experiments like that of figure 1 is a two-stage process, as follows.

The initial stage is one of rapid propagation across the sample of rotated layers inside the band. This is assumed to take place at the orientation *β* that, at the point of maximum dilation of figure 3*b*, exactly eliminates the initial transverse squashing: as propagation occurs, all work against friction is thereby avoided. Once the band has formed, further rotation within it can take place only through the configuration shown as figure 3*c* and beyond. Frictional effects are now brought into play, and, under further imposed end-shortening, the system would eventually lock up (figure 3*d*), with deformation most probably progressing to the next band. It should be noted that the apparently conflicting demands of constant thickness (parallel folding) and identical layers (similar folding) are mutually reconcilable when the deformation comprises straight sections and sharp corners as in a kink band, but only when *α*=2*β* as in figure 3*c* [12].

The paper starts with a brief outline of an earlier model of kink banding [2], based on flexurally rigid blocks of material sliding under Coulomb friction. Particular attention is paid to through-thickness effects, and good comparisons are found with experiments such as figure 1. This is then complemented by a purely elastic model of similar folding comprising layers that are offset one to another but otherwise identical. A nonlinear differential equation is developed based on the standard equation for an elastic beam on a linear elastic foundation [13]. A collocation method is introduced to solve the full nonlinear problem, giving a minimum energy solution that emerges from the flat fundamental state in the form of a kink band, orientated at angle *β* across the sample, at a finite compressive in-plane critical load. Finally, we take a similar elastic model for a constrained, two-dimensional, regularized granular medium [14]. Solutions are traced, using the bifurcation code Auto [15]. The linearized problem can be seen as a four-dimensional map, with eigenvalues that describe a discrete bifurcation problem. Shear band solutions, closely related to those of kink bands, are found, represented again by a stable energy minimum. As extra cells join in, and the shear band grows, the load-deflection response follows a well-known ‘snakes-and-ladders’ scenario centred around the classical Maxwell load [16,17].

## 2. Block model

A simple experiment on constrained and compressed layers of paper [2] was shown earlier in figure 1. About 800 half-sheets of A4 cut longitudinally were initially squashed together between two steel plates by a transverse load. Subsequently, an in-plane load was introduced by slowly imposing an axial displacement, whereas the transverse displacement was kept constant. Load cells recorded both loads at 1 s intervals, and an in-plane transducer registered the axial end-displacement. To aid visualization, a black-edged sheet was introduced about every 25 layers.

Output from a typical experiment is seen in figure 1*b*. Kink bands form from the loaded end, starting with B and followed by D, F and G. Each band appears with a clear-cut instability, marked by a drop in applied axial load, and is followed by subsequent re-stabilization and eventual lock-up. There are many approaches that could be taken to modelling such phenomena and here we introduce two, a block model and a continuous beam model.

Block models of kink banding, based on assumptions of straight sections punctuated by sharp corners, have been compared against experiments such those of figure 1 with good agreement [2]. A typical such model is shown in figure 4. Kink bands form in blocks that are rigid against flexure, with bending energy localized into hinges where it is stored in rotational elastic springs of stiffness *c*, the value of which can be estimated from the constraining overburden pressure *q*, and known properties of paper [12]. The blocks slide over each other resisted by Coulomb friction, with *N* representing the normal reaction between layers in the band. Work is also done by the external load *P*. The blocks have finite stiffness in both the axial and transverse senses. Springs *k* represent initial in-plane stiffness and give the system the ability to store energy, later to be released in a snapback process as seen in the experiment of figure 1. Transverse stiffness orientates the band at an angle *β* across the sample, so that *N* can drop to zero. A single elastic spring of stiffness *k*_{f} per layer represents the penalty for entering the stiff supporting foundation. Estimates for all parameter values, via a range of modelling techniques, are given by Wadee *et al.* [2], and enable direct comparisons with experiments such as those of figure 1.

Looking at the kink bands of figure 1 (also figure 3*b*), it can be seen that, within a band, the same number of layers is stretched at oblique angle *β* across the sample as is directly (transversely) compressed outside it. This means that, within the band, as *α* grows from zero the layers initially dilate, with maximum dilation occurring when they are positioned perpendicular to the orientation of the band across the sample (*α*=*β* in figures 3 and 4). If this is taken as the unstressed state, it is a simple matter [2] to show that the transverse (through-thickness) compressive strain is
2.1

The dilation reduces the transverse compression *N* within the band, and consequently the frictional forces that resist the shearing action. This gives the system the opportunity to reduce drastically the effects of friction, by propagating a dynamic jump at exactly the *β* value that drops *ϵ*_{t} to zero according to equation (2.1). This would take place in the state of maximum dilation, *α*=*β*, and would mean of course that *N* would also drop to zero.

Wadee *et al.* [2] propose that, under controlled end-shortening, it is this mechanism that selects the orientation angle *β*. Analysis of the development of a band then has two distinct stages:

Orientation

*β*is first determined from the imposed transverse squash, such that the state of maximum dilation corresponds precisely to*N*=0 but with the layers just touching. Thus, the greater the initial applied transverse strain, the more obliquely the band is orientated across the sample.Further evolution in

*α*then takes place at this fixed angle of orientation. This is based on the premise that it would be energetically very expensive for*β*to vary and hinge lines to move. As*α*grows such that*α*>*β*, transverse compression and friction again come into play, and the band eventually locks up.

We note that both *α* and bandwidth *b* are considered free variables, so the system can at all stages select its favoured bandwidth [18]. This would seem reasonable for the unstable part of the response, but fails to penalize hinge migration once the system has re-stabilized and lock-up is approached; however, the assumption may not be as crucial as it might at first appear, because bandwidth growth is also restricted by the increase in deflection into the foundation stiffness *k*_{f}.

Figure 5 shows a direct comparison between the predictions of this single-band block model and the experiment of figure 1, as described by Wadee *et al.* [2]. A fundamental equilibrium solution comprising the flat state is nominally stable at all loads, but in practice would be increasingly vulnerable to small disturbances as load increases into a state of extreme metastability. Under increasing load, the flat equilibrium path is approached asymptotically by a post-buckled path for which *α*≠0. As load drops on this path, *α* grows until eventually a minimum is reached in load/end-shortening space, whereupon the load starts to increase with *α* continuing to grow. The experiments show a sudden drop in load to this stable post-buckling state at a constant value of end-shortening; note that the 1 s interval between recorded readings gives a slightly misleading impression here, in that the instability appears to take place over a finite time step, whereas it is, to all intents and purposes, instantaneous.

It can be seen that, as end-shortening drops to a minimum and reverses direction with increasing *α*, *b* also continues to grow. Experimentally, however, this process would be interrupted by the development of further kink bands, not allowed for in the present analysis.

## 3. Continuous model

The block model performs well when compared with experiment, and reflects much of the geometry of kink banding. But it is based entirely on observation: there is no rationale to suggest why it should have the characteristics it does—straight but compressible limbs and sharp corners for instance—apart from the fact that it is close to what is seen experimentally. It can, however, be usefully complemented by a continuous model derived from the underlying well-known differential equation of beam bending. As with the blocks, it is first assumed that all (internal) layers take the same shape (similar, figure 2*b*, as opposed to parallel, figure 2*a*) folding. This means that the underlying model can be reduced to that of a typical layer in a stack of identical counterparts, all supported on an elastic foundation. The layer is given bending energy, a share of the foundation energy, and strain energy of squashing according to equation (2.1), along with work done by the end load. Calculus of variations then leads to a fourth-order nonlinear ordinary differential equation (ODE), which when solved numerically gives kink band solutions.

### (a) Total potential energy

#### (i) Bending energy

The bending energy of a typical layer of bending stiffness *EI* is given by
3.1

#### (ii) Work done into foundation

For a linear elastic Winkler foundation of equivalent stiffness *k* per layer, the work done into the foundation by a single layer is taken to be
3.2

#### (iii) Squash energy

For the geometry of figure 3 and a linear elastic filling material of transverse stiffness *k*_{V}, from equation (2.1), we can write the energy stored in the squashing of the beam as
3.3For displacement from the flat fundamental state, this is most usefully written in increasing powers of *α*. Remembering that *β* is a constant, for simplicity we assume that *β*^{2} can be ignored in comparison with unity, while noting that corrections to accommodate finite *β* could readily be applied to the coefficients determined below if required. Writing *α*=*u*_{x}, ignoring the constant term in *β*^{4}, and truncating after quartic terms in *u*_{x}, we have
3.4

#### (iv) Work done by load

If load *P* remains constant as it moves through displacement *Δ*, we must consider also its first-order work contribution [19],
3.5

### (b) Reversibility

The total potential energy, comprising total strain energy minus work done by load, can be written
3.6where is the Lagrangian [19]. The theory of dynamical systems can be applied to variational problems on an interval; within this representation, by replacing the normal role of time with the spatial coordinate *x*, the close linking of a Hamiltonian structure with reversibility can be exploited to the full [20,21]. Here, we simply note that *V* is invariant under the transformation
3.7suggesting the possibility of an antisymmetric homoclinic mode as a stationary solution of *V* .

#### (i) Euler–Lagrange equation

Calculating the derivative *δV* of *V* with respect to a vertical perturbation *φ*, it follows that
3.8If we consider homoclinic boundary conditions such that as , the boundary terms in (3.8) disappear. If *δV* =0 is to hold for all choices of *φ*, minimizers of *V* must satisfy the Euler–Lagrange equation
3.9leading to the fourth-order nonlinear differential equation
3.10Rescaling the variables such that , this Euler–Lagrange equation can be written in the dimensionless form
3.11where , *C*_{1}=(*k*/*EI*)^{3/4}, *C*_{2}=*k*/*EI* and primes denote differentiation with respect to the scaled parameter . From this point on, we work exclusively with this dimensionless form, dropping the tilde.

### (c) Linearized system

For small values of *u*, i.e. close to the flat state *u*=0, the nonlinear terms can be neglected and the Euler–Lagrange equation (3.11) linearized to become
3.12Substituting *u*=e^{ξx}, the characteristic equation is
3.13and the eigenvalues are
3.14The linearized behaviour can then be split into three distinct cases:

If

*γ*<−2,*ξ*is purely real. This leads to homoclinic behaviour with non-wavy tails, as might be expected for a kink band.If −2<

*γ*<2,*ξ*is complex conjugate. This leads to homoclinic behaviour with wavy tails [13].If

*γ*>2,*ξ*is purely imaginary. This leads to periodic behaviour [22].

Movement of these eigenvalues under changing *γ* is summarized in figure 6.

#### (i) Boundary conditions

The general solution to (3.12) has the form
3.15If we set *A*_{1}=*A*_{3}=0, then the response is purely that of the convergent inset. Differentiation with respect to *x* then gives expressions for *u*^{′}, *u*^{′′} and *u*^{′′′} in terms of *A*_{2} and *A*_{4}, which can be eliminated for the following *inset condition*:
3.16The corresponding outset condition is then obtained by reversing the direction of *x*, replacing the earlier-mentioned conditions by their negative counterparts. These conditions apply for small *u* only, but can be used as appropriate boundary conditions for a kink band at large but finite value, *x*=±*L*. As *u* grows, the inset and outset are influenced by the nonlinearity, linking in the form of a *homoclinic connection* [23]; fully localized solutions are represented by paths following the route, outset to inset. We thus refer to (3.16) as *homoclinic boundary conditions*.

The origin *x*=0 can also usefully be treated as a boundary, and clearly *u*_{0}=0 supplies one boundary condition. A second obvious choice is , as might be expected from the antisymmetric form of a kink band. However, this is of no use, if the second derivative is allowed to take a sudden jump at the origin, as in the following section: we return to this point in §3*d*.

#### (ii) Wrench instability

We have particular interest in the case where *γ*≤−2 and the roots are real. For *γ*<−2, the general solution involves two unequal values of *ξ*, meaning that any deflected shape satisfying the linearized equation (3.12) must grow exponentially in either the positive or the negative sense of *x*. There is a possible exception to this however, when *γ*=−2, and the real roots coalesce with *ξ*=1. This suggests that this point could act as a bifurcation, from which nonlinear solutions could emerge into the regime *γ*<−2, where the linear equation (3.12) need only be satisfied in the tails as . A possible mode shape can be written as
3.17and is shown in figure 7, along with its continuously varying first derivative *u*′, and discontinuous second derivative *u*^{′′}.

The jump in *u*^{′′} at *x*=0 implies a corresponding jump in bending moment at this point, as equal and opposite forces, infinitesimally misaligned, act as a dipole to create a critical moment or wrench and overcome the natural resistance of the material bending stiffness. The instability would appear as a bifurcation at critical load *γ*=−2, and run across the sample at angle *β*; note that the linear equation has a second solution at −*β*, but, to follow this, the corresponding nonlinear equation would require a sign change on terms involving *β*. We refer to the phenomenon as a *wrench instability*, and see it as a likely candidate for nucleation of a kink band. The associated bifurcation point at *γ*=−2 appears at the same parameter values as a Belyakov–Devaney bifurcation [20], but differs from the latter significantly in that it sits on the fundamental equilibrium path and so has a solution found in the roots of the linear equation.

Although the mode shape appears instantaneously centred about a point in *x*, it is significant that it retains continuity in both *u* and *u*′; discontinuities in *u*^{′′} are of course allowable in normal beam theory, at points of applied external moment for example.

The dipole can be found from Fourier transforms. Consider the fourth-order ODE
3.18Because the equation is defined on the whole real line, we can take its Fourier transform so that
3.19We note that the Fourier transform of the proposed solution *u*(*x*)=*Ax* e^{−|x|} is given by the standard calculation
3.20This implies that and therefore *f*=*Aδ*^{′}, where *δ*(*x*) is the Dirac delta function, so that *δ*^{′} represents a dipole or a point-wise moment at the origin. Thus, *u*(*x*)=*Ax* e^{−|x|} solves the following fourth-order ODE:
3.21and shows that if *A* is infinitesimally small, so too are *u*, the derivatives of *u* and the dipole.

### (d) Numerical solution

We now look for numerical solutions to the nonlinear fourth-order boundary problem (3.11) using pseudo-arclength continuation [24]. We seek solutions with the reversibility of (3.7), over a large but finite domain *x*∈[0,*L*]. To provide numerical stability the Euler–Lagrange equation is solved subject to the end-shortening constraint , so *γ* becomes a free parameter that is found as part of the solution. Physically, this can be related to rigid loading [19], and allows continuation in the space (*Δ*,*γ*). Representing the Euler–Lagrange equation (3.11) and the constraint by a system of five first-order equations, the boundary value problem is open to solution by standard collocation methods (e.g. Matlab's solver `bvp4c`). Antisymmetric boundary conditions, , are used at *x*=0, and homoclinic conditions (3.16) at *x*=*L*.

Figure 8 shows corresponding numerical solutions, emerging from a transcritical bifurcation point at *γ*=2. By using complex variables, Matlab automatically deals with the switch in homoclinic boundary conditions from real to complex conjugate eigenvalues at *γ*=−2. For good results, as the fundamental path is approached at *γ*=2, longer and longer lengths of strut are needed computationally. Although the bifurcation of figure 8 is clearly transcritical, we have not attempted to follow the solution as it enters the stable part of the post-buckling regime where *γ*>2. Periodic solutions would then be expected, and entirely different questions would arise over the appropriate selection of sample length *L*, together with other related issues [22].

The kink band behaviour of figure 8 initiates from a bifurcation point at *γ*=2, and starts life as a periodic solution. As *γ* drops, it localizes and develops wavy tails, before going through *γ*=−2 and losing the waviness. The question arises: can there be an equilibrium solution emerging directly from *γ*=−2 in the form of the wrench instability? At present, we are unable to provide any compelling numerical evidence to answer this. Difficulties lie in the choice of a fourth boundary condition to replace the antisymmetry condition ; for the linear solution of the wrench, this ought to be replaced by a finite, amplitude-dependent value, representing one half of the jump in at the origin. Our current numerical scheme therefore does not hold for the linear equation at *γ*=−2, and the search continues for a suitable fourth boundary condition.

## 4. Granular model

Recently, attention has been drawn to similarities and differences between models of discrete struts on elastic foundations, and regularized chains of constrained granular media [25]. This has led to the idea that layered and granular media may show similar behaviour under the same form of loading. More recently still, the single force chain model has been extended to embrace a number of parallel chains, constrained to act together much as the layers of the previous section [14]. For certain chosen parameter values, equilibrium solutions found using Auto [15] are seen to have the appearance of shear bands, similar in form to the kink bands of figure 1. In the following, we explore similarities between layered and granular systems, comparing the movements of the eigenvalues of §3*c*(i) with those for the equivalent linear map of the latter in §4*c*.

### (a) Outline formulation

A basic single chain model is given in figure 9. Following an initial formulation by Tordesillas & Muthuswamy [26], a set of *N* rigid discs is supported laterally by springs of stiffness *k*^{s} and loaded by a compressive force *F*. Resistance is provided against the discs rolling and sliding relative to one another by elastic springs of stiffnesses *k*^{r} and *k*^{t}, respectively [25]. Extra end discs are allowed to rotate but constrained from lateral movement. Sideways displacements are measured by a set of *N* generalized displacements *q*_{i}*R*, and there is a corresponding set of rotations *ω*_{i} as shown in figure 10. The *simple support* boundary conditions shown here allow the system to be fully represented by the 2*N* degrees of freedom *q*_{i} and *ω*_{i} [14]. In what follows, to model more realistically what might be expected in a shear band, these boundary conditions are replaced by the discrete map equivalents of the *homoclinic boundary conditions* of §3*a*(i). Figure 10 shows the general relationship between two adjacent cells, and demonstrates geometrically how rotation *ω*_{i}, rolling angle *ψ*_{i}=*ω*_{i+1}−*ω*_{i}, slope , slip displacement *d*_{i}=*R*(*α*−*β*)=*R*(2*θ*_{i}−*ω*_{i+1}−*ω*_{i}), and contribution to overall end-shortening relate one to another.

With the addition of an extra row of discs, the precise position of each disc in the second row can be fully determined with respect to a counterpart in the first [14]. Consider, for example, a pair of parallel force chains that begin life both flat (*q*_{i}=0) and close-packed. As well as contacting with two immediate neighbours in its own force chain, each disc is also then in contact with two in the parallel chain, in a circumstance sometimes referred to as a ‘force cycle’ [27]. As the system displaces from the flat state, kinematics dictates that the discs move away from being close-packed. If the slope remains within the bounds −60^{°}≤*θ*_{i}≤60^{°}, the shape of the second chain simply mirrors that of the first but displaced by 2*R*, with each disc in the second row sitting at either 11 o'clock or 1 o'clock with respect to a counterpart in the first [14]. The change in area of the void space from the close-packed configuration then provides a slope-dependent energy penalty, much like *U*_{V} of equation (3.3) for the continuous model.

### (b) Total potential energy

The total nonlinear potential function for the single chain model becomes
4.1with the measures *q*_{i},*d*_{i},*ψ*_{i} and *Δ*_{i} being described according to figure 10. The first term represents the energy stored in the foundation, the second is the slip energy, the third the rolling energy and the final term is the work done by the load *F*. With an extra row of discs, this contribution to the total energy per row remains unchanged, but there is an added term representing an energy penalty against movement away from being close-packed. For our present purposes, where *θ*_{i} is small, this can be taken as
4.2where *k*^{v} is a bulk stiffness measuring the resistance against an increase in void area. For more details, see Hunt & Hammond [14].

The total potential energy per row is then
4.3which can be linearized by dropping all terms above quadratic order in *q*_{i} and *ω*_{i} to give
4.4Differentiation with respect to the two fundamental degrees of freedom *q*_{i} and *ω*_{i} then leads to the two second-order linearized equilibrium equations,
4.5and
4.6valid if both *q*_{i} and *ω*_{i} are small. Note that as the (*i*−1)th, *i*th and (*i*+1)th terms all appear, this is effectively a pair of second-order difference equations.

### (c) Eigenvalues of the linear map

As with the continuous problem of the previous section, true homoclinic solutions can only be found with application of correct homoclinic boundary conditions at the ends of our finite length system to replace the simple supports of figure 9. These can be found from the eigenvalues and eigenvectors of the linear map, which is developable from equations (4.5) and (4.6). This involves substituting *q*_{i+1}=2*θ*_{i}+*q*_{i} and *ω*_{i+1}=*ψ*_{i}+*ω*_{i} to remove the (*i*+1)th terms, and then reworking to obtain *q*_{i} and *ω*_{i} in terms of the four discrete variables *q*_{i−1},*θ*_{i−1},*ω*_{i−1} and *ψ*_{i−1}. The map can then be written in matrix form as follows
4.7where A is a 4×4 non-symmetric matrix. Unfortunately, the elements of A are too complicated to be given explicitly here, so, from this point, we resort to a numerical description, taking typical parameter values from earlier work, *R*=0.00114 m, *k*^{s}=0.000113*k*^{n},*k*^{t}=0.2*k*^{n} and *k*^{r}=0.5*k*^{n}, where *k*^{n}=10^{5} N m^{−1} [25,26], along with *k*^{v}=16.1 GPa m^{−1} [14]. This gives, correct to four significant figures,
4.8The four eigenvalues *λ*_{i} of A will indicate the behaviour in homoclinic tails, where the four variables *q*_{i},*θ*_{i},*ω*_{i} and *ψ*_{i} are small. With changes in *F*, the movement of these eigenvalues follows the same trends as those of the continuous model of §3*c*, with three separate cases as follows:

For

*F*<45.98 (to four significant figures), all four eigenvalues are real. These appear in symplectic pairs so, numbering from largest to the smallest,*λ*_{1}*λ*_{4}=1 and*λ*_{2}*λ*_{3}=1. Together, they describe rates of growth or decay in non-wavy homoclinic tails.For 45.98<

*F*<49.41 (four significant figures), all four eigenvalues are complex and sit off the unit circle. This implies growth or decay in wavy homoclinic tails.For 49.41<

*F*, all four eigenvalues are complex and sit on the unit circle. This implies periodic solutions to the linear equation.

#### (i) Homoclinic boundary conditions

Here, we are mostly interested in case 1, which might be expected in the tails of a shear band. Following §3*c*(i), we would like to apply boundary conditions that represent either pure divergence from, or convergence to, the flat fundamental state as it is approached exponentially. For a second-order system, this is straightforward as there is a single eigenvector to cover each circumstance [28], but now there are two candidates for each with different rates of growth or decay. Remembering that *λ*_{1}>*λ*_{2}>1>*λ*_{3}>*λ*_{4}, the pure inset condition can be found by first setting the combination *ηλ*_{4}+(1−*η*)*λ*_{3} to zero and solving for *η*. The symplectic counterparts to *λ*_{4} and *λ*_{3} (*λ*_{1} and *λ*_{2}, respectively) can then be applied in the same combination (*ηλ*_{1}+(1−*η*)*λ*_{2}) to give an expected rate of growth for all four variables *q*_{i},*θ*_{i},*ω*_{i} and *ψ*_{i}. This means that *q*_{0} and *ω*_{0} (figure 9) can be written in terms of the degrees of freedom *q*_{1} and *ω*_{1}, and leads to an appropriate outset boundary condition to use within the continuation package Auto [15]. The opposite can be done for the inset condition at the other end.

### (d) Nonlinear solutions

A full set of nonlinear equilibrium equations can be generated by setting the first derivative of the total energy (4.3) with respect to each degree of freedom *q*_{i} or *ω*_{i} to zero, while applying the earlier-mentioned homoclinic boundary conditions. These can be solved in Auto as described by Hunt & Hammond [14], and solutions plotted in load/end-shortening space to mirror those of figure 8 for the continuous model. Our earlier paper [14] focused on the development of shear bands through a ‘snakes-and-ladders’ sequence of rising and falling equilibrium paths as extra discs are added to a shear band. Here, although the advanced post-buckling behaviour is much the same, we choose to concentrate on solutions relatively close to the flat state, to help to understand more fully the nucleation process. Figure 11 shows equilibrium solutions on a plot of load against end-shortening, and figure 12 six examples of deflected shapes on these paths, as the shear band develops. These have the same characteristics as the kink band plots of figure 8, i.e. rapidly falling load for very small end-shortening values.

Unlike for the continuous strut, we cannot yet proceed into the region of complex eigenvalues (case 2 above) as the programme is able to accept only real eigenvalues. However, it does seem likely that we have a situation like that of figure 8, with the mode shape emerging from the uppermost available value of *F*=45.98, at a small but finite maximum displacement of, in this case, about 0.4*R*. It is interesting to observe that, having adopted homoclinic boundary conditions, finite length effects are effectively removed; if extended beyond the ends of the sample in either the positive or negative sense, then the solution would continue to approach the flat state exponentially, with longer chains producing exactly the same solution over the given range. This suggests that, provided the sample is long enough, shear banding is typically free from boundary effects. Moreover, as the system gets longer the tails get flatter, and the relative distance of point (*a*), for example, in figure 11 would get closer to the fundamental path in comparison with the overall length. By this stage, the length scale of interest is more likely to be grain size than sample size, and a useful numerical test to run within Auto might be to vary *R* parametrically.

## 5. Concluding remarks

This paper does perhaps raise more questions than it answers. By extending the well-known beam equation to include through-thickness squashing, we have found good numerical evidence that a kink band can initiate from a finite load (*γ*=2) position on the fundamental (flat) equilibrium path. The fact that a finite critical load state can be identified perhaps represents a step forwards from earlier studies [12,2], although the post-buckling path that emerges is so thoroughly unstable, with the post-bifurcation load falling off so rapidly as the band width grows, that a computed critical load may have little bearing on real situations. We have also identified a similar scenario for our numerical granular model, with a shear band forming at an equivalent critical load of *F*=49.41. In each case, it is suggested that the band starts life with a bifurcation into a periodic waveform, which immediately starts to localize under rapidly falling load. These buckle patterns continue to retain wavy tails until reaching a load level of *γ*=−2 for the layered model or *F*=45.98 for the granular one, whereupon the tails lose their waviness. As the corresponding load changes take place over relatively small changes in applied end-shortening, the instability would be expected to be characterized by severe imperfection sensitivity.

From a practical point of view, this seems perhaps a rather unlikely scenario, and so some efforts have been expended looking also for a post-buckling path emerging in the form of a wrench instability, at *γ*=−2 for the continuous model and *F*=45.98 for the granular one. In each case, no immediate evidence is found, although the results of both studies are clearly of enough interest to warrant further examination. With nucleation occurring at a point in *x*, the wrench seems to have something in common with the recently identified phenomenon of a creasing instability [29]. Quite what this means for a discrete system such as the granular model, where grain size may be an issue, is at present an open question.

## Acknowledgements

We are grateful to Mark Peletier, Chris Budd and Alan Champneys for many useful discussions.

## Footnotes

One contribution of 17 to a Theme Issue ‘A celebration of mechanics: from nano to macro’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.