## Abstract

When twisting a strip of paper or acetate under high longitudinal tension, one observes, at some critical load, a buckling of the strip into a regular triangular pattern. Very similar triangular facets have recently been found in solutions to a new set of geometrically exact equations describing the equilibrium shape of thin inextensible elastic strips. Here, we formulate a modified boundary-value problem for these equations and construct post-buckling solutions in good agreement with the observed pattern in twisted strips. We also study the force–extension and moment–twist behaviour of these strips by varying the mode number *n* of triangular facets and find critical loads with jumps to higher modes.

## 1. Introduction

When twisting a strip of paper or acetate under high longitudinal tension, one observes, at a critical load, a buckling of the strip into a regular triangular pattern (figure 1*a*). The deformation is reversible. Sheets of paper or acetate are, for practical purposes, inextensible, and the observed pattern, consisting of helically stacked nearly flat triangular facets, appears to be Nature’s way of achieving global twisting by means of local bending and minimal stretch. This mode of buckling does not appear to have been reported in the literature. Perhaps this is because the phenomenon occurs at relatively large twisting angles and under relatively high tension (in order to suppress the more common looping instability). The buckling pattern observed has ridges running at roughly 45^{°} angles to the centreline of the strip. The ridges radiate out from vertices on the edge of the strip where stress concentration occurs. There is superficial similarity with the well-known Yoshimura or diamond buckling pattern of thin cylindrical shells (Yoshimura 1951), but this pattern requires compressive rather than tensile loading.

The phenomenon of buckling of an *extensible* elastic strip under tension and twisting moment has been studied before. In the 1930s, Green (1936, 1937) considered buckling of twisted strips under constant tensile force, treating both the case of zero and non-zero end force. Buckling of a twisted, helical, orthotropic plate into a sinusoidal buckling pattern in the longitudinal direction was studied by Crispino & Benson (1986). A numerical investigation of wrinkling of twisted plates was carried out by Mockensturm (2001), considering both the case of constant end-to-end distance and the case of constant end force, and found different results in the two cases. In the latter case, it was found that buckling may occur in both the lateral and longitudinal directions.

A perturbation method was used recently to further explore the wrinkling instability under small twist (Coman & Bassom 2008). As the twist increases from zero, the surface of a strip first becomes helicoidal, which causes extensional forces near the edges, while the core domain is in compression. When a critical twist and tension are reached, the core domain of the strip buckles into an oscillatory pattern. The surface near the edges remains helicoidal; no strain localization occurs.

The problem of the deformation of an *inextensible* flat plate was discussed at length by Clebsch (1862). Such a plate deforms isometrically and its surface is therefore developable (i.e. has zero Gaussian curvature), making the analysis more geometrical. Sadowsky (1931) developed a large-deformation theory of narrow elastic strips as early as in 1930. Approximate equations for wide strips were derived in the mid-1950s by Mansfield (1989). These equations predict the distribution of generators of the developable surface while ignoring the actual three-dimensional geometry. This work was followed up by Ashwell (1962), where localization of stresses at two diagonally opposite corners was found for a strip in its first buckling mode. The actual shape of the strip was not computed.

The constraint of zero Gaussian curvature causes the buckling patterns of inextensible strips to differ strongly from the oscillatory buckling patterns of extensible strips. The helicoidal shape of the edges of the strip found in Coman & Bassom (2008) is not a solution for inextensible strips, not even for infinitesimally small twist. Rather, we observe a sequence of relatively flat triangular domains that are not restricted to the core of the strip. The edges show a sequence of points with high curvature.

It is worth noting that both responses, the smooth sinusoidal one and the localized one, can be observed on the same paper-strip model, depending on the environmental conditions. When the humidity is low and the paper is dry, it behaves like an inextensible material. When it is slightly wet, it becomes noticeably stretchable. Note also that the solution for a slightly extensible strip is well described by the inextensible model almost everywhere except for small domains where the stress concentrates (Witten 2007). In this paper, we compute geometrically exact developable solutions of inextensible strips.

A geometrically exact set of equilibrium equations for the large deformation of thin inextensible plates of finite width has recently been derived (Starostin & van der Heijden 2007). The equations are ordinary differential equations (ODEs) and are obtained by using the inextensibility constraint to reduce stresses and strains to the centreline of the strip. They are much easier to analyse than the usual partial differential equations of elastic-plate theory. The new set of equations was used to solve a classical problem in mechanics, namely to find the equilibrium shape of a developable Möbius strip (Starostin & van der Heijden 2007). Numerical solutions revealed the existence, for any aspect ratio of the strip, of a nearly flat triangular region associated with the (unique) inflection point of the strip (figure 1*b*). The triangular facet of the Möbius-strip solution clearly resembles the facets of the buckling pattern of the twisted strip in figure 1*a*, and in this paper, we use the new system of equations to construct post-buckling solutions in good agreement with experiment. The central idea is to modify the boundary-value problem (BVP) for the Möbius strip such as to ‘cut out’ the triangular (more precisely, trapezoidal) region (figure 1*c*) and to use symmetry to reflect and multiply the elementary triangular facet into a periodic triangular pattern. This procedure avoids having to integrate numerically through the bending energy singularity associated with the vertex of the triangle on the edge of the strip, as found in Starostin & van der Heijden (2007) and analysed further in Hornung (2009). The buckling pattern is then obtained without having to compute bifurcations from an initially flat state, which in any case occurs at zero load and is a highly degenerate state.

