Nonlinear Dynamics of Filaments IV:
Spontaneous Looping of Twisted Elastic Rods
Alain Goriely∗† and Michael Tabor∗
∗
University of Arizona, Program in Applied Mathematics,
Building #89, Tucson, AZ85721, USA
†
Universit´e Libre de Bruxelles, D´epartement de Math´ematique, CP218/1
1050 Brussels, Belgium, e-mail: agoriel@ ulb.ac.be
Submitted to Proc. Roy. Soc. A
Abstract
Everyday experience shows that twisted elastic filaments spontaneously form loops. We model
the dynamics of this looping process as a sequence of bifurcations of the solutions to the Kirchhoff
equation describing the evolution of thin elastic filament. The control parameter is taken to be
the initial twist density in a straight rod. The first bifurcation occurs when the twisted straight
rod deforms into a helix. This helix is an exact solution of the Kirchhoff equations whose stability
can be studied. The secondary bifurcation is reached when the helix itself becomes unstable
and the localization of the post-bifurcation modes is demonstrated for these solutions. Finally,
the tertiary bifurcation takes place when a loop forms at the middle of the rod and the looping
becomes ineluctable. Emphasis is put on the dynamical character of the phenomena by studying
the dispersion relation and deriving amplitude equations for the different configurations.
1
Introduction
Everybody has been faced at one point or another with the impossible tangles formed by coiled
telephone chords. The mechanism by which a telephone chord entangles with itself is a generic
phenomenon encountered in many contexts. It is one of the most basic forms of instability encountered
with elastic filaments and is usually referred to as the writhing instability, i.e., a change in spatial
configuration of the filament to reduce the overall twist of the unstable structure. Consider the
following experiment: a straight elastic filament is held between one’s fingers and twist is injected
by rotating one end while holding the other one fixed. After a small amount of twist has been
added (small rotation), the filament is no longer straight but assumes a helical form (with very small
radius). As the twist is increased the deformation of the filament tends to localize in the middle of the
rod and eventually a loop forms (see Fig. 1). Experimental studies of filament twisting (Thompson
& Champneys, 1996) have shown that this sequence of bifurcations is qualitatively correct. It is the
purpose of this paper to study and describe this looping process as a sequence of dynamical instabilities
within the framework of elasticity theory.
The starting point of almost all analyses of twisted filaments is the Kirchhoff equations. This
set of partial differential equations describes the time and space evolution of a filament subjected to
external stresses (induced by applied forces and moments at the ends). The classical analysis (Love,
1892; Timoshenko & Gere, 1961) of the stationary equations suggests that a straight, twisted elastic
filament deforms into a helix. However, these analyses are limited in many respects. First, they do
not include or explain the dynamical phenomena triggering such an instability. Second, while the
period of the helix can be easily computed by linear theory, its radius cannot be obtained. Third, it is
not known from this analysis whether the new helix is itself stable. Other analyses of the stationary
equation (Coyne, 1990) propose a completely different scenario for the looping by exhibiting a family
of stationary solutions (valid only for a restricted class of boundary conditions) deforming continuously
1
γ : twist
a
b
γ increases
c
d
e
f
Figure 1: A schematic description of the sequence of bifurcation of a looping process
2
from a straight rod to a loop. We will come back to these solutions and their relationship with our
model in the final discussion.
The analysis presented here departs radically from previous studies. We consider the dynamical
instabilities as the main mechanism triggering the changes in filament configuration. The looping
process is then explained as a sequence of bifurcations of the solutions to the Kirchhoff equations in
which the twist density is taken as a control parameter. The first bifurcation occurs when the straight
twisted filament becomes unstable. It deforms to a helix whose radius can be computed by a nonlinear
analysis (Goriely & Tabor, 1996; Goriely & Tabor, 1997a; Goriely & Tabor, 1997b). It is then argued
that this new helix is itself an exact solution to the Kirchhoff model and that it’s stability can in turn
be studied by recently developed methods (Goriely & Tabor, 1997c). The new physical parameters
associated with this helix (which are required for the associated stability analysis), such as the twist
density, can be deduced from energetic considerations.
Another view of looping, which does not utilize Kirchhoff dynamics, has been proposed by
Ricca (Ricca, 1995) and is based on a consideration of the elastic energy of various space curve
configurations.
The secondary bifurcation occurs when the helix itself becomes unstable. The unstable mode of the
helix tends to localize the deformation at one point of the rod and forms a loop. Mode localization is a
well-known experimental and theoretical phenomena for plates or filaments under stress (Tvergaard &
Needleman, 1980; Pomeau, 1981; Damil & Pottier-Ferry, 1986; Champneys & Thompson, 1997). Our
analysis provides a simple dynamical model explaining the tendency of stressed filamentary structures
to localize the deformations. Finally, the nonlinear analysis is used to compute the amplitude of the
deformation and the tertiary bifurcation is reached when the loop collapses onto itself.
As already mentioned, Thompson and Champneys have performed detailed experiments on the
buckling, localization and looping of straight twisted rod (Thompson & Champneys, 1996). They
show that the straight rod first bifurcates into a helix. Then, the helical solution further localizes
until a “small dynamic jump” occurs and a loop is formed. The analysis presented here is only
slightly different than the one experimentally observed. The discrepancies between the sequence of
events they describe in the paper and the one we use in ours are of two types: First a quantitative
difference, the first bifurcation does not occur with the same wavelength (the famous one-turn-perwave). The authors have indeed later shown that this puzzling discrepancy from the classical buckling
theory results from the influence of small initial curvature in the rod (Champneys, van der Heyden
& Thompson, 1997). Therefore this discrepancy cannot appear in our analysis (since we never take
this effect into account and work only with the simplest possible hypothesis). The second difference
come from the observed “continuous localization” of the helical rod whereas we talk of a secondary
bifurcation. One of the conclusion that we draw from our analysis is that for long enough rods, the
first and secondary bifurcation are so close that they would actually be indistinguishable. After the
secondary bifurcation, our theory also predicts continuous localization until the tertiary bifurcation
point is reached (equivalent to the “dynamic jump” of Thompson and Champneys). Therefore, we
believe that the sequence of events described in our paper is qualitatively consistent with the one
described in Thompson and Champneys.
The paper is organized as follows. In Section 2, we describe the general setup, the Kirchhoff model
and briefly review the linear and nonlinear methods that will be used. Since our proposed model is
quite elaborate and uses ideas from both linear and nonlinear stability theory we give a summary of
the proposed three step process in Section 3 and then each of these steps is discussed in detail in the
sections that follow. In Section 4 we study the primary bifurcation of the straight rod and obtain the
post-bifurcation helix. In Section 5, the linear and nonlinear analysis of the helix is performed and
we show the tendency to localization. In Section 6 we find the tertiary bifurcation where the looping
first occurs.
2
General setup
The derivation of the Kirchhoff model has been explained in detail in earlier paper (Goriely & Tabor,
1997a; Goriely & Tabor, 1997b). Here, we review the most important aspects of the theory relevant to
the problem. We first explain the kinematics of space-curves by introducing a director basis attached
to a general space-curve and then the dynamics of rods is described within the approximation of linear
elasticity theory. The section ends with a brief review of the linear and nonlinear stability methods.
3
2.1
Spin and twist
Let X = X(s, t) : R × R → R be a time dependent space curve, parameterized by the arc length s.
A director basis {d1 , d2 , d3 } can be attached to the curve as follows: The vector d3 (s, t) = X 0 (s, t) is
the tangent vector of X at s (the prime denotes the s-derivative). The vectors {d1 (s, t), d2 (s, t)} are
chosen such that {d1 , d2 , d3 } forms, a right-handed orthonormal triad (d1 × d2 = d3 , d2 × d3 = d1 ). If
d1 is along d03 , the director basis specializes to the well-know Frenet triad for which d1 is the normal
vector and d2 the Rbi-normal vector. The curve X = X(s, t) can be obtained by integrating the tangent
s
d3 (s, t)ds.
vector: X(s, t) =
The space and time evolution of the director basis specifies the kinematics of the space-curve X :
3
d0i =
3
X
Kij dj
i = 1, 2, 3,
(1.a)
Wij dj
i = 1, 2, 3,
(1.b)
j=1
d˙i =
3
X
j=1
where (˙) stands for the time derivative. The antisymmetry of W and K:




