Triangular buckling patterns of twisted inextensible strips

Proc. R. Soc. A doi:10.1098/rspa.2010.0200 Published online Triangular buckling patterns of twisted inextensible strips BY A. P. KORTE, E. L. STAROSTIN AND G. H. M. VAN DER HEIJDEN* Centre for Nonlinear Dynamics, University College London, London WC1E 6BT, UK 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. Keywords: twisted inextensible strip; triangular buckling pattern; calculus of variations; boundary-value problem; load–displacement behaviour 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 1a). 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 *Author for correspondence ([email protected]). Received 13 April 2010 Accepted 26 May 2010 1 This journal is © 2010 The Royal Society 2 A. P. Korte et al. (a) (b) (c) n s=L b s=0 Figure 1. (a) Twisted acetate model strip under tension. (b) Möbius developable structure of aspect ratio 2p with generators shown. (c) Trapezoidal segment taken from the Möbius structure to construct the periodic strip, with normal, binormal and centreline shown, h′ (0) = 0, h(L) = 0. The colouring changes according to the local bending energy density, from violet for regions of low bending to red for regions of high bending. Note the singularity on the edge of the strip. 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 Proc. R. Soc. A Buckling of inextensible strips 3 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 1b). The triangular facet of the Möbiusstrip solution clearly resembles the facets of the buckling pattern of the twisted strip in figure 1a, 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 1c) 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 Proc. R. Soc. A 4 A. P. Korte et al. 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, then x(s, t) = r(s) + t[b(s) + h(s)t(s)], t(s) = h(s)k(s),  s ∈ [0, l], t ∈ [−w, w], (2.1) is a parametrization of an embedded developable strip of length l and width 2w (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 k and t 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 arctan(1/h) 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 is M =− k 1 + h2 , 2 1 + th′ (2.2) where the prime denotes differentiation with respect to arclength s. Introducing rectangular coordinates (u1 , u2 ) by developing the surface into a rectangle, u1 = s + th(s) and u2 = t, (2.3) the bending energy of a strip of thickness 2h can be written as the following integral over the surface of the strip (Love 1927): U = 2D  M 2 du1 du2 , (2.4) where D = 2Eh 3 /[3(1 − n2 )] is the flexural rigidity, E is Young’s modulus and n is Poisson’s ratio. When M (given by equation (2.2)) is substituted into Proc. R. Soc. A Buckling of inextensible strips 5 equation (2.4), and the coordinates changed to s and t, the t integration can be carried out (Wunderlich 1962), giving L U = Dw g(k, h, h′ )ds, (2.5) 0 with ′ 2  g(k, h, h ) = k 1 + h  2 2   1 1 + wh′ log . 2wh′ 1 − wh′ (2.6) In the zero-width limit, w → 0, this reduces to the Sadowsky result g = k2 (1 + h2 )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 T ∗ J ∞ of the infinite jet bundle J ∞ of a fibred manifold, which induces a bigrading of the differential forms on J ∞ 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 = (Ft , Fn , Fb )T , and moment, M = (Mt , Mn , Mb )T , in the directions of the Frenet frame and two scalar equations (Starostin & van der Heijden 2007, 2009) F ′ + u × F = 0, vk g + hMt + Mb = 0 M ′ + u × M + t × F = 0, and ′ (vh′ g) − vh g − kMt = 0, (2.7) (2.8) where u = k(h, 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 · F and M · F are first integrals of the equations. Note that the first of equations (2.8) is algebraic in the variables (k, h, h′ ), while the second is a second-order ODE in h. 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 k, t, k′ and t′ the examples in Capovilla et al. (2002) (which covers functional dependence up to k, t only). This is done in the appendix, where it is shown that it is straightforward to accommodate the additional functional dependence on k′ , t′ 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 k′ and t′ are easily managed. It is straightforward to extend this method to a functional involving any number of derivatives of k and t, 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 Proc. R. Soc. A 6 A. P. Korte et al. Euler angles (q, j, f) (Love 1927), the derivatives of the angles are related to the curvature and torsion as q′ = k cos f, sin f sin q ′ f = −k sin f cot q + kh, j′ = k and (2.9) while the centreline equation in component form gives x ′ = sin q cos j, y ′ = sin q sin j and (2.10) ′ z = cos q. Equation (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 j and f interchanged (Goldstein 1980). Alternatively, it can be obtained from the Frenet–Serret equations. The convention is such that the polar singularity (here at q = 0) usually associated with Euler angles does not cause any problems for the solutions we are interested in (the flat strip will have q = p/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 1a and the Möbius-strip solution in figure 1b. The Möbius-strip solution has special points where either h or h′ is zero. Points where h′ = 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 1b, this corresponds to a generator of constant colour. Points where h = 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 h′ = 0 between any two points with h = 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 k = 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, |h′ | → 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 1b, coming out of this singular point is a nearly flat (violet) triangular region. Now, turning to figure 1a, 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 h = 0) and the nearest cylindrical point (where h′ = 0). A cut-out of this section is shown in figure 1c. Note that both bounding generators are of constant colour, illustrating that the section can be reflected about both end generators. Proc. R. Soc. A 7 Buckling of inextensible strips b6 n1 b2 r2 r1 P4 P3 b8 P2 r0 P1 b4 b0 Figure 2. Mode n = 4 symmetry operations. The first half period is rotated about n1 to form the first period P1 . Rotating P1 about b2 produces the second period P2 . Rotating P2 about b4 , then P3 about b6 , produces the final two periods of the figure. A two-step symmetry operation is therefore used to construct a strip of length 2nL from a trapezoid of length L. Let the arclength parameter of the centreline be s = 0 at the cylindrical point (h′ = 0) and s = L at the inflection point (h = 0) (figure 1c). The first operation is a rotation through 180◦ of the trapezoid about the normal n1 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 (P1 in figure 2), of length 2L. The binormal to the strip at the inflection point is denoted b0 . 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 2nL. In particular, let r(0) ≡ r0 , r(L) ≡ r1 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 p around the unit vector g, Rg (a) = 2g(g · a) − a, then the first and subsequent periods are obtained by the iteration  r2 = r1 + Rn1 (r0 − r1 ), b2 = Rn1 (b0 ), r2i = r2i−2 + Rb2i−2 (r2i−4 − r2i−2 ) and b2i = Rb2i−2 (b2i−4 ), (2.11) (2.12) where i = 2, . . . , n. The b2i 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 Proc. R. Soc. A 8 A. P. Korte et al. one boundary aligned with b0 ) is rotated about n1 giving the first period P1 . This first period is then rotated about b2 to give the second period P2 , which itself is rotated about b4 to give the third period P3 . Finally, the third period is rotated about b6 to give the fourth period P4 . It remains to formulate the boundary conditions for the initial trapezoid. At s = 0, these are  Fn (0) = 0, Mn (0) = 0, h′ (0) = 0, (2.13) x(0) = 0, y(0) = 0 and z(0) = 0, where the values for x, y and z fix an arbitrary position in space, and Mb (0) = −h(0)Mt (0) − 2k(0)(1 + h2 (0))2 , (2.14) which can be read off as the s = 0 limit of the first of equations (2.8), where h′ (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 Mt (0) = (2/3)(1 + h2 (0))k(0)(−6h(0) + w 2 (1 + h2 (0))h′′ (0)).) The boundary conditions at s = L are ⎫ k(L) = 0, Fb (L) = 0, Mb (L) = 0, ⎬ (2.15) p q(L) = , j(L) = 0 and f(L) = p,⎭ 2 where 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 n1 and b0 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/|e|, e = r2n − r0 , F(L) · ê = F¯ and M(L) · ê = M̄ , (2.16) where F¯ and M̄ 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 Fn (0) and Mn (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 k(L) = 0 enforces an inflection point at x = L. At this point, the torsion t must also be zero and, in fact, it must go to zero, in s, faster than the curvature k (Randrup & Røgen 1996). Therefore, we also have h(L) = 0. In practice, to compute a starting solution, we first compute an approximate solution with k(L) ≃ 0.1 to stay away from the Proc. R. Soc. A Buckling of inextensible strips 9 singularity at x = L. When all other boundary conditions are satisfied we ‘pull’ the solution into the singularity by continuing k(L) to zero as far as possible, typically reaching values of 0.001. At this point, we typically have h(L) ≃ 0.002, while 1/w − |h′ (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 h′ = 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 |h′ | = 1/w. As a check on the numerical results, the first integrals F · F and M · F are typically found to be constant to within 10−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 a as representative independent variables, which might also be used as controls in an experiment. As the component of b0 perpendicular to e is a0 ≡ b0 − (b0 · ê)ê, with a2n similarly defined, a 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. cos a = a0 · a2n /(|a0 ||a2n |). Thus, |e| and a 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 M̄ ) of applied force F¯ against |e| and response curves (at constant F¯ ) of applied moment M̄ against a. 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 F¯ for any solution, continuing in M̄ to another fixed value of M̄ , then continuing in F¯ to obtain a new force–extension curve at the new fixed value of M̄ . 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 (M scales as k, cf. equation (2.14)), k scales as k → k/k, whereas h 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 2nL, width 2w and period n = 2m, selecting n/2 consecutive periods gives a strip of length nL, continuing this solution to w and rescaling s → 2s gives a period n/2 strip with the same aspect ratio as the original strip. All strips are scaled to an aspect ratio of 2nL/2w = 10.5300 unless otherwise stated. Proc. R. Soc. A 10 A. P. Korte et al. (a) (b) (c) Figure 3. Shapes of strips for n = 8: (a) F¯ = 13.42, M̄ = 0.162; (b) F¯ = 6.80, M̄ = 2.71; and (c) F¯ = 0.197, M̄ = 12.15. (a) 14 (b) – M = –4.15 – M = –3.19 – M = –1.80 – M = –1.19 – M = –0.50 – M = 0.04 – M = 0.52 – M = 1.04 – M = 1.53 – M = 2.01 – M = 2.48 – M = 2.99 – M = 7.25 12 10 8 – F 6 4 – 30 25 20 – M 15 10 2 F= –0.9119 – F = –0.3986 – F = –0.3126 – F = –0.1874 – F = –0.0009 – F = 0.3931 – F = 0.7893 – F = 1.1856 – F = 1.5813 – F = 1.9604 – F = 2.3568 5 0 0 −2 0 2 4 6 |e| 8 10 0 1 2 3 α 4 5 6 Figure 4. n = 2 mode. (a) Force F¯ versus scaled end-to-end distance |e|. The curves of constant M̄ form a nested sequence on the plot, with M̄ increasing from the top middle dashed curve (M̄ = −4.15) to the dotted curve, then from the dot-dashed curve to the dashed curve on the right. (b) Moment M̄ versus twist angle a (in radians). The curves of constant F¯ form a nested sequence on the plot, with F¯ increasing from the dashed curve (F¯ = −0.9119) to the dotted curve. 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 1a. 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 2nL/2w = 10.5300, except for figure 6a, which has an aspect ratio 42.30). On the left of each figure is the force response for the sequence of axial moments (M̄ ) shown on the legend. The length scale is arbitrary and is such that L = 0.6581. Proc. R. Soc. A 11 Buckling of inextensible strips (a) (b) – M = –25.1 – M = –15.1 – M = –2.52 – M = –0.24 – M = 0.075 – M = 1.07 – M = 2.08 – M = 3.08 – M = 4.07 – M = 5.08 – M = 6.08 – M = 9.23 – M = 12.1 – M = 15.1 – M = 22.6 – M = 30.0 – M = 39.0 – M = 54.9 20 – F 25 15 20 10 15 5 M 10 0 5 −5 0 −10 F = –4.06 – F = –3.67 – F = –3.26 – F = –1.64 – F = 0.01 – F = 1.67 – F = 3.32 – F = 4.98 – F = 6.63 – F = 8.15 – F = 9.73 – F = 11.4 – 0 2 4 6 8 −5 0 10 2 4 α |e| 6 8 Figure 5. n = 4 mode. (a) Force F¯ versus scaled end-to-end distance |e|. The curves of constant M̄ form a sequence on the plot, with M̄ increasing from the top left dashed curve (M̄ = −25.1) to the dotted curve, then from the dot-dashed curve to the top right dashed curve. (b) Moment M̄ versus twist angle a (in radians). The curves of constant F¯ form a sequence on the plot, with F¯ increasing from the dashed curve (F¯ = −4.06) to the dotted curve (F¯ = −3.26), then to the top dotted curve. (a) (b) 50 50 40 – M = –1.34 – M = –0.7596 – M = –0.56 – M = –0.513 – M = –0.2568 – M = 0.3955 – M = –20.816 – M = –10.13 30 20 – F 10 0 −10 −20 – F = –12.87 – F = –6.683 – F = –0.01 – F = 6.80 – F = 13.42 – F = 20.04 – F = 26.67 – F = 33.3 – F = 39.93 – F = 46.56 40 30 – M 20 10 0 −30 −10 −40 0 1 2 |e| 3 0 2 4 6 8 10 α Figure 6. n = 8 mode. (a) Force F¯ versus scaled end-to-end distance |e|. The curves of constant M̄ form a nested sequence on the plot, with M̄ increasing from the left dashed loop (M̄ = −1.34) to the dotted curve, then from the dot-dashed curve to the top dashed curve. (b) Moment M̄ versus twist angle a (in radians). The curves of constant F¯ form a nested sequence on the plot, with F¯ increasing from the dashed curve (F¯ = −12.87) to the dot-dashed curve. On the right of each figure is shown the moment response for the sequence of axial forces (F¯ ) shown on the legend. These curves were obtained by continuation, where F¯ was varied, keeping M̄ fixed or vice versa. Exceptions to this are some curves in figures 5a,b and 6a, where both F¯ and M̄ were allowed to vary, only Proc. R. Soc. A 12 A. P. Korte et al. so that a continuation could be started at an arbitrary point on the graph (at which point one of F¯ or M̄ would then be held fixed). This accounts for the curve-crossing in these figures. In figure 4a, 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 4b 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 5a, 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 5b 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 a. The three lowest negative applied force continuations terminate at a = 2p, again forming a closed strip, this time the T4,1 elastic torus ribbon knot. In figure 6a, 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 (M̄ = −20.816) terminates in the double cover of the T4,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 6b 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 a. 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 (M̄ = −35.13) terminates in the T8,1 elastic torus ribbon knot. Both ends of the dot-dashed curve (M̄ = 13.67) terminate in a double covering of the T4,1 elastic torus ribbon knot, whereas the T8,3 structure (shown at the minimum end-toend 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 T8,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 8b 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 M̄ . To obtain this plot, strips for different modes have been scaled to the same aspect ratio 2nL/2w = 10.5300. Parts of these curves with positive Proc. R. Soc. A 13 Buckling of inextensible strips – F 150 100 T8,1 T4,1 50 T8,3 0 −50 0 2 4 6 |e| Figure 7. n = 8 mode. Force F¯ versus scaled end-to-end distance, |e|, showing strip solutions. Solid line, M̄ = −9.19; dashed-dotted line, M̄ = 13.67; dashed line, M̄ = −35.13. (a) (b) Figure 8. Shapes of strips: (a) n = 2, figure-of-eight and (b) n = 8, F¯ = −19.975, M̄ = 0.1424. 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 Proc. R. Soc. A 14 A. P. Korte et al. 8 6 |e| 4 2 0 −5 0 5 10 15 20 – F Figure 9. Fold structure of scaled end-to-end distance versus applied force F¯ , n = 2, . . . , 8 at constant M̄ = 14.49. Grey dash-dotted line, n = 2; grey dashed line, n = 3; grey solid line, n = 4; black solid line, n = 5; black small dashed line, n = 6; black dash-dotted line, n = 7; black long dashed line, n = 8. 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 a) 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 t, or h, equal to zero from the start (in equation (2.2)) leading to the Euler elastica equation. In this case, no singularities will occur. Proc. R. Soc. A Buckling of inextensible strips 15 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. This work was supported by the UK’s Engineering and Physical Sciences Research Council (EPSRC) under grant no. EP/F023383/1. Appendix A. Derivation of the equations using standard variational techniques  We wish to calculate the variation of the functional H = f (k, t, k′ , t′ ) ds. 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. Proc. R. Soc. A 16 A. P. Korte et al. 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 H = f1 (k, t) ds, stating that it can 1   be written down from the variations of H2 = f2 (k) ds and H3 = f3 (t) ds, which are previously calculated. The resulting force and moment corresponding to H1 can then be written down from the expressions for the force and moment for H2 and H3 . This follows from the  chain rule. Similarly, one can calculate the variation of the functional H = f (k, t, k′ , t′ ) ds using dH1 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 H1 , 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 dH expressed in terms of dk and dt only. This method is easily generalized to functionals involving higher derivatives of k and t. Following the notation of Capovilla et al. (2002), the infinitesimal deformation of a space curve parts is dr = J|| t + J1 n + J2 b, and denoting d|| , d⊥ as the tangential and normal  of the deformation, respectively, the variation of any functional H = f ds is  dH = d0 H + d⊥ f ds, (A 1)   where d0 H ≡ d|| H + f d⊥ ds = ((f J|| )′ − f kJ1 ) ds. We therefore need to calculate the second term on the right-hand side of equation (A 1), where by the chain rule, d⊥ f = fk d⊥ k + ft d⊥ t + fk′ d⊥ k′ + ft′ d⊥ t′ , using the notation fb ≡ vf /vb. Since for any scalar h, d⊥ (h ′ ) = kh ′ J1 + (d⊥ h)′ , then d⊥ f = (fk − fk′′ )d⊥ k + (ft − ft′′ )d⊥ t + k(k′ fk′ + t′ ft′ )J1 + (fk′ d⊥ k + ft′ d⊥ t)′ , (A 2) where d⊥ k and d⊥ t are given in Capovilla et al. (2002). Equation (A 1) along with equation (A 2) is of the form   (A 3) dH = ds Ei Ji + ds Q′ , with Q the Noether charge, so that the Euler–Lagrange equations Ei = 0 can be immediately written down, since these are just the coefficients of Ji , which can be read off from equation (A 1). These are two coupled differential equations in the unknowns k and t, 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 dr = e, where F is defined by Q = −e · F. Similarly, the conserved moment T is obtained by specializing the deformation to a constant infinitesimal rotation dr = U × r, where T is defined by Q = −U · T 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 ad⊥ b, where b = k, t, has to Proc. R. Soc. A 17 Buckling of inextensible strips 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 (ad⊥ b)′ , one just needs the contribution from d⊥ b. Thus, for an infinitesimal rotation, dr′ = U × r′ = U × t and dr′′ = kU × n. Using the compact expressions of Langer & Perline (1991, eqn (3.7b,c)), dk = dr′′ · n − 2kdr′ · t gives dk = 0, i.e. d⊥ k = −d|| k = −J|| k′ . But J|| ≡ dr · t = (r × t) · U, giving d⊥ k = −U · (k′ r × t). As this is of the form −U · (r × F), Q ≡ ad⊥ k in equation (A 2) does not contribute to the moment M. Similarly, dt = dr′′ · b/k′ + dr′ · (kb − tt) = 0, giving d⊥ t = −U · (r × t′ t). As this is of the form −U · (r × F), Q ≡ ad⊥ t 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)), ad⊥ k contributes −ab to the moment, whereas ad⊥ t contributes −at − a ′ n/k. Thus, from equations (A 1) and (A 2), one can immediately write down the moment as (f ′′′ − ft′ )n (A 4) + (fk′′ − fk )b. M = (ft′′ − ft )t + t k In order to calculate the contribution to the force of Q ≡ fk′ d⊥ k + ft′ d⊥ t, one notes that since dr = e is a constant, it follows from the formulation of Langer & Perline (1991, eqn (3.7b,c)) that dk = 0, i.e. d⊥ k = −d|| k = −J|| k′ = −e · (k′ t). Thus, the contribution of Q ≡ fk′ d⊥ k in (A 2) to F is k′ fk′ t. In a similar manner, for this constant translation, dt = 0, giving d⊥ t = −e · (t′ t), so that Q ≡ ft′ d⊥ t in equation (A 2) contributes t′ ft′ t to the force. Now from Capovilla et al. (2002) (eqns (48) and (64)), ad⊥ k contributes kat + a ′ n + tab to the force, whereas ad⊥ t contributes tat + ta ′ n/k − ((a ′ /k)′ + ka)b. Also, d|| H contributes −f t to the force. Thus, from equations (A 1) and (A 2), one can immediately write down the force as F = (−f + k′ fk′ + t′ ft′ + k(fk − fk′′ ) + t(ft − ft′′ ))t t + fk′ − fk′′′ + (ft′ − ft′′′ ) n k + t(fk − fk′′ ) − k(ft − ft′′ ) −  (ft′ − ft′′′ ) k ′  b. (A 5) From these components of F and M, the differential equations for the moment components and the differential equation for Ft in equations (2.7) are easily obtained. For instance, from equation (A 4), it is seen that Mt′ = kMn , while equations for Mn , Mb and Ft can be extracted similarly. The remaining two equations for Fn and Fb are obtained from equation (A 5) and the Euler–Lagrange equations Ei = 0.  In summary, the variation of H = f (k, t, k′ , t′ ) ds can be easily obtained from the variations dk and dt. 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 k′ and t′ are therefore easily managed. Proc. R. Soc. A 18 A. P. Korte et al. One way of regarding equations (A 4) and (A 5) is that they prescribe F and M once k and t for the given density f are known. The two Euler–Lagrange equations for k and t are given by Ei = 0, i.e. by setting the coefficients of Ji 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 M , and, for the functional (2.6), to change variables from (k, t) to (k, h). By considering functionals of the form g(k, h, h′ ) = f (k, t, k′ , t′ ), one can show that (vg/vk)h,h′ = fk + hft + h′ ft′ , (vg/vh)k,h′ = kft + k′ ft′ and (vg/vh′ )k,h = kft′ . 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 = (Ft , Fn , Fb )T , M = (Mt , Mn , Mb )T , one can show, using the tangent and binormal components of M from equation (A 4), that  vk g + hMt + Mb = 0 (A 6) and (vh′ g)′ − vh g − kMt = 0. The two scalar equations (A 6), along with the six balance equations (equations (2.7)) for the components of the internal force F and moment M (see Starostin & van der Heijden 2007, 2009), are then the equations satisfied by the centreline of the developable strip for the unknowns F, M , k, h, F′ + u × F = 0 vk g + hMt + Mb = 0 and and M ′ + u × M + t × F = 0, ′ (vh′ g) − vh g − kMt = 0, (A 7) (A 8) where u = k(h, 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). References Anderson, I. M. 1989 The variational bicomplex. Technical report, Utah State University, Logan. See http://www.math.usu.edu/∼fg_mp/Publications/VB/vb.pdf. Ashwell D. G. 1962 The inextensional twisting of a rectangular plate. Q. J. Mech. Appl. Math. 15, 91–107. (doi:10.1093/qjmam%2F15.1.91) Capovilla, R., Chryssomalakos, C. & Guven, J. 2002 Hamiltonians for curves. J. Phys. A: Math. Gen. 35, 6571–6587. (doi:10.1088/0305-4470/35/31/304) Clebsch, A. 1862 Theorie der Elasticität fester Körper. Leipzig, Germany: Teubner. Coman, C. D. & Bassom, A. P. 2008 An asymptotic description of the elastic instability of twisted thin elastic plates. Acta Mech. 200, 59–68. (doi:10.1007/s00707-007-0572-3) Crispino, D. J. & Benson, R. C. 1986 Stability of twisted orthotropic plates. Int. J. Mech. Sci. 28, 371–379. (doi:10.1016/0020-7403(86)90056-1) Doedel, E. J., Champneys, A. R., Fairgrieve, T. R., Kuznetsov, Y. A., Sandstede, B. & Wang, X. J. 2007 AUTO-07p: continuation and bifurcation software for ordinary differential equations. See http://indy.cs.concordia.ca/auto/. Goldstein, H. 1980 Classical mechanics. Reading, MA: Addison-Wesley Publishing Company, Inc. Graustein, W. C. 1966 Differential geometry. New York, NY: Dover. Green, A. E. 1936 The equilibrium and elastic stability of a thin twisted strip. Proc. R. Soc. Lond. A 154, 430–455. (doi:10.1098/rspa.1936.0061) Green, A. E. 1937 The elastic stability of a thin twisted strip. II. Proc. R. Soc. Lond. A 161, 197–220. (doi:10.1098/rspa.1937.0141) Proc. R. Soc. A Buckling of inextensible strips 19 Hangan, T. 2005 Elastic strips and differential geometry. Rend. Sem. Mat. Univ. Pol. Torino 63, 179–186. Hornung, P. 2009 Minimizers of Kirchhoff’s plate functional: Euler–Lagrange equations and regularity. C. R. Acad. Sci. Paris 347, 647–650. (doi:10.1016/j.crma.2009.03.031) Kogan, I. A. & Olver, P. J. 2003 Invariant Euler–Lagrange equations and the invariant variational bicomplex. Acta Applicandae Mathematicae 76, 137–193. (doi:10.1023/A:1022993616247) Langer, J. & Perline, R. 1991 Poisson geometry of the filament equation. J. Nonlinear Sci. 285, 131–144. (doi:10.1007/BF01209148) Love, A. E. H. 1927 A treatise on the mathematical theory of elasticity. Cambridge, UK: Cambridge University Press. Mansfield, E. H. 1989 The bending and stretching of plates, 2nd edn. Cambridge, UK: Cambridge University Press. Mockensturm, E. M. 2001 The elastic stability of twisted plates. J. Appl. Mech. Trans. ASME 68, 561–567. (doi:10.1115/1.1357517) Randrup, T. & Røgen, P. 1996 Sides of the Möbius strip. Arch. Math., 66, 511–521. Sadowsky, M. 1931 Theorie der elastisch biegsamen undehnbaren Bänder mit Anwendungen auf das Möbius’sche Band. Verhandl. des 3. Intern. Kongr. f. Techn. Mechanik, 1930, Teil II, 444–451. Starostin, E. L. & van der Heijden, G. H. M. 2007 The shape of a Möbius strip. Nat. Mater. 6, 563–567. (doi:10.1038/nmat1929) Starostin, E. L. & van der Heijden, G. H. M. 2009 Force and moment balance equations for geometric variational problems on curves. Phys. Rev. E 79, 066602. (doi:10.1103/PhysRevE. 79.066602) Thamwattana, N., McCoy, J. A. & Hill, J. M. 2008 Energy density functions for protein structures. Q. J. Mech. Appl. Math. 61, 431–451. (doi:10.1093/qjmam%2Fhbn012) van der Heijden, G. H. M. & Thompson, J. M. T. 2000 Helical and localised buckling in twisted rods: a unified analysis of the symmetric case. Nonlinear Dynm. 21, 71–99. (doi:10.1023/A:1008310425967) Wakamatsu, H., Yamasaki, T., Arai, E. & Hirai S. 2007 Modeling of flexible belt objects toward their manipulation. In Proc. 2007 IEEE Int. Symp. Assembly and Manufacturing, Ann Arbor, MI, 22–25 July 2007, pp. 1–6. (doi:10.1109/ISAM.2007.4288440) Witten, T. A. 2007 Stress focusing in elastic sheets. Rev. Mod. Phys. 79, 643–675. (doi:10.1103/ RevModPhys.79.643) Wunderlich, W. 1962 Über ein abwickelbares Möbiusband. Monatsh. Math. 66, 276–289. (doi:10.1007/BF01299052) Yoshimura, Y. 1951 On the mechanism of buckling of a circular cylindrical shell under axial compression. Reports of the Institute of Science and Technology of the University of Tokyo, vol. 5(5), 179–198. [English translation: Technical Memorandum 1390 of the National Advisory Committee for Aeronautics, Washington DC, 1955.] Proc. R. Soc. A