Having obtained a post-buckling strip solution with a certain number, say *n*, of triangular facets, we then study the strip’s force–extension and moment–twist behaviour for various mode numbers *n* and infer some jump instabilities.

The paper is organized as follows. In §2, we formulate the BVP for the centreline-reduced equations for a developable strip. We also use symmetry properties of the solution to construct triangular buckling patterns for the strip by concatenating elementary facets. In §3, we numerically solve the BVP and compute load–displacement curves for the post-buckling strip solutions for various mode numbers *n*. In §4, we discuss our results and draw some conclusions.

## 2. The boundary-value problem

A sufficiently thin elastic surface will deform by bending only (Witten 2007), and will therefore deform isometrically. If such a surface is flat in its unstressed state, it will remain so under deformation and therefore have zero Gaussian curvature (Graustein 1966). It is said to be developable.

Now consider a rectangular sheet or strip. If **r**(*s*) is a parametrization of the centreline of the strip, with *s* being the arclength, then2.1is a parametrization of an embedded developable strip of length *l* and width 2*w* (Randrup & Røgen 1996). Here, **t** and **b** are two unit vectors of the Frenet frame {**t**,**n**,**b**}, tangent, principal normal and binormal to the centreline, while *κ* and *τ* are, respectively, the curvature and torsion of the centreline, which uniquely specify (up to Euclidean motions) the centreline of the strip (Graustein 1966). By equation (2.1), the surface, in turn, is completely determined by the centreline of the structure. The straight lines *s*=const. are the generators of the surface, which make an angle with the positive tangent direction. Given this parametrization, the mean curvature *M* can be easily calculated, e.g. by using the coefficients of the first and second fundamental forms of the surface, themselves calculated from partial derivatives of **x** with respect to *s* and *t* (Graustein 1966). The result is2.2where the prime denotes differentiation with respect to arclength *s*.

Introducing rectangular coordinates (*u*_{1},*u*_{2}) by developing the surface into a rectangle,2.3the bending energy of a strip of thickness 2*h* can be written as the following integral over the surface of the strip (Love 1927):2.4where *D*=2*Eh*^{3}/[3(1−*ν*^{2})] is the flexural rigidity, *E* is Young’s modulus and *ν* is Poisson’s ratio. When *M* (given by equation (2.2)) is substituted into equation (2.4), and the coordinates changed to *s* and *t*, the *t* integration can be carried out (Wunderlich 1962), giving2.5with2.6In the zero-width limit, *w*→0, this reduces to the Sadowsky result *g*=*κ*^{2}(1+*η*^{2})^{2} (Sadowsky 1931).

Minimization of this elastic energy functional is a one-dimensional variational problem cast in a form that is invariant under Euclidean motions. In Anderson (1989), Euler–Lagrange equations for such geometric variational problems are derived in a general context via a splitting of the cotangent bundle of the infinite jet bundle of a fibred manifold, which induces a bigrading of the differential forms on known as the variational bicomplex (see also Kogan & Olver 2003). In more physical terms, they can be written in the form of six balance equations for the (normalized) components of the internal force, * F*=(

*F*

_{t},

*F*

_{n},

*F*

_{b})

^{T}, and moment,

*=(*

**M***M*

_{t},

*M*

_{n},

*M*

_{b})

^{T}, in the directions of the Frenet frame and two scalar equations (Starostin & van der Heijden 2007, 2009)2.7 2.8where

**ω**=

*κ*(

*η*,0,1)

^{T}is the Darboux vector. Equations (2.7) are just the vectorial fixed-frame force and moment-balance equations

**F**′=

**0**,

**M**′+

**r**′×

**F**=

**0**, written out in the Frenet frame. (Here, we adopt the convention that bold roman symbols are used for vectors, while bold italic symbols are used for triples of components of these vectors in the Frenet frame.) It follows immediately that

*⋅*

**F***and*

**F***⋅*

**M***are first integrals of the equations. Note that the first of equations (2.8) is algebraic in the variables (*

**F***κ*,

*η*,

*η*′), while the second is a second-order ODE in