0
ω3 −ω2
0
κ3 −κ2
0
κ1  ,
0
ω1  ,
W =  −ω3
K =  −κ3
κ2 −κ1
0
ω2 −ω1
0
(2)
is a consequence of the orthonormality of the
The elements
P3of K and W make up the components
Pbasis.
3
of the twist and spin vectors; namely κ = i=1 κi di and ω = i=1 ωi di .
2.2
The Kirchhoff model
The Kirchhoff model describes the space and time evolution of thin filaments, i.e. filaments whose
length is much greater than their radius and whose curvature is sufficiently large relative to the small
length scales in the problem. In this approximation all physically relevant quantities characterizing
the three-dimensional elastic body are averaged over the cross sections attached to the central axis of
the filament. This results in a one-dimensional theory.
The total force, F = F (s, t), and the total moment, M = M (s, t), are expanded in terms of
P3
P3
the director basis; namely, F = i=1 fi di , M = i=1 Mi di . The Kirchhoff model expresses the
conservation of linear and angular momentum together with the constitutive relationship of linear
elasticity relating the moments to the strains as characterized by the twist vector. For a naturally
straight rod (i.e., no intrinsic curvature or twist) with circular cross-section the scaled equations
read (Dill, 1992; Coleman et al., 1993):
F 00 = d¨3 ,
M 0 + d3 × F = d1 × d¨1 + d2 × d¨2 ,
M = κ1 d1 + κ2 d2 + Γκ3 d3 ,
(3.a)
(3.b)
(3.c)
where Γ = 1/(1 + σ) (where σ is the Poisson ratio) measures the ratio between bending and twisting
coefficients of the rod. These equations, together with (1) can be reduced to a set of 9 equations,
second order in space and time, for the 9 unknowns (κ, ω, f ).
A constant twist γ of an elastic rod about its axis can be conveniently introduced by using the
extra degree of freedom provided by the director basis. The vectors d1 and d2 can be chosen such
that they rotate around d3 in the normal plane at a constant rate (with respect to the arc-length), γ,
independent of the space curve torsion. This axial twist is conveniently visualized as the twist of a
“ribbon” about the rod axis and, where appropriate, it is referred to as the ribbon twist. An equivalent
way of introducing the twist γ is to define a new set of variables rotating with the basis. That is,
rather than using the variables {κ, ω, f }, a new set of variables {ˆ
κ, ω
ˆ , fˆ} is obtained by introducing a
constant rotation along d3 and a translation of the vector κ. The rotated vectors are now:
4
fˆ = R.f,
(4.a)
(Rij − δj3 γ) κj
(4.b)
Rij .dj ,
(4.c)
ω
ˆ = R.ω,
κ
ˆi =
3
X
j=1
dˆi =
3
X
i = 1, 2, 3,
j=1
where

cos(γs) − sin(γs)
R =  sin(γs) cos(γs)
0
0

