# Instability and self-contact phenomena in the writhing of clamped rods

G.H.M. van der Heijden, Sébastien Neukirch, V. G. A. Goss, J. M. T. Thompson

► To cite this version:

G.H.M. van der Heijden, Sébastien Neukirch, V. G. A. Goss, J. M. T. Thompson. Instability and self-contact phenomena in the writhing of clamped rods. International Journal of Mechanical Sciences, 2003, 45 (1), pp.161-196. 10.1016/S0020-7403(02)00183-2. hal-05384962

HAL Id: hal-05384962

https://hal.science/hal-05384962v1

Submitted on 27 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Instability and self-contact phenomena in the writhing of clamped rods

G.H.M. van der Heijden, S. Neukirch, V.G.A. Goss & J.M.T. Thompson

Centre for Nonlinear Dynamics
University College London
Gower Street, London WC1E 6BT, UK

(To appear in the *Int. J. Mech. Sciences*)

## Abstract

We use the Cosserat rod theory to present a unified picture of jump phenomena, associated with looping, snap-through, pop-out, etc., in twisted clamped rods undergoing large deflections. Both contact-free rods and rods with isolated points of self-contact are considered. Taking proper account of the symmetries of the problem we find that an arbitrary contact-free solution is fully characterised by four parameters; each point contact adds another two. A shooting method is used for solving the boundary value problem. An intricate bifurcation picture emerges with a strong interplay between planar and spatial rod configurations. We find new jump phenomena by treating the ratio of torsional to bending stiffness of the rod as a bifurcation parameter. Load-deflection curves are computed and compared with results from carefully conducted experiments on contact-free as well as self-contacting metal-alloy rods.

**Keywords:** Cosserat rod, clamped boundary conditions, end rotation, (out-of-plane) buckling, instability, bifurcation, self-contact, looping, snap-through, pop-out, experiments

## 1 Introduction

We study the global behaviour of clamped rods, that is, rods whose ends are fixed into parallel end-fittings that can be moved inwards towards each other and rotated (see Fig. 1). We especially focus on dynamic jump phenomena that may occur when the end-to-end displacement or the end-to-end rotation is varied. Examples of such jump phenomena are looping, snap-through, pop-out, etc., where by ‘looping’ we mean a jump of a twisted rod into self-contact, ‘snap-through’ is used for a jump not leading to self-contact, and ‘pop-out’ is used for a jump out of self-contact. All these phenomena, and their associated hysteresis cycles, are readily observed in simple hand-held experiments, and we encourage the interested reader to keep a piece of rod within reach (a silicone rubber specimen will do nicely). More sophisticated laboratory tests with metal-alloy rods (see Fig. 1) have been found to give results, to be presented later, which are also in reasonable quantitative agreement with the theory. Mathematically, most of these jump phenomena can be characterised as fold or pitchfork bifurcations.
<figure><img src="image_2.png" /><figcaption>Figure 1: The testing rig showing the force transducer on the right and the torque transducer on the left of the clamps holding the rod.</figcaption></figure>

Clamped end conditions tend to impose stronger constraints on the rod's configuration in space and to produce more stable solutions and a richer bifurcation set than for instance found in a rod held in a sleeve (so the end is free to rotate about its axis, and no twisting moment can be transmitted), a rod supported by a Cardan (or Hookean, or universal) joint (so the end is free to rotate about two orthogonal axes, one fixed in the rod and one fixed in space), or in a rod supported by a ball-and-socket joint (so the end is free to rotate about three independent axes, and the total moment is zero at the end). Clamping is also the natural way of supporting a rod in many engineering applications.

Jump phenomena in end-loaded rods have been studied in scattered places in the literature. On the one hand there are studies of looping and pop-out of twisted cables (e.g., [55], [52], [9]),

<table><caption>Table 1: List of frequently used symbols.</caption><tr><td><i>L</i>, <i>ρ</i></td><td>length and dimensionless radius of the rod</td></tr><tr><td><i>B</i>, <i>C</i>, <i>γ</i></td><td>bending and torsional stiffness, and stiffness ratio <i>C</i>/<i>B</i></td></tr><tr><td><b><i>F</i></b>, <i>F</i>, <i>f</i></td><td>(dimensionless) internal force</td></tr><tr><td><i>T</i>, <i>t</i></td><td>(dimensionless) end tension</td></tr><tr><td>Δ<i>F</i></td><td>magnitude of point contact force</td></tr><tr><td><b><i>M</i></b>, <b><i>m</i></b></td><td>(dimensionless) internal moment</td></tr><tr><td><b><i>M</i></b><sub>0</sub>, <b><i>m</i></b><sub>0</sub>, <b><i>m</i></b><sub>0x</sub>, <b><i>m</i></b><sub>z</sub></td><td>(dimensionless) reference moments at the midpoint of the rod</td></tr><tr><td><i>M</i></td><td>end twisting moment</td></tr><tr><td><i>D</i>, <i>d</i></td><td>(dimensionless) end shortening</td></tr><tr><td><i>R</i></td><td>end rotation</td></tr></table>
or compressed rods (e.g., [25], [33], [41]). On the other hand there are more mathematical works concerned mainly with the existence of points of instability (e.g., [1], [29]). In most of these studies, however, only one particular phenomenon is considered ([33] and [29] are more inclusive). Here we present a more complete picture of how pre- and post-jump solutions are connected in the solution space. Along the way we derive several exact analytical expressions for critical buckling loads not seen previously in the literature. New jump phenomena are found by exploring a range of values for the ratio of torsional to bending stiffness of the rod. We treat both free rods and rods with one or multiple points of self-contact for which we model the rod as an impenetrable elastic tube of finite radius.

Our starting point is formed by the Kirchhoff equations for a perfect, inextensible, isotropic rod undergoing arbitrary deformations under the action of end loads. A singularity-free (in fact, polynomial) ODE is obtained by avoiding the usual Euler angles and working with the position and tangent vectors as variables instead. Although the equations for a perfect, isotropic rod are known to be solvable in terms of elliptic integrals and functions, this is not particularly helpful if one is interested in numerical results, especially when considering rigid, i.e., displacement-controlled, loading, as the implicit transcendental equations have to be solved numerically. Also, working with closed-form solutions leads to annoying discontinuities when parameters are varied. Therefore we use numerical solution methods throughout. The results give a coherent and unified picture of the buckling and post-buckling of relatively short pieces of rod and allow full loading sequences to be followed and hysteresis loops to be identified. Since most of the instabilities encountered pose a threat to the proper operation of pipes and cables in industrial environments, this approach may be useful in designing operating conditions with minimum risk of looping or snap-through, or conditions which favour pop-out in order to avoid kinking (i.e., permanent plastic deformation).

We mostly study open rods as these are more conventionally found in engineering applications. Nevertheless, the planar ring is found to be an important organising centre for deformations of both the planar and the spatial elastica. The ring is an example of a closed rod and for such rods we can define the so-called link as the topological linking number of two lines drawn on opposite sides of the unstressed rod [14]. This link is a topological invariant, i.e., it is unaffected by continuous deformations. This and the fact that the (clamped!) boundary conditions for a closed rod are naturally prescribed make closed rods an attractive object for studies from a more mathematical viewpoint. Indeed, their symmetry and bifurcation properties have been studied extensively by Maddocks and co-workers (e.g., [31]). Closed rods with self-contact have more recently been studied by Coleman and co-workers [46, 8] as models for DNA plasmids.

Clamped rods have several attractive features. Firstly, as we shall show, clamped boundary conditions applied to the (reversible) equilibrium equations imply that all open-rod solutions have reflection (anti-)symmetry about the midpoint. So one need only compute half solutions; full solutions follow by suitable reflection. Secondly, because in clamped rods the end tangents are aligned, the natural end displacements most readily accessible in an experiment, namely the end shortening (or slack) and the end rotation, are also the ones through which the end tension and twisting moment do work. Thus, by using end shortening and end rotation in load-deflection diagrams one can infer *static* stability properties of solutions by making use of the so-called fold rule [43, 30]. The end rotation, in a sense to be made precise, can be seen as a generalisation of the link, well-defined for closed rods, to arbitrary clamped configurations. Although an important quantity, it has not been considered in the literature so much ([33] is a notable exception).

For self-contacting configurations we find a bifurcation sequence in which the rod jumps to states of increasing number of self-contacts as a control parameter (end rotation or end shortening).
is increased, resulting in a ply structure known, for instance, from DNA supercoiling [8, 44] and textile yarn twisting [13]. This process of ply formation, although remarkably delicate in its theoretical detail, can easily be observed in a rubber rod. The beginning of ply formation can be seen in Fig. 12.

The organisation of the paper is as follows. In Section 2 we formulate the problem of a clamped contact-free rod. We introduce the end rotation as a measure of torsional deflection and discuss its relationship to link and writhe defined for closed rods. Taking into account the symmetries of the problem we show that for a given rod a clamped solution is fully characterised by four parameters. In Section 3 we present our numerical results focussing on the various jump phenomena found. We also use the fold rule to ascertain the (in)stability of solution branches. Section 4 sets up the equations for a rod with an arbitrary number of isolated point contacts. For a symmetric solution each point contact requires an additional two parameters to be specified. Numerical results, in the form of load-deflection curves, for self-contacting rods are given in Section 5. Section 6 presents some results from our experiments, and Section 7 closes this study with a discussion. In an appendix we spell out a proof that deformations of clamped rods are symmetric.

# 2 Formulation of a clamped rod

## 2.1 Equilibrium equations

Consider a thin elastic rod of length $L$ held by end forces and moments. If we denote the position of the rod's centreline by the vector function $R$, then the internal force, $F$, and moment, $M$, in an element of the rod are governed by the balance equations [2]:

$$F' = 0, \qquad (1)$$

$$M' + R' \times F = 0, \qquad (2)$$

where $(\cdot)' = d/dS$, $S$ being the arclength parameter measured along the centreline of the rod. We shall take the origin of arclength at the middle of the rod, so that $-L/2 \le S \le L/2$. We assume the rod to be perfect (i.e., intrinsically straight and untwisted), inextensible and unshearable. Thus the only distortions which the rod suffers are on account of bending and twisting moments, and not of shear forces or tension. Let $\{d_1, d_2, d_3\}$ be a right-handed rod-centred orthonormal co-ordinate frame (the *directors*) with $d_3$ the local tangent to the rod and $d_1$ and $d_2$ two vectors in the normal cross-section coincident with the principal bending axes. Thus we have

$$R' = d_3, \qquad (3)$$

while the evolution of the co-ordinate frame along the rod is governed by the equation

$$d_i' = u \times d_i \quad (i = 1, 2, 3). \qquad (4)$$

Here $u$ is the strain vector whose components in the moving frame are the curvatures and the twist. We introduce linear constitutive relations between the moments and the strains:

$$M \cdot d_1 = B_1 u_1, \quad M \cdot d_2 = B_2 u_2, \quad M \cdot d_3 = C u_3. \qquad (5)$$

$B_1$ and $B_2$ are the principal bending stiffnesses about $d_1$ and $d_2$, respectively, while $C$ is the torsional stiffness about $d_3$.
Clearly, **F** is a constant vector in space. Assuming this vector to be non-zero we can define a fixed orthonormal co-ordinate frame {$e_x, e_y, e_z$} with $e_z$ taken in the direction of the internal force:

$$ \mathbf{F} = F \mathbf{e}_z, \quad F > 0, \qquad (6) $$

and $e_x$ and $e_y$ such that at $S = 0$ the tangent to the rod lies in the ($e_x, e_z$) plane (i.e., $d_3(0) \cdot e_y = 0$). This can be done without loss of generality. (The force **F** will only be zero at isolated points along solution curves to be presented in the following sections. These points correspond to straight rods or twistless rings, and since these can be dealt with separately (see Table 2), the above definition gives no problems.) Eq. (2) can now also be integrated to give

$$ \mathbf{M} = F \mathbf{e}_z \times \mathbf{R} + \mathbf{M}_0, \qquad (7) $$

where $\mathbf{M}_0$ is an integration constant which has to be determined through the boundary conditions as part of the problem. It easily follows that the equilibrium equations have the two first integrals (constants or conserved quantities, not varying with $S$)

$$ \mathbf{M} \cdot \mathbf{F} = F \mathbf{M} \cdot \mathbf{e}_z \quad \text{and} \quad \frac{1}{2} \mathbf{M} \cdot \mathbf{u} + \mathbf{F} \cdot \mathbf{d}_3 = \frac{1}{2} (B_1 u_1^2 + B_2 u_2^2 + C u_3^2) + F \mathbf{d}_3 \cdot \mathbf{e}_z. \qquad (8) $$

The first of these, after division by $F$, is the moment about the internal force vector ('the wrench moment'), while the second, interpreted properly, is the Hamiltonian [23]. The fact that the resultant of the vector force and vector moment at any cross-section is statically equivalent to a wrench fixed in space is an important concept. As we vary the loading of the rod, and move along an equilibrium path, the orientation of the wrench axis will normally vary, typically becoming displaced and rotated relative to the fixed axis of the loading rig. We shall return to this in Section 5.2.