*η*.

It will be of interest to a more general audience how the same equations can be obtained from first principles using standard variational methods, extending to a function of *κ*,*τ*, *κ*′ and *τ*′ the examples in Capovilla *et al.* (2002) (which covers functional dependence up to *κ*,*τ* only). This is done in the appendix, where it is shown that it is straightforward to accommodate the additional functional dependence on *κ*′,*τ*′ without computing any new variations, by simply reusing the results already given in Langer & Perline (1991) and Capovilla *et al.* (2002). The additional terms generated by the dependence of *g* on *κ*′ and *τ*′ are easily managed. It is straightforward to extend this method to a functional involving any number of derivatives of *κ* and *τ*, and hence obtain the closed-form expressions given in Starostin & van der Heijden (2009).

The shape of the strip’s centreline is found by differentiating the first of equations (2.8) in order to turn it into a differential equation and then numerically solving equations (2.7) and (2.8) as a BVP in conjunction with three Euler-angle equations describing the evolution of the Frenet frame relative to a fixed frame, and the centreline equation **r**′=**t**. Taking Love’s convention for the Euler angles (*θ*,*ψ*,*ϕ*) (Love 1927), the derivatives of the angles are related to the curvature and torsion as2.9while the centreline equation in component form gives2.10Equation (2.9) can be written down directly from the Darboux vector by noting that Love’s convention is the usual *y*-convention (van der Heijden & Thompson 2000) in classical mechanics with *ψ* and *ϕ* interchanged (Goldstein 1980). Alternatively, it can be obtained from the Frenet–Serret equations. The convention is such that the polar singularity (here at *θ*=0) usually associated with Euler angles does not cause any problems for the solutions we are interested in (the flat strip will have *θ*=π/2). Alternatively a four-parameter quaternion representation can be used to avoid the singularity at the expense of a norm condition.