0
0 
1
(5)
The
spin and
vectorPare form invariant with respect to the change of variables:
P3 force,
P3 twist
3
κ
ˆ = i=1 κ
ˆ i dˆi , ω
ˆ = i=1 ω
ˆ i dˆi , Fˆ = i=1 fˆi dˆi . The transformed Kirchhoff equations read:
¨
Fˆ 00 = dˆ3 ,
(6.a)
ˆ 0 + dˆ3 × Fˆ = dˆ1 × dˆ¨1 + dˆ2 × dˆ¨2 ,
M
ˆ =κ
ˆ 2 dˆ2 + Γ(ˆ
κ3 + γ)dˆ3 ,
M
ˆ 1 dˆ1 + κ
(6.b)
(6.c)
Note also that due to the translation on the last component of the curvature vector, we now have:
ˆ3 + γ
κ
ˆ .dˆ3 = κ
The main advantage of this representation is that simple stationary filament configurations can
be conveniently expressed as constant solutions of the Kirchhoff equations. For instance, the straight
rod with constant twist, γ, and axial force, f3 , has:
ω = (0, 0, 0),
κ = (0, 0, γ) ,
f = (0, 0, f3 ) .
(7)
which transform under (4) to:
κ
ˆ=ω
ˆ = (0, 0, 0) ,
fˆ = (0, 0, f3 ) .
(8)
In the same way, a helical filament with twist γ = γH and Frenet curvature and torsion κF , τF which
takes the form (Goriely & Tabor, 1997c):
´
³τ s κ
κF
F
F
,
cos(δs),
sin(δs) , δ 2 = κ2F + τF2 ,
(9)
XH =
δ
δ
δ
with
κ = (κF sin(γs), κF cos(γs), τF + γ) ,
(10.a)
ω = (0, 0, 0) ,
(10.b)
τF
f0 )
f = (f0 sin(γs), f0 cos(γs),
κF
(10.c)
κ, fˆ), i.e.
where f0 = −τF κF + ΓκF τF , can also be written as a constant solution in terms of (ˆ
κ
ˆ = (0, κF , τF ) , ω = (0, 0, 0) ,
fˆ = (0, κF γΓ + κF τF (Γ − 1), τF (τF (Γ − 1) + Γγ)) .
We now drop the hats and consider system (6) as the main system of equations studied here.
5
(11.a)
(11.b)
2.3
Perturbation expansions
To study the stability of stationary solutions, we developed a perturbation method at the level of the
director basis (Goriely & Tabor, 1996; Goriely & Tabor, 1997a; Goriely & Tabor, 1997b; Goriely &
Tabor, 1997c) . The main idea is to expand all local quantities around the local stationary solutions
and close the system to each order in the perturbation parameter by demanding that the basis remain
orthonormal up to a given order:
(0)
di = di
(1)
(2)
+ ²2 di
+ ²di
+ ...
i = 1, 2, 3,
(12)
The orthonormality condition di .dj = δij gives rise to a system of constraints which can be solved to
each order:
(1)
di
=
3
X
(1) (0)
Aij dj ,
(13.a)
j=1
(k)
di
=
3 ³
X
(k)
(k)
Aij + Sij
´
(0)
dj ,
k > 1,
(13.b)
j=1
(j)
where S (k) is a symmetric matrix whose entries depend only on αi
antisymmetric matrix:


(k)
(k)
−α2
0
α3

(k) 
A(k) =  −α3(k)
,
0
α1
(k)
(k)
α2
−α1
0
with j < k and A(k) is the
(14)
Once the vector α(k) is known, it is an easy matter to obtain the perturbed solution by integrating
the tangent vector up to any given order:


Z s X
3
k
X
(k)
(k)
(0)
δ3i + ²A3i +
ds
²k (A3i + S3i ) di + O(²k+1 ).
(15)
X(s, t) =
i=1
j=2
P3
Any local vector V =
i=1 vi di can be expanded in terms of the perturbed basis; namely
V = V (0) + ²V (1) + ²2 V (2) + ... :
´
X ³ (1)
(0)
vi + (A(1) .v (0) )i di .
(16)
V (1) =
i
Hence we can express the first order perturbation of the twist and spin matrix, i.e. K = K (0) +
²K (1) + ..., W = W (0) + ²W (1) + ..., where
∂A(1) h (1) (0) i
+ A ,K
,
∂s
i
∂A(1) h (1)
=
+ A , W (0) .
∂t
K (1) =
(17.a)
W (1)
(17.b)
where [A, B] = A.B − B.A. Higher order perturbations can be obtained in terms of the lower order
terms.
Using these equations, one can write the k-th order perturbation of Newton’s equation (6.a)
and moment’s equation (6.b) in terms of {α(k) , f (k) }. The first of this system (k = 1), referred
to as the dynamical variational equations, controls the stability of the stationary solutions with
respect to linear time-dependent modes. Higher-order modes provide better approximations of the
perturbed solutions and are used to perform a nonlinear analysis. To emphasize the linear character
of these equations we rewrite them as a linear system of 6 equations for the 6-dimensional vector
µ(0) = {κ(0) , f (0) }, µ(k) = {α(k) , f (k) }, k > 0:
6
LE (µ(0) ).µ(1) = 0,
LE (µ(0) ).µ(k) = Hk (µ(0) , µ(1) , . . . , µ(k−1) ),
(18.a)
(18.b)
k>1
where LE is a linear differential operator in s and t whose coefficients may depend on s through
the unperturbed solution µ(0) = {κ(0) , f (0) } and Hk are vector valued functions depending on the
variables (µ(0) , . . . , µ(k−1) ).
2.4
Linear stability analysis
The linear stability of stationary solutions is determined through the use of (18.a). It has a set of
fundamental solutions which we label by the spatial mode number n, i.e.
ns
σt+i L
,
µ(1)
n = ξn e
(19)
where ξn ∈ C . The growth rate of this mode, σ = σ(n), is determined by the dispersion relation:
∆(σ, n) = det(LE ) = 0 obtained by substituting (19) into (18.a). Typically ∆ is a very complicated
expression which is best derived by symbolic manipulation. The threshold of instability is heralded
by a change in sign of (the real part of) σ and can be determined by examining the neutral curves
corresponding to the parameter values, for given n, for which σ = 0. These curves are thus solutions
of ∆(0, n) = 0. In what follows all our statements concerning critical values of twist (or tension) at
which a given configuration becomes unstable are based on their determination from the dispersion
relations.
Of course the linear analysis can only identify the initial instabilities as a function of the parameters.
As these instabilities grow (exponentially) in time the linear approximation breaks down and any
further description of the bifurcation requires a nonlinear analysis.
6
2.5
Nonlinear analysis
The techniques of nonlinear analysis enable one to derive equations for the amplitude of a solution
close to bifurcation. The distance from the bifurcation point is considered to be of the order of the
perturbation itself and this relationship is used to introduce new, longer, space and time scales over
which the solution varies. For the problems considered here, the twist γ is taken as the the control
(or “stress”) parameter and one sets
²2 = γ − γc .
(20)
where γc is the critical value at which the bifurcation occurs.
Stretched time and space scales appropriate to the problem at hand are introduced:
t0 = t,
t1 = ²t,
s0 = s,
s1 = ²s.
(21)
(22)
and taking into account the expansion in the bifurcation parameter and the new scales, one seeks
solutions of the full system (c.f. equations (18) order by order in ²:
0(²0 ) :
E(µ(0) ; s0 , t0 ) = 0
(23.a)
0(²1 ) :
(0)
LE (µ(0) ).µ(1)
(0)
LE (µ(0) ).µ(2)
(0)
LE (µ(0) ).µ(3)
(23.b)
0(²2 ) :
3
0(² ) :
=0
+
+
(1)
LE (µ(0) ).µ(1)
(1)
LE (µ(0) ).µ(2)
..
.
7
= H2 (µ(1) )
+
(1)
LE (µ(0) ).µ(1)
(23.c)
(1)
= H3 (µ
(2)
,µ
)
(23.d)
where now the µ(i) are functions of the stretched variables, i.e. µ(i) = µ(i) (s0 , s1 , ..., t0 , t1 , ..) and
P i (i)
² LE is the expansion of the linear operator in terms of the new variables (si , ti ). These
LE =
linear solutions (c.f. (19)) are written as a superposition of the neutral modes: for example, to O(²),
nc s0
µ(1) will involve terms of the form Xn (s1 , t1 )ξn ei L where Xn (s1 , t1 ) represents the slowly varying
amplitude of the unstable mode n = nc . A nonlinear amplitude equation for this slowly varying
amplitude arises as a condition (the Fredholm alternative) that the solutions µ(i) remain bounded.
Typically one finds this condition at O(²3 ) and its derivation involves, at least in the cases considered
here, massive symbolic manipulation.
In the case of a straight rod, discussed in the next section, the nonlinear amplitude equation takes
the form of a (system of) nonlinear partial differential equations and the amplitude(s) can be thought
of as the envelope of a small packet of unstable wave numbers centered around the unstable mode. In
other cases, such as the helical problem considered later, we are only interested in the amplitude of an
isolated mode and in this case, as will be shown explicitly, one derives a simpler nonlinear ordinary
differential equation.
3
Summary of looping mechanism
The sequence of bifurcations that constitutes our model of spontaneous looping are, in summary, as
follows:
• Primary Bifurcation
(a) Starting with a straight rod, with arbitrary twist, γS , and tension φ2 , linear stability analysis
is used to show that it becomes unstable (for a given φ) at a critical twist value γS = γ1 and then
bifurcates into a helix.
(b) The precise geometrical form of the helix is identified by using nonlinear analysis to determine it’s
amplitude.
(c) Using energy considerations the redistribution of the straight rod twist, γS , into twist and torsion
of the new helix is determined. The helix is then specified by its curvature κF , torsion τF and twist
γH .
• Secondary Bifurcation
(a) A linear stability analysis of the helix obtained in step (a) is carried out and the critical twist
value, γH = γ2 , at which it becomes unstable is determined.
(b) The (linear) post-bifurcation solutions are constructed and a one-loop solution, of arbitrary amplitude B, is identified.
(c) A nonlinear analysis of the unstable helical mode found in step (b) is used to determine the amplitude of the loop.
• Tertiary Bifurcation
(a) A simple criteria is developed to determine the critical B value, Bc , at which the loop will flip (
i.e. the onset of looping).
(b) The twist value, γ3 , at which this loop amplitude equals Bc is found.
The various special values of the twist, γH , γ2 , γ3 , can all be related back to the initial twist density,
γS , injected into the straight rod. Thus, in what follows, the term “control” parameter refers to the
value of γS .
8
Process
Initial straight
twisted rod
(Fig 1a)
Primary bifurcation
Straight rod → helix
(Fig 1b)
Secondary bifurcation
Helix becomes linearly
unstable (Fig 1c)
Tertiary bifurcation
Looping occurs
(Fig 1d)
Parameters
γS
φ2
γ1
κF , τF
γH
A
γ2
B
γ3
Bc
Definition of parameters
Twist density in straight rod
Tension in straight rod
Value of γS where the rod is unstable
Frenet curvature and torsion of helix
Axial twist
Amplitude of the helix
Value of γS where the helix is unstable
Amplitude of the first mode of deformation
of the helix
Value of γS where the loop collapses
Critical value of the amplitude for looping
Table 1: The different bifurcations and the new parameters introduced at each step
4
Primary bifurcation: Instability of a straight twisted rod
We consider a stationary rod of length L, simply supported, with tension T = φ2 along the x-axis and
twist γS . Rather than considering general boundary conditions, we just consider the case in which the
tangents at the ends can assume any values while the tension along the x-axis is kept constant. We
assume that these boundary conditions are maintained throughout the sequence of bifurcation and
that the unique control parameter is the twist γS . The effect of other type of boundary conditions
will be further discussed in the final section. The stationary solution is given by:
¡
¢
µ(0) = 0, 0, 0, 0, 0, φ2 .
4.1
(24)
Linear analysis
Using the techniques of linear stability analysis we are able to derive the dispersion relations for this
case (Goriely & Tabor, 1996; Goriely & Tabor, 1997b) and obtain the neutral curve from the relation
∆(0, n; γS ) = 0, where
i
¡
¢ h¡ 2 2
¢2
L γS (Γ − 1) − L2 φ2 − n2 − L2 γS 2 (Γ − 2)2 n2
(25)
∆(n; γS ) = L2 γS 2 − n2
From this one can deduce that the first instability occurs for the critical mode
nc =
φL(2 − Γ)
Γ
(26)
with critical twist value
2φ
,
(27)
Γ
This condition for the first instability is the classical buckling condition that can be found
in (Timoshenko & Gere, 1961) but written in re-scaled variables. For, γS bigger than, but close
to, γ1 , the straight rod bifurcates to the helix:
¶
µ
A
A
(28)
XH (s) = s, cos φs, sin φs .
φ
φ
γ1 =
4.2
Nonlinear analysis
The linear analysis does not specify the amplitude A of the new solution. However, a nonlinear
analysis can be performed (Goriely & Tabor, 1997b) yielding an amplitude equation consisting of a
system of nonlinear partial differential equations coupling the amplitude of the unstable mode to the
filament twist density. The stationary solution of this amplitude equation gives the amplitude A as a
function of γS , namely:
9
A2 = 2
4.3
(γS − γ1 )
.
φ
(29)
Resummation
The remarkable feature of the first bifurcation is that the new solution (28) is an exact solution of
the Kirchhoff equations. Indeed, as was shown explicitly in the first section, it is well-know that
the Kirchhoff system sustains helicoidal solutions. This property can also be seen at the level of
the perturbation expansion. Thus, if we carry on the perturbation analysis to higher-order in ² and
solve the subsequent linear system (18), it can be observed that the effect of higher-order corrections
µ(2) , µ(3) , . . . is to provide correction to the radius of the helix without changing the overall form of the
solution (that is, higher-orders do not introduce higher-order “harmonics”). Since the second order
in ² provides correction of the radius of order A2 and, as shown below, the radius A is small at the
second bifurcation, the first order approximation of the helix is sufficient for our analysis. Therefore,
to study the stability of the perturbed helicoidal solution, one can study the stability of the new exact
stationary solutions by taking the zeroth-order solution to be the helix itself with parameters given
by the nonlinear analysis of the straight rod.
In order to study the stability of the helix (28), we first re-write it in the standard form (9):
¶
µ
κF
τF s˜ κF
,
cos(δ˜
s),
sin(δ˜
s) , δ 2 = κ2F + τF2 ,
(30)
XH =
δ
δ
δ
where s˜ is the arc-length along XH (as opposed to s which is the arc length along the original straight
rod). The identification of (30) with (28) gives:
φ
Aφ
, κF =
,
2
1+A
1 + A2
φ
s
δ=√
, s˜ = √
,
1 + A2
1 + A2
τF =
(31.a)
(31.b)
(H)
The tangential force is determined from the relation f3 = φ2 cos θ, where θ is the pitch angle of
the helix (tan θ = κτFF ). Hence
φ2
(H)
.
(32)
f3 = √
1 + A2
So far, we have computed the geometrical parameters of the new helical (space) curve. However,
our elastic filament also has an imposed twist. In the deformation from straight rod to helical filament,
the twist density changes. In order to find the exact new twist density corresponding to the helical
filament, we compute the energy for both structures.
The total energy of the system is:
Z
H=
L
£
¤
ds κ21 + κ22 + Γ(κ3 + γ)2 + φ2 cos θ ,
(33)
0
where, the two first terms in the integral represent the elastic energy due to curvature effects, the
third one the elastic energy due to torsion and twist and the last one, the potential energy due to
external constraints.
For the straight rod and the helices, we have, respectively:
γ = γS ,
κS = (0, 0, 0),
κH = (0, κF , τF ), γ = γH ,
(34.a)
(34.b)
where γH is the unknown ribbon twist superimposed on the helical curve of Frenet torsion τF and
curvature κF . The ribbon twist is the actual twist of the rod: it represents the rotation in the crosssection of the director basis with respect to the Frenet basis. The energies of the straight rod and
helix are, respectively:
10
HS = L(ΓγS2 + φ2 ),
£
¤
HH = L κ2F + Γ(τF + γH )2 + φδ .
(35.a)
(35.b)
Assuming conservation of energy and equating (35.a) and (35.b), the helical ribbon twist is found to
be
r
ΓγS2 + φ(φ − δ) − κ2F
,
(36)
γH = −τF ±
Γ
where the choice of sign is yet to be determined.
There is a degeneracy in the limit where the radius of a helix shrinks to zero. Indeed, in the limit
A → 0, both torsion τF and twist γH contribute to the total twist of the straight rod γS . Formally,
a twisted straight rod can be written as a helix with zero radius, zero twist and finite torsion or a
helix with zero radius, zero torsion and finite twist or any combination of the two. Therefore, in order
to relate the ribbon twist γH to the twist γS as A → 0 we introduce the pseudo-torsion τ0 and the
residual twist γ0 of a straight rod as the limits:
τ0 = lim τF ,
(37.a)
γ0 = lim γH .
(37.b)
A→0
A→0
The twist of the straight rod, viewed as the limit of the helix with vanishing radius, is simply
γS = τ0 + γ0 , and in general, we have:
κ.d3 = τF + γH ,
(38)
for all A ≥ 0.
In the limit A → 0, the twist density goes to γ1 and the residual twist is γ0 = γ1 − τ0 = φ (2−Γ)
Γ .
The sign determination in (36) can now be found by demanding that γH → γ0 as the radius vanishes.
Thus
r
ΓγS2 + φ(φ − δ) − κ2F
(39)
γH = −τF +
Γ
We now have all the parameters of the helix as a function of the initial twist density γS . Indeed,
the radius A is a function of γS (see (29)), and the Frenet curvature and torsion are function of A
through (31.a). Finally the ribbon twist is expressed in terms of all the other parameters and γS .
Therefore, the parameters (κF , τF , γH ) describing the helical filament are know in terms of γS and
the stability of the helix as a function of the control parameter γS and the number of helical turn
N = δL
2π can now be analyzed.
5
Secondary bifurcation: Instability of the helix
An extensive study of the stability of helices has been given in (Goriely & Tabor, 1997c). We base
the following analysis on this work and specialize it to the family of helices, obtained in the previous
section, with one free parameter (the control parameter γS ). The stationary solution is written as
(c.f. equations (11)):
©
ª
µ(0) = 0, κF , τF , 0, κF γH Γ + κF τF (Γ − 1), τF γH Γ + τF2 (Γ − 1) .
(40)
where the helical ribbon twist, γH , the curvature and the torsion should be thought of in terms of it’s
relationship to the control parameter γS .
11
σ: growth rate ( x 10
-6
)
1
n: mode
1
0
2
γ = γ1
5
γ = γ2
-3
-4
4
γ > γ2
-1
-2
3
γ1 < γ < γ2
-5
-6
Figure 2: Solutions of the dispersion relation for σ 2 as a function of n with increasing values of
γS = 0.2267, 0.2675, 0.2687, 0.2695. Γ=3/4, φ = 1/10, N = 5
5.1
Linear analysis
The fundamental mode solutions to the variational equations (18.a) are of the form:
ns
σn t+iδ N
,
µ(1)
n = ξn e
(41)
where ξn ∈ C and n denotes the mode number which, due to the choice of boundary conditions, is an
integer between 1 and N . The explicit form of the linear operator LE as a function of the curvature,
torsion, twist and tension is given in (Goriely & Tabor, 1997c) and from this the dispersion relations,
∆(σ, n; γH ) = 0 can be obtained in the usual way. A typical plot of them is shown on Fig. 2 for
increasing value of γH .
The neutral curve is determined from ∆(0, n; γH ) = 0, where
6
n2
)
(42)
N2
For given n, one can then read off the the value of γH at which the instability occurs. In general,
different modes n can become unstable; however, within the family of helices parameterized by γS ,
the mode n = 1 is always the first unstable mode as the control parameter is increased. It is this case
that corresponds to our “secondary” bifurcation.
Let γ2 be the critical value at which new solutions appear, i.e. ∆(0, 1; γ2 ) = 0. This last
relation is transcendental in γS ; therefore, the exact value ofγ2 as a function of γS cannot be
obtained. Nevertheless, by expanding all the parameters to first order in (γS − γ1 ), a remarkably
good approximation of γ2 is found:
2
∆(0, n; γH ) = (Γ − 1)(Γ − 2)τF2 + 2Γ(Γ − 2)γH τF + Γ2 γH
+ δ 2 (1 −
γ2 = γ1 +
4φ
4N 2 + (3Γ + 16)N + 4
(43)
The approximation obtained by expanding the parameters to second-order in (γS −γ1 ) only slightly
improves this last result. The difference γ2 − γ1 is always very small and decreases to zero as the
12
a
b
Figure 3: The deformation of a helix due to the unstable mode n = 1, a) to first-order, b) to second
order. The parameters are: Γ = 3/4, φ = 1/10, N = 5, K = 20 and s runs from 80 to 80+100π
number of loops N increases. In the limit of an infinite long rod, γ1 → γ2 . Therefore, the delay of the
bifurcation is only due to the discretization of modes (induced by the boundary conditions).
At this point, it is of interest to reflect on the classical results of elasticity theory. The instability
of the twisted infinite (or finite) straight rod is a well-known result obtained by many authors (see
(Love, 1892; Timoshenko & Gere, 1961) for instance). However,the stability of the new helix obtained
after bifurcation has, to the best of our knowledge, never been investigated (or even questioned). The
condition (43) gives the instability threshold of the twisted helix obtained after the first bifurcation of
a straight rod. It is different from the delayed bifurcation condition obtained for a finite straight rod.
The derivation of this new condition stems from our determination of stability properties directly from
the dynamical Kirchhoff equations. However, it is likely that in real, well-controlled, experiments on
the bifurcation of rods (where the number of loops N may be large) the twisted rod could jump past
the secondary bifurcation (see Thompson and Champneys (1996) for recent experimental data) and
the helical solution might then not be observed.
Since the mode n = 1 is the first unstable mode, the effect of the instability is to localize the
solution at one point. Indeed, to first-order the explicit form of the bifurcated solution is:
2N KRν1
nδs
),
cos(
nτF
N
·
K ν2 − ν1
n−N
sin(
δs) +
x2 (s, t) = R cos(δs) −
δ n−N
N
·
n−N
K ν2 − ν1
cos(
δs) −
x3 (s, t) = R sin(δs) −
δ n−N
N
x1 (s, t) = P δs −
(44.a)
¸
ν2 + ν1
n+N
sin(
δs) ,
n+N
N
¸
ν2 + ν1
n+N
cos(
δs)
n+N
N
(44.b)
(44.c)
where K = ²Beσt and
£
¤
ν1 = n3 τF κ2F δ 5 Γ (n2 − N 2 )2 δ 4 n2 + (n2 − N 2 )2 N 2 δ 2 (σ 2 − 2κF ) + N 4 σ 2 (N 2 + n2 )
£
¤
ν2 = −N n4 δ 4 κ2F τF Γ δ 4 (n2 − N 2 )(2 − Γ) + 2N 5 n4 σ 2
(45)
To second-order the approximate solution is much closer to the exact solution since the arc-length
is conserved to order O(²4 ). The shape of the bifurcated solution is shown on Fig. 3 to first and second
order in ².
5.2
Nonlinear analysis
The unstable modes of the helix are discretized due to the boundary conditions. Therefore, a nonlinear
analysis can only take into account the temporal evolution of these discrete modes. At the first
instability, there is only one unstable mode, namely n = 1. Thus, we seek to derive an equation
describing the temporal evolution of the mode amplitude B. The corresponding stationary solution of
13
this amplitude equation will give us a relation between the amplitude and the control parameter. As
discussed in Section 2, the main idea behind a (weakly) nonlinear analysis is to look at the solution
near threshold and introduce new scales (in this case, a new time scale) proportional to the distance
to the bifurcation point and let the arbitrary amplitude, B, vary on this time scale. If the system is
close enough to the secondary bifurcation, the difference between the twist density and the critical
twist γ2 is proportional to the perturbation parameter itself, namely:
²2 = γS − γ2
(46)
The new, longer, time scale is t1 = ²t. The choice of this new scale can be justified by expanding the
dispersion relation in power of ² (see (Goriely & Tabor, 1997b)). Taking into account the possibility of
the solutions varying on these (independent) different scales, one can now solve the linear system (18).
To first order one obtain the linear solution:
δs
µ(1) = B(t1 )ξ1 ei N + c.c.
(47)
where c.c. stands for the complex conjugate and ξ1 ∈ C is specified such that ξ1 .ξ1 = 1.
The second-order solution can be found in the same way by solving (18.b) with k = 2:
6
δs
δs
µ(2) = χ0 + χ1 ei N + χ2 e2i N + c.c.,
(48)
In order to find a condition on the amplitude B(t1 ) we demand that the solutions to the thirdorder system (18.b) with k = 3 remains bounded. Thus, we apply the Fredholm alternative to the
system (18.b): here this consists of integrating H3 against all neutral solutions of the adjoint operator
L†E :
Z
L
ζ1 .H3 (µ(0) , µ(1) , µ(2) )ds = 0,
(49)
0
where ζ1 is the adjoint solution to µ1 (i.e. L†E .ζ1 = 0) . This compatibility condition gives rise to a
differential equation for B as a function of t1 :
(1)
∂2B
= B(c1 − c3 |B|2 )
∂t21
(50)
where c1 and c3 are
h
i
2
c1 = 2Γ (Γ − 1) τF 2 + 2 Γ γc (Γ − 1) τF + Γ2 γc 2 [(2 − Γ) τF − Γ γc ]
¢
¡
¢
¤
£¡
× Γ2 − 4 Γ + 3 τF 2 + −4 Γ γc + 2 Γ2 γc τF + Γ2 γc 2 + δ 2
¢
¡
¢
£¡
× 2 Γ4 + 5 − 12 Γ3 + 23 Γ2 − 18 Γ τF 4 + 46 Γ2 γc + 8 Γ4 γc − 36 γc Γ3 − 18 Γ γc τF 3
¡
¢
+ 3 δ 2 Γ2 − 8 Γ + 23 Γ2 γc 2 + 3 δ 2 + 12 Γ4 γc 2 − 36 Γ3 γc 2 + 6 + 2 Γ2 − 6 Γ δ 2 τF 2
¡
¢
+ Γγc 6 Γδ 2 − 12 Γ2 γc 2 + 8 Γ3 γc 2 − 8 + 4 Γ − 6δ 2 τF
¤−1
+2 Γ2 γc 2 + 2 Γ4 γc 4 + 3 Γ2 γc 2 δ 2 + 2 δ 2
(51)
c3 =
£ ¡ 2
¢
¡
¢
¤
2 Γ − 4 Γ + 3 τF 2 + 2 −4 Γ γc + 2 Γ2 γc τF + 2 Γ2 γ2 2 + 2 δ 2
¢
¡
¢
¤
£¡
× 2 Γ3 − 9 + 17 Γ − 10 Γ2 τF 2 + Γγ2 10 − 12 Γ + 4Γ2 τF − 3 δ 2 − 2 Γ2 γc 2 + 2 Γ3 γc 2 + Γ δ 2
i
h
3
2
× (Γ − 1) τF 3 + 3 Γ γc (Γ − 1) τF 2 + 3 Γ2 γc 2 (Γ − 1) τF + Γ3 γc 3 [(Γ − 2) τF + Γ γc ]
¢
¡
¢
£¡
× 2 Γ4 + 5 − 12 Γ3 + 23 Γ2 − 18 Γ τF 4 + 46 Γ2 γc + 8 Γ4 γc − 36 γc Γ3 − 18 Γ γc τF 3
¡
¢
+ 3 δ 2 Γ2 − 8 Γ + 23 Γ2 γc 2 + 3 δ 2 + 12 Γ4 γc 2 − 36 Γ3 γc 2 + 6 + 2 Γ2 − 6 Γ δ 2 τF 2
¡
¢
+ Γγ2 6 Γδ 2 − 12 Γ2 γc 2 + 8 Γ3 γc 2 − 8 + 4 Γ − 6δ 2 τF
¤−1
+2 Γ2 γc 2 + 2 Γ4 γc 4 + 3 Γ2 γc 2 δ 2 + 2 δ 2
(52)
14
0.285
0.28
γ
γ3
0.275
γ2
0.27
γ1
0.265
N
0.26 2
4
6
8
10
Figure 4: The primary, secondary and tertiary bifurcation as a function of the parameters for N = 2
to 10 with Γ = 3/4, φ = 1/10
where γc = γH (γ2 ), that is the critical value of the helical ribbon twist (as given by (39) with γS = γ2 )
at which the helix becomes unstable with respect to the first unstable mode n = 1.
All the parameters involve in the coefficients of the amplitude equation depends on the control
2
= c1 /c3 gives a relation between the amplitude
parameter γS . The stationary solution of (50) Bstat
and the initial twist density.
6
Tertiary bifurcation
The effect of the instability is to locate the perturbation at one point of the rod (chosen here to be
the middle point). Eventually, as the amplitude B is varied, a critical point is reached where the
projection of the solution in the x − y plane becomes multivalued. The situation is schematically
depicted on Fig. 1.d. The critical point at which the loop becomes perpendicular to the x-axis is a
new bifurcation point where the loop is about to flip over on itself. The critical value of the amplitude
Bc at which the loop centered at s = sc will flip corresponds to the condition:
x01 (sc ; Bc ) = 0,
x001 (sc ; Bc ) = 0,
(53)
This “flipping” condition is somewhat similar to the one used in Coyne’s analysis (Coyne, 1990) of
the static Kirchhoff model with the main difference that the parameter used to characterize the
Flipping in Coyne’s paper is the end displacements. The flipping condition (53) can be solved
analytically using the second order solution. However, the explicit form is too involved to be useful
and is not given here.
It is now possible to find the value of the tertiary bifurcation by finding the value of γS = γ3 such
that Bstat = Bc . A plot of the sequence of bifurcation is shown on Fig. 4 for different values of N .
Fig. 5 shows the filament as γ3 is varied.
15
a
b
γ increases
c
d
e
f
g
Figure 5: A sequence of helical filaments for varying values of γS , a) γS < γ1 , b) γ1 < γS < γ2 , c)
γ2 < γS > γ3 , d) γS = γ3 , e) γS > γ3 . The parameters are Γ = 3/4, φ = 1/10, N = 5, γ1 = 0.2667,
γ2 = 0, 2684, γ3 = 0.27376,and γS = γ2 + µ with µ × 103 =0.1, 1.5, 3, 4.9, 5.36, 5.4.
16
7
7.1
Discussion and conclusions
The dynamical picture
In the previous section, we studied the looping process as a sequence of bifurcations as the control
parameter is varied from γS < γ1 to γS > γ3 . This sequence of configurations can be obtained by
considering that the system attains its stationary state for every small change of the control parameter.
In order to obtain these configurations we rely on two assumptions: First, that in a real system, the
presence of damping will allow the system to relax in time (otherwise, none of the stationary solution
could be obtained as the system is conservative). Second, that the change in the control parameter
occurs on a much smaller time scale than the relaxation time; thereby giving rise to a sequence of
quasi-stationary configuration. However, these assumptions are not necessary for the basic phenomena
of looping. Indeed, we chose to present the sequence of bifurcations as quasi-stationary for the sake of
simplicity. The real problem of looping consists in bringing a system to a region of instability where
looping takes place, that is (in our setting) by suddenly varying the twist density of a stable straight
rod to a value γS > γ3 . Then, the dynamics of the system will directly take the system from a straight
rod to a loop without passing through the intermediary steps (this situation is actually closer to the
everyday experience of playing with telephone chords). The sequence of configurations of Fig. 5 can
then be seen as a dynamical sequence.
In this analysis, we have used the twist density as a control parameter. This is only one of the
possible choices. We could have used, equivalently, the tension (think of suddenly decreasing the
tension at the end of a straight twisted rod to produce the same dynamics) or the intrinsic twist
density. This last choice might be relevant to some problems occurring in biology where filaments
may grow with no twist but a large intrinsic twist deficit: effectively creating a twist density resulting
in the formation of a loop (for a remarkable example of this biological phenomena in the growth of
bacterial filament, see (Mendelson, 1978; Mendelson, 1990)). In any cases, the dynamics obtained by
using different parameters as control parameters are equivalent and it is likely that in any real system
changes in all parameters will actually occur simultaneously.
7.2
Why do loops forms at the middle of the rod?
One commonly observes that the loops created by twisting the ends of a rod usually form in the middle
of the rod. This is due to the choice of boundary conditions. Indeed for the sake of simplicity we have
consider that the ends are fixed in space but the tangents are free. Doing so, all the points of the
strings are actually equivalent and the loop can form at any point. However, a simple argument can
show why loops form at the middle for different boundary conditions. If the ends are held clamped
(that is the tangents at the ends are constrained along the axis), the first bifurcation does not lead
to a helix but rather (as shown in (Goriely & Tabor, 1997b)) to a helix modulated by a envelope-like
shape (so that the radius of the modulated helix goes to zero at the ends). The maximal amplitude
of this solution is reached at the middle. We have seen here that the secondary instability is triggered
by increasing the radius of the helix (controlled by the initial twist density). Therefore, the secondary
instability for a twisted rods with clamped ends will be first triggered at the middle of the string
where the amplitude of the deformation is maximal.
7.3
Conclusions
We have now completed our picture of the looping as a dynamical process: a phenomena triggered
by a sequence of instabilities for different configurations of the filament. The first bifurcation occurs
when a pulled, twisted straight rod becomes unstable. The linear analysis predicts that it deforms to
a helix. The radius of such a helix can be computed via a nonlinear analysis (Goriely & Tabor, 1996;
Goriely & Tabor, 1997b) while its twist density can be obtained by energetics consideration. The
important feature of the bifurcated solution is that it is an exact solution of the model. Therefore,
its stability can be easily studied by the same method used to study the stability of the straight rod.
The linear stability of the helix shows that it rapidly becomes unstable, given rise to a secondary
bifurcation. Moreover, the instability tends to localize the deformation of the rod at one point. The
amplitude of this localization can be obtained by performing a one-mode amplitude expansion. This
nonlinear analysis provides a relationship between the amplitude of the deformation and the control
17
parameters. The tertiary bifurcation is reached when a loop is formed at the middle and looping takes
place.
The dynamical mechanism proposed here differs in many respects from previous analysis and
exploits a series of results on the linear and nonlinear stability of filaments. These techniques are
quite general and should be applicable to other problems involving dynamical changes in form.
Acknowledgments This work is supported by DOE grant DE-FG03-93-ER25174 and Flinn
Foundation’s Biomathematics and Dynamics Initiative program (ID#048-1000206-94).
References
Champneys, A. R., & Thompson, J. M. T. 1997. A multiplicity of localised buckling modes for
twisted rod equations. Proc. Roy. Soc. London A, 452, 2467–2491.
Champneys, A. R., van der Heiden, G. H. M. & Thompson, J. M. T. 1997. Spatially complex
localisation after one-twist-per-wave equilibria in twisted rod circular rods with initial curvature.
To be published in Phil. Trans. Roy. Soc. A.
Coleman, B. D., Dill, E. H., Lembo, M., Lu, Z., & Tobias, I. 1993. On the dynamics of rods
in the theory of Kirchhoff and Clebsch. Arch. Rational Mech. Anal., 121, 339–359.
Coyne, J. 1990. Analysis of the formation and elimination of loops in twisted cable. IEE Journal of
Oceanic Engineering, 15, 72–83.
Damil, N., & Pottier-Ferry, M. 1986. Wavelength selection in the postbuckling of a long
rectangular plate. Int. J. Solids Structures, 22, 511–526.
Dill, E. H. 1992. Kirchhoff’s theory of rods. Arch. Hist. Exact. Sci., 44, 2–23.
Goriely, A., & Tabor, M. 1996. New amplitude equations for thin elastic rods. Phys. Rev. Lett.,
77, 3537–3540.
Goriely, A., & Tabor, M. 1997a. Nonlinear dynamics of filaments I: Dynamical instabilities.
Physica D, 105, 20–44.
Goriely, A., & Tabor, M. 1997c. Nonlinear dynamics of filaments II: Nonlinear Analysis. Physica
D, 105, 45–61.
Goriely, A., & Tabor, M. 1997b. Nonlinear dynamics of filaments III: Instabilities of helical rods.
Proc. Roy. Soc. London (A).
Love, A. E. H. 1892. A treatise on the mathematical theory of elasticity. Cambridge University
Press, Cambridge.
Mendelson, N. H. 1978. Helical Bacillus subtilis macrofibers: morphogenesis of a bacterial
multicellular macroorganism. Proc. Natl. Acad. Sci. USA, 75, 2472–2482.
Mendelson, N. H. 1990. Bacterial macrofibers: the morphogenesis of complex multicellular bacterial
forms. Sci. Progress Oxford, 74, 425–441.
Pomeau, Y. 1981. Nonlinear pattern selection in a problem of elasticity. J. Physique Lettres, 42,
L1–L4.
Ricca, R. L. 1995. The energy spectrum of a twisted flexible string under elastic relaxation. J. Phys.
A, 28, 2335–2352.
Thompson, J. M. T., & Champneys, A. R. 1996. From helix to localized writhing in the torsional
post-buckling of elastic rods. Proc. Roy. Soc. London A, 452, 117–138.
Timoshenko, S. P., & Gere, J. M. 1961. Theory of elastic stability. Mc Graw-Hill, New York.
Tvergaard, V., & Needleman, A. 1980. On the localization of buckling patterns. J. Appl. Mech.,
47, 613–619.
18