In this paper we shall only be concerned with an isotropic rod, i.e., one for which $B_1 = B_2 = B$, say. Then the twisting moment $\mathbf{M} \cdot \mathbf{d}_3$ is a further first integral and we have $B \mathbf{d}_3' = \mathbf{M} \times \mathbf{d}_3$. Therefore

$$ B \mathbf{d}_3' = (F \mathbf{e}_z \times \mathbf{R} + \mathbf{M}_0) \times \mathbf{d}_3. \qquad (9) $$

As we shall be interested in rigid loading sequences in which the end rotation of the rod is controlled we must follow the evolution along the rod of a material vector such as $\mathbf{d}_1$. Noting that

$$ \mathbf{d}_3 \times \mathbf{d}_3' = \mathbf{u} - u_3 \mathbf{d}_3, \qquad (10) $$

we can use (4) to obtain

$$ B \mathbf{d}_1' = (F \mathbf{e}_z \times \mathbf{R} + \mathbf{M}_0) \times \mathbf{d}_1 + (B - C) u_3 \mathbf{d}_3 \times \mathbf{d}_1, \qquad (11) $$

where $u_3$ is the (constant) twist.

Thus we have arrived at a system of three equations (3), (9) and (11) for the vector functions $\mathbf{R}$, $\mathbf{d}_3$ and $\mathbf{d}_1$, with parameters $F$, $\mathbf{M}_0$ and $u_3$. The equations satisfy the constraints

$$ \mathbf{d}_3 \cdot \mathbf{d}_3 = 1, \quad \mathbf{d}_1 \cdot \mathbf{d}_1 = 1, \quad \mathbf{d}_1 \cdot \mathbf{d}_3 = 0. \qquad (12) $$

Since the left-hand sides in (12) are first integrals of the equations, these constraints can be ensured by suitable boundary conditions.
## 2.2 Symmetries and boundary conditions

For reasons explained in the Introduction we consider clamped boundary conditions, which can be formulated as

$$ \mathbf{d}_3(-L/2) = \mathbf{d}_3(L/2), \qquad (13) $$

$$ \mathbf{R}(L/2) - \mathbf{R}(-L/2) = \lambda \mathbf{d}_3(L/2), \quad \text{for some } \lambda \in \mathbb{R}. \qquad (14) $$

Now the rod equations are well-known to have (reversing) symmetry properties [48]. These symmetries imply that solutions are either themselves reversible or come in reversible pairs. In our case of clamped rods one can prove that solutions selected by the boundary conditions are in fact reversible. Since the proof is most easily given using Euler angles we defer it to the Appendix. A consequence of the symmetry is that $\mathbf{M} \cdot \mathbf{e}_y$ is an odd function of arclength, so provided we take the initial conditions

$$ \mathbf{R}(0) = (X(0), Y(0), Z(0)) = (0, 0, 0), \qquad (15) $$

so that by (7) we have $\mathbf{M}_0 = \mathbf{M}(0) = (M_{0x}, 0, M_{0z})$, the reversing symmetry manifests itself as invariance of the equations (3), (9) and (11) under the following transformation:

$$ S \to -S, \qquad \begin{array}{l} X \to -X, \\ d_{3x} \to d_{3x}, \\ d_{1x} \to d_{1x}, \end{array} \qquad \begin{array}{l} Y \to Y, \\ d_{3y} \to -d_{3y}, \\ d_{1y} \to -d_{1y}, \end{array} \qquad \begin{array}{l} Z \to -Z, \\ d_{3z} \to d_{3z}, \\ d_{1z} \to d_{1z}, \end{array} \qquad (16) $$

where we have introduced components of vectors with respect to the fixed co-ordinate frame $\{\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z\}$. (Since $M_z(s)$ is constant, we shall henceforth refer to $M_{0z}$ simply as $M_z$.)

This invariance implies that we can restrict ourselves to the arclength interval $[0, L/2]$, and specify boundary conditions at $S = 0$ and $S = L/2$. Then (13) and (14) reduce to two clamp conditions,

$$ d_{3y}(L/2) = 0, \qquad X(L/2) d_{3z}(L/2) = Z(L/2) d_{3x}(L/2), \qquad (17) $$

and a formula for the span length $\lambda$,

$$ \lambda = 2 \{X(L/2) d_{3x}(L/2) + Z(L/2) d_{3z}(L/2)\}. \qquad (18) $$

Further, since $\mathbf{d}_3$ is a unit vector and $d_{3y}(0) = 0$ we can specify

$$ \mathbf{d}_3(0) = (d_{3x}(0), d_{3y}(0), d_{3z}(0)) = (\sin \theta, 0, \cos \theta), \qquad (19) $$

where $\theta \in [0, 2\pi)$ is the angle the rod at the midpoint makes with $\mathbf{e}_z$, and

$$ \mathbf{d}_1(0) = (\cos \theta, 0, -\sin \theta). \qquad (20) $$

Note that this choice satisfies the constraints (12).

Since $F$ (the ‘wrench force’) by definition is positive, the sign of the tangent vector $\mathbf{d}_3$ determines whether the rod is locally in tension or compression. In particular, the component of the end force in the direction of the rig axis (the ‘rig force’) is given by

$$ T = \mathbf{F} \cdot \mathbf{d}_3(L/2) = F d_{3z}(L/2), \qquad (21) $$

positive for tension, negative for compression. Meanwhile, the end twisting moment (or ‘rig moment’) is, from (5),

$$ M = \mathbf{M} \cdot \mathbf{d}_3(L/2) = C u_3. \qquad (22) $$
These $T$ and $M$ will be used to present our results since they are the quantities most readily measured in an experiment.

The end shortening, or slack, is defined as

$$D = L - \lambda = L - 2 \{X(L/2) d_{3x}(L/2) + Z(L/2) d_{3z}(L/2)\}. \quad (23)$$

It measures the difference in the end-to-end distance of the buckled and the straight rod. Note that $D$ takes on values between 0 and $2L$, and that $D = L$ for a closed rod. For $D > L$ the rod leaves the clamps on the outside (see, e.g., Fig. 2(a)). We shall call such solutions *inverted*. Notice that we are here ignoring the self-intersection at the clamps that would arise in a loading sequence in which $D$ passed through the value $L$.

Since the end tangents are aligned (by (13)) we can define the unique reduced end rotation $\hat{R} \in [0, 2\pi)$ by

$$\cos \hat{R} = \mathbf{d}_1(L/2) \cdot \mathbf{d}_1(-L/2), \quad \sin \hat{R} = \mathbf{d}_2(L/2) \cdot \mathbf{d}_2(-L/2), \quad (24)$$

where $\mathbf{d}_2 = \mathbf{d}_3 \times \mathbf{d}_1$. The full end rotation $R$ can be written as $R = \hat{R} + 2k\pi$, where $k$ is an integer that can be determined by using continuity. More precisely, we define the end rotation of a clamped rod in arbitrary deformation to be the end rotation that rod has ‘when completely pulled straight’ without (further) rotation of the ends (i.e., in the limit $D \to 0$, where simply $R = ML/C$). For this definition to be applicable any arbitrary clamped solution must be connected, by a continuous path of 3D solutions, to a unique (twisted) straight solution. Although we do not have a proof of this it is our firm conjecture, based on experience in our extensive numerical computations, that this is indeed the case, i.e., there appear to be no isolated (either open or closed) branches of solutions. This means that numerical continuation, in the right direction, can be employed to ‘trace back’ a given solution to $D = 0$ where $k$ (or $R$) can be read off immediately.

The (*numerical*) end rotation $R$ as defined above has the advantages that it is well defined for any arbitrary configuration and that it permits us to distinguish between different configurations that have the same reduced end rotation $\hat{R}$. However, it has to be realised that during parameter continuation in, for instance, $D$ (as in Fig. 3), the rod configuration may go through a self-intersection (which would not affect $R$). In consequence, $R$ does not necessarily agree with the *physical* end rotation (which a real rod would have that could not intersect itself). The problem with this physical end rotation, however, is that it is not straightforward to define generally.

One approach would be to ‘pull the rod straight’ by monotonically decreasing $D$ without (further) rotation of the ends. In this process, which mimics the actual straightening of an arbitrarily deformed physical rod, the rod would be allowed to jump dynamically to another solution branch (possibly one of self-contacting solutions, as computed in Section 5) when a limit point is encountered. This new branch may be disconnected from the initial branch (if the mode of the solution is different). The problem is that if multiple stable branches are available, then it is not a priori clear to which branch the rod will jump. This would require a dynamical theory. Even worse, if, along the loading path, the solution has gone through a self-intersection, then the rod may end up in a knotted configuration, in which case $D = 0$ cannot be reached (for a rod of finite thickness) and the approach to defining a physical end rotation fails utterly.

A more promising route would be to generalise the link, a topological quantity defined for a closed rod, to open rods. For a closed rod the link, $Lk$, is given by the linking number of the two curves $\mathbf{r} + \epsilon\mathbf{d}_1$ and $\mathbf{r} - \epsilon\mathbf{d}_1$ (which can be thought of as two curves drawn on opposite sides of the rod). For sufficiently small $\epsilon$ this is well defined and we have the celebrated formula [7, 14]

$$Lk = Tw + Wr. \quad (25)$$
Here $Tw = (1/2\pi) \int_L u_3 dS$ the total twist, and $Wr$ the writhe. The writhe is a property of the centreline of the rod and is a measure of the out-of-plane deformation (it is not a topological invariant). It can be computed as the double integral