Before we specify the boundary conditions, let us have a closer look at the buckling pattern in figure 1*a* and the Möbius-strip solution in figure 1*b*. The Möbius-strip solution has special points where either *η* or *η*′ is zero. Points where *η*′=0 are called cylindrical points, as the surface is locally a cylinder, i.e. the mean curvature in equation (2.2) is constant along the local generator. In figure 1*b*, this corresponds to a generator of constant colour. Points where *η*=0, by contrast, are called conical because the edge of regression, on which nearby generators intersect each other, has a cusp. At these conical points, the generator is perpendicular to the centreline, as follows from equation (2.1). Clearly, there must be at least one point where *η*′=0 between any two points with *η*=0. In fact, the Möbius strip has three cylindrical points and three conical points, one of the latter being special because it corresponds to the only inflection point of the centreline (where *κ*=0). At this point, the binormal component of both the force and the moment are zero. Furthermore, it is found that at the inflection point, |*η*′|→1/*w* (i.e. the edge of regression reaches the edge of the strip) and the bending energy density *g* diverges, i.e. we have stress concentration. As is seen in figure 1*b*, coming out of this singular point is a nearly flat (violet) triangular region.

Now, turning to figure 1*a*, we observe that the buckling pattern consists of points of high stress located alternatingly on both edges of the strip, while locally cylindrical ridges bound flat triangular (more precisely, trapezoidal) regions similar to those found in the Möbius-strip solution. This suggests that we can describe the buckling pattern by a solution of the equations built up of alternating copies of the trapezoidal section between the inflection point (where *η*=0) and the nearest cylindrical point (where *η*′=0). A cut-out of this section is shown in figure 1*c*. Note that both bounding generators are of constant colour, illustrating that the section can be reflected about both end generators.

A two-step symmetry operation is therefore used to construct a strip of length 2*nL* from a trapezoid of length *L*. Let the arclength parameter of the centreline be *s*=0 at the cylindrical point (*η*′=0) and *s*=*L* at the inflection point (*η*=0) (figure 1*c*). The first operation is a rotation through 180^{°} of the trapezoid about the normal **n**_{1} at *s*=0 (see also figure 2). The original and the rotated trapezoid together make a continuous surface, which forms one period of the full strip (*P*_{1} in figure 2), of length 2*L*. The binormal to the strip at the inflection point is denoted **b**_{0}.

In the second step, a symmetrical strip of *n*+1 periods is obtained from a parent strip of *n* periods by rotating an end period of the parent strip 180^{°} about the end binormal. This binormal is located at the inflection point of that end period and therefore aligned with the end generator. Thus, strips of any period can be built up in this way by successively reflecting an end period about its end binormal. The length of the centreline of a symmetrical strip of *n* periods made in this way is therefore 2*nL*. In particular, let **r**(0)≡**r**_{0}, **r**(*L*)≡**r**_{1} be, respectively, the position vectors of the centreline at the inflection and cylindrical points of the trapezoid solution of the BVP. Define a rotation of π around the unit vector **g**,2.11then the first and subsequent periods are obtained by the iteration2.12where *i*=2,…,*n*. The **b**_{2i} are then the unit binormals at the inflection points of the periods of the strip. The reflection rules in equation (2.12) are for the centreline, but can be easily generalized to the whole surface, as was done in Starostin & van der Heijden (2007) to obtain closed-strip solutions from the trapezoid solution of the BVP. These reflection rules are shown in figure 2 for the example *n*=4 with the relevant binormals shown. The initial trapezoid (with one boundary aligned with **b**_{0}) is rotated about **n**_{1} giving the first period *P*_{1}. This first period is then rotated about **b**_{2} to give the second period *P*_{2}, which itself is rotated about **b**_{4} to give the third period *P*_{3}. Finally, the third period is rotated about **b**_{6} to give the fourth period *P*_{4}.

It remains to formulate the boundary conditions for the initial trapezoid. At *s*=0, these are2.13where the values for *x*, *y* and *z* fix an arbitrary position in space, and2.14which can be read off as the *s*=0 limit of the first of equations (2.8), where *η*′(0)=0. Equation (2.14) is required to fix the integration constant for the ODE obtained by differentiating the first of equations (2.8). (Note that by Taylor expanding the second of equations (2.8) around *s*=0, one can also show that *M*_{t}(0)=(2/3)(1+*η*^{2}(0))*κ*(0)(−6*η*(0)+*w*^{2}(1+*η*^{2}(0))*η*′′(0)).)

The boundary conditions at *s*=*L* are2.15where the angles fix an arbitrary orientation of the strip in space.

The force and moment boundary conditions in equations (2.13) and (2.15) are enforced by the rotational symmetry. Since **n**_{1} and **b**_{0} are local axes of rotational (or reflection) symmetry, any component of force and moment along those axes must vanish. This guarantees continuity of forces and moments in the strip and therefore yields a valid solution of the mechanical problem.

The remaining boundary conditions for the system of ODEs are the projections of the end force and moment of the strip on the end-to-end unit vector , **e**=**r**_{2n}−**r**_{0},2.16where and are also continuation (control) parameters. Equations (2.13)–(2.16) give a final count of 15 boundary conditions for the same number of first-order ODEs.

The BVP is solved by continuation of each of the boundary conditions in equation (2.16) in turn. To obtain the starting solution for the continuation, the numerical solution of the Möbius strip between the inflection point and nearest cylindrical point is continued in *F*_{n}(0) and *M*_{n}(0) until equations (2.13) are satisfied.

There are significant numerical difficulties solving this BVP, as both ends of the integration interval have a singularity. The boundary condition *κ*(*L*)=0 enforces an inflection point at *x*=*L*. At this point, the torsion *τ* must also be zero and, in fact, it must go to zero, in *s*, faster than the curvature *κ* (Randrup & Røgen 1996). Therefore, we also have *η*(*L*)=0. In practice, to compute a starting solution, we first compute an approximate solution with *κ*(*L*)≃0.1 to stay away from the singularity at *x*=*L*. When all other boundary conditions are satisfied we ‘pull’ the solution into the singularity by continuing *κ*(*L*) to zero as far as possible, typically reaching values of 0.001. At this point, we typically have *η*(*L*)≃0.002, while 1/*w*−|*η*′(*L*)|, the distance from the singularity, is typically as small as 10^{−6}.

At the other singularity, at *s*=0, numerical convergence requires Taylor expansions (up to fourth order in our case) of the left-hand sides of equations (2.8) about *η*′=0 to be used for a small interval at *s*=0. In addition, to improve convergence, the absolute value of the argument of logarithm is taken in equation (2.6). We note that it has not proved possible (for instance by similar Taylor expansions) to remove the singularity in equations (2.8) at |*η*′|=1/*w*.

As a check on the numerical results, the first integrals * F*⋅

*and*

**F***⋅*

**M***are typically found to be constant to within 10*

**F**^{−9}.

Once a solution to the BVP is obtained for the centreline of the first half period, the surface of the full strip is found from equation (2.1) and the reflection rules (2.12). We note that no measures in the above formulation are taken to prevent self-intersection of the strip.

## 3. Numerical results: response of the strip to applied loads

We study the force–extension and moment–twist behaviour of strips with the same aspect ratio but differing number of periods, *n* (i.e. different modes). Since it is the response of the full strip that is of interest, we choose the end-to-end distance |**e**| and the accumulated twist angle *α* as representative independent variables, which might also be used as controls in an experiment. As the component of **b**_{0} perpendicular to **e** is , with **a**_{2n} similarly defined, *α* is defined as the angle between the components of the two end-of-strip binormals perpendicular to the end-to-end vector of the strip, i.e. . Thus, |**e**| and *α* are both calculated from the reflection rules (2.12).

The AUTO continuation software (Doedel *et al.* 2007) was used to compute response curves (at constant ) of applied force against |**e**| and response curves (at constant ) of applied moment against *α*. Solutions for strips under high tension are therefore obtained by continuation from the starting solution, which is under low tension. In the case of the force–extension behaviour, further response curves are obtained by fixing for any solution, continuing in to another fixed value of , then continuing in to obtain a new force–extension curve at the new fixed value of . A similar scheme is used for the moment–twist behaviour. In this way, the force–extension and moment–twist behaviour for the full strip is obtained from the solution to the BVP.

Under a distance rescaling *s*→*s*′=*sk*, the force and moment scale as ** F**→

**/**

*F**k*

^{2}and

**→**

*M***/**

*M**k*(

**scales as**

*M**κ*, cf. equation (2.14)),

*κ*scales as

*κ*→

*κ*/

*k*, whereas

*η*and the Euler angles are scale invariant. Rescaling can be used to obtain strip solutions of a given period from solutions of another different period, with the same aspect ratio, via a

*w*continuation. For example, starting from a strip of length 2

*nL*, width 2

*w*and period

*n*=2

*m*, selecting

*n*/2 consecutive periods gives a strip of length

*nL*, continuing this solution to

*w*and rescaling

*s*→2

*s*gives a period

*n*/2 strip with the same aspect ratio as the original strip. All strips are scaled to an aspect ratio of 2

*nL*/2

*w*=10.5300 unless otherwise stated.

Many different strip shapes are obtained depending on the boundary conditions (2.16). Figure 3 shows three shapes of a strip with *n*=8, one under a high axial tension, one under a relatively high axial moment and an intermediate one that is in reasonably good agreement with the experiment in figure 1*a*. Here, and in all surface plots in this paper, the colouring varies according to the bending energy density, from violet for regions of low bending to red for regions of high bending (scales are individually adjusted).

Continuation results for modes *n*=2, 4, 8 are shown in figures 4–6 (for aspect ratio 2*nL*/2*w*=10.5300, except for figure 6*a*, which has an aspect ratio 42.30). On the left of each figure is the force response for the sequence of axial moments () shown on the legend. The length scale is arbitrary and is such that *L*=0.6581. On the right of each figure is shown the moment response for the sequence of axial forces () shown on the legend. These curves were obtained by continuation, where was varied, keeping fixed or vice versa. Exceptions to this are some curves in figures 5*a*,*b* and 6*a*, where both and were allowed to vary, only so that a continuation could be started at an arbitrary point on the graph (at which point one of or would then be held fixed). This accounts for the curve-crossing in these figures.

In figure 4*a*, for the *n*=2 mode, the sequence of curves at constant negative applied end moment is the group at the bottom left, with the moment increasing from top (large negative) to bottom (near zero). The group of curves at the bottom right of the figure is for constant positive applied end moment, with the moment increasing from bottom (near zero) to top (large positive). In figure 4*b* is a sequence of curves at constant applied end force, increasing from bottom (negative) to top (positive). The curves for smallest (negative) constant force start at zero moment and have a maximum twist angle before returning to zero twist angle where the solution is the figure-of-eight shown below.

In figure 5*a*, for the *n*=4 mode, the sequence of curves at constant negative applied end moment is the group at the bottom left, with the moment increasing from top (large negative) to bottom (near zero). The group of curves at the bottom right of the figure is for constant positive applied end moment, with the moment increasing from bottom (near zero) to top (large positive). In figure 5*b* is a sequence of curves at constant applied end force, increasing from bottom (negative) to top (positive), showing the response curve of applied end moment to the twist angle *α*. The three lowest negative applied force continuations terminate at *α*=2π, again forming a closed strip, this time the *T*_{4,1} elastic torus ribbon knot.

In figure 6*a*, for the *n*=8 mode, is a sequence of curves at constant applied end moment (positive and negative). This picture is more complicated than the other modes. Some of the curves terminate on the vertical axis, i.e. when the end-to-end distance vanishes. These closed-ribbon solutions are invariably torus ribbon knots. For example, the upper dashed curve () terminates in the double cover of the *T*_{4,1} elastic torus ribbon knot. The series of nested closed curves have a quadruple cover of the figure-of-eight solution where the curves approach vanishing end-to-end distance. In figure 6*b* is a sequence of curves at constant applied end force, increasing from bottom (negative) to top (positive), showing the response curve of applied end moment to the twist angle *α*.