$$Wr = \frac{1}{4\pi} \int_L \int_L \frac{[\mathbf{R}'(s) \times \mathbf{R}'(s')] \cdot [\mathbf{R}(s) - \mathbf{R}(s')]}{|\mathbf{R}(s) - \mathbf{R}(s')|^3} ds ds', \quad (26)$$

and is also equal to the signed crossing number averaged over planar projections from all possible directions [14]. An elementary exposition of link and writhe can be found in [51].

The definitions of link and writhe can be applied to open rods by imagining a closing piece of material connecting the ends of the deforming rod (such a closure was introduced earlier in [50]). This closure can be assumed to be twistless (by imagining it to be infinitely stiff in torsion) so that all the twist of the closed structure is in the rod. Then the total structure (rod plus closure) is closed and $Lk$ and $Wr$ are well-defined. They will however depend on the shape of the closure. For instance, if the rod is bent into a twistless circle then the closure will also form a closed curve and the total structure has a self-crossing. If both loops of this self-crossing structure are in the same plane then the writhe of the total structure will be +1 (for a suitable sign convention) and, $Tw$ being 0, by (25), $Lk = 1$. Since for a twistless circular rod $R = 2\pi$, we conclude that $Lk = R/(2\pi)$, i.e., the end rotation is essentially the link of the closed structure. For this relationship to hold generally, self-intersections in either the rod or the closure must be disallowed during any given loading sequence. We have already concluded that self-intersections do not affect the numerical end rotation $R$, but the writhe $Wr$ and the link $Lk$ will jump by $\pm 2$. (In case two (or any even number $n$ of) symmetrically placed self-intersections occur simultaneously, then the total change in $Wr$ and $Lk$ will be $\pm 2n$.)

This shows the significant problems with trying to generalise the link to open rods and to relate it to the physical end rotation, i.e., the measure of torsional deflection that is relevant for energy purposes. The matter will be pursued in detail elsewhere. In the remainder of this paper we shall only use the numerical end rotation $R$ as defined above. As long as no self-intersections occur during a loading sequence this $R$ agrees with the physical end rotation.

## 2.3 Solving the equations

Throughout this paper we present results in dimensionless form by using the length $L$ and the bending stiffness $B$ as characteristic scales. Occasionally we shall find it convenient to use the following abbreviations:

$$s = \frac{S}{L}, \quad \mathbf{r} = \frac{\mathbf{R}}{L}, \quad d = \frac{D}{L}, \quad f = \frac{FL^2}{4\pi^2 B}, \quad t = \frac{TL^2}{4\pi^2 B}, \quad \mathbf{m} = \frac{ML}{2\pi B}, \quad \gamma = \frac{C}{B}, \quad (27)$$

and to write $\mathbf{r} = (x, y, z)$ and $\mathbf{m}(0) = \mathbf{m}_0 = (m_{0x}, m_{0y}, m_{0z})$. We then have a system of nine equations (3), (9), (11) with nine initial conditions (15), (19), (20), a solution of which is fully specified by the constant $\gamma$ (fixed for a given rod) and the four dimensionless parameters

$$m_{0x}, m_z, f, \theta, \quad (28)$$

from which the dimensionless twist can be computed as

$$\frac{Lu_3}{2\pi} = \frac{1}{\gamma} (m_{0x} \sin \theta + m_{0z} \cos \theta). \quad (29)$$
<table><caption>Table 2: Summary of special solutions for a clamped rod (* means any value).</caption><thead><tr><th></th><th>m<sub>0x</sub></th><th>m<sub>z</sub></th><th>f</th><th>θ</th><th>type of solution</th></tr></thead><tbody><tr><td>1</td><td>0</td><td>0</td><td>0</td><td>*</td><td>straight and prismatic, unstressed</td></tr><tr><td>2</td><td>≠ 0<br>m<sub>0x</sub><sup>2</sup> + m<sub>z</sub><sup>2</sup> = n<sup>2</sup></td><td>0</td><td>tan θ = -m<sub>z</sub>/m<sub>0x</sub><br>(limit π/2)</td><td></td><td>untwisted n-covered planar ring<br>(M = const.)</td></tr><tr><td>3</td><td>0</td><td>0</td><td>≠ 0</td><td>0, π</td><td>straight and prismatic,<br>in tension (θ = 0) or compression (θ = π)</td></tr><tr><td>4</td><td>0</td><td>≠ 0</td><td>≠ 0</td><td>0, π</td><td>straight and twisted,<br>in tension (θ = 0) or compression (θ = π)</td></tr><tr><td>5</td><td>±nF</td><td>±n</td><td>≠ 0</td><td>π/2, 3π/2</td><td>twisted n-covered planar ring</td></tr><tr><td>6</td><td>≠ 0</td><td>0</td><td>≠ 0</td><td>0, π</td><td>odd-mode inflectional (θ = π) or<br>any-mode non-inflectional (θ = 0 or π)<br>elastica in y-z plane</td></tr><tr><td>7</td><td>0</td><td>0</td><td>≠ 0</td><td>*</td><td>even-mode inflectional elastica in x-z plane</td></tr><tr><td>8</td><td>≠ 0</td><td>≠ 0</td><td>≠ 0</td><td>tan θ = -m<sub>z</sub>/m<sub>0x</sub></td><td>untwisted 3D solution</td></tr><tr><td>9</td><td>≠ 0</td><td>0</td><td>≠ 0</td><td>*</td><td>3D solution with M ⋅ F = 0</td></tr></tbody></table>

The end displacements *d* and *R* are hereby determined through (23) and (24). The parameters (28) are found by applying the boundary conditions. In a rigid-load situation the end displacements *d* and *R* will be prescribed, and the two clamp conditions (17) together with (23) and one of the equations in (24) provide a well-posed problem of four equations for the four unknowns (28). In a typical loading sequence *d* (or *R*) is then varied and a curve of solutions is traced out.

It should be remarked that the equations (3) and (9), with the first integrals as mentioned in Section 2.1, are completely integrable and solvable by quadrature. The integrations are most conveniently carried out in a cylindrical co-ordinate system parallel to *e*<sub>*z*</sub> [23] and yield expressions in terms of elliptic integrals and functions (see, e.g., [38]). However, in the displacement-controlled loading case, as opposed to the dead (i.e., traction-controlled) loading case, *m*<sub>0</sub> and *f* are unknown and one ends up with a set of implicit transcendental equations which have to be solved numerically. (This is even more so in the contact problem considered in Section 4.) Therefore we shall use a numerical approach throughout (although some analytical expressions for planar solutions and critical parameter values will be given in Section 3).

Thus the final boundary value problem is solved over the interval [0, 1/2] using a shooting method in which the parameters (28) are varied. The other half of the solution is obtained by suitable reflection (cf. (16)). For the numerical integration we use the highly accurate 8th-order Runge-Kutta code DOP853 as described in [18], with local error tolerance 10<sup>-12</sup>. For single-parameter continuation of solutions we use the code AUTO [10] which uses orthogonal collocation allied to pseudo-arclength continuation. Equation (11) is slaved to (3) and (9), and is only integrated in case one requires (control of) the end rotation.

Special solutions, such as planar ones, are obtained for specific parameter values. Table 2 gives a summary. A more extensive classification of 2D as well as 3D clamped rod configurations will be given in [35]. As an example we now demonstrate that *m*<sub>*z*</sub> = 0 = *u*<sub>3</sub> implies a planar solution (the converse is not true). From (29) it follows that there are two cases: (i) θ = 0, π (case 6 in Table 2) and (ii) *m*<sub>0x</sub> = 0 (case 7 in Table 2). In case (i) *d*<sub>3x</sub>(0) = 0 by (19), and we have
the equations $d'_{3x} = m_y d_{3z}$, $m'_y = f d_{3x}$. Since $m_y(0) = 0$ by reversibility, it follows that all the derivatives of $d_{3x}$ are zero at $s = 0$, and hence that $d_{3x} \equiv 0$ (and thus $m_y \equiv 0$). Since $x(0) = 0$ this means that $x \equiv 0$ and the solution lies in the $y-z$ plane. If we introduce the angle $\phi$ such that $d_{3y} = \sin \phi$, $d_{3z} = \cos \phi$, then the remaining equations $d'_{3z} = m_x d_{3y}$, $m'_x = -f d_{3y}$ yield the familiar equations for the Euler elastica, $\phi'' - f \sin \phi = 0$, $m_x = -\phi'$. In case (ii) we have the equations $d'_{3y} = -m_x d_{3z}$, $m'_x = -f d_{3y}$, and since $d_{3y}(0) = 0$ we similarly find that $d_{3y} \equiv 0$ (and thus $m_x \equiv 0$). This, with $y(0) = 0$, gives $y \equiv 0$, so the solution is seen to lie in the $x-z$ plane.

### 3 Results of numerical continuation

Initially we take the value 5/7 for the stiffness ratio $C/B$. This is the value for the nitinol rods which we used in experiments (see Section 6). According to the theory of linear elasticity for rods of solid homogeneous circular cross-section, as the ones we used, this corresponds to a Poisson ratio $\nu = B/C - 1 = 0.4$. Later $C/B$ will be varied.

Where possible we use the fold rule to infer (static) stability properties of the solution branches computed. Standard bifurcation theory [22] predicts that folds in bifurcation diagrams mark exchanges of stability. From elastic stability theory it is well-known that in diagrams in which pairs of conjugate variables (one of which doing work through the other) are plotted, one can not only infer a change of stability but also the direction of change as a solution curve goes through a fold [43]. Examples of such distinguished diagrams in the present context are load-deflection diagrams of rig tension $T$ against end shortening $D$ and rig moment $M$ against end rotation $R$.

This fold rule is an instance of the application of the Morse index as defined for a certain class of boundary value problems having an underlying variational structure [30]. The index gives the maximum dimension of the space on which the second variation can be made negative. Provided some genericity conditions are satisfied the index changes by one each time a fold is traversed, and a solution can be stable only if it has index zero. Since only changes of index are recorded one needs stability information of one solution branch by other means in order to start the method off. Once that information is available (for instance, by physical means) the shape of solution branches in distinguished diagrams determines which solutions have non-zero index, and which therefore are unstable. Recent extensions of the use of an index (for instance to constrained problems and to cases where bifurcations occur) can be found in [36, 32, 20].

In the present case the total potential energy $V$ in the rod is given by [50]:

$$V = \frac{1}{2} \int (Bu_1^2 + Bu_2^2 + Cu_3^2) dS + TD - MR. \quad (30)$$

It then follows that in $D-T$ diagrams presented below the lower branch of solutions in a fold opening to the right, and the upper branch in a fold opening to the left, cannot represent local minima of the functional $V$ considered as a function of the variable $T$ and bifurcation parameter $D$. More generally, the index in the former case is decreased by 1, and in the latter case increased by 1, as one traverses the fold in upward direction (at constant $R$). Because of the different signs of the work terms in (30) the situation is the opposite in the $R-M$ diagrams presented below: in a fold opening to the right the upper branch is unstable, while in a fold opening to the left the lower branch is unstable. These rules are applied in Figs 3 and 7 below.
<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Labelled solution branches for the first four inflectional and non-inflectional modes of the planar elastica as given by (31) and (32). The inflectional ones bifurcate from the straight rod (vertical axis) and have $R = 0$, the non-inflectional ones all pass through the ring at $D/L = 1$ (marked by the solid circle) and have $R = 2n\pi$, where $n$ ($n \ge 1$) is the mode number. Insets show true-view, and fixed-length, shapes of the first and second mode. (b) Buckling curves for the first four modes of a clamped twisted rod as given by (33). Insets of the first and second mode of the planar elastica, at $M = 0$, illustrate the difference between odd and even modes, with the dotted lines indicating the wrench axes.</figcaption></figure>
## 3.1 Primary buckling

A clamped straight non-twisted rod buckles into a planar inflectional Euler elastica at the critical load $t = TL^2/(4\pi^2B) = -1$. The load-deflection curve for a general n-mode (i.e., a clamped planar elastica with $n + 1$ inflection points) for $n$ odd is given parametrically in terms of the elliptic modulus $k$ ($0 \le k \le 1$) by the expressions (see, e.g., [45])

$$ \frac{TL^2}{B} = -4(n + 1)^2 K^2(k), \qquad \frac{D}{L} = 2 \left( 1 - \frac{E(k)}{K(k)} \right), \qquad (31) $$

where $D$ is the end shortening and $K(k)$ and $E(k)$ are the Legendre complete elliptic integrals of the first and second kind, respectively. For $n$ even no such explicit parametrisation exists, but one can solve the equation for the planar elastica and apply the boundary conditions (13) and (14) to obtain a set of transcendental equations. For a uniform rod all buckling modes but the primary one ($n = 1$) are unstable, both in the plane and spatially, as first shown in [29].

For an n-mode (any $n$) non-inflectional clamped elastica the expressions are:

$$ \frac{TL^2}{B} = 4n^2k^2K^2(k), \qquad \begin{array}{l} \displaystyle \frac{D}{L} = 2 - \frac{2}{k^2} \left( 1 - \frac{E(k)}{K(k)} \right) \quad \left( \frac{D}{L} < 1 \right), \\ \displaystyle \frac{D}{L} = \frac{2}{k^2} \left( 1 - \frac{E(k)}{K(k)} \right) \quad \left( \frac{D}{L} > 1 \right). \end{array} \qquad (32) $$

Fig. 2(a) shows curves for the first four inflectional as well as non-inflectional modes. Note that the ‘figure-of-eight’ appears on both the first- and second-mode inflectional branch. On the former it is held at the edge and requires a compression, while on the latter it is held at the centre and requires a tension. (Despite the superficial similarity this figure-of-eight is not Bernoulli’s lemniscate.) In [11] a branch of (non-reversible) solutions is computed connecting the two figures-of-eight.

A non-trivial planar elastica (i.e., not a straight rod or a ring) always has $m_z = 0 = u_3$ (cf. Table 2). From the scaling in (27), and (31) and (32), it follows that if the parameters ($m_{0x}, t, \theta$) give a first-mode non-trivial planar elastica with end shortening $d$, then ($nm_{0x}, n^2t, \theta$) give an n-mode non-trivial planar elastica with end shortening $d$ for an odd-mode inflectional or an odd-mode non-inflectional elastica, and end shortening $2-d$ for an even-mode non-inflectional elastica. The scaling behaviour of the even-mode inflectional elastica, for which the rig and wrench axes make an angle (see the shapes in the insets of Fig. 2), is less straightforward.

If a clamped straight rod is subject to both compression and twist then it buckles at a critical load given by

$$ \cos \pi \sqrt{m_z^2 - 4t} - \cos \pi m_z = \left( 2\pi t \sin \pi \sqrt{m_z^2 - 4t} \right) / \sqrt{m_z^2 - 4t}, \qquad (33) $$

as follows from a linearisation of the equations (3) and (9) about the trivial solution. This condition defines a countable infinity of disjoint buckling curves in the $t-m_z$ load plane, the first four of which are drawn in Fig. 2(b). We are not aware of the appearance of (33) anywhere in classical works such as [45] and [56], although the latter deals with the special case of pure torsional buckling (i.e., $t = 0$), when (33) reduces to $\tan \pi m_z = \pi m_z$. Condition (33) is, however, derived in [54], where it is mentioned that the result agrees with results in [4] and [47]. (Greenhill [17], who gave the correct condition $m_z^2 = 1+4t$ for pinned ends, guessed the formula $m_z^2 = 4+4t$ for the first clamped buckling mode, but this is incorrect except for pure-thrust buckling for which $m_z = 0$. Incidentally,
pinned boundary conditions for a twisted rod in space are here taken to mean: $X(-L/2) = 0 = X(L/2), Y(-L/2) = 0 = Y(L/2), M_x(-L/2) = 0 = M_x(L/2), M_y(-L/2) = 0 = M_y(L/2).$

Although, as in the planar case, higher-order modes are unstable [20] and will not be observed near buckling, their solution curves may go through the requisite number of folds to produce stable solutions. For instance, Figs 3 and 6 show second-mode curves which give rise to stable branches going through the ring. Also, different modes interact with each other away from buckling, so they have to be considered if one is to get a global picture of the solution set. Indeed, it is not straightforward to define the ‘mode’ of a solution away from the linear regime near primary buckling. However, since a given solution can be traced back to a unique nearly straight solution near buckling we can define the mode of a solution by reference to this ‘linear’ solution. Although not always satisfactory (for instance, it makes the mode of the planar ring undefined), this definition is sufficient for our present purposes. A more careful definition is given in [35] by using the representation in terms of elliptic integrals and functions. It is then found that the mode changes when and only when $D/L$ goes through one (again rendering the mode of the planar ring ambiguous).

## 3.2 Out-of-plane buckling and beyond

Fig. 3 shows $D-T$ solution curves for several values of $R$. Two different sets of curves can be distinguished: those that have $R > 2\pi$, emerge from first-order modes at buckling (on the vertical axis), then go through a fold, and end with $T \to \infty$, and those that go through the planar ring (marked by the solid circle). The former will be further discussed in the next section; the latter, which show more complicated behaviour, form the topic of this section.

To start with the planar solutions, we observe that end shortening of the planar rod ($R = 0$) initially takes place under a rising compressive load. At a critical compression $-T$ and end shortening $D$, given by a critical modulus $k_{\text{out}}$ in (31), the planar rod becomes unstable and buckles out of the plane. The critical point, described by a pitchfork bifurcation, is indicated by the triangle in Fig. 3. At this point the axial compressive load starts to drop with increasing end shortening. The load becomes slightly tensile for a brief period before the shape of the rod reaches the circular ring with one turn of twist. Upon further increase of $D$ the configuration becomes inverted and slightly non-planar (it was shown in [49] that the only planar configurations of a twisted rod that can be supported by end loads only are the straight line and the circle).

The location of the out-of-plane bifurcation can be calculated analytically by linearising the 3D equilibrium equations about the planar inflectional elastica and solving the resulting equations subject to our clamped boundary conditions to derive a condition for the existence of neighbouring non-planar solutions. The calculations are most conveniently carried out using Euler angles and involve the evaluation of numerous integrals of combinations of Jacobian elliptic functions. We omit the details here and refer to Zajac [55] for a similar calculation for the non-inflectional elastica (a result we shall discuss shortly). The result can be presented in the following parametrised form:

$$ \frac{D}{L} = 2 \left( 1 - \frac{E(k)}{K(k)} \right), \quad \frac{C}{B} = \frac{K(k) [2(1 - k^2)K(k) + (-3 + 4k^2)E(k)]}{2(1 - k^2)K^2(k) + (-5 + 4k^2)K(k)E(k) + 3E^2(k)}. \quad (34) $$

The corresponding critical buckling load (for an arbitrary odd mode!) is obtained from (31).

Although this out-of-plane buckling instability has a sizable literature, (34) appears to be new, although Miyazaki & Kondo [33] give a highly implicit expression for the dependence of the critical load on the torsional stiffness. Previous studies of this secondary bifurcation further include: Kovári [25], who considers the case of an isotropic rod free to twist at one end, Antman
<figure><img src="image_4.png" /><figcaption>Figure 3: Load-deflection diagram with solution curves for $R = 0$ (a), $\pi$ (b), $2\pi$ (c), $3\pi$ (d), $4\pi$ (e) and $R = R_c = 8.9527$ (f). Solid lines represent stable branches, dashed lines unstable branches under rigid loading. Short-dashed curves are for the planar elastica, either looped and non-inflectional (upper curve), or inflectional (lower curve, retraced by the $R = 0$ curves, but emanating from the Euler compressive buckling load $-1$). The triangle indicates the out-of-plane instability for $R = 0$ at $D/L = 0.5590$, the diamond the into-the-plane instability for $R = 2\pi$ at $D/L = 0.6767$. The solid circle marks the ring. ($\gamma = 5/7$.)</figcaption></figure>

[1], who studies the linearisation of the nonlinear equations about a planar solution in order to obtain existence results, and Maddocks [29], who considers bifurcation under several types of boundary conditions in dead load using variational methods. The instability was studied in a perturbative approach allowing extension to anisotropic and imperfect rods in [41] in the context of birdcaging in wound cables.

The solid curve in Fig. 4 gives a graphical representation of (34). Notice that for a large enough ratio of torsional to bending stiffness the onset of lateral instability may occur for configurations beyond the figure-of-eight, i.e., for inverted solutions, a fact readily demonstrated with a vacuum cleaner tube.

Moving on to the non-planar solutions we first remark that for $R$ somewhere between $\pi$ and $2\pi$ the primary bifurcation from the trivial state (on the vertical axis) changes from being supercritical to being subcritical, with the bifurcating branch of solutions rising with increasing $D$. Although this doesn't mean much for a rod under rigid loading, for dead loading it means that the rod jumps, rather than gradually evolves, out of its straight configuration. This critical point was studied in [5] through a Lyapunov-Schmidt reduction for the slightly different case of a rod with prescribed parallel end tangents but with one end free to undergo side deflection.

The curves for $0 < R < 2\pi$ that emerge from first-order buckling modes are stable for all $D$. However, there are additional curves for the same $R$ values which emanate from (unstable) higher-order buckling modes, and undergo the right number and orientation of folds to have a
<figure><img src="image_5.png" /><figcaption>Figure 4: Loci of the out-of-plane (solid) and into-the-plane (dashed) bifurcations occurring along the $R = 0$ and $R = 2\pi$ curves, respectively. Insets show true views of the rod's shape at several of the points of lateral instability. The interval $2/3 < C/B < 1$ ($0 < \nu < 1/2$) delimited by the diamonds is the region of prime, though not sole, physical interest.</figcaption></figure>

stable branch of solutions through the ring. This stable branch shrinks as $R$ is increased and disappears at $R = R_c := 2\pi(\sqrt{3}/\gamma - 1)$ at which point the planar ring becomes unstable (see curve (f) in Fig. 3 with vertical tangent at the ring). None of the curves going through the ring having $R > R_c$ appear to contain stable branches.

The critical point at $R_c$ corresponds to the well-known twisted-ring instability criterion derived by Zajac [55]. More generally, it can be shown [27] that the nth buckling mode ($n \ge 1$) of the m-times covered ring ($m \ge 1$) bifurcates from the primary branch at a critical twisting moment given by

$$ \frac{ML}{2\pi B} = \sqrt{(n+m)^2 - m^2}. \qquad (35) $$

(Apart from this unpublished work [27] to our knowledge this result is only quoted, albeit incorrectly, in [15].) As shown in [36], this bifurcation is subcritical (with a jump into self-contact) for small $\gamma$ and supercritical (with a gradual deformation out of the plane) for large $\gamma$, the transition occurring at $\gamma = \gamma_c$, where

$$ \gamma_c = 1 + \frac{3m^2}{2(n+m)^2}. \qquad (36) $$

An early reference to stability studies of twisted elastic rings in the context of DNA is [26].

Notice that the curve for $R = R_c$ in Fig. 3 ends in a third-order buckling mode at $D = 0$. What is the highest mode with a stable branch through the ring depends on $\gamma$. Numerically we find this to be the 2-mode for $0.803 < \gamma < \sqrt{3}$ (consistent with Fig. 6) and the 3-mode for $0.11 < \gamma < 0.803$. An arbitrarily high mode may have a (small) stable branch through the ring by taking $\gamma$ sufficiently small. As follows from (35), for $\gamma > \sqrt{3}$ the ring will be unstable even for $R = 0$, i.e., with one full turn of twist put in, so no stable branches exist through the
<figure><img src="image_6.png" /><figcaption>Figure 5: Solution curves for $R = 0$ and varying $\gamma$. Shown are curves for $\gamma = 1.7$ (a), $1.8$ (b), $3.2$ (c) and $12.0$ (d). Short-dashed curves are for the planar elastica. All curves are going through the planar ring at $D/L = 1$ indicated by the solid circle. This ring, with one full turn of twist, is stable for $\gamma < \sqrt{3}$ and unstable for $\gamma > \sqrt{3}$. The diamond indicates the figure-of-eight for which $-TL^2/4\pi^2 B = 4K(k_8)^2/\pi^2 = 2.1833...$, where $k_8 = 0.9089...$ solves $K(k) = 2E(k)$, while the triangle marks the self-touching elastica at $k = 0.8551...$.</figcaption></figure>
ring. Incidentally, returning to the problems with defining a satisfactory end rotation discussed in Section 2.2, along the $R_c$ curve in Fig. 3 two simultaneous symmetric self-intersections of the rod axis occur somewhere between $D/L = 0.83$ and $D/L = 0.84$.

The behaviour for even larger values of $\gamma$ is depicted in Fig. 5 showing how the solution curve for $R = 0$ gets more and more drawn into a spiralling structure around the ring point. Solution curves for $R$ between 0 and $2\pi$ follow this spiralling behaviour but with fewer and fewer oscillations occurring as $2\pi$ is approached. Thus for these first-mode solutions increasing $R$ has a stabilising effect (the converse is true for second-mode solutions). In view of the fold rule for stability change, the spiralling does not create stable branches. Notice that the ring becomes unstable before the figure-of-eight gains stability but after the point where the planar elastica first touches itself (at $k = 0.8551...$ in (31)), so that the rod would have to intersect itself in order to follow the path in Fig. 5(b) with the dynamic jump across the unstable circle. Also notice that at $\gamma = 12.0$ (Fig. 5(d)) the out-of-plane pitchfork bifurcation has become subcritical. The transition from super- to subcritical occurs around $\gamma = 9$.

At first sight the critical value $\gamma = \sqrt{3}$, let alone the higher values in Fig. 5, might appear to be unrealistically high given that for solid circular rods of normal material the ratio of torsional to bending stiffness satisfies $2/3 < \gamma < 1$ (corresponding to $0 < \nu < 1/2$). However, Kehrbaum & Maddocks [24] have shown that a rod whose cross-section is anisotropic but which is rapidly twisting in the unstressed state can be approximated by an effectively isotropic rod in which the ratio of the effective torsional and bending stiffnesses can be much larger than 1. Also, direct experiments to estimate the bending and torsional stiffnesses of DNA molecules have led to 0.7–1.5 as the accepted range of values for the ratio $\gamma$ [21, 39, 19, 6].

### 3.3 Jump into self-contact (looping)

Under increasing $D$ the $R = 2\pi$ solution emerging from the first-order buckling mode in Fig. 3 becomes planar at the ‘into-the-plane’ critical point indicated by the diamond and given by a critical modulus $k_{in}$ in (32). At that point the solution has lost exactly one turn of twist and is twistless. The dependence of this critical point on the torsional stiffness of the rod is illustrated by the dashed curve in Fig. 4 which was obtained by Zajac [55] in parametrised form as

$$ \frac{D}{L} = 2 - \frac{2}{k^2} \left(1 - \frac{E(k)}{K(k)}\right), \quad \frac{C}{B} = \frac{K(k) \left[(4 - 3k^2)K(k) + (-4 + k^2)E(k)\right]}{(1 - k^2)K^2(k) + (2 - k^2)K(k)E(k) - 3E^2(k)} \quad (37) $$

from analysis of the 3D variational equations about the clamped planar looped elastica. The corresponding critical load (for an arbitrary $n$-mode!) is found by referring to (32). As $C/B$ tends to infinity $D/L$ tends to 0 (and $T$ to infinity), i.e., the planar elastica retains the loop for all end loads. As Zajac observes, the insensitivity of this pop-out point for the twistless planar elastica to $C/B$ means that one cannot greatly mitigate the formation of kinks by designing cables with favourable $C/B$ ratios. But this assumes the cable has zero thickness and no twist; we shall consider pop-out in twisted rods of finite thickness in Section 5.1. This Zajac pop-out of a zero-radius rod, associated with a loss of stability of the untwisted planar elastica under controlled decreasing $D$, can be seen as a super-critical pitchfork bifurcation.

The $R = 2\pi$ solution curve separates different post-buckling behaviours. Curves emerging from first-order modes with $R > 2\pi$ form limit points in $D$, where, under rigid loading, the rod dynamically jumps into a self-contacting looped solution (to be considered in Section 5). If
we could artificially follow through the entire (stable and unstable) post-buckling path for these
solutions, to zero *D* and infinite *T*, we would remove exactly one turn worth of twist.

Looping of rods with residual twist, i.e., those with more than one turn of end rotation put in, was studied by Coyne [9]. He took a 'semi-infinite' approach in which the rod is taken to be of finite length but of the shape of the homoclinic solution of the equations (i.e., the solution which is asymptotically straight towards both ends and strictly needs infinite length). Since an infinite rod always effectively experiences dead load he found for the looping instability the infinite-rod Greenhill [17] condition $m_z^2 = 4t$. See also [34] for a more systematic study of the semi-infinite approximation and how it relates to the finite-length case.

Besides the out-of-plane bifurcation (for $R=0$) and the into-the-plane bifurcation (for $R = 2\pi$) there is a third bifurcation associated with the planar elastica, again for $R=0$. It occurs along the second-mode branch of the inflectional planar elastica (see Fig. 2(a)), and defines the point where the $R=0$ curve connects to that planar elastica curve. (So the spatial elastica curve for $R=0$ is seen to connect the first-mode inflectional planar elastica to the second-mode inflectional planar elastica, both of which are independent of $\gamma$.) For the present value of $\gamma$ this second-mode bifurcation occurs too close to $D/L = 2$ to be visible along the (a) curves in Fig. 3. It can be observed in Fig. 6 for $\gamma = 1.3$ and 1.7. This value 1.7 is larger than $\gamma_c$ (for $n=m=1$) and generates different behaviour of solution curves through the ring. For instance, additional limit points develop along these curves, leading, for $R=0.3$, to two stable branches of solutions.

## 3.4 Snap-through

In Section 3.2 we saw a pitchfork-type secondary bifurcation where the rod buckles out of the plane. Put into energy terms: it is favourable for the rod to release some bending energy at the cost of some torsional energy. It is well-known from bifurcation theory that if a second parameter in the system is varied, a pitchfork bifurcation gets unfolded (generically) into solution curves with a simple fold (see the curves labelled (b) in Fig. 3). In structural engineering this is known from systems in which imperfection parameters are present. In our rod the end rotation plays this role as secondary bifurcation parameter, and the fold gives rise to snap-through behaviour in which the rod dynamically jumps from one configuration to another, symmetrically related, configuration. Associated with this are possibilities of hysteresis loops in which the rod jumps back and forth between two paths.

The phenomenon is illustrated in the *M*-*R* load-deflection diagram of Fig. 7. All curves are running from (*R*, *M*) = (-2π, 0) to (*R*, *M*) = (2π, 0). These end points have to be identified and the curves imagined closed: at (*R*, *M*) = (2π, 0) the rod is a planar non-inflectional elastica with no twist and with a self-intersection (for *D*/*L* > 1 the hypothetical closing piece connecting the ends of the rod intersects itself), and when going through this self-intersection *R* jumps from 2π to -2π (and the writhe of the hypothetical closed centreline from 1 to -1). At the origin the rod is a first-mode inflectional planar elastica (for instance, a figure-of-eight for *D*/*L* = 1).

Snap-through is possible for $d_{\text{out}} = 0.5590 < D/L < 0.936$ and $D/L > 1.109$ when there is a branch of stable solutions available as one jumps off the fold. The limits 0.936 and 1.109 are determined numerically and correspond to values of $D/L$ for which the fold occurs precisely at $R = 2\pi$. For $0.936 < D/L < 1.109$ the rod will jump into a self-contacting solution at some $|R| > 2\pi$. However, in Section 5 we shall see that even within the range of parameters where snap-through is possible in principle, jumps to self-contact may occur because of self-contacting solution branches interfering with the curves in Fig. 7.
<figure><img src="image_7.png" /><figcaption>Figure 6: *D-T* diagrams showing the second-mode bifurcation: (a) $\gamma = 1.3$ and curves are for $R = 0$ (p), 0.628 (q), 1.885 (r) and 3 (s); (b) $\gamma = 1.7$ and curves are for $R = 0$ (p), 0.1 (q), 0.3 (r), 0.6 (s) and 1 (t). Solid branches are stable, dashed branches unstable. The squares indicate the second-mode bifurcation, the triangles the out-of-plane bifurcation, while the circles mark the ring.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Load-deflection diagram with solution curves for $D/L = 0.5, 0.62, 0.8, 0.95$ and $1.2$. These curves have to be imagined closed by identifying the end points at $R = \pm 2\pi$, where the rod is a planar self-intersecting non-inflectional elastica. At the origin the rod is an inflectional elastica. Solid lines represent stable branches, dashed lines unstable branches. ($\gamma = 5/7$.)</figcaption></figure>

Note that for values of $D/L$ smaller than that of the into-the-plane instability ($D/L = 0.6767$), the curves in Fig. 7 also show the fold instability with jump into self-contact discussed in the previous section.

As in Fig. 3, Fig. 7 shows only curves with stable branches. Other curves exist at the same value of $D/L$. For example, look ahead to the dashed curve in Fig. 10 drawn there because a (stable) branch of self-contacting solutions bifurcates from it. This dashed curve, for $D/L = 1.2$, runs from $(R, M) = (-4\pi, 0)$ to $(R, M) = (4\pi, 0)$, at which points the solution is a second-mode non-inflectional (i.e., 2-looped) elastica.

## 3.5 A scaling law for $R$ and $\gamma$

In Fig. 5 we varied $\gamma$ to study what happens when the ring goes unstable. In fact, as pointed out in [8], there exists a scaling law which allows one to obtain the behaviour for any $\gamma$ from knowledge about solution curves for any particular $\gamma$. The rule, explained in Fig. 8, is based on the observation that the difference $W := R - ML/C = R - 2\pi m/\gamma$, $m := ML/(2\pi B)$, like the writhe $Wr$ to which it is related, is only a function of the centreline of the rod and therefore independent of $\gamma$. Thus, if for $\gamma = \gamma_1$ the solution curve in the $R-m$ plane is known, then the solution curve for $\gamma = \gamma_2$ can be found by translating $R$, at fixed $m$, by an amount equal to the horizontal distance between the straight lines of slope $\gamma_1/(2\pi)$ and $\gamma_2/(2\pi)$. It may be convenient in calculations to take for $\gamma$ the value 1, so that the equations for $d_3$ and $d_1$, (9) and (11), become identical, with just different initial conditions.

In view of this scaling law, the complicated spiralling behaviour of Fig. 5 for fixed $R$ and
<figure><img src="image_9.png" /><figcaption>Figure 8: Scaling behaviour under a change of $\gamma$ of solution curves in the R-m diagram, where $m = ML/(2\pi B)$. The straight lines are given by $m = R\gamma/(2\pi)$. The two double arrows have equal lengths.</figcaption></figure>

varying $\gamma$ must have an analogue in the case that $\gamma$ is kept fixed and $R$ is varied.

# 4 Formulation of a rod with point contacts

The thickness of the rod, and the fact that a real rod cannot intersect itself, can be taken into account on the level of 1D rod theory by considering an impenetrable elastic tube of certain radius about the centreline of the rod. This technique was applied by Shi & Hearst [38] and Miyazaki & Kondo [33] to configurations with a single point of contact, and by Coleman and his co-workers (see [8], [42], [46] and references therein) to more general situations. Most of these studies, but not [33], are concerned with supercoiling of DNA plasmids, i.e., closed pieces of DNA. In one study in [42] open rods are considered, but subject to pinned boundary conditions as defined in Section 3.1, reflecting the situation in recent experiments on segments of individual DNA molecules placed in a magnetic trap [40]. (The author actually calls it Cardan-joint boundary conditions but these are more commonly taken to mean that no moment can be transmitted about one fixed axis and one body axis; see, e.g., Section VIII.12 of [2].) Both pinned and Cardan-joint boundary conditions are rarely encountered in engineering applications. In [33] a clamped strut is considered.

Here we apply the same technique to our clamped rod. We only consider configurations with isolated points of contact. Coleman & Swigon [8] also studied rods with intervals of continuous self-contact for the special case of a straight line of contacts. Such configurations are encountered in ply structures seen in DNA molecules and twisted textile yarns, and in [8] were found to be the natural state reached after going through a number (two or three) of point contacts. We have something to say about rods contacting along an interval in a separate publication [44].

We follow the analysis of Section 2 and choose the origin of arclength in the middle of the rod. Although in case of one or several points of self-contact there may be non-reversible solutions satisfying clamped boundary conditions (there were for pinned boundary conditions in [8]), we shall again only look at reversible solutions having $m_{0y} = 0$. At a point of self-contact (assumed to be
frictionless) a reactive concentrated force will act between two non-neighbouring cross-sections of the rod. As usual, this force will be assumed to be normal to the surface of the rod. Consequently, the constant internal force **F** will undergo a discontinuous change. The moment **M**, on the other hand, will be assumed to be continuous (i.e., no concentrated moments), as will be **R** and the directors **d**<sub>i</sub> (and hence u<sub>3</sub>). From equation (7) we see that in order to ensure continuity of **M** across a point of contact, **M**<sub>0</sub> will have to change. Thus we consider equations (3), (9) and (11) with subscripts and write

$$
\mathbf{R}' = \mathbf{d}_3, \tag{38}
$$

$$
B \mathbf{d}'_3 = (\mathbf{F}_i \times \mathbf{R} + \mathbf{M}_{i0}) \times \mathbf{d}_3, \quad (39)
$$

$$
B \mathbf{d}'_1 = (\mathbf{F}_i \times \mathbf{R} + \mathbf{M}_{i0}) \times \mathbf{d}_1 + (B - C) u_3 \mathbf{d}_3 \times \mathbf{d}_1. \quad (40)
$$

Here **F**<sub>i</sub> and **M**<sub>i0</sub> are the as yet unknown constant internal force and matching moment in the section [S<sub>i</sub>, S<sub>i+1</sub>], where S<sub>0</sub> = 0, S<sub>n+1</sub> = L/2, and n is the number of contacts. (The property of having normal reactive forces and no concentrated moments has been rigorously established for a wide class of contact problems by Schuricht [37].)

We shall only consider solutions with symmetric contacts, i.e., contacts at *S*-symmetric points **R**(*S*) and **R**(*−*S). These include the configurations the rod jumps into at the limit points in Fig. 3, as well as the writhing configurations obtained when subsequently moving the ends in. At the first point of contact, *S*<sub>1</sub>, we have the geometrical relations

$$
|\mathbf{R}(S_1) - \mathbf{R}(-S_1)| = 2\rho L, \quad (41)
$$

$$
(\mathbf{R}(S_1) - \mathbf{R}(-S_1)) \cdot \mathbf{d}_3(\pm S_1) = 0, \quad (42)
$$

where ρ is the dimensionless radius of the rod. Since **F** is defined to be the force from the element at *S'* acting on the element at *S*, where *S'* > *S*, we have

$$
\mathbf{F}_0 = \mathbf{F}_1 + \Delta \mathbf{F}_1, \quad (43)
$$

with

$$
\Delta \mathbf{F}_1 = \frac{\Delta F_1}{2\rho L} (\mathbf{R}(S_1) - \mathbf{R}(-S_1)) \quad (44)
$$

the force, of magnitude Δ*F*<sub>1</sub>, pointing from *S* = −*S*<sub>1</sub> to *S* = *S*<sub>1</sub>. For **M** to be continuous we need

$$
\mathbf{M}(S_1) = \mathbf{F}_0 \times \mathbf{R}(S_1) + \mathbf{M}_{00} = \mathbf{F}_1 \times \mathbf{R}(S_1) + \mathbf{M}_{10}, \quad (45)
$$

which yields

$$
\mathbf{M}_{10} = \mathbf{M}_{00} + (\mathbf{F}_0 - \mathbf{F}_1) \times \mathbf{R}(S_1) = \mathbf{M}_{00} + \Delta \mathbf{F}_1 \times \mathbf{R}(S_1). \quad (46)
$$

So we can write

$$
\mathbf{F}_1 = (0, 0, F) - \frac{\Delta F_1}{\rho L} (X(S_1), 0, Z(S_1)), \quad (47)
$$

$$
\mathbf{M}_{10} = (M_{00x}, 0, M_z) + \frac{\Delta F_1}{\rho L} (-Y(S_1) Z(S_1), 0, X(S_1) Y(S_1)). \quad (48)
$$

The same can be done at each contact and we obtain

$$
\mathbf{F}_i = \mathbf{F}_{i-1} - \frac{\Delta F_i}{\rho L} (X(S_i), 0, Z(S_i)), \quad (49)
$$

$$
\mathbf{M}_{i0} = \mathbf{M}_{(i-1)0} + \frac{\Delta F_i}{\rho L} (-Y(S_i) Z(S_i), 0, X(S_i) Y(S_i)) \quad (50)
$$
as the constants to be used in (39), (40).

Thus in our shooting method the jump conditions (49) and (50) give two extra equations for the two unknowns $\Delta F_i$ and $S_i$ for each contact-free section of rod. Together with the two clamp conditions (17) applied to the terminal section, and the equations (23) and (one of) (24) for the end shortening and end rotation this gives a well-posed problem of $2n+4$ equations with the same number of unknowns for a (symmetric) solution with $n$ points of self-contact. This count may be compared with the $12n+6$ according to the method used in [8] for solving full rather than half solutions.

Since at contact points $M$ is continuous while $F$ has a jump in a direction orthogonal to $d_3$, we see that of the first integrals identified in Section 2.1 the twist $u_3$ and the Hamiltonian $M \cdot u/2 + F \cdot d_3$ will be conserved along the entire self-contacting rod, while $f$ and $m_z$ will jump at contact points.

# 5 Numerical results for self-contact

As before, we take $\gamma = 5/7$, and we further fix $\rho = 0.03/(2\pi)$. Results for fixed end rotation are shown in Fig. 9. It is found that for $R$ less than about $4.8\pi$ the rod delocalises under increasing end shortening $D$ with the 1-contact solution approaching the ring. For $R$ larger than about $4.8\pi$ curves in the $D-T$ diagram form limit points in $D$. Physically, the loop starts to rotate and the rod writhes up into a ply. Note that the 1-contact solution goes into compression before it becomes unstable and jumps into a state with more contacts. Diamonds along the unstable branches indicate where self-penetration starts to occur. No bifurcating branch of solutions could be detected, but the numerical approach of this self-penetration suggests that a branch of solutions with an interval of self-contacts might bifurcate, which we can't pick up with our present techniques. This is different from the constant-$D$ case discussed below where a branch of 2-contact solutions bifurcates off the 1-contact branch. Fig. 9(b) shows how the force $\Delta F$ varies with $D$.

Results for fixed end shortening are shown in Figs 10, 11 and 12. In Fig. 10(a), for $D/L = 0.5$, we find solutions with up to 3 points of self-contact. The transition from the contact-free to the 1-contact, and the transition from the 1-contact to the 2-contact solution, are through dynamic jumps at folds, while the transition from 2-contact to 3-contact is smooth. True-view 3D shapes of the rod at these events are shown in Fig. 12. Extensive numerical searches did not give any hint of the existence of a branch of 4-contact solutions as $R$ is increased further. This is consistent with the work of Coleman & Swigon [8] who find, for closed rods, that at some point the third self-contact grows into an interval of self-contacts, still flanked by two points of self-contact.

Fig. 10(b) shows that for $D/L = 0.86$ multistability occurs: at the fold along the (solid) contact-free curve the rod can either snap-through and remain without self-contacts, or jump into a configuration with one point of self-contact. (Of course, we cannot rule out the existence of further stable solutions.) Experiments reveal that the contact-free option is normally taken. In fact, the 1-contact curve goes round a fold and intersects the free rod curve a second time at $R = -1.2841$ where again $\Delta F$ drops to zero. In view of the fold rule this 'pop-out' only involves unstable branches. Upon increase of $D/L$ the two pop-out points coalesce at $D/L = 0.8946$, leaving a contact curve disconnected from the contact-free curve.

Multistability is also seen in Fig. 10(c), for $D/L = 1.2$ (i.e., for inverted solutions). Experiments show that this time the contact-free route is normally taken. The branch of 1-contact solutions does not bifurcate from the usual contact-free curve but from an entirely unstable
<figure><img src="image_10.png" /><figcaption>Figure 9: Writhing a rod or not, by varying $D/L$ and fixing $R$ at $4\pi$, $5\pi$, $6\pi$ and $8\pi$. For $R = 4\pi$ the rod does not writhe and the curve of 1-point contacts runs up to the ring at $D/L = 1$ (indicated by the circle). In contrast, the curves for $R = 5\pi$, $6\pi$ and $8\pi$ form a limit point, signifying that the rod writhes up into a ply. The transition to writhing takes place between $R = 4.8\pi$ and $R = 5\pi$. Triangles indicate pop-out points where $\Delta F$ drops to zero (see (b)) and self-contact is lost. Diamonds indicate points where, presumably, an interval of self-contact is formed. ($\gamma = 5/7$, $\rho = 0.03/(2\pi)$.)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>Figure 10: Ply formation under varying $R$ and (a) $D/L = 0.5$, (b) $D/L = 0.86$, (c) $D/L = 1.2$. All the lower branches coming in at the various folds are unstable. Triangles indicate pop-out points where one self-contact is lost/gained. In (a) curves for up to three points of self-contact are included. Labels refer to the 3D shapes shown in Fig. 12, while the inset shows that the 2-contact curve forms a limit point to the left. In (b) the dashed curve is for 1-contact solutions. Insets give true $z-y$ views of the rod's centreline at the indicated points. In (c) the 1-contact curve bifurcates from the dashed curve of non-contacting solutions running from $(R, M) = (-4\pi, 0)$ to $(R, M) = (4\pi, 0)$, at which end points the solution is an (inverted) second-mode non-inflectional elastica. ($\gamma = 5/7$, $\rho = 0.03/(2\pi)$.)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11: Formation of the third point of contact along the 2-contact branch in Fig. 10(a). $\Delta = 2\sqrt{x^2 + z^2}$ is the symmetric self-distance of the rod; $s$ is arclength. The minimum possible value of $\Delta$ is $2\rho$, attained at self-contact, and is indicated by the dotted horizontal line in the left figures. The right-most self-contact turns out to have the largest pressure force $\Delta F$. The projections of the shape of the rod (centre and right) are true views (from infinity) and show the rotation of the end loop. ($\gamma = 5/7$, $\rho = 0.03/(2\pi)$.)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 12: True 3D views of the rod along stable solution paths in Fig. 10(a). Labels correspond to the labels in that figure. The starting point (a) is at the planar elastica. The first jump, into 1-point contact (from (b) to (c)), occurs at $R = 11.6315$; the second jump, into 2-point contact (from (d) to (e)), occurs at $R = 18.4027$. The transition to 3-point contact at (f) (corresponding to (c) in Fig. 11) is smooth.</figcaption></figure>
contact-free curve (dashed in the figure) which at $(R, M) = (\pm 4\pi, 0)$ connects to a second-mode non-inflectional planar elastica. The solution at $R = 0$ along this dashed curve is, geometrically, close to a figure-of-eight held at the edge (for $D/L = 1$ it would be an exact figure-of-eight). The dashed curve has an additional intermediate zero of $M$ where the solution is non-planar.

Fig. 11 illustrates the formation of the third point of contact along the 2-contact branch in Fig. 10(a) using the symmetric self-distance

$$ \Delta(s) = \sqrt{[x(s) - x(-s)]^2 + [y(s) - y(-s)]^2 + [z(s) - z(-s)]^2} = 2\sqrt{x(s)^2 + z(s)^2} \quad (51) $$

to monitor the process. Notice the rotation of the end loop relative to the rig as $R$ is increased.

## 5.1 Pop-out

Yabuta et al. [52, 53] studied the stability of thin looped self-contacting rods by looking at the second variation of the total potential energy. By assuming one turn of a circular helix as the shape of the loop, and a perfectly straight rod elsewhere, a stability criterion could be derived. A simplified criterion was obtained by further assuming that all the twist in the rod is removed by looping (i.e., no residual twist) and taking zero rod thickness, but this gives non-physical results in sharp contrast with the exact result (Eq. (37) above) by Zajac (such as the prediction of a critical $\gamma$ above which no pop-out occurs in the planar elastica, which therefore would kink, or ‘hockle’).

If the looped rod has residual twist, i.e., if $R > 2\pi$, then at the point of self-contact a force will act to maintain the contact. If the end shortening is decreased by pulling the ends (keeping $R$ fixed) the loop will be tightened up to a point where the loop pops open. This pop-out only occurs in rods of finite thickness; an infinitely thin rod with $R > 2\pi$ would not go unstable and could be pulled out straight. In many practical cases a cable or pipe will kink (i.e., suffer plastic deformation) before pop-out is reached, damaging the structure.

Zajac [55] studied pop-out of a clamped rod of zero thickness by analysing the 3D variational equations about the clamped planar looped elastica. The result is the into-the-plane critical curve we have already seen in Fig. 4. Indeed, Coyne [9] points out that the Zajac condition (37) is best interpreted as a criterion for looping rather than pop-out. Our analysis shows precisely how Zajac’s result (for zero-radius rods) relates to both, since it applies to $R = 2\pi$ only, at which point looping and pop-out become the same phenomenon evolving in opposite directions.

Pop-out may also occur in rods with multiple points of self-contact when one contact point is lost (as in Fig. 10(a)). At the point of pop-out (out of a single contact) the shape of the rod is again governed by the contact-free rod equations of Section 2.1. Therefore, in the numerical continuations which produced Fig. 9 we can infer pop-out when the self-distance $\Delta$ reaches $2\rho$. The locus of pop-out points thus obtained is shown by the solid curve in Fig. 13 on which some $R$-values are highlighted. Along a 1-point contact curve in Fig. 9 pop-out is defined by the intersection point of this curve with the free rod curve at the same value of $R$. At this point the force jump $\Delta F$ drops to zero (see Fig. 9(b)). For a rod of zero radius pop-out is described mathematically by a pitchfork bifurcation; at finite radius it becomes a type of discontinuous bifurcation.

A semi-analytical estimate of the conditions at pop-out can be obtained by following (and slightly improving upon) Coyne [9] in using the homoclinic solution of the rod equations already mentioned in Section 3.3. For this homoclinic solution the shape of the rod is given by (see, e.g., [38]):

$$ x(s) = \frac{1}{2\pi t} \sqrt{4t - m_z^2} \operatorname{sech} \left( \pi s \sqrt{4t - m_z^2} \right) \sin(m_z \pi s), $$
<figure><img src="image_16.png" /><figcaption>Figure 13: Pop-out locus parametrised by $R$. Markers are for $R = 2.5\pi$, $3\pi$, $3.5\pi$, $4\pi$, $5\pi$, $6\pi$, $7\pi$ and $8\pi$ (in increasing order of tension). The solid curve is obtained from our criterion $\Delta = 2\rho$, the long dashed curve represents results based on the homoclinic orbit, while the short dashed curve shows Yabuta's approximation. The inset shows further results with markers for $R = 2.25\pi$ and $R = 2.1\pi$. For $R = 2.25\pi$ and higher the homoclinic results are very accurate. ($\gamma = 5/7$, $\rho = 0.03/(2\pi)$.)</figcaption></figure>

$$y(s) = -\frac{1}{2\pi t} \sqrt{4t - m_z^2} \operatorname{sech} \left( \pi s \sqrt{4t - m_z^2} \right) \cos(m_z \pi s), \quad (52)$$

$$z(s) = s - \frac{1}{2\pi t} \sqrt{4t - m_z^2} \tanh \left( \pi s \sqrt{4t - m_z^2} \right),$$

with the self-contact occurring at arclength $s = \bar{s}$ such that

$$z(\bar{s}) = 0, \quad x(\bar{s}) = \rho. \quad (53)$$

Further, in [50] the expressions

$$D = (2/T) \sqrt{4BT - M^2} \quad \text{and} \quad R = \frac{ML}{C} + 4 \arccos \left( \frac{M}{2\sqrt{BT}} \right) \quad (54)$$

are derived for the end shortening and end rotation of the homoclinic solution. From (52)–(54) we can find the tension and end shortening at pop-out as a function of the radius of the rod, $\rho$, and the amount of end rotation, $R$, applied. The result is the long dashed curve in Fig. 13. Also displayed in this figure are results (having $d = 1/\sqrt{2t}$) following from Yabuta's analysis taking $u_3 = (R - 2\pi)/L$ for the residual twist. (For the Yabuta curve, Eq. (34) in [53] is used, which gives better agreement than the corresponding but conflicting result quoted in [52].) It is seen that for $R > 2.25\pi$ the results based on the homoclinic approximation are virtually indistinguishable from the 'exact' numerical results. Notice that the pop-out curve in the $D-T$ diagram (in fact any of
the three curves in Fig. 13) considered as a whole is independent of $\gamma$; its parametrisation by $R$, of course, does depend on $\gamma$.

Depending on whether the 1-contact branch emanates from an unstable or stable branch of contact-free solutions, at pop-out the rod undergoes either a dynamic jump or a smooth transition out of contact. The former is the case in Fig. 9(a) where the rod jumps to a contact-free state at much lower tension (or perhaps compression) $T$ described by the point vertically down from the pop-out point on the same constant-$R$ curve. Note from Fig. 9(b) that $\Delta F$ goes through a maximum before it rapidly drops to zero and pop-out occurs. An alternative way of jumping out of contact is provided by the small-scale limit point seen along the 2-contact branch in Fig. 10: if one were to decrease $R$ along this branch the rod would jump out of 2-point contact before the force $\Delta F$ dropped to zero.

An example of smooth pop-out is seen in Fig. 10(b).

## 5.2 Rig-wrench behaviour

As we have seen in our formulation, the forces and moments acting across any cut through a finite-length rod will be equivalent to a wrench acting about a ‘wrench axis’ in space. It is for this reason that any rod of finite length, deforming in three dimensions under end forces and moments, will lie conceptually within the field of solutions (the spatial elastica) that can be adopted by an infinite rod loaded by equal and opposite wrenches. This makes the spatial elastica a natural framework within which to explore and understand the behaviour of finite rods. We do, however, have to take cognisance of the fact that during loading the wrench of the elastica solution, which before buckling lies along the axis of the testing rig, will move spatially relative to this ‘rig axis’ during post-buckling.

The straight lines in Fig. 14 indicate the equivalent wrench lines for a contact-free rod as well as for (contact-free pieces of) a rod with two points of self-contact. The wrench line undergoes a step change at a contact point on account of the step change in the internal force $f$. The sequence of plots in Fig. 11 also nicely illustrate the movement of the wrench relative to the rig during the loading process (recall that $e_z$ was chosen along the wrench axis). This movement is composed of a translation $h$ in the $y$ direction and a rotation $\phi$ in the $x-z$ plane given by

$$h = \frac{M_{0x}}{F}, \quad \cos \phi = d_{3z}(L/2). \qquad (55)$$

This may be compared with the situation for the clamped planar elastica which, as shown by the insets in Fig. 2(b), has a pure translation for the odd modes and a pure rotation for the even modes.

For a well-localised solution in a long rod the wrench line is close to the line of the rig.

## 6 Experiments

In this section we report on some laboratory experiments we carried out to supplement our theoretical investigations of the various jump phenomena. Surprisingly few experimental studies exist in the literature. We can only mention Yabuta’s study [53] of pop-out for which jacketed optical fibres were used because of their high flexibility. However, the author mentions neither the exact boundary conditions nor the applied loading sequence.
<figure><img src="image_17.png" /></figure>

Figure 14: Examples of a contact-free (a) and a 2-contact (b) rod configuration with the equivalent wrench lines included. Also shown are the projections of these lines on the plane parallel to the x-z plane and going through both ends of the rod. The wrench line of a (piece of) free rod is the line about which the internal force and the equivalent moment act.

The choice of material for our experiments was constrained by the requirements that the rod must be highly flexible in bending and twisting but must not suffer from sag due to self-weight. Also, for a quantitative study the loads involved must be measurable (which effectively excludes rubber rods). The superelastic nickel titanium alloy nitinol proved to be a suitable material. Nitinol has found extensive use in a wide range of applications. For instance, in the medical industry it is used in orthopaedic devices and orthodontic arches. It is also employed as spectacle frames. The properties of nitinol are strongly dependent on processing history and a wide range of values for *B* are given by different manufacturers. Therefore, the elastic rigidities, *B* and *C*, were estimated directly by means of simple cantilever and torque tests. These gave linear constitutive relations and a ratio of $\gamma = 5/7$. Experiments were limited to the range $0 \le R \le 5\pi$. Within these limits the specimens remained in the linear elastic range.

Circular-cross-section nitinol rods with a radius of 0.5 mm and ranging in length from 350 to 500 mm (from 1000 to 1200 mm for inverted experiments) were clamped in a specially constructed rig (see Fig. 1) in which their ends could be rotated and displaced independently. By means of suitably calibrated transducers fixed in line at the ends of the rod the rig moment *M* and rig force *T* could be measured to an accuracy of 0.0001 Nm and 0.0001 N, with values ranging from -0.2 to 0.2 Nm and -0.5 to 1 N, respectively. Some experiments involved the recording of loads on the limit of the sensitivity of the transducers. However, with care in selecting lengths of rod (short enough to avoid important sagging effects, yet long enough to have an appreciable linear elastic range) good repeatability of results was obtained. Each of the data points shown in Figs 15(a) and (b) represents an average over at least three and at most five measurements.

Starting with a straight untwisted rod, two sets of experiments were performed. In one set the end rotation, *R*, was input and kept fixed while the end shortening, *D*, was controlled, i.e., the aligned ends were moved towards each other in small increments. For each adjustment of *D* the end tension, *T*, and twisting moment, *M*, were recorded. Some results are shown in Fig. 15 where it can be seen that reasonably good correlation between the mathematical model and the experiment is obtained. For *R* = 0 and *R* = $2\pi$ the pitchfork bifurcations to and from a planar state occur as predicted theoretically. For *R* > $2\pi$ the values of *D* at which the rod jumps into a loop are also in good agreement with the theory. After loop formation experiments were continued
<figure><img src="image_18.png" /><figcaption>Figure 15: Experimental results for a clamped nitinol rod superimposed on theoretical curves: (a) tests at fixed R for R = 0, π, 2π, 3π, 4π and 5π (including tests on inverted rods) with curves of Fig. 3; (b,c,d) tests at fixed D for D/L = 0.5, 0.677 and 0.8, respectively, with curves of Fig. 7. Notice the hysteresis loops (indicated by the arrows) associated with both snap-through and pop-out. ($ρ = 1/900$.)</figcaption></figure>
under both increasing and decreasing *D*. It is found that by increasing *D* the loop remains intact up to values of *R* ≈ 4.5π with the shape of the loop evolving into a more circular form as *D* approaches *L*. However, for *R* > 4.5π the loop itself rotates to form a ply, because the loop does not release enough torsional strain energy during loop formation. Theoretically the transition to writhing was located near 4.8π for a slightly different slenderness ratio *ρ* (see Section 5).

It is also possible to induce the rod to pop out of the loop by decreasing *D*. In such an experiment the loop dimensions diminish causing the curvature, and therefore the bending moment around the loop, to increase. This can cause permanent plastic damage (i.e., kinking). In fact, it was found that a completely successful (i.e., kink-free) pop-out experiment is only possible for a small amount of end rotation (up to about 0.1 radian above *R* = 2π).

In the second set of experiments the initially straight rod was first bent by an input of *D*, and the experiment was continued under control of *R* keeping *D* fixed. For *D* less than that at the out-of-plane bifurcation (*D*/*L* = 0.5990) the only jump encountered is into a loop at high *R*. This jump corresponds to the fold for *R* > 2π in the previous set of experiments. However, once *D* reaches the critical value at which the out-of-plane bifurcation occurs, configurations in the vicinity of *R* = 0 become unstable and snap-through occurs. The rod jumps through the unstable region, and if *R* is reversed it does not jump back again at the same value, i.e., an hysteresis loop emerges. As the end shortening is increased, the extent of this hysteresis cycle spreads, and with it the magnitude of the snap buckling. However, this development is accompanied by a corresponding diminution of the jump into (and out of) a self-contacted loop (which all takes place at much higher values of *R* than snap buckling).

At *D*/*L* = 0.6767 the jump into a loop finally disappears and is replaced by a smooth path both in and out of self-contact. Some results are shown in Fig. 15(b,c,d) where the changing size of these jumps for four different settings of *D*/*L* can be clearly seen. Also in Fig. 15 the hysteresis cycles associated with the snap-through and the jumping in and out of a loop are clearly illustrated.

In the vicinity of *D*/*L* = 1 the rod behaves like a ring. A twisted ring becomes unstable at a critical twisting moment given by the Zajac condition (35) and jumps into a figure-of-eight with two points of self-contact. Reversing the experiment causes the rod to jump to a stable figure-of-eight with one point of self-contact, creating a small hysteresis cycle, before jumping out of contact back into a ring. This behaviour is in agreement with the theoretical results in [8] (although these are for a slightly different value of *γ*). Unfortunately, the rig does not allow load measurements in this region. As *D* is increased beyond it the loads are measurable again and under control of *R* the rod jumps from one side of *R* = 0 to the other, but without self-contact.

In conclusion, most of the jump phenomena and associated hysteresis loops found theoretically can be observed, with reasonably good quantitative agreement, in relatively simple laboratory tests. Deviations from the theory are likely to be the result of the effect of gravity and sensitivity to small errors in satisfying exact clamped boundary conditions. For instance, experiments with anisotropic rods (i.e., ones for which $B_1 \neq B_2$) suggest that the vertical offset between the theoretical curve and the experimental data at primary buckling (on the vertical axis in Fig. 15(a)) might be caused by gravity. A strongly anisotropic rod will tend to bend only about its weakest principal axis, either up (hogging) or down (sagging), and experiments reveal a difference in behaviour depending on which direction the rod buckles in. Details of these and other experiments can be found in [16].
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

Figure 16: Diagrams summarising the various jump phenomena for (a) $\gamma = 5/7$, (b) $\gamma = 1.7$ and (c) $\gamma = 2.1$. Triangles mark the out-of-plane bifurcation, diamonds the into-the-plane bifurcation and squares the second-mode bifurcation.
# 7 Discussion

We have made a comprehensive study of jump phenomena in clamped rods with and without self-contact. A summary of the results is given in Fig. 16 where D-R displacement diagrams are shown for three typical values of $\gamma$. The curves are loci of fold and pop-out points obtained by numerical continuation (pop-out being defined by the condition $\Delta = 2\rho$). Bifurcation points are indicated by the same markers as used in the main text. Diagram (a) is for $\gamma = 5/7$, the value for nitinol used throughout. Diagram (b) is for $\gamma = 1.7$ just larger than $\gamma_c$ in (36) (for $n = m = 1$). Notice that the fold lines near $D/L = 1$ have formed some cusps leading to more complicated behaviour. Diagram (c) is for $\gamma = 2.1$. For this value of $\gamma$ the out-of-plane bifurcation occurs for $D/L > 1$ (cf. Fig. 4), i.e., for inverted solutions, and the out-of-plane and second-mode bifurcations have changed order. The diagrams can be used to locate regions in the displacement plane where stable solutions exist.

We have only looked at clamped boundary conditions. Although they occur widely in engineering applications they are by no means the only ones found, and it is clear that the bifurcation behaviour observed depends strongly on the type of boundary conditions considered. For instance, for clamped rods we found the circular ring to be an organising centre of bifurcation behaviour; yet for instance for ball-and-socket boundary conditions (i.e., $M(L/2) = 0$) the ring cannot even be supported. Meanwhile, the solutions for a rod whose ends are held in a sleeve form a subset of the solutions for a clamped rod, namely those that have no twist, i.e., $u_3 = 0$ (cf. case 8 in Table 2).

Our analysis does not depend on the length scale $L$ of the problem and therefore applies to DNA molecules, carbon nanotubes [28] and submarine communication cables [9] alike. However, gravity effects, ignored in this study, do depend on the length scale. Consider for example the relative deflection (deflection/length) due to self-weight of a horizontal rod of circular cross-section of radius $r$. Engineering bending theory shows that this scales as $(\rho g/E)(L/r)^2 L$, where $\rho$ and $E$ are the density and Young's modulus of the material, $g$ is the gravitational acceleration and $L$ is the length of the rod. To minimise gravitational effects in an experiment we therefore want a material with low $\rho/E$ and a low aspect ratio, $L/r$, consistent with still having a 'thin' rod. More significantly, we want an absolute low value for $L$. Relative gravitational deflections are small for physically small objects, large for physically large objects. This is why we cannot have land animals much larger than elephants or much taller than giraffes.

# A Appendix: clamped solutions are reversible

We show that clamped solutions to the equilibrium equations (3), (9) are reversible. This generalises a result by Domokos [12] stating the same for closed rods. Introducing Euler angles $\theta$ and $\psi$ to parametrise the unit vector $d_3$ as

$$d_{3x} = \sin \theta \cos \psi, \quad d_{3y} = \sin \theta \sin \psi, \quad d_{3z} = \cos \theta, \qquad (56)$$

and using the conserved quantities discussed in Section 2.1, we can reduce (3), (9) to the following system of (dimensionless) equations:

$$\theta'' + \frac{dV(\theta)}{d\theta} = 0, \quad V(\theta) = \frac{(m_z - m_3 \cos \theta)^2}{2 \sin^2 \theta} + f \cos \theta, \qquad (57)$$

$$\psi' = \frac{m_z - m_3 \cos \theta}{\sin^2 \theta}, \qquad (58)$$
<figure><img src="image_22.png" /><figcaption>Figure 17: Phase portraits for the planar oscillator (57): (a) typical case, ($m_z, m_3, t$) = ($-0.52, -0.76, 0.8$), (b) special case for $m_z = m_3 \neq 0$ and $m_z^2/f \in (0, 4)$ with a pair of homoclinic orbits, ($m_z, m_3, t$) = ($-0.6, -0.6, 0.8$).</figcaption></figure>

where $m_3 = m \cdot d_3$ (details can be found in [50]). The equation for $\theta$ is the well-known equation for the inclination angle of the Lagrange top. Equation (58) for $\psi$ can be solved once $\theta$ is known. This system of equations is invariant under the transformation $s \to -s$, $\theta \to \theta$, $\theta' \to -\theta'$, $\psi \to -\psi$, and is therefore reversible.

Phase portraits of (57) have been studied in [3]. By reversibility they have reflection symmetry about the $\theta$ axis. In addition they are symmetric about the vertical axis. The generic phase portrait is depicted in Fig. 17(a). It has two sets of closed level curves of the Hamiltonian $H = \theta'^2/2 + V(\theta)$ around two fixed points. (These fixed points correspond to helical solutions.) Degenerate cases occur for special values of the parameters $m_z, m_3$ and $f$: for $m_z = m_3 = 0$ (57) reduces to the equation for the planar elastica, for $m_z = m_3 \neq 0$ and $m_z^2/f \in (0, 4)$ we obtain the phase portrait with homoclinic orbits shown in Fig. 17(b), while for $m_z = m_3 \neq 0$ and $m_z^2/f \notin (0, 4)$ the phase portrait consists solely of closed level curves around the origin, much like those outside the homoclinic orbit in Fig. 17(b). In the following we shall only deal with the generic case and consider the closed level curve not crossing the vertical axis highlighted in Fig. 17(a). The proof can easily be adapted to deal with a level curve encircling the origin.

For $\theta$ the boundary condition (13) implies that

$$ \theta(-1/2) = \pm\theta(1/2) \pmod{2\pi}. \qquad (59) $$

By symmetry we can restrict our attention to the positive sign, so that (13) further yields (since $\theta \neq 0$)

$$ \psi(-1/2) = \psi(1/2) \pmod{2\pi}. \qquad (60) $$

In order to be compatible with our choice of co-ordinate frames in Section 2.1 and to ensure that $d_{3y}(0) = 0$ we have to take as initial condition $\psi(0) = 0$, since $\theta(0) \neq 0 \pmod{\pi}$.

Now, on our level curve consider an orbit with initial conditions at O. To obtain a full solution we have to integrate from this point, where $s = 0$, backward to $s = -1/2$ (point A, say) and forward to $s = 1/2$ (point B, say). Let point C be such that $s_C = T_p/2$, where $T_p$ is the period of the orbit, and let D and E be the reflections of C and A about the $\theta$ axis. S and T are the
intersection points of the level curve with the $\theta$ axis and lie in the fixed point set of the reversing transformation.

We now apply the boundary conditions (59) and (60). There are two cases:

(i) The orbit under consideration has exactly period $T_p = 1$. Then the points A, B, and C coincide, and we shall show below that the boundary conditions imply that the centreline of the rod forms a closed curve. Hence the solution is congruent with the solution that has $s = 0$ at S (just 'slide the clamp along the rod'), and is thus reversible.

(ii) By (59), points A and B must have the same $\theta$ and hence B coincides with E. But since the arclength from A to C, $s_{AC}$, is equal to the arclength from C to B, $s_{CB}$, we have $s_{AC} = s_{CE}$. Further, we have $s_{CE} = s_{CD} + s_{DE}$, and since reversibility of (57) implies that $s_{AC} = s_{DE}$ we conclude that $s_{CD} = 0$. Hence C, T and D coincide, and by the definition of C this means that O coincides with S, i.e., the starting point must lie on the symmetry axis, making the solution reversible.

It remains to show that $T_p = 1$ implies a closed rod (case (i)). To see this, form the vector product of (9) with $\mathbf{d}_3$ and use (56) and (58) to obtain

$$fx + m_{0y} = \frac{m_3 - m_z \cos \theta}{\sin \theta} \sin \psi + \theta' \cos \psi, \quad (61)$$

$$-fy + m_{0x} = \frac{m_3 - m_z \cos \theta}{\sin \theta} \cos \psi - \theta' \sin \psi. \quad (62)$$

Then

$$x(1/2) - x(-1/2) = \frac{1}{f} \{\theta'(1/2) - \theta'(-1/2)\} \cos \psi(1/2), \quad (63)$$

$$y(1/2) - y(-1/2) = \frac{1}{f} \{\theta'(1/2) - \theta'(-1/2)\} \sin \psi(1/2), \quad (64)$$

and since A and B coincide (hence $\theta'(1/2) = \theta'(-1/2)$) we find that the $x$ and $y$ co-ordinates at both end points are equal. The second clamp condition (14) then gives that $\lambda = 0$, i.e., the rod has a closed centreline.

These symmetric solutions with the clamps at a non-symmetric position are not found in our numerical work. In our shooting method we use one of the clamps as a reference position and therefore only solutions with symmetrically placed clamps are obtained. Closed solutions with non-symmetric clamps lie on straight lines ($d = 1$) in the $d-t$ diagram, as on the line connecting the two figures-of-eight in Fig. 2a.

## Acknowledgments

GH was supported by a Royal Society Research Fellowship, SN by a Royal Society Visiting Fellowship.
References

[1] S.S. Antman, Large lateral buckling of nonlinearly elastic rods, *Arch. Rational Mech. Anal.* **84**, 293-305 (1984).

[2] S.S. Antman, *Nonlinear Problems of Elasticity* (Springer-Verlag, Berlin, 1995).

[3] S.S. Antman & K.B. Jordan, Qualitative aspects of the spatial deformation of non-linearly elastic rods, *Proc. R. Soc. Edinb. 73A*, 85-105 (1974).

[4] M. Beck, Knickung gerader Stäbe durch Druck under konservative Torsion, *Ing.-Arch. 23*, 231-253 (1955).

[5] P.B. Béda, A. Steindl & H. Troger, Postbuckling of a twisted prismatic rod under terminal thrust, *Dyn. Stab. Systems* **7**, 219-232 (1992).

[6] C. Bouchiat & M. Mezard, Elasticity model of supercoiled DNA molecule, *Phys. Rev. Lett.* **80**, 1556-1559 (1998).

[7] G. Călugăreanu, Sur les classes d'isotopie des nœuds tridimensionnels et leurs invariants, *Czechoslovak Mathematical Journal* **11**, 588-625 (1961).

[8] B.D. Coleman & D. Swigon, Theory of supercoiled elastic rings with self-contact and its application to DNA plasmids, *J. Elasticity* **60**, 173-221 (2000).

[9] J. Coyne, Analysis of the formation and elimination of loops in twisted cable, *IEEE J. Ocean. Eng.* **15**, 72-83 (1990).

[10] E.J. Doedel, A.R. Champneys, T.R. Fairgrieve, Yu.A. Kuznetsov, B. Sandstede & X.J. Wang, AUTO97: Continuation and bifurcation software for ordinary differential equations (1997) (available by anonymous ftp from ftp.cs.concordia.ca/pub/doedel/auto).

[11] G. Domokos, Global description of elastic bars, *Z. Angew. Math. Mech.* **74**, T289-T291 (1994).

[12] G. Domokos, A group-theoretic approach to the geometry of elastic rings, *J. Nonlinear Sci.* **5**, 453-478 (1995).

[13] W.B. Fraser & D.M. Stump, The equilibrium of the convergence point in two-strand yarn plying, *Int. J. Solids Struct.* **35**, 285-298 (1998).

[14] F.B. Fuller, The writhing number of a space curve, *Proc. Nat. Acad. Sci. USA* **68**, 815-819 (1971).

[15] A. Goriely & M. Tabor, Nonlinear dynamics of filaments I. Dynamical instabilities, *Physica D* **105**, 20-44 (1997).

[16] V.G.A. Goss, G.H.M. van der Heijden, J.M.T. Thompson & S. Neukirch, Experiments on snap buckling, hysteresis and loop formation in twisted rods (submitted).

[17] A.G. Greenhill, On the strength of shafting when exposed both to torsion and to end thrust, *Proc. Instit. Mech. Eng.*, 182-225 (April 1883).
[18] E. Hairer, S.P. Nørsett & G. Wanner, *Solving Ordinary Differential Equations I: Nonstiff Problems*, 2nd ed. (Springer-Verlag, Berlin, 1993).

[19] P.J. Heath, J.B. Clendinning, B.S. Fujimoto & J.M. Schurr, Effect of bending strain on the torsion elastic constant of DNA, *J. Mol. Biol.* **260**, 718-730 (1996).

[20] K.A. Hoffman, R.S. Manning & R.C. Paffenroth, Calculation of the stability index in parameter-dependent calculus of variations problems: buckling of a twisted elastic strut, *SIAM J. Appl. Dyn. Systems* **1**, 115-145 (2002).

[21] D.S. Horowitz & J.C. Wang, Torsional rigidity of DNA and length dependence of the free energy of DNA supercoiling, *J. Mol. Biol.* **173**, 75-91 (1984).

[22] G. Iooss & D.D. Joseph, *Elementary Stability and Bifurcation Theory* (Springer-Verlag, New York, 1980).

[23] S. Kehrbaum & J.H. Maddocks, Elastic rods, rigid bodies, quaternions, and the last quadrature, *Phil. Trans. R. Soc. Lond. A* **355**, 2117-2136 (1997).

[24] S. Kehrbaum & J.H. Maddocks, Effective properties of elastic rods with high intrinsic twist, in *Proceedings of the 16th IMACS World Congress, Lausanne, 2000* (to appear).

[25] K. Kovári, Räumliche Verzweigungsprobleme des dünnen elastischen Stabes mit endlichen Verformungen, *Ing.-Arch.* **37**, 393-416 (1969).

[26] M. Le Bret, Catastrophic variation of twist and writhing of circular DNAs with constraint?, *Biopolymers* **18**, 1709-1725 (1979).

[27] Y. Li & J.H. Maddocks, On the computation of equilibria of elastic rods. Part I: Integrals, symmetry and a Hamiltonian formulation, preprint, Department of Mathematics, University of Maryland, College Park (1996) (available at http://lcvmwww.epfl.ch/publis.html).

[28] O. Lourie, D.M. Cox, H.D. Wagner, Buckling and collapse of embedded carbon nanotubes, *Phys. Rev. Lett.* **81**, 1638-1641 (1998).

[29] J.H. Maddocks, Stability of nonlinearly elastic rods, *Arch. Rational Mech. Anal.* **85**, 311-354 (1984).

[30] J.H. Maddocks, Stability and folds, *Arch. Rational Mech. Anal.* **99**, 301-328 (1987).

[31] R.S. Manning & J.H. Maddocks, Symmetry breaking and the twisted ring, *Comput. Methods Appl. Mech. Engrg* **170**, 313-330 (1999).

[32] R.S. Manning, K.A. Rogers & J.H. Maddocks, Isoperimetric conjugate points with application to the stability of DNA minicircles, *Proc. R. Soc. Lond. A* **454**, 3047-3074 (1998).

[33] Y. Miyazaki & K. Kondo, Analytical solution of spatial elastica and its application to kinking problem, *Int. J. Solids Struct.* **34**, 3619-3636 (1997).

[34] S. Neukirch, G.H.M. van der Heijden & J.M.T. Thompson, Writhing instabilities of twisted rods: from infinite to finite length, *J. Mech. Phys. Solids* **50**, 1175-1191 (2002).
[35] S. Neukirch & M.E. Henderson, Classification of the spatial clamped elastica, part I: symmetries and zoology of solutions (in preparation).

[36] K.A. Rogers, *Stability exchange in parameter-dependent constrained variational principles with applications to elastic rod models of DNA minicircles*, PhD thesis, University of Maryland, 1997.

[37] F. Schuricht, A variational approach to obstacle problems for shearable nonlinearly elastic rods, *Arch. Rat. Mech. Anal.* **140**, 103-159 (1997).

[38] Y. Shi & J.E. Hearst, The Kirchhoff elastic rod, the nonlinear Schrödinger equation, and DNA supercoiling, *J. Chem. Phys.* **101**, 5186-5200 (1994).

[39] T.R. Strick, J.-F. Allemand, D. Bensimon, A. Bensimon & V. Croquette, The elasticity of a single supercoiled DNA molecule, *Science* **271**, 1835-1837 (1996).

[40] T.R. Strick, J.-F. Allemand, D. Bensimon & V. Croquette, Behavior of supercoiled DNA, *Biophys. J.* **74**, 2016-2028 (1998).

[41] D.M. Stump & G.H.M. van der Heijden, Birdcaging and the collapse of rods and cables in fixed-grip compression, *Int. J. Solids Struct.* **38**, 4265-4278 (2001).

[42] D. Swigon, *Configurations with self-contact in the theory of the elastic rod model for DNA*, PhD thesis, Rutgers, State University of New Jersey, 1999.

[43] J.M.T. Thompson, Stability predictions through a succession of folds, *Phil. Trans. R. Soc. Lond. A* **292**, 1-23 (1979).

[44] J.M.T. Thompson, G.H.M. van der Heijden & S. Neukirch, Supercoiling of DNA plasmids: mechanics of the generalized ply, *Proc. R. Soc. Lond. A* **458**, 959-985 (2002).

[45] S. Timoshenko & J. Gere, *Theory of Elastic Stability* (McGraw-Hill, New York, 1961).

[46] I. Tobias, D. Swigon & B.D. Coleman, Elastic stability of DNA configurations. I. General theory, *Phys. Rev. E* **61**, 747-758 (2000).

[47] A. Trosch, *Ing.-Arch.* **20**, 265 (1952).

[48] G.H.M. van der Heijden & J.M.T. Thompson, Lock-on to tape-like behaviour in the torsional buckling of anisotropic rods, *Physica D* **112**, 201-224 (1998).

[49] G.H.M. van der Heijden, A.R. Champneys & J.M.T. Thompson, Spatially complex localisation in twisted elastic rods constrained to lie in the plane, *J. Mech. Phys. Solids* **47**, 59-79 (1999).

[50] G.H.M. van der Heijden & J.M.T. Thompson, Helical and localised buckling in twisted rods: a unified analysis of the symmetric case, *Nonlinear Dynamics* **21**, 71-99 (2000).

[51] J.H. White, An introduction to the geometry and topology of DNA structure, in *Mathematical Methods for DNA Sequences* (CRC, Boca Raton FL, 1989), pp. 225-253.
[52] T. Yabuta, N. Yoshizawa & N. Kojima, Cable kink analysis: cable loop stability under tension, *J. Appl. Mech.* **49**, 584-588 (1982).

[53] T. Yabuta, Submarine cable kink analysis, *Bull. Jap. Soc. Mech. Eng.* **27**, 1821-1828 (1984).

[54] D.W. Zachmann, Nonlinear analysis of a twisted axially loaded elastic rod, *Quart. Appl. Math.*, 67-72 (April 1979).

[55] E.E. Zajac, Stability of two planar loop elasticas, *Trans ASME, Ser. E., J. Appl. Mech.* **29**, 136-142 (1962).

[56] H. Ziegler, *Principles of Structural Stability* (Blaisdell Publishing Co., London, 1968).