In figure 7, for the *n*=8 mode, is a sequence of curves at constant applied end moment, with some of the resulting structures shown. The upper curve () terminates in the *T*_{8,1} elastic torus ribbon knot. Both ends of the dot-dashed curve () terminate in a double covering of the *T*_{4,1} elastic torus ribbon knot, whereas the *T*_{8,3} structure (shown at the minimum end-to-end distance of the loop) does not quite close, as can be seen by the fact that the curve from which it originates does not meet the vertical axis. The other structure shown on the same loop in the line diagram (shown at the maximum end-to-end distance of the loop on the figure) is related to the *T*_{8,3} torus knot, but it is also not closed, and moreover the surface of the strip is self-intersecting.

As an example of further strip solutions, in figure 8*b* is shown a solution, with *n*=8 and under a compressive axial force, that is not located on any of the computed response curves.

Figure 9 displays force–extension curves for varying mode number, *n*=2,…,8, at fixed moment . To obtain this plot, strips for different modes have been scaled to the same aspect ratio 2*nL*/2*w*=10.5300. Parts of these curves with positive slope are expected to correspond to stable solutions. The curves predict that, for *n*≥5, under increasing tension, solutions jump to higher mode (down in |**e**|) as the force is increased beyond the folds seen in the diagram.

## 4. Discussion

When twisted, and pulled, an acetate model strip buckles into a regular pattern of triangular facets. We have computed periodic solutions describing this buckling pattern by formulating and numerically solving a geometrically exact BVP for the large deformation of a thin, wide, inextensible strip. We have also obtained response curves of force against end-to-end distance and twisting moment against end-to-end angle for mode numbers *n* up to 8. The computations can be straightforwardly extended to higher modes. Our results predict critical forces and jumps into higher buckling modes that would be interesting to explore experimentally.

We note that post-buckling solutions have been obtained without considering bifurcation, and branch switching, from the trivial flat state. Indeed, the flat state is highly degenerate in our formulation, as it has zero curvature everywhere, and therefore the Frenet frame and the torsion are undefined. Our solutions by construction have an inflection point at one of the ends of the elementary period, and this is always accompanied by a divergence of the bending energy density. Our numerical results suggest that periodic solutions with these singularities approach the flat state arbitrarily closely, suggesting that these solutions do bifurcate from the flat state. This bifurcation would have to occur at zero load, as any applied twisting moment (or twist angle *α*) would deflect the strip from its flat state.

Other solutions may exist. These may or may not have singularities, but experiments suggest that solutions with singularities (stress localization) dominate if a relatively large tensile force is applied. If we consider boundary conditions without an applied axial moment, then twistless solutions with a planar centreline may be obtained. Such solutions could be studied by setting *τ*, or *η*, equal to zero from the start (in equation (2.2)) leading to the Euler elastica equation. In this case, no singularities will occur.

By construction, our solutions are periodic, which tends to be what one observes in experiments. However, non-periodic solutions can be constructed in a similar way by matching different trapezoidal segments. One would keep the first symmetry operation (reflection about the normal at the cylindrical point), but instead of reflecting about the binormal at the inflection point in the second step, one could match the segment to a (suitably rescaled) trapezoidal segment of different length *L*. The resulting solution would no longer be symmetric about the inflection points.

Our analysis is made possible by the fact that the inextensibility constraint can be used to reduce the problem to finding the centreline of the strip, which is formulated as a system of ODEs. This explicit reduction may not be possible if effects not considered here, such as gravity or other distributed forces acting on the strip, are included. For extensible strips, no such reduction would be possible, and one would have to solve the partial differential equations of plate theory. For small amounts of extensibility, one would expect that solutions for twisted pulled strips would still show stress concentration on the edge of the strip, but without the strain energy density actually going infinite.

In Ashwell (1962) and Mansfield (1989), analyses were performed similar to ours in that there too the double strain energy integral was reduced by integrating along the generator. However, the integrand expression was subsequently simplified by using an approximate moment balance equation that ignores the non-planarity of the strip. This approach gives a realistic evolution of the generator with arclength, but fails to predict the appearance of a cylindrical point at finite distance from the clamped end of the strip.

Triangular patterns are known to occur in a variety of problems of elastic sheets, including fabric draping and paper crumpling (Witten 2007). A sheet crumples (forming sharp points and straight creases) when it is forced into a constrained domain. The sheet is predominantly under compression. It is interesting to note that, by constrast, in our triangular buckling pattern, the strip is in (relatively high) tension. In both cases, we observe a focusing of the strain energy that may lead to fracture of the material. Strain-energy localization thus appears to be a generic response of a thin elastic sheet to an external constraint.

Our results may be relevant for paper, fabric and sheet-metal processing. They may also be of importance in the robotic manipulation of flexible belts, e.g. film circuit boards (Wakamatsu *et al.* 2007). In all these sheet manipulations, it is desirable to avoid shapes with high concentrated stresses that may lead to tearing. Our formulation may help to choose boundary conditions that avoid unwanted configurations.

## Acknowledgements

This work was supported by the UK’s Engineering and Physical Sciences Research Council (EPSRC) under grant no. EP/F023383/1.

- Received April 13, 2010.
- Accepted May 26, 2010.

- This journal is © 2010 The Royal Society

## Appendix A. Derivation of the equations using standard variational techniques

We wish to calculate the variation of the functional . This calculation can be done using the formulation of Capovilla *et al.* (2002), substituting explicit expressions for quantities wherever they arise and grouping similar terms, which gives fully simplified expressions for all quantities of interest.

Alternatively, one can re-use the variations already calculated in Capovilla *et al.* (2002) to yield the same results, without having to compute new variations. There, for example, the variation is given for , stating that it can be written down from the variations of and , which are previously calculated. The resulting force and moment corresponding to *H*_{1} can then be written down from the expressions for the force and moment for *H*_{2} and *H*_{3}. This follows from the chain rule. Similarly, one can calculate the variation of the functional using *δH*_{1} without computing new variations. This follows as a consequence of the chain and product rules in differentiation. The resulting force and moment corresponding to *H* can then be written down from the expressions for the force and moment for *H*_{1}, with a small number of additional terms arising from the chain and product rules.

Instead of following these approaches, we here use a third and more economical method by leaving *δH* expressed in terms of *δκ* and *δτ* only. This method is easily generalized to functionals involving higher derivatives of *κ* and *τ*. Following the notation of Capovilla *et al.* (2002), the infinitesimal deformation of a space curve is *δ***r**=*Ψ*_{‖}**t**+*Ψ*_{1}**n**+*Ψ*_{2}**b**, and denoting *δ*_{‖},*δ*_{⊥} as the tangential and normal parts of the deformation, respectively, the variation of any functional isA1where . We therefore need to calculate the second term on the right-hand side of equation (A 1), where by the chain rule, *δ*_{⊥}*f*=*f*_{κ}*δ*_{⊥}*κ*+*f*_{τ}*δ*_{⊥}*τ*+*f*_{κ′}*δ*_{⊥}*κ*′+*f*_{τ′}*δ*_{⊥}*τ*′, using the notation *f*_{β}≡∂*f*/∂*β*.

Since for any scalar *h*, *δ*_{⊥}(*h*′)=*κh*′*Ψ*_{1}+(*δ*_{⊥}*h*)′, thenA2where *δ*_{⊥}*κ* and *δ*_{⊥}*τ* are given in Capovilla *et al.* (2002). Equation (A 1) along with equation (A 2) is of the formA3with the Noether charge, so that the Euler–Lagrange equations can be immediately written down, since these are just the coefficients of *Ψ*_{i}, which can be read off from equation (A1). These are two coupled differential equations in the unknowns *κ* and *τ*, and the known density *f*, which are quoted in general form in Thamwattana *et al.* (2008) and in Hangan (2005) for the Sadowsky functional.

Apart from the Euler–Lagrange equations, it is of interest to find expressions for the conserved force **F** and the moment **M**. The force **F** is obtained by specializing the deformation to a constant infinitesimal translation *δ***r**=**e**, where **F** is defined by . Similarly, the conserved moment **T** is obtained by specializing the deformation to a constant infinitesimal rotation *δ***r**=**Ω**×**r**, where **T** is defined by and decomposed into **T**=**r**×**F**+**M**. Thus, only the total derivative parts of equation (A 1) contribute to **F** and **M**. But these contributions can be read off from the results in Capovilla *et al.* (2002) by noting that a term in equation (A 2) of the form *aδ*_{⊥}*b*, where *b*=*κ*,*τ*, has to be re-expressed in the form of equation (A 3) to isolate its total derivative part. This has already been done in Capovilla *et al.* (2002). For a term in equation (A 2) of the form (*aδ*_{⊥}*b*)′, one just needs the contribution from *δ*_{⊥}*b*.

Thus, for an infinitesimal rotation, *δ***r**′=**Ω**×**r**′=**Ω**×**t** and *δ***r**′′=*κ***Ω**×**n**. Using the compact expressions of Langer & Perline (1991, eqn (3.7b,c)), *δκ*=*δ***r**′′⋅**n**−2*κδ***r**′⋅**t** gives *δκ*=0, i.e. *δ*_{⊥}*κ*=−*δ*_{‖}*κ*=−*Ψ*_{‖}*κ*′. But *Ψ*_{‖}≡*δ***r**⋅**t** = (**r**×**t**)⋅**Ω**, giving *δ*_{⊥}*κ*=−**Ω**⋅(*κ*′**r**×**t**). As this is of the form −**Ω**⋅(**r**×**F**), in equation (A 2) does not contribute to the moment **M**. Similarly, *δτ*=*δ***r**′′⋅**b**/*κ*′+*δ***r**′⋅(*κ***b**−*τ***t**)=0, giving *δ*_{⊥}*τ*=−**Ω**⋅(**r**×*τ*′**t**). As this is of the form −**Ω**⋅(**r**×**F**), in equation (A 2) also does not contribute to the moment **M**. (Note that by retaining these terms, one would obtain the force **F**, but in the next paragraph we will, instead, calculate **F** by considering a constant infinitesimal translation.) Now, from Capovilla *et al.* (2002, eqns (49) and (67)), *aδ*_{⊥}*κ* contributes −*a***b** to the moment, whereas *aδ*_{⊥}*τ* contributes −*a***t**−*a*′**n**/*κ*. Thus, from equations (A 1) and (A 2), one can immediately write down the moment asA4In order to calculate the contribution to the force of , one notes that since *δ***r**=**e** is a constant, it follows from the formulation of Langer & Perline (1991, eqn (3.7b,c)) that *δκ*=0, i.e. *δ*_{⊥}*κ*=−*δ*_{‖}*κ*=−*Ψ*_{‖}*κ*′=−**e**⋅(*κ*′**t**). Thus, the contribution of in (A 2) to **F** is *κ*′*f*_{κ′}**t**. In a similar manner, for this constant translation, *δτ*=0, giving *δ*_{⊥}*τ*=−**e**⋅(*τ*′**t**), so that in equation (A 2) contributes *τ*′*f*_{τ′}**t** to the force. Now from Capovilla *et al.* (2002) (eqns (48) and (64)), *aδ*_{⊥}*κ* contributes *κa***t**+*a*′**n**+*τa***b** to the force, whereas *aδ*_{⊥}*τ* contributes *τa***t**+*τa*′**n**/*κ*−((*a*′/*κ*)′+*κa*)**b**. Also, *δ*_{‖}*H* contributes −*f***t** to the force. Thus, from equations (A 1) and (A 2), one can immediately write down the force asA5

From these components of **F** and **M**, the differential equations for the moment components and the differential equation for *F*_{t} in equations (2.7) are easily obtained. For instance, from equation (A 4), it is seen that *M*_{t}^{′}=*κM*_{n}, while equations for *M*_{n}, *M*_{b} and *F*_{t} can be extracted similarly. The remaining two equations for *F*_{n} and *F*_{b} are obtained from equation (A 5) and the Euler–Lagrange equations .

In summary, the variation of can be easily obtained from the variations *δκ* and *δτ*. Their explicit expressions are not required; their contributions can be read off from Capovilla *et al.* (2002). The additional terms generated by the additional dependence of *f* on *κ*′ and *τ*′ are therefore easily managed.

One way of regarding equations (A 4) and (A 5) is that they prescribe **F** and **M** once *κ* and *τ* for the given density *f* are known. The two Euler–Lagrange equations for *κ* and *τ* are given by , i.e. by setting the coefficients of *Ψ*_{i} to zero in equation (A 3). Instead of solving the problem this way, however, it is preferable to set up an equivalent system of coupled one-dimensional ODEs as in Starostin & van der Heijden (2009). One way to do this is to use the natural variables * F* and

*, and, for the functional (2.6), to change variables from (*

**M***κ*,

*τ*) to (

*κ*,

*η*).

By considering functionals of the form *g*(*κ*,*η*,*η*′)=*f*(*κ*,*τ*,*κ*′,*τ*′), one can show that (∂*g*/∂*κ*)_{η,η′}=*f*_{κ}+*ηf*_{τ}+*η*′*f*_{τ′}, (∂*g*/∂*η*)_{κ,η′}=*κf*_{τ}+*κ*′*f*_{τ′} and (∂*g*/∂*η*′)_{κ,η}=*κf*_{τ′}. Using these identities, and writing the components of the internal force and moment in the directions of the Frenet frame of tangent, principal normal and binormal as * F*=(

*F*

_{t},

*F*

_{n},

*F*

_{b})

^{T},

*=(*

**M***M*

_{t},

*M*

_{n},

*M*

_{b})

^{T}, one can show, using the tangent and binormal components of

*from equation (A 4), thatA6The two scalar equations (A 6), along with the six balance equations (equations (2.7)) for the components of the internal force*

**M***and moment*

**F***(see Starostin & van der Heijden 2007, 2009), are then the equations satisfied by the centreline of the developable strip for the unknowns*

**M***,*

**F***,*

**M***κ*,

*η*,A7 A8where

**ω**=

*κ*(

*η*,0,1)

^{T}is the Darboux vector. Equations (A 7) and (A 8) constitute a system of differential-algebraic equations that are turned into a system of ODEs by differentiation of the first of equations (A 8).