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

HAL
open science

# Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method

Arnaud Lazarus, J.T. Miller, P.M. Reis

► To cite this version:

Arnaud Lazarus, J.T. Miller, P.M. Reis. Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method. Journal of the Mechanics and Physics of Solids, 2013, 61 (8), pp.1712-1736. 10.1016/j.jmps.2013.04.002. hal-01447361

HAL Id: hal-01447361

https://hal.science/hal-01447361

Submitted on 24 May 2017

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.

Public Domain
# Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method

A. Lazarus<sup>a</sup>, J.T. Miller<sup>b</sup>, P.M. Reis<sup>a,b</sup>

We present a theoretical and numerical framework to compute bifurcations of equilibria and stability of slender elastic rods. The 3D kinematics of the rod is treated in a geometrically exact way by parameterizing the position of the centerline and making use of quaternions to represent the orientation of the material frame. The equilibrium equations and the stability of their solutions are derived from the mechanical energy which takes into account the contributions due to internal moments (bending and twist), external forces and torques. Our use of quaternions allows for the equilibrium equations to be written in a quadratic form and solved efficiently with an asymptotic numerical continuation method. This finite element perturbation method gives interactive access to semi-analytical equilibrium branches, in contrast with the individual solution points obtained from classical minimization or predictor-corrector techniques. By way of example, we apply our numerics to address the specific problem of a naturally curved and heavy rod under extreme twisting and perform a detailed comparison against our own precision model experiments of this system. Excellent quantitative agreement is found between experiments and simulations for the underlying 3D buckling instabilities and the characterization of the resulting complex configurations. We believe that our framework is a powerful alternative to other methods for the computation of nonlinear equilibrium 3D shapes of rods in practical scenarios.

**Keywords:**

Elastic rods

Quaternions

Path-following techniques

Equilibrium

Stability

## 1. Introduction

Filaments, rods and cables are encountered over a wide range of length-scales, both in nature and technology, providing outstanding kinematic freedom for practical applications. Given their slender geometry, they can undergo large deformations and exhibit complex mechanical behavior including buckling, snap-through and localization. A predictive understanding of the mechanics of thin rods has therefore long motivated a large body of theoretical and computational work, from Euler's elastica in 1744 (Levien, 2008) and Kirchhoff's kinetic analogy in 1859 (Dill, 1992) to the burgeoning of numerical approaches such as finite element-based methods in the late 20th century (Zienkiewicz et al., 2005), and the more recent algorithms based on discrete differential geometry (Bergou et al., 2008). Today, these advances in modeling of the mechanics of slender elastic rods are helping to tackle many cutting-edge research problems. To name just a few, these

<footnote><sup>a</sup> Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote><sup>b</sup> Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>
range from the supercoiling of DNA (Coleman and Swigon, 2000; Marko and Neukirch 2012), self-assembly of rod-coil block copolymers (Wang et al., 2012), design of nano-electromechanical resonators (Lazarus et al., 2010b, 2010a), development of stretchable electronics (Sun et al., 2006), computed animation of hairs (Bertails et al., 2006) and coiled tubing operations in the oil-gas industries (Wicks et al., 2008).

An ongoing challenge in addressing these various problems involves the capability to numerically capture their intrinsic geometric nonlinearities in a predictive and efficient way. These nonlinear kinematic effects arise from the large displacements and rotations of the slender structure, even if its material properties remain linear throughout the process (Audoly and Pomeau, 2010). As a slender elastic rod is progressively deformed, the nonlinearities of the underlying equilibrium equations become increasingly stronger leading to higher densities in the landscape of possible solutions for a particular set of control parameters. When multiple stable states coexist, classic step-by-step algorithms such as Newton-Raphson methods (Crisfield, 1991) or standard minimization techniques (Luenberger, 1973) are often inappropriate since, depending on the initial guess, they may not converge toward the desired solution, or any solution. Addressing these computational difficulties calls for alternative numerical techniques, such as well-known continuation methods (Riks, 1979; Koiter, 1970). Continuation techniques are based on coupling nonlinear algorithms (e.g. predictor-corrector (Riks, 1979) or perturbation methods (Koiter, 1970)) with an arc-length description to numerically follow the fixed points of the equilibrium equations as a function of a control parameter, that is often a mechanical or geometrical variable of the problem. With the goal of determining the complete bifurcation diagram of the system, these methods enable the computation of all of the equilibrium solution branches, as well as their local stability.

Two main approaches can be distinguished for continuing the numerical solutions of geometrically nonlinear problems. The first includes predictor-corrector methods whose principle is to follow the nonlinear solution branch in a stepwise manner, via a succession of linearizations and iterations to achieve equilibrium (Crisfield, 1991). These methods are now widely used, particularly for the numerical investigation of solutions of conservative dynamical systems, with the free path-following mathematical software AUTO being an archetypal example (Doedel, 1981). Quasi-static deformations of slender elastic rods have been intensively studied using this software (Thompson and Champneys, 1996; Furrer et al., 2000; Healey and Mehta, 2005), mostly due to the analogy between the rod's equilibrium equations with the spinning top's dynamic equations (Davies and Moon, 1993). Although popular and widely used, the main difficulty with these algorithms involves the determination of an appropriate arc-length step size, which is fixed *a priori* by the user, but may be intricately dependent on the system's nonlinearities along the bifurcation diagram. A smaller step size will favor the computation of the highly nonlinear part of the equilibrium branch, such as bifurcation points, but may also impractically increase the overall computational time. On the other hand, a larger step size may significantly compromise the accuracy and resolution of the results.

The second class of continuation algorithms, which have received less attention, is a perturbation technique called the *Asymptotic Numerical Method* (ANM), which was first introduced in the early 1990s (Damil and Potier-Ferry, 1990; Cochelin, 1994). The underlying principle is to follow a nonlinear solution branch by applying the ANM in a stepwise manner and represent the solution by a succession of local polynomial approximations. This numerical method is a combination of asymptotic expansions and finite element calculations which allows for the determination of an extended portion of a nonlinear branch at each step, by inverting a unique stiffness matrix. This continuation technique is significantly more efficient than classical predictor-corrector schemes. Moreover, by taking advantage of the analytical representation of the branch within each step, it is highly robust and can be made fully automatic. Unlike incremental-iterative techniques, the arc-length step size in ANM is adaptative since it is determined *a posteriori* by the algorithm. As a result, bifurcation diagrams can be naturally computed in an optimal number of iterations. The method has been successively applied to nonlinear elastic structures such as beams, plates and shells but the geometrical formulations were limited to the early post-buckling regime and to date, no stability analyses were performed with ANM (Cochelin et al., 1994; Zahrouni et al., 1999; Vannucci et al., 1998).

In this paper, we develop a novel implementation of the semi-analytical ANM algorithm to follow the equilibrium branches and local stability of slender elastic rods with a geometrically-exact 3D kinematics. In Section 2, we first describe the 3D kinematics where the rod is represented by the position of its centerline and a set of unit quaternions to represent the orientation of the material frame. In Section 3, we then derive the closed form of the rod's cubic nonlinear equilibrium equations. To do this, we minimize the geometrically-constrained mechanical energy including internal bending and twisting energy, as well as the work of external forces and moments. Introducing the flexural and torsional internal moments in the vector of unknowns yields differential equilibrium equations that are quadratic. In Section 4, we proceed by presenting the numerical method developed to compute the equilibrium solutions. Using a finite-difference scheme, the discretized system of equilibrium equations can be solved with the ANM algorithm, which is particularly efficient for computing our algebraic quadratic form. The local stability of the computed equilibrium branches is assessed by a second order condition on the constrained energy. Finally, we describe how to implement this numerical method in the open source software MANlab; a user-friendly, interactive and Matlab-based path-following and bifurcation analysis program (Arquier, 2007; Karkar et al., 2010). In Section 5, we develop our own precision model experiment for the fundamental problem of the writhing of a clamped, heavy and naturally curved elastic rod. Because the writhing configuration has been studied extensively in the past (Thompson and Champneys, 1996; Goriely and Tabor, 1998a; VanderHeijden and Thompson, 2000; Goyal et al., 2008), this is an ideal scenario in which to challenge our numerical model against experimental results. Our simulations are robust, computationally time-efficient and exhibit excellent quantitative agreements with our experiments, demonstrating the predictive power of our framework.
## 2. Kinematics

In this section, we present the formulation for the geometry and 3D kinematics of the slender elastic rod that we will use in our study. Assuming no shear strains and inextensibility, the mechanical deformations are represented by the rate of change of the orientation along the rod, characterized by a set of geometrically constrained unit quaternions.

### 2.1. Cosserat theory of elastic rods

An elastic rod is a slender elastic body which has a length along one spatial direction that is much larger than its dimensions in the two other perpendicular directions, that define the cross section (Fig. 1(a)). We denote the typical size of the cross-section by $h$ and the other length scale by $L$. At large scales, the rod can be regarded as an adapted material curve: its centerline. If $s$ denotes the curvilinear coordinate along the centerline of the undeformed rod, we can represent this line by a position vector function (with respect to some fixed origin) of the material point originally at $s$ in the reference configuration

$$ \mathbf{r}(s) = [r_x(s) \ r_y(s) \ r_z(s)]^T = [x(s) \ y(s) \ z(s)]^T. \qquad (1) $$

We consider unstretchable rods whose centerline remains inextensible upon deformation. As explained in detail in Audoly and Pomeau (2010), this assumption is physically justified for a wide range of loading conditions, provided that the aspect ratio of the rod, $h/L$, is small. Under this assumption, the variable $s$ is also the curvilinear coordinate along the centerline in the actual configuration. The configuration of the rod is not only characterized by the path of its centerline but also by how much it twists around this line. We consider this twist by introducing the material frame ($d_1(s) \ d_2(s) \ d_3(s)$) in the deformed configuration. At each particular location $s$, we associate an orthonormal basis $d_k$, ($k = 1, 2, 3$) attached to the centerline. The centerline, together with this set of material frames, form what is called a Cosserat curve. We choose the orientation of these material frames in a way such that the directors $d_1$ and $d_2$ lie in the plane of the cross-section, while the third director $d_3$ is always parallel to the tangent of the curve (see Fig. 1(a)). Considering the case of small strains, the triad ($d_1(s) \ d_2(s) \ d_3(s)$) remains approximately orthonormal upon deformation. This is known as the Euler-Bernoulli kinematical hypothesis (assumption of no shear deformations).

Before we are able to establish the constitutive relation, we have to quantify the rate of change of position and orientation along the rod's centerline. The rate of change in the position of the centerline is a strain vector $v(s) = [v_1(s) \ v_2(s) \ v_3(s)]^T$ that vanishes since shearing in both transverse directions and stretching are neglected. Therefore, the strains arise from the orientational rate of change of the cross-sections alone. In the framework of differential geometry of curves in 3D space (Audoly and Pomeau, 2010), this quantity is called the Darboux vector $\Omega(s)$

$$ \Omega(s) = \kappa_1(s)d_1(s) + \kappa_2(s)d_2(s) + \kappa_3(s)d_3(s). \qquad (2) $$

The physical interpretation of the Darboux vector is that the material frame rotates in the fixed frame with a rotation velocity, $\Omega(s)$, when following the centerline at unit speed. The quantities $\kappa_1$ and $\kappa_2$ in Eq. (2), called the material curvatures, illustrated in Fig. 1(b1) and (b2), represent the extent of rotation of the material frame, with respect to the directions $d_1$ and

<figure><img src="image_2.png" /><figcaption>Fig. 1. Kinematics of the Cosserat rod in the global cartesian frame $(x, y, z)$. (a) The configuration of the rod is defined by its centerline $r(s)$. The orientation of each mass point of the rod is represented by an orthonormal basis ($d_1(s) \ d_2(s) \ d_3(s)$), called the directors, where $d_3(s)$ is constrained to be tangent to $r(s)$ and (b) the three local modes of deformation of the elastic rod, associated with the change of (b1) material curvature $\kappa_1$ related to the direction $d_1$ of the cross-section, (b2) material curvature $\kappa_2$ related to $d_2$, and (b3) twist.</figcaption></figure>
$d_2$ of the cross-section. The quantity $κ_3$ quantifies the rotation of the material frame with respect to the tangent $d_3$, and is called the material twist of the rod (see Fig. 1(b3)). In order to write the material curvature and twist in an explicit form, the Darboux vector has to be rotated into the local frame. Using the condensed notation, $κ(s) = [κ_1(s) κ_2(s) κ_3(s)]^T$, and the rotation matrix of the Euclidean 3D space $R(s) \in \mathbb{R}^{3 \times 3}$, this rotated Darboux vector is

$$κ(s) = R^T(s)Ω(s), \qquad (3)$$

or, in terms of the directors $d_k(s)$,

$$κ_k(s) = d_k(s) \cdot Ω(s), \qquad (4)$$

since the directors $d_k(s)$ constitute the columns of the rotation matrix $R(s) = [d_1(s) d_2(s) d_3(s)]$.

The 3D kinematics formulation of our inextensible and unshearable elastic rod is not yet complete because we won't be able to derive the equilibrium equations directly from the material curvatures. In fact, a difficulty arises when trying to compute the infinitesimal work of the external forces using the variables $κ_1, κ_2$ and $κ_3$. A perturbation of these quantities yields a non-local perturbation to the centerline and attached material frame so that the work of the external forces cannot be written in a straightforward manner (Audoly and Pomeau, 2010). Instead, the classic approach is to choose as degrees of freedom the orientation of the material frame characterized, in this paper, by a set of quaternions. We shall now explain how to represent the rotation matrix $R(s)$ or the directors, $d_k(s)$, and the strain rate vector, $κ(s)$, in the framework of quaternions.

## 2.2. Quaternion representation

Quaternions are a number system that extends the complex number representation of geometry in a plane to the three-dimensional space (Altmann, 1986). They were first described by Hamilton in 1843 (Hamilton, 1847, 1853) and were extensively used in many physics and geometry problems before losing prominence in the late 19th century following the development of numerical analysis. Quaternions were then revived in the late 20th century, primarily due to their power and simplicity in describing spatial rotations, and have since been revived in a wide range of fields: applied mathematics (Kuipers, 1999), computer graphics (Dam et al., 1998; Hanson, 2005), optics (Horn, 1987; Tweed et al., 1990), robotics (Chou and Kamel, 1991), orbital mechanics (Arribas et al., 2006; Waldvogel, 2008) or the mechanics of slender elastic rods (Healey and Mehta, 2005; Dichmann et al., 1996; Kehrbaum and Maddocks, 1997; Balaeff et al., 2006; Spillmann and Teschner, 2009; Spillmann and Harders, 2010). It is beyond the scope of this article to discuss a detailed evaluation of the advantages and disadvantages of using quaternions over other rotation parameterizations. However, we highlight that quaternions are a non-singular representation of rotation, unlike Euler angles for instance, even if they are less intuitive than direct angles. Moreover, we favor quaternions over trigonometric approaches because of their remarkably compact quadratic polynomial form. We will show that one striking outcome of using quaternions is that the equilibrium equations we shall derive are, at most, cubic in terms of the degrees of freedom. This property is at the heart of the numerical continuation method presented in Section 4.

The fundamental relation of the algebra of quaternions, denoted by $H$, is

$$i^2 = j^2 = k^2 = ijk = -1, \qquad (5)$$

where $i, j$, and $k$ are the basis elements of $H$. A quaternion number in $H$ is written in the form $q_1 i + q_2 j + q_3 k + q_4$ where the imaginary part $q_1 i + q_2 j + q_3 k$ is an element of the vector space $\mathbb{R}^3$ and the real part $q_4$ is a scalar. Using the basis $i, j, k, 1$ of $H$ makes it possible to write a quaternion as a set of quadruples, usually expressed as a vector in $\mathbb{R}^4$

$$q = [q_1 \ q_2 \ q_3 \ q_4]^T. \qquad (6)$$

Quaternions of norm one, or unit quaternions, are a particularly convenient mathematical notation for representing orientations of objects in three dimensions. Using Euler's rotation theorem which states that a general re-orientation of a rigid-body can be accomplished by a single rotation about some fixed axis, one can represent a rotation by a set of quaternions, known as Euler parameters

$$q = [b_x \sin(\Phi/2) \ b_y \sin(\Phi/2) \ b_z \sin(\Phi/2) \cos(\Phi/2)]^T, \qquad (7)$$

where $\Phi$ is the Euler principal angle and $b = [b_x \ b_y \ b_z]^T$ is the unit length principal vector such that $b_x^2 + b_y^2 + b_z^2 = 1$ (see Fig. 2). Given that four Euler parameters are needed to define a three-dimensional rotation, a natural constraint equation

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Rotation of a rigid body using Euler's rotation theorem and a set of unit quaternions. Knowing the rotation angle $\phi$ around the unit vector $b$, we can associate a rotation matrix $R(q)$ and three directors $d_1(q), d_2(q)$ and $d_3(q)$ expressed exclusively in terms of quaternions according to Eqs. (7)–(10).</figcaption></figure>
prescribing that $\mathbf{q}$ is indeed a unit quaternion follows from Eq. (7)

$$q_1^2 + q_2^2 + q_3^2 + q_4^2 = 1. \qquad (8)$$

The orthogonal matrix representation corresponding to a rotation by the quaternion $q = q_1 i + q_2 j + q_3 k + q_4$ with $\|q\| = 1$ is

$$\mathbf{R}(\mathbf{q}) = \begin{bmatrix} q_1^2 - q_2^2 - q_3^2 + q_4^2 & 2(q_1 q_2 - q_3 q_4) & 2(q_1 q_3 + q_2 q_4) \\ 2(q_1 q_2 + q_3 q_4) & -q_1^2 + q_2^2 - q_3^2 + q_4^2 & 2(q_2 q_3 - q_1 q_4) \\ 2(q_1 q_3 - q_2 q_4) & 2(q_2 q_3 + q_1 q_4) & -q_1^2 - q_2^2 + q_3^2 + q_4^2 \end{bmatrix}. \qquad (9)$$

Returning to the context of a thin elastic rod discussed above, its material frame ($\mathbf{d}_1(s), \mathbf{d}_2(s), \mathbf{d}_3(s)$) remains orthonormal upon deformations and its rigid-body re-orientation can be expressed by the rotation matrix given in Eq. (9) such that

$$\mathbf{R}(\mathbf{q}(s)) = [\mathbf{d}_1(\mathbf{q}(s)) \ \mathbf{d}_2(\mathbf{q}(s)) \ \mathbf{d}_3(\mathbf{q}(s))]. \qquad (10)$$

The local frame is now parameterized in terms of the curvilinear unit quaternion coordinates vector $\mathbf{q}(s) = [q_1(s) \ q_2(s) \ q_3(s) \ q_4(s)]^T$ along the slender rod. Following the classical derivation given in Dichmann et al. (1996) and Kehrbaum and Maddocks (1997), the material curvatures, $\kappa_1(s)$ and $\kappa_2(s)$, and the twist, $\kappa_3(s)$, given in Eq. (3), can be derived solely in terms of the Euler parameters

$$\kappa_k(s) = 2\mathbf{B}_k \mathbf{q}(s) \mathbf{q}'(s) \quad \text{for } k = 1, 2, 3, \qquad (11)$$

where $(\cdot)'$ denotes differentiation with respect to $s$ and the skew-symmetric matrices $\mathbf{B}_k$ read

$$\mathbf{B}_1 = \begin{bmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & -1 & 0 & 0 \\ -1 & 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{B}_2 = \begin{bmatrix} 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \end{bmatrix}, \quad \mathbf{B}_3 = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -1 & 0 \end{bmatrix}. \qquad (12)$$

With the new expression of $\kappa_k(s)$, we have been able to write the total strains in terms of the locally perturbable variables $\mathbf{q}(s)$, which will be used in the derivation of the equations of equilibrium by variation of elastic energy presented below in Section 3.

It is important to note, however, that the kinematic formulation is not yet complete since the four quaternions $q_1(s), q_2(s), q_3(s)$ and $q_4(s)$ are not geometrically independent. First, to represent a three-dimensional rotation with four coordinates, the unit quaternion assumption $\|\mathbf{q}(s)\| = 1$ given in Eq. (8), must be verified. Secondly, whereas thus far we have treated the centerline position $\mathbf{r}(s)$ and the orientations $\mathbf{q}(s)$ as separate entities, the positions and the orientations cannot be considered independently. Indeed, the material frames parameterized by $\mathbf{r}(s)$ and $\mathbf{q}(s)$ are coupled by the constraint that the third director $\mathbf{d}_3(\mathbf{q}(s))$ is always parallel to the tangent $\mathbf{r}'(s)$

$$\mathbf{r}'(s) = \mathbf{d}_3(\mathbf{q}(s)), \qquad (13)$$

where $\mathbf{r}'(s)$ is the unit tangent vector to the Cosserat curve and $\|\mathbf{r}'(s)\| = 1$ along the centerline since we assumed inextensibility. The three constraints set by Eq. (13) assure that the directors are adapted to the Cosserat curve (see Fig. 1).

The three-dimensional kinematics of our inextensible and unshearable rod (including bending and twist) is represented by Eq. (11), which links the strain rates to the local orientation of the material frame, together with the four geometrical constraints given in Eqs. (8) and (13). For the remainder of this article, the three positions $r_x(s), r_y(s), r_z(s)$ and four quaternion coordinates $q_1(s), q_2(s), q_3(s), q_4(s)$ constitute the seven degrees of freedom of our slender elastic rod (Spillmann and Teschner, 2009; Spillmann and Harders, 2010). After taking into account the four constraint equations, only three of the DOFs are, in fact, geometrically independent. Their values are determined by the three-dimensional equilibrium equations, which we now address in the following section.

### 3. Mechanical equilibrium

Having formulated the kinematics of our system, we proceed by analyzing the energetics of an arbitrary configuration of the slender elastic rod. We will then derive the equations for equilibrium obtained under the assumption that this energy is stationary under small deformations for the given boundary conditions and geometrical constraints introduced above. We highlight the fact that the equilibrium equations are highly nonlinear due to geometry, rather than the material response.

#### 3.1. Energy formulation

For simplicity, and to avoid loss of generality, we shall adopt the framework of Hookean elasticity and consider linear isotropic constitutive laws. For practical purposes, this hypothesis is usually appropriate since, for slender elastic rods, the strains at the material level are typically small. Under this assumption, the total elastic energy of the slender elastic rod can be written as the uncoupled sum of bending and twisting contributions (Audoly and Pomeau, 2010). Although the reference configuration of the rod is assumed to be stress-free, we can readily account for rods with intrinsic natural curvature and
twist. Doing so, the elastic energy of a rod with length L and a constant cross-section reads

$$ \mathcal{E}_e = \frac{EI_1}{2} \int_0^L (\kappa_1(s) - \hat{\kappa}_1(s))^2 ds + \frac{EI_2}{2} \int_0^L (\kappa_2(s) - \hat{\kappa}_2(s))^2 ds + \frac{GJ}{2} \int_0^L (\kappa_3(s) - \hat{\kappa}_3(s))^2 ds, \quad (14) $$

where we used the previously defined rotational strain rate vector $\kappa(s) = [\kappa_1(s) \ \kappa_2(s) \ \kappa_3(s)]^T$, and where the quantities $\hat{\kappa}_1(s)$, $\hat{\kappa}_2(s)$ and $\hat{\kappa}_3(s)$ are the intrinsic natural curvature and twist of the rod along the directors $\mathbf{d}_1$, $\mathbf{d}_2$ and $\mathbf{d}_3$, respectively. In this expression, E is the Young's modulus of the material and $G = E/2(1 + \nu)$ is the shear modulus of the material with Poisson's ratio $\nu$. The constants $I_1$ and $I_2$ are the moments of inertia along the principal directions of curvature in the plane of the cross-section $\mathbf{d}_1$ and $\mathbf{d}_2$ and J is the moment of twist which, similarly to $I_1$ and $I_2$ for the bending energy, depends only of the geometry of the cross-section. Replacing the material curvatures $\kappa_1(s)$, $\kappa_2(s)$ and twist $\kappa_3(s)$ by their expression given in Eq. (11) allows us to write the elastic energy $\mathcal{E}_e$ in a more compact form, in term of the rotational degrees of freedom $\mathbf{q}(s)$ alone

$$ \mathcal{E}_e(\mathbf{q}(s)) = \sum_{k=1}^{3} \frac{E_k I_k}{2} \int_{0}^{L} (2\mathbf{B}_k \mathbf{q}(s) \mathbf{q}'(s) - \hat{\kappa}_k(s))^2 ds, \quad (15) $$

where $E_1 = E_2 = E$, $E_3 = G$ and $I_3 = J$.

## 3.2. Variation of the energy

We now follow a variational approach for the elastic energy in Eq. (15), and consider an infinitesimal perturbation from an arbitrary configuration of the rod. The perturbed quantities are preceded by $\delta$. Carrying out the first variation of Eq. (15), the corresponding variation of the energy $\mathcal{E}_e$ is

$$ \delta \mathcal{E}_e = \sum_{k=1}^{3} E_k I_k \int_{0}^{L} (2\mathbf{B}_k \mathbf{q}\mathbf{q}' - \hat{\kappa}_k)(2\mathbf{B}_k \delta \mathbf{q}\mathbf{q}' + 2\mathbf{B}_k \mathbf{q}\delta \mathbf{q}') ds, \quad (16) $$

where $\delta \mathbf{q} = [\delta q_1 \ \delta q_2 \ \delta q_3 \ \delta q_4]^T$ is the vector of the arbitrary perturbations of the rotational degrees of freedom $\mathbf{q}$. Upon integration by parts, we transform Eq. (16) into an integral that depends on $\delta \mathbf{q}$ alone to arrive at

$$ \delta \mathcal{E}_e = \left[ \sum_{k=1}^{3} E_k I_k (2\mathbf{B}_k \mathbf{q}\mathbf{q}' - \hat{\kappa}_k) 2\mathbf{B}_k \delta \mathbf{q} \right]_0^L - \int_0^L \sum_{k=1}^{3} E_k I_k [(2\mathbf{B}_k \mathbf{q}\mathbf{q}' - \hat{\kappa}_k) 4\mathbf{B}_k \mathbf{q}' + (2\mathbf{B}_k \mathbf{q}\mathbf{q}' + 2\mathbf{B}_k \mathbf{q}\mathbf{q}' - \hat{\kappa}_k') 2\mathbf{B}_k \mathbf{q}] \delta \mathbf{q} ds, \quad (17) $$

where the first term stands for the variation of elastic energy over the entire interval and is the boundary term from the integration by parts assuming that the rod is parameterized from $s=0$ to $s=L$. Physically, this first term represents the work done by the operator upon a change of orientation applied to the ends of the rod. We can rewrite this term in the concise form $[\mathbf{T}(s)\delta\mathbf{q}]$ where

$$ \mathbf{T}(s) = G_1(s)2\mathbf{B}_1\mathbf{q} + G_2(s)2\mathbf{B}_2\mathbf{q} + G_3(s)2\mathbf{B}_3\mathbf{q}. \quad (18) $$

The vector $\mathbf{T}(s) = [T_1(s) \ T_2(s) \ T_3(s) \ T_4(s)]^T$ is the internal moment projected in the quaternion basis defined as a linear superposition of the internal moments due to elementary modes of deformation. The functionals $G_1(s)$, $G_2(s)$ and $G_3(s)$ given by

$$ G_k(s) = E_k I_k (\kappa_k(s) - \hat{\kappa}_k(s)) = E_k I_k (2\mathbf{B}_k \mathbf{q}(s) \mathbf{q}'(s) - \hat{\kappa}_k(s)) \quad (19) $$

are respectively the two flexural and torsional moments, defined as the components of $\mathbf{T}(s)$ in the local material frame. The second term in Eq. (17) is the work done by the operator upon a change of orientation applied along the rod. The elementary contribution to the integral can be rewritten $\int_0^L \tau(s)\delta\mathbf{q}ds$ where $\tau(s) = [\tau_1(s) \ \tau_2(s) \ \tau_3(s) \ \tau_4(s)]^T$ as a four-dimensional vector written in the quaternion basis that reads

$$ \tau(s) = \sum_{k=1}^{3} (G_k(s)4\mathbf{B}_k \mathbf{q}'(s) + G'_k(s)2\mathbf{B}_k \mathbf{q}(s)), \quad (20) $$

where $G'_k(s) = E_k I_k (2\mathbf{B}_k \mathbf{q}'(s) \mathbf{q}'(s) + 2\mathbf{B}_k \mathbf{q}(s) \mathbf{q}''(s) - \hat{\kappa}_k(s))$ is the differential of $G_k(s)$ with respect to $s$. The quantity $\tau(s) ds$ is the net moment applied on an infinitesimal element of the rod located between the cross-sections at $s$ and $s+ds$.

Before arriving to the equilibrium equations from this variation, we need to consider the external loads that are applied to the rod, and whose work must balance the variation of energy at equilibrium. Here, we consider two types of external loads: point forces ($\mathbf{P}(0)$, $\mathbf{P}(L)$) and torques ($\mathbf{M}(0)$, $\mathbf{M}(L)$) that are applied at the two ends $s=0$ and $s=L$, and distributed forces and torques that are applied along the length of the rod, with linear densities $\mathbf{p}(s)$ and $\mathbf{m}(s)$, respectively. The density of forces, $\mathbf{p}(s)$, can represent, for instance, the weight of the rod, and the density of moments, $\mathbf{m}(s)$, hydrostatic loadings such as the result of viscous stresses due to a swirling flow around the rod. The total work done by these external forces upon an infinitesimal perturbation of the rod's configuration is

$$ \delta W = \mathbf{P}(0)\delta\mathbf{r}(0) + \mathbf{M}(0)\delta\mathbf{q}(0) + \mathbf{P}(L)\delta\mathbf{r}(L) + \mathbf{M}(L)\delta\mathbf{q}(L) + \int_{0}^{L} (\mathbf{p}(s)\delta\mathbf{r}(s) + \mathbf{m}(s)\delta\mathbf{q}(s)) ds, \quad (21) $$

where $\delta\mathbf{r} = [\delta r_1 \ \delta r_2 \ \delta r_3]^T$ is the vector of the small arbitrary perturbations of the translational degrees of freedom $\mathbf{r}$. According to Eq. (21), the external forces $\mathbf{P}(s) = [P_x(s) \ P_y(s) \ P_z(s)]^T$ and $\mathbf{p}(s) = [p_x(s) \ p_y(s) \ p_z(s)]^T$ are defined in terms of the
global directions x, y, z whereas the external moments **M**(s) = [M₁(s) M₂(s) M₃(s) M₄(s)]ᵀ and **m**(s) = [m₁(s) m₂(s) m₃(s) m₄(s)]ᵀ are expressed in the quaternion basis i, j, k, 1 of H, defined by Eq. (5). In Section 5, we will show through the specific example of the writhing of a rod how to express physical rotational quantities (e.g. boundary conditions or external moments) in terms of quaternions.

### 3.3. Equilibrium equations

Thus far, we have implicitly assumed that the perturbations [δ**r** δ**q**]ᵀ can be chosen freely. This is, however, not the case since our rod is subject to the kinematical constraints introduced previously in Eqs. (8) and (13). These constraints are imposed in the derivation of the equations of equilibrium by adding a number of Lagrange multipliers into the variation of the elastic energy εe and external loads δW. In this Lagrangian formalism, the enforcement of the unicity of quaternion in Eq. (8), translates as the continuous functional constraint

$$C_{\alpha}[\mathbf{q}(s)] = \mathbf{q}(s)\mathbf{q}(s)-1 = 0, \qquad (22)$$

where the brackets indicate that C_α depends on the function q(s), globally. Moreover, Eq. (13), which ensures that the directors are adapted to the Cosserat curve, translates to three conditions on the continuous vector-valued function

$$C_{\mu}[\mathbf{r}(s), \mathbf{q}(s)] = \mathbf{r}'(s) - \mathbf{d}_3(\mathbf{q}(s)) = \mathbf{0}. \qquad (23)$$

With the expressions for the energy of an arbitrary configuration of the rod in Eqs. (17) and (21) in hand, the equations of equilibrium are now obtained by assuming that the energy is stationary under small deformations for a given set of boundary conditions and geometrical constraints; Eqs. (8) and (13). This is equivalent to requiring that the first order variation of the functionals δεe and δW, combined linearly with the variation of the constraints δC_α and δC_μ over the interval from s=0 to s=L (i.e. the Lagrangian) vanish

$$\delta \mathcal{E}_{\text{Lag}} = \delta \mathcal{E}_e - \delta W + \int_0^L \alpha(s) \delta C_\alpha(s) \, ds + \int_0^L \mu(s) \delta C_\mu(s) \, ds = 0. \qquad (24)$$

In this equation, the variation of the constraint C_α(s) given in Eq. (22) takes the form

$$\int_{0}^{L} \alpha(s) \delta C_{\alpha}(s) \, ds = \int_{0}^{L} 2\alpha q \delta q \, ds, \qquad (25)$$

where the scalar function α(s) is the Lagrange multiplier that imposes the norm of the quaternions to be one. The variation of the constraints C_μ(s) given in Eq. (23) reads, after integration by parts

$$\int_{0}^{L} \mu(s) \delta C_{\mu}(s) \, ds = [\mu \delta \mathbf{r}]_{0}^{L} - \int_{0}^{L} \mu' \delta \mathbf{r} + 2 \mathbf{D}(\mathbf{q}) \mu \delta \mathbf{q} \, ds, \qquad (26)$$

where the terms of the vector valued function μ(s) = [μₓ(s) μᵧ(s) μ_z(s)]ᵀ are the Lagrange multipliers ensuring the condition of inextensibility of the slender elastic rods and the operator D(q) reads

$$\mathbf{D}(\mathbf{q}(s)) = \begin{bmatrix} q_3(s) & -q_4(s) & -q_1(s) \\ q_4(s) & q_3(s) & -q_2(s) \\ q_1(s) & q_2(s) & -q_3(s) \\ q_2(s) & -q_1(s) & -q_4(s) \end{bmatrix}. \qquad (27)$$

Now, substituting Eqs. (17), (21), (25) and (26) into the Lagrangian of Eq. (24), we arrive at the first variation of the geometrically constraint elastic energy of the slender elastic rod

$$\begin{aligned} \delta \mathcal{E}_{\text{Lag}} = & [\mathbf{T}(s)\delta\mathbf{q}(s) + \mu(s)\delta\mathbf{r}(s)]_{0}^{L} - \mathbf{M}(0)\delta\mathbf{q}(0) - \mathbf{M}(L)\delta\mathbf{q}(L) \\ & - \mathbf{P}(0)\delta\mathbf{r}(0) - \mathbf{P}(L)\delta\mathbf{r}(L) - \int_{0}^{L} (\mathbf{p}(s) + \mu'(s))\delta\mathbf{r}(s) \, ds \\ & - \int_{0}^{L} (\tau(s) + \mathbf{m}(s) - 2\alpha(s)\mathbf{q}(s) + 2\mathbf{D}(\mathbf{q}(s))\mu(s))\delta\mathbf{q}(s) \, ds. \end{aligned} \qquad (28)$$

The condition that the variation in Eq. (28) must vanish for an arbitrary perturbations δ**r**(s) and δ**q**(s) yields the strong form of the equilibrium equations for our elastic rod as second-order differential equations

$$\mathbf{0} = \mathbf{p}(s) + \mu'(s) \qquad (29a)$$

$$\mathbf{0} = \tau(s) - \mathbf{m}(s) + 2\alpha(s)\mathbf{q}(s) - 2\mathbf{D}(\mathbf{q}(s))\mu(s). \qquad (29b)$$

When projected along the three directions of the global cartesian frame (x, y, z), the vector equation Eq. (29a) yields a set of three differential equations that can be interpreted as the balance of forces. The vector of Lagrange multiplier μ(s) measures the resultant of the contact forces transmitted through the rod's cross-section. Indeed, calculating the forces acting on a small element of the rod of length ds, we find that the element is submitted to the contact forces μ(s + ds) and −μ(s) from the
neighboring elements, and to the external force $\mathbf{p}$ ds. At equilibrium, the total forces ($\mu'$ ds + $\mathbf{p}$ ds) is zero as described by Eq. (29a).

When projected along the four elements (i, j, k, 1) of the quaternion basis $\mathcal{H}$, the vector equation Eq. (29b) yields a set of four differential equations that can be interpreted as the balance of moments. Working in the quaternion basis, it is, however, not straightforward to find an obvious physical interpretation for each of the terms but it suffices to say that they are related to the internal moments acting on a small element of the rod ds.

For the equilibrium equations written in Eq. (29) to be complete and well-posed, one must add the geometrical constraints given by Eqs. (22) and Eq. (23), which, in their projected and developed form, read as

$$0 = q_1^2(s) + q_2^2(s) + q_3^2(s) + q_4^2(s) - 1 \qquad (30a)$$

$$0 = r'_{x}(s) - 2q_{1}(s)q_{3}(s) - 2q_{2}(s)q_{4}(s) \qquad (30b)$$

$$0 = r'_{y}(s) - 2q_{2}(s)q_{3}(s) + 2q_{1}(s)q_{4}(s) \qquad (30c)$$

$$0 = r'_{z}(s) + q_{1}^{2}(s) + q_{2}^{2}(s) - q_{3}^{2}(s) - q_{4}^{2}(s). \qquad (30d)$$

In the seven differential equilibrium equations of (29) plus the four differential equations in Eq. (30), the eleven unknowns are the four Lagrange multipliers $\alpha(s)$, $\mu_x(s)$, $\mu_y(s)$ and $\mu_z(s)$, the four rotational degrees of freedom $q_1(s)$, $q_2(s)$, $q_3(s)$ and $q_4(s)$ and the three translational degrees of freedom $r_x(s)$, $r_y(s)$ and $r_z(s)$. Thanks to the use of quaternions, the kinematics is geometrically-exact and the resultant equilibrium equations are simply polynomial since the highest geometric nonlinearity comes from the vector $\tau(s)$ given in Eq. (18), which is cubic in $\mathbf{q}(s)$. In Section 4, while developing the numerical implementation, we will make extensive use of this smooth and regular nonlinearity to efficiently compute the numerical solutions of these equations.

So far, in Eq. (28), we have only considered the vanishing of the integral term. Likewise, boundary terms should also vanish since this equation is also to be satisfied for perturbations localized at its extremities. The first boundary terms, associated with rotations $\delta\mathbf{q}(L)$ and $\delta\mathbf{r}(0)$ yield

$$0 = (\mathbf{T}(0) + \mathbf{M}(0))\delta\mathbf{q}(0), \qquad (31a)$$

$$0 = (\mathbf{T}(L) - \mathbf{M}(L))\delta\mathbf{q}(L). \qquad (31b)$$

The remaining boundary terms associated with displacements $\delta\mathbf{r}(0)$ and $\delta\mathbf{r}(L)$ of the ends $s=0$ and $s=L$, respectively, yield

$$0 = (\mu(0) + \mathbf{P}(0))\delta\mathbf{r}(0), \qquad (31c)$$

$$0 = (\mu(L) - \mathbf{P}(L))\delta\mathbf{r}(L). \qquad (31d)$$

To provide a physical interpretation of the behavior at the boundary conditions, we first consider Eq. (31a). If the endpoint $s=0$ is free to rotate, the vector $\delta\mathbf{q}(0)$ is arbitrary and one is led to the boundary condition $\mathbf{T}(0) + \mathbf{M}(0) = \mathbf{0}$. This is the total torque applied on the section $s=0$, which is the sum of the internal moments $\mathbf{T}(0)$ transmitted by the downstream part of the rod, $s > 0$, and of the moment $\mathbf{M}(0)$ applied by the operator. At equilibrium, the total torque should vanish when the end is free to rotate. If the endpoint $s=0$ is fixed, the perturbations that are consistent with the kinematics are such that $\delta\mathbf{q}(0) = \mathbf{0}$ and the equation is automatically satisfied. The boundary condition is then the one imposing the rotation of the fixed end, which leaves the total number of boundary conditions unchanged. The same reasoning holds for Eq. (31b) near the opposite end, $s=L$, although the total torque is now $\mathbf{T}(L) - \mathbf{M}(L)$, since, in this case the internal moment is applied by the downstream part of the rod, $s < L$.

The two other boundary conditions written in Eqs. (31c) and (31d) can be handled in a similar fashion. Near an end where the displacement is unconstrained, the total force should be zero. This total force is $\mu(0) + \mathbf{P}(0)$ near the end $s=0$, by a similar reasoning as above. However, the total force is $\mu(L) - \mathbf{P}(L)$ near the opposite end, $s=L$, given that the internal forces $\mu(L)$ are now applied by the downstream part of the rod, $s < L$. This remark validates our previous interpretation as for the physical interpretation of the Lagrange multipliers $\mu_x(s)$, $\mu_y(s)$ and $\mu_z(s)$; they are the internal forces along the three directions of the global frame that constrain the directors to be adapted to the Cosserat curve.

Together, Eqs. (29)–(31) constitute the set of geometrically-exact cubic differential equations that describe the mechanical behavior of the slender elastic rod represented in Fig. 1(a). These nonlinear differential equations could be solved with classic boundary value problem algorithms upon knowing the boundary conditions in terms of external forces or kinematics. Moreover, coupled with traditional predictor-corrector methods, one should be able to continue, step-by-step, the solutions of this nonlinear elastic problem in terms of given geometric or mechanical control parameters (Crisfield, 1991; Doedel, 1981). In the following section, in an alternative point of departure, we use a continuation method based on the Asymptotic Numerical Method (ANM) developed in the early 1990's to solve elastic structural problems in the early post-buckled regime (Damil and Potier-Ferry, 1990; Cochelin, 1994). Taking advantage of the particular cubic form of the geometrically-exact equilibrium Eqs. (29)–(31), this path-following perturbation technique will enable the determination of semi-analytical nonlinear solution branches by inverting a simple stiffness matrix at each step of the continuation. This outstanding numerical property makes the ANM algorithm highly robust and computationally efficient at determining the various equilibria of our slender elastic rod.
## 4. Numerical method

In this section, we solve the differential equilibrium equations Eqs. (29)–(31) using a finite element based semi-analytical path-following method. We first approximate the continuous degrees of freedom using finite differences approximation to interpolate the mechanical and geometrical variables at each nodes and elements. Thanks to the quaternion formalism introduced above, the equilibrium equations can be reduced to an algebraic set of quadratic equations by considering the flexural and torsional internal moments as unknowns. This quadratic form is particularly well suited to ANM which is a semi-analytical continuation algorithm to compute the branches of solution of a set of nonlinear polynomial equations. To follow all the bifurcated branches, we show how the local stability of the computed equilibria can be assessed by using the second order conditions of constrained minimization problems. Finally, we describe the implementation of our algorithm into MANlab, a free and interactive bifurcation analysis software based in MATLAB.

### 4.1. Discretization

In order to compute the equilibrium equations, Eqs. (29)–(31), we first explain how to discretize the main function unknowns such as strain rate vector $κ(s)$, material frame ($d_1(s)$, $d_2(s)$, $d_3(s)$), positional and rotational degrees of freedom $r(s)$ and $q(s)$ or Lagrange multipliers $α(s)$ and $μ(s)$.

The position of the rod is represented by discretizing its centerline into $N$ elements separated by $N + 1$ spatial control points, $r'(s^i) = r^i = [r_x^i r_y^i r_z^i]^T$ in $R^3$, located by the discrete curvilinear coordinate $s^i$ as illustrated in Fig. 3(a). The spatial derivative of the positional degrees of freedom is approximated by the forward finite difference between two successive nodes

$$ \mathbf{r}'(s^i) \approx \frac{\mathbf{r}(s^i + ds^i) - \mathbf{r}(s^i)}{ds^i} = \frac{\mathbf{r}^{i+1} - \mathbf{r}^i}{ds^i} \quad (32) $$

where $ds^i = \|\mathbf{r}^{i+1} - \mathbf{r}^i\| = L/N$ is the length of the ith element and L is the total length of the rod. To ensure the inextensibility condition of our rods, $ds^i$ is constant upon deformation and the stretch along the centerline is forced to verify $\|\mathbf{r}'(s^i)\| = 1$.

The orientations of the centerline elements are represented by employing $N$ material frames $R(\mathbf{q}^i) = [\mathbf{d}_1^i \mathbf{d}_2^i \mathbf{d}_3^i]$ in $R^{3 \times 3}$ where $\mathbf{q}^i = [q_1^i q_2^i q_3^i q_4^i]^T$ is the set of quaternions associated with each element j. According to Eq. (9), the directors of the jth element are vectors in $R^3$ represented at the midpoints on the centerline segments (see Fig. 3(a)) such that

$$ \mathbf{d}_1^i = \begin{bmatrix} q_1^i q_1^i - q_2^i q_2^i - q_3^i q_3^i + q_4^i q_4^i \\ 2(q_1^i q_2^i + q_3^i q_4^i) \\ 2(q_1^i q_3^i - q_2^i q_4^i) \end{bmatrix}, \quad \mathbf{d}_2^i = \begin{bmatrix} 2(q_1^i q_2^i - q_3^i q_4^i) \\ -q_1^i q_1^i + q_2^i q_2^i - q_3^i q_3^i + q_4^i q_4^i \\ 2(q_2^i q_3^i + q_1^i q_4^i) \end{bmatrix}, \quad \mathbf{d}_3^i = \begin{bmatrix} 2(q_1^i q_3^i + q_2^i q_4^i) \\ 2(q_2^i q_3^i - q_1^i q_4^i) \\ -q_1^i q_1^i - q_2^i q_2^i + q_3^i q_3^i + q_4^i q_4^i \end{bmatrix}. \quad (33) $$

Replacing the quaternions functions $q(s)$ by their discrete counterparts $\mathbf{q}^i$ in the expression of strain rates given in Eq. (11), we can write the discrete material curvatures $κ_1^j, κ_2^j$ and the twist $κ_3^j$ expressing the extent of rotation around the directors, $\mathbf{d}_1^i, \mathbf{d}_2^i$ and $\mathbf{d}_3^i$, between two successive elements (see Fig. 3(b)) in the form

$$ κ_k^j = 2B_k \bar{q}^j \bar{q}^j \quad \text{for } k = 1, 2, 3. \quad (34) $$

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Finite-difference method. (a) The centerline of the rod is discretized into $N$ elements separated by $N+1$ nodes, $r^i$. We also consider $N$ material frames $R(\mathbf{q}^i)$ to express the orientation of the jth element and (b) the change of orientation between two successive elements $j$ and $j+1$ is expressed by the $N-1$ discrete material curvatures at the interconnected nodes: (b1) $κ_1^j$ around the director ($\mathbf{d}_1^j + \mathbf{d}_1^{j+1})/2$, (b2) $κ_2^j$ around the director ($\mathbf{d}_2^j + \mathbf{d}_2^{j+1})/2$, and (b3) $κ_3^j$ around the director ($\mathbf{d}_3^j + \mathbf{d}_3^{j+1})/2$.</figcaption></figure>
In Eq. (34), we introduced the average and the spatial derivative of the rotational degrees of freedom of the Pjth element $\mathbf{q}^j$ as

$$ \bar{\mathbf{q}}^j = \frac{\mathbf{q}^{j+1} + \mathbf{q}^j}{2} \quad \text{and} \quad \mathbf{q}^j = \frac{\mathbf{q}^{j+1} - \mathbf{q}^j}{ds^j}, \qquad (35) $$

where $ds^j = 1/2 (\|\mathbf{r}^{j+2} - \mathbf{r}^{j+1}\| + \|\mathbf{r}^{j+2} - \mathbf{r}^{j+1}\|) = L/N$, taking into consideration the rod's inextensibility condition.

In a similar fashion, replacing the continuous function $\mathbf{q}(s)$ and its derivative $\mathbf{q}'(s)$ by their discretized counterparts, $\bar{\mathbf{q}}^j$ and $\mathbf{q}^j$, respectively, given in Eq. (35), the first variation of elastic energy previously given in Eq. (16) can be approximated by a Riemann sum over the elements from $j=1$ to $N$

$$ \delta\mathcal{E}_e \approx \sum_{k=1}^{3} \sum_{j=1}^{N-1} 2G_k^j(\mathbf{B}_k\mathbf{q}^j\delta\mathbf{q}^{j+1} - \mathbf{B}_k\mathbf{q}^{j+1}\delta\mathbf{q}^j) ds^j. \qquad (36) $$

In this equation, the $N$ vectors $\delta\mathbf{q}^j = [\delta\bar{\mathbf{q}}_1^j \delta\bar{\mathbf{q}}_2^j \delta\bar{\mathbf{q}}_3^j \delta\bar{\mathbf{q}}_4^j]^T$ are the discrete version of the perturbed rotational degrees of freedom $\delta\mathbf{q}(s)$ and are associated with each element $j$. The $3(N-1)$ constants $G_k^j$ are an approximation of the flexural and torsional internal torques $G_k(s)$ given in Eq. (19), which read

$$ G_k^j = E_k I_k \left( \mathbf{B}_k(\mathbf{q}^j + \mathbf{q}^{j+1}) \frac{1}{ds^j} (\mathbf{q}^{j+1} - \mathbf{q}^j) - \bar{\mathcal{C}}_k^j \right), \qquad (37) $$

and are defined between two successive elements for $1 \le j \le N$.

Replacing $\mathbf{q}(s)$ and $\mathbf{r}(s)$ by their discrete counterparts at each element $\mathbf{q}^j$ and each node $\mathbf{r}^j$, respectively, the variation of the work done by the external forces and torques previously given in Eq. (21), can be approximate by its discrete version

$$ \delta W \approx \mathbf{P}^0 \delta \mathbf{r}^1 + \mathbf{M}^0 \delta \mathbf{q}^1 + \mathbf{P}^L \delta \mathbf{r}^{N+1} + \mathbf{M}^L \delta \mathbf{q}^N + \sum_{i=2}^{N} \mathbf{p}^i \delta \mathbf{r}^i ds^i + \sum_{j=2}^{N-1} \mathbf{m}^j \delta \mathbf{q}^j ds^j, \qquad (38) $$

where $\delta \mathbf{r}^i = [\delta \bar{\mathbf{r}}_x^i \delta \bar{\mathbf{r}}_y^i \delta \bar{\mathbf{r}}_z^i]^T$ is the vector of the perturbed displacement at each node. In this equation, we have introduced the vector of point forces at the two ends $\mathbf{P}^0 = [P_x^0 P_y^0 P_z^0]^T$ and $\mathbf{P}^L = [P_x^L P_y^L P_z^L]^T$ and the vector of density of external forces at each node $\mathbf{p}^i = [p_x^i p_y^i p_z^i]^T$ for $2 \le i \le N$ defined in term of the global directions x, y, z. In the discrete version of $\delta W$, Eq. (38), we have also introduced the torques applied at the two ends $\mathbf{M}^0 = [M_1^0 M_2^0 M_3^0 M_4^0]^T$ and $\mathbf{M}^L = [M_1^L M_2^L M_3^L M_4^L]^T$ and the vector of density of external moment at each element $\mathbf{m}^i = [m_1^i m_2^i m_3^i m_4^i]^T$ for $2 \le j \le N-1$, also expressed in the quaternion basis.

Before we derive the algebraic system of equilibrium equations, we still need to write the discrete form of the variation of work due to the geometrical constraints in Eqs. (22) and (23). Replacing $\mathbf{q}(s)$ by its discrete counterpart, $\mathbf{q}^j$, we can expand Eq. (25) in the form of a Riemann sum

$$ \int_{0}^{L} \alpha(s) \delta C_{\alpha}(s) ds \approx \sum_{j=1}^{N} 2 \alpha^{j} \mathbf{q}^{j} \delta \mathbf{q}^{j} ds^{j}, \qquad (39) $$

where $\alpha^j$ is a discrete scalar at each element $j$, approximating the continuous Lagrange parameter given in Eq. (25). Introducing the $N$ vectors of Lagrange parameters $\mu^j = [\mu_x^j \mu_y^j \mu_z^j]^T$ which prescribe that each element is parallel to the tangent $\mathbf{r}'(s^i)$ given in (32), we can rewrite Eq. (26) in its discrete form

$$ \int_{0}^{L} \mu(s) \delta C_{\mu}(s) ds \approx \mu^{N} \delta \mathbf{r}^{N+1} - \mu^{1} \delta \mathbf{r}^{1} - \sum_{i=1}^{N-1} (\mu^{i+1} - \mu^{i}) \delta \mathbf{r}^{i+1} + \sum_{j=1}^{N} \mathbf{D}(\mathbf{q}^{j}) \mu^{j} \delta \mathbf{q}^{j} ds^{j}, \qquad (40) $$

where the vector of Lagrange parameters, or internal contact forces, $\mu(s)$, has been approximated by the forward finite difference between each successive elements

$$ \mu^j = \frac{\mu^{j+1} - \mu^j}{ds^j}, \qquad (41) $$

and $\mathbf{D}(\mathbf{q}^j)$ is the discrete counterpart of $\mathbf{D}(\mathbf{q}(s))$ introduced in Eq. (27) which reads, at each element

$$ \mathbf{D}(\mathbf{q}^j) = \begin{bmatrix} q_3^j & -q_4^j & -q_1^j \\ q_4^j & q_3^j & -q_2^j \\ q_1^j & q_2^j & -q_3^j \\ q_2^j & -q_1^j & -q_4^j \end{bmatrix}. \qquad (42) $$

As we did previously in the continuum case, we now require that the discrete variation of the Lagrangian $\delta\mathcal{E}_{\text{Lag}}$ (given in Eq. (24) as the sum of Eqs. (36), (38), (39), and (40)) vanishes for any arbitrary perturbations $\delta\mathbf{r}^i$ and $\delta\mathbf{q}^j$. This condition yields the set of algebraic equilibrium equations of the discrete unshearable and inextensible slender elastic rod. The condition that this variation is zero for any perturbed displacements $\delta\mathbf{r}^i$ leads to the balance of forces as a set of algebraic equations

$$ \mathbf{P}^0 + \mu^i = \mathbf{0} \quad \text{for } i=1 \qquad (43a) $$

$$ \mathbf{p}^i L/N + \mu^i - \mu^{i-1} = \mathbf{0} \quad \text{for } 2 \le i \le N \qquad (43b) $$
$$ \mathbf{P}^L - \mu^{L-1} = \mathbf{0} \quad \text{for } i=N+1. \qquad (43c) $$

Projected along the three directions of the global cartesian frame $(x, y, z)$, Eq. (43) yield $3(N+1)$ linear equations for the $3N$ unknowns $\mu^j$. In the limit $N$ very large, these equations converge to the continuous differential equations Eq. (29a). The condition that the variation of the Lagrangian is zero for any arbitrary perturbations $\delta q^i$ leads to the balance of moments as a set of discrete algebraic equations

$$ \tau^j + \mathbf{M}^0 + 2\mathbf{D}(\mathbf{q}^j)\mu^j - 2\alpha^j \mathbf{q}^j = \mathbf{0} \quad \text{for } j=1 \qquad (44a) $$

$$ \tau^j + \mathbf{m}^j + 2\mathbf{D}(\mathbf{q}^j)\mu^j - 2\alpha^j \mathbf{q}^j = \mathbf{0} \quad \text{for } 2 \le j \le N-1 \qquad (44b) $$

$$ \tau^j - \mathbf{M}^L + 2\mathbf{D}(\mathbf{q}^j)\mu^j - 2\alpha^j \mathbf{q}^j = \mathbf{0} \quad \text{for } j=N. \qquad (44c) $$

In these equations, $\tau^j = [\tau_1^j \tau_2^j \tau_3^j \tau_4^j]^T$ is the vector of net internal moment applied on the element $j$ written in the quaternion basis

$$ \tau^j = \sum_{k=1}^{3} G_k^j 2\mathbf{B}_k \mathbf{q}^{j+1} \quad \text{for } j=1 \qquad (44d) $$

$$ \tau^j = \sum_{k=1}^{3} (G_k^{j-1} 2\mathbf{B}_k \mathbf{q}^{j-1} - G_k^j 2\mathbf{B}_k \mathbf{q}^{j+1}) \quad \text{for } 2 \le j \le N-1 \qquad (44e) $$

$$ \tau^j = \sum_{k=1}^{3} G_k^{j-1} 2\mathbf{B}_k \mathbf{q}^{j-1} \quad \text{for } j=N. \qquad (44f) $$

and $2\mathbf{D}(\mathbf{q}^j)\mu^j$ is the moment resultant from the internal contact forces $\mu^j$ applied on the element $j$. In the limit of large $N$, Eq. (44) converge to the continuous differential equation Eq. (29a). Projected along the four elements $(i, j, k, 1)$ of the quaternion basis $\mathcal{H}$, Eq. (44) yield $4N$ nonlinear equations for the $4N$ unknowns $\mathbf{q}^j$ and $N$ unknowns $\alpha^j$. The missing equations required to compute all the unknowns are given by the geometrical constraints Eq. (30) which can be rewritten in the algebraic form

$$ (\mathbf{q}_1^j)^2 + (\mathbf{q}_2^j)^2 + (\mathbf{q}_3^j)^2 + (\mathbf{q}_4^j)^2 - 1 = 0 \quad \text{for } 1 \le j \le N \qquad (45a) $$

$$ \mathbf{r}^{j+1} - \mathbf{r}^j - \mathbf{d}_3^j \, ds^j = \mathbf{0} \quad \text{for } 1 \le j \le N, \qquad (45b) $$

where we used the forward finite difference in Eq. (32) to approximate $\mathbf{r}'(s)$.

The $4N$ geometrical constraints given in Eq. (45), together with the $7N+3$ equilibrium equations Eqs. (43) and (44) form the set of algebraic equations describing the constrained equilibrium configuration of the rod represented by the $7N+3$ degrees of freedom $\mathbf{r}^j$ and $\mathbf{q}^j$ and the $4N$ Lagrange parameters $\alpha^j$ and $\mu^j$.

We highlight the fact that the only approximations made in the above equations arise from the finite difference discretization since the initial continuous formulation is geometrically-exact due to the use of quaternions. Furthermore, it is remarkable to notice that the equilibrium configurations of the extremely twisted and bended elastic rod can be represented by the smooth polynomial equations Eqs. (43)-(45). In the next section, we exploit the particularly smooth nonlinearities of the equilibrium equations by using Asymptotic Numerical Methods (ANM) (Damil and Potier-Ferry, 1990; Cochelin, 1994; Cochelin et al., 1994) which are efficient path-following techniques that give access to semi-analytical solution branches of polynomial nonlinear algebraic systems.

## 4.2. Asymptotic Numerical Method

We now explain and adapt the particular ANM introduced in Cochelin (1994) for solving the equilibrium equations of slender elastic rods described above. This ANM is a perturbation technique allowing for the computation of a large part of a solution branch of quadratic algebraic system of equations with only one stiffness inversion. Applied in a step-by-step manner, one can compute a complex nonlinear branch by a succession of local asymptotic expansions and thus determine a semi-analytical bifurcation diagram. Because of the local analytical representation of the branch within each step, this continuation technique has a number of important advantages when compared to classical predictor-corrector schemes (Cochelin, 1994). In particular, the algorithm is fully automatic, remarkably robust, and faster than incremental-iterative methods.

To apply the asymptotic numerical method to the mechanics of elastic rods, we first rewrite the algebraic nonlinear systems of equilibrium equations Eqs. (43)-(45) in the compact form

$$ \mathbf{f}(\mathbf{w}, \lambda) = \mathbf{0}. \qquad (46) $$

where $\mathbf{f}$ is a smooth nonlinear vector valued function in $\mathbb{R}^{11N+3}$ with $N$ the number of elements of the discretized rod, $\lambda$ is a scalar control parameter (usually a mechanical or geometrical parameter of the physical problem such as the rotation or
displacement at one end of the rod) and **w** is the vector of unknowns which, in our case reads

$$ \mathbf{w} = [\mathbf{r}^1 \mathbf{q}^1 \dots \mathbf{r}^N \mathbf{q}^N \mathbf{r}^{N+1} \dots \alpha^1 \mu^1 \dots \alpha^N \mu^N]^T. \quad (47) $$

According to the discretization presented in Section 4.1, **w** is a vector of size 11N + 3 which includes the 7N+3 mechanical degrees of freedom separated into positions **r**^i and quaternions **q**^i. Moreover, the 4N Lagrange parameters are required to impose the geometrical constraints.

In what follows, in a process that we refer as recasting, we now transform Eq. (46) into a quadratic form, which is a particular framework of the ANM that allows us to formally and systematically write a large class of physical problems including rods (Cochelin, 1994; Cochelin et al., 1994). Given the original cubic form of Eq. (46), this quadratic recast is achieved introducing a new vector of unknowns **u** of size 14N+1,

$$ \mathbf{u} = [\mathbf{w} G_1^1 G_2^1 G_3^1 \dots G_N^1 G_N^2 G_N^3 \lambda]^T, \quad (48) $$

which includes the initial vector of unknowns **w** given in Eq. (47), the control parameter λ and where we added the 3(N-1) flexural and torsional internal torques G_k^i(q^i) introduced in Eq. (37). Using the new vector **u** instead of **w**, we can recast the cubic nonlinear vector valued function **f**(**w**, λ) given in Eq. (46) into the quadratic form

$$ \mathbf{f}(\mathbf{u}) = \mathbf{L}_0 + \mathbf{L}(\mathbf{u}) + \mathbf{Q}(\mathbf{u}, \mathbf{u}) = \mathbf{0}, \quad (49) $$

where **f**(**u**) is a vector in $\mathbb{R}^{14N}$ since we added the 3(N-1) nonlinear quadratic equations Eq. (37) to **f**(**w**), $\mathbf{L}_0$ is a constant vector and $\mathbf{L}(\bullet, \bullet)$ and $\mathbf{Q}(\bullet, \bullet, \bullet)$ are a linear and bilinear vector valued operators, respectively. The expression of **f**(**u**), representing the equilibrium equations of our inextensible and unshearable elastic rod given in Eqs. (43)-(45), is provided in Appendix A. In Section 5 where we will apply our method to the quasi-static writhing of a double-clamped elastic rod, we will illustrate, by way of example, how to introduce the boundary conditions and the control parameter in the vector **f**(**u**) of Eq. (49).

We can now proceed and compute the solutions **u** of the set of quadratic equations Eq. (49) with the asymptotic numerical method. This technique is based on the perturbation of the vector of unknowns **u** in terms of a path-parameter **a** in the form of the asymptotic expansion

$$ \mathbf{u}(a) = \mathbf{u}_0 + \sum_{p=1}^{m} a^p \mathbf{u}_p, \quad (50) $$

where $\mathbf{u}_0$ is the starting fixed point, solution of Eq. (49), **m** is the truncation order of the power series and **a** is the path-parameter which will be formally defined below. Replacing $\mathbf{u}(a)$ by its asymptotic expansion Eq. (50) in the quadratic form Eq. (49), we obtain the quadratic Taylor series in the neighborhood of $\mathbf{u}_0$,

$$ \mathbf{f}(\mathbf{u}(a)) = \mathbf{L}_0 + \mathbf{L}(\mathbf{u}_0) + \mathbf{Q}(\mathbf{u}_0, \mathbf{u}_0) + a[\mathbf{L}(\mathbf{u}_1) + 2\mathbf{Q}(\mathbf{u}_0, \mathbf{u}_1)] + \sum_{p=1}^{m} a^p \left[ \mathbf{L}(\mathbf{u}_p) + 2\mathbf{Q}(\mathbf{u}_0, \mathbf{u}_p) + \sum_{i=1}^{p-1} \mathbf{Q}(\mathbf{u}_i, \mathbf{u}_{p-i}) \right]. \quad (51) $$

Recalling that $\mathbf{u}_0$ is a solution of Eq. (49), we can rewrite Eq. (51) in the form of a power series of a quadratic vector valued function $\mathbf{f}_p$ of size 14N,

$$ \mathbf{f}(\mathbf{u}(a)) = a\mathbf{f}_1 + a^2\mathbf{f}_2 + \dots + a^m\mathbf{f}_m = \mathbf{0}. \quad (52) $$

Since Eq. (52) has to be verified for every value of **a**, we need $\mathbf{f}_p = \mathbf{0}$ for every order $p \le m$. This leads to **m** linear systems in $\mathbf{u}_p$ in the form

$$ \forall p \in [1 \dots m], \quad \left. \frac{\partial \mathbf{f}}{\partial \mathbf{u}} \right|_{\mathbf{u}^0} \mathbf{u}_p = \mathbf{f}_p^{nl}, \quad (53a) $$

where, due to the particular quadratic form of Eq. (49), the Jacobian matrix of $\mathbf{f}(\mathbf{u})$ evaluated at the initial solution vector $\mathbf{u}_0$ reads

$$ \left. \frac{\partial \mathbf{f}}{\partial \mathbf{u}} \right|_{\mathbf{u}^0} \mathbf{u}_p = \mathbf{L}(\mathbf{u}_p) + 2\mathbf{Q}(\mathbf{u}_0, \mathbf{u}_p), \quad (53b) $$

and the nonlinear vector $\mathbf{f}_p^{nl}$ on the right-hand side of Eq. (53a) consists of a quadratic sum that only depend on the previous order

$$ \mathbf{f}_p^{nl} = 0 \quad \text{for } p=1 \quad (53c) $$

$$ \mathbf{f}_p^{nl} = - \sum_{i=1}^{p-1} \mathbf{Q}(\mathbf{u}_i, \mathbf{u}_{p-i}) \quad \text{for } 1 < p \le m. \quad (53d) $$

The original nonlinear problem in Eq. (49) has thereby been reduced to a definite set of **m** linear systems given in Eq. (53a) where the matrix on the left-hand side is identical for each order.

However, each linear system in Eq. (53) is, so far, under-determined since the dimension of **f** is 14N whereas the size of the vector of state variable at order **p**, $\mathbf{u}_p$, is 14N+1. The remaining equation is provided by the definition of the path
parameter $a$ as defined in Cochelin (1994). We consider a measure that includes the entire set of physical unknowns and that is also robust toward limit and bifurcation points, i.e. an arc-length measure. Mathematically, we identify the path parameter $a$ as the projection of the vector of state variables increment $\mathbf{u}-\mathbf{u}_0$ on the normalized tangent vector $\mathbf{u}_1$ (see Fig. 4(a)),

$$a = (\mathbf{u}(a) - \mathbf{u}_0)^T \mathbf{u}_1. \qquad (54)$$

Replacing $\mathbf{u}(a)$ by its asymptotic expansion, Eq. (50), in Eq. (54), we obtain

$$a \mathbf{u}_1^T \mathbf{u}_1 - a + a^2 \mathbf{u}_1^T \mathbf{u}_2 + \cdots + a^m \mathbf{u}_1^T \mathbf{u}_m = 0. \qquad (55)$$

Verifying Eq. (55) at every power of $a$ provides us with the supplementary equations at every $p$

$$\mathbf{u}_1^T \mathbf{u}_p = \delta_{p1} \quad \text{for } 1 \le p \le m, \qquad (56)$$

where $\delta_{p1}$ is the Kronecker delta, $\delta_{11} = 1$ at the first order $p=1$ and it is zero otherwise.

Finally, the original nonlinear problem in Eq. (49) has now been transformed in the well-posed $m$ linear systems in $\mathbb{R}^{14N+1}$

$$\left[ \begin{matrix} \partial f \\ \partial \mathbf{u} \end{matrix} \right]_{\mathbf{u}0} \mathbf{u}_1^T \mathbf{u}_p = \left\{ \begin{matrix} \mathbf{f}_p^{nl} \\ \delta_{p1} \end{matrix} \right\} \quad \text{for } 1 \le p \le m, \qquad (57)$$

with a unique solution $\mathbf{u}_p$ that we can solve iteratively since each of these vectors is defined with the solution of the previous order according to the definition of $\mathbf{f}_p^{nl}$ given in Eq. (53). In this linearized numerical problem, the only matrix to inverse is the one on the left-hand side of Eq. (57) since it is the same at every order $p$. This is in striking contrast with classical predictor-corrector methods where one needs to actualize the Jacobian for every linear systems (Crisfield, 1991).

For practical purposes, one will inverse the matrix and compute the unknown $\mathbf{u}_1$ at first order independently and then compute the higher orders $\mathbf{u}_p$ for $1 < p \le m$ from the well-defined systems in Eq. (57).

Once each $\mathbf{u}_p$ has been found, we still have to estimate the validity domain of the asymptotic expansion since Eq. (49) can only be true for values of the perturbation parameters $a$ inside the radius of convergence of the power series given in Eq. (50) (see Fig. 4(b)). A simple, robust and accurate way of calculating an approximation of the convergence radius $a_{max}$, explained in detail in Cochelin (1994), is to assume that a solution branch is acceptable as long as the norm of the nonlinear $(14N+1)$-dimensional vector field $\mathbf{f}(\mathbf{u}(a))$ is less than a tolerance criterion $\epsilon$,

$$\forall a \in [0, a_{max}], \quad \|\mathbf{f}(a)\| < \epsilon, \qquad (58)$$

where $\epsilon$ determines the accuracy of our numerical results. We have computed the $\mathbf{u}_p$ according to the power series expansion of $\mathbf{f}(a)$ given in Eq. (52) so that the norm of $\mathbf{f}$ is zero up to the truncation order $m$. Consequently, the residue of this series is given by the norm of $\mathbf{f}(a)$ for $p > m$. Assuming that the order $m+1$ dominates in the residue, we obtain the relation between the norm of $\mathbf{f}(a)$ and the vector at the order $m+1$,

$$\|\mathbf{f}(a)\| \approx a^{m+1} \|\mathbf{f}_{m+1}\| \qquad (59)$$

<figure><img src="image_5.png" /><figcaption>Fig. 4. (a) The path parameter $a$ is identified as the projection of the solution branch on the tangent vector. This projection is normalized by the length of the tangent vector $[u(1) \lambda(1)]$ which is set to unity and (b) behavior of a power series close to the radius of convergence. The considered nonlinear unidimensional equation is $\mathbf{f}(u(a)) = u(1+a) - 1 = 0$. Its unique solution $u(a) = 1/(1+a)$, represented as solid line, can be expanded asymptotically as $u(a) \approx 1 - a + a^2 + \cdots + (-1)^m a^m + (-1)^{m+1} a^{m+1}$. The asymptotic expansions for $m=5, 10, 15$ and $20$ are represented as different lines. We see that $u(a)$ can be approximated up to the radius of convergence $a=1$. Applying our ANM method, the step length calculated through Eqs. (59)-(61) would give $a_{max} = \epsilon^{1/(m+1)} = 0.4642$ with $\epsilon = 1 \times 10^{-7}$ and $m=20$.</figcaption></figure>
where $\mathbf{f}_{m+1} = \partial \mathbf{f} / \partial \mathbf{u}|_{\mathbf{u}_0} \mathbf{u}_{m+1} - \mathbf{f}_{m+1}^{nl}$ according to Eqs. (52) and (53). Replacing $\|\mathbf{f}(a)\|$ by its definition Eq. (59) in Eq. (58) leads to

$$ \|\mathbf{f}(a)\| < \varepsilon \Leftrightarrow a < \left( \frac{\varepsilon}{\|\mathbf{f}_{m+1}\|} \right)^{1/m+1} \qquad (60) $$

which sets an upper limit to the path-parameter $a$. While truncating the asymptotic series $\mathbf{u}(a)$ given in Eq. (50) at the order $m$, we implicitly assumed that $\mathbf{u}_{m+1} = \mathbf{0}$ so that $\mathbf{f}_{m+1} = \mathbf{f}_{m+1}^{nl}$. Replacing $\mathbf{f}_{m+1}$ by the nonlinear term $\mathbf{f}_{m+1}^{nl}$ in Eq. (60), we obtain an estimation of the maximum step length,

$$ a_{max} = \frac{\varepsilon}{\|\mathbf{f}_{m+1}^{nl}\|} \qquad (61) $$

For practical purposes, when applying Eq. (61) to any quadratic vector valued function $\mathbf{f}(\mathbf{u}(a))$, we found that $\mathbf{f}(\mathbf{u}(a_{max})) \approx \varepsilon$. In general, the power series in Eq. (50) converges slowly, close to the radius of convergence (Cochelin et al., 1994). Decreasing $\varepsilon$ leads to a diminishing of $a_{max}$ but, more importantly, results in an increased accuracy of the computed series. An optimal value of $m$ and $\varepsilon$ in terms of convergence of the asymptotic series, accuracy of the solution given by Eq. (58) and size of $a_{max}$ is found empirically for the following range of parameters (Cochelin et al., 1994): $\varepsilon = 1 \times 10^{-7}$ and $15 \le m \le 20$.

The power series expansion given in Eq. (50) and computed with the $m$ linear systems in Eq. (57), together with the maximum step size $a_{max}$ given by Eq. (61) define a portion of the nonlinear equilibrium branches of the slender elastic rod in terms of a given control parameter $\lambda$. The next step of our calculation, the continuation of the solution branch, is now computed by applying the present asymptotic numerical method taking $\mathbf{u}(a_{max})$ as the new starting equilibrium $\mathbf{u}_0$ of the new portion. A complete solution branch is therefore constructed as a succession of semi-analytical portions in the form of Eq. (50), whose length is automatically determined through the estimation of the convergence radius of each power series as sketched in Fig. 5(a). Unlike classical predictor-corrector methods (Riks, 1979; Doedel, 1981), our step length is adaptive; it is naturally large for weakly nonlinear solutions and becomes shorter when strong nonlinearities occur. As a consequence, the automatization of the continuation method is significantly easier and more robust than with standard predictor-corrector methods.

One would expect that the residue $\|\mathbf{f}(\mathbf{u}(a))\|$ would increase progressively at every continuation step so that the accuracy of the new starting equilibrium $\mathbf{u}_0$ would gradually decrease (Cochelin et al., 1994). In practice however, it is rare to see the residue increase up to $10\varepsilon$, especially given the smooth nature of the nonlinearities of our equilibrium equations Eqs. (43)–(45) for thin elastic rods. In Section 5, where we implement this continuation method to a series of specific test-case problems, all the bifurcation diagrams are computed with a residue smaller than $\varepsilon = 1 \times 10^{-7}$, with no correction step (note that a correction step may be necessary in the general ANM framework in the case of non-polynomial nonlinearities (Karkar et al., 2013)).

When the accuracy $\varepsilon$ is set to be small enough by the user, the ANM is able to follow the branch whenever bifurcation are encountered (Baguet and Cochelin, 2003). This is a remarkably robust property for a path-following algorithm, especially when compared to predictor-corrector techniques which typically would systematically bifurcate because of the discrete nature of their continuation steps. Nevertheless, we now need a special procedure to switch branches in order to determine the full bifurcation diagram. A classic strategy is to slightly modify the original equilibrium equations Eq. (49) by adding a low-norm perturbation vector

$$ \mathbf{f}_p(\mathbf{u}) = \mathbf{L}_0 + \mathbf{L}(\mathbf{u}) + \mathbf{Q}(\mathbf{u}, \mathbf{u}) + c\mathbf{P}, \qquad (62) $$

where $\mathbf{f}_p(\mathbf{u})$ is the perturbed problem, $\mathbf{P}$ a normalized vector of constant random numbers and $c$ is the intensity of the perturbation. This additional perturbation procedure transforms the exact bifurcation into a perturbed bifurcation (see Fig. 5(b)). The idea is to use the perturbed branch to bifurcate on the non-crossing branch (Allgower and Georg, 2003). Changing the sign of $c$ allows us to explore a symmetrical quasi-bifurcation as represented in Fig. 5(b) with the perturbed branches $P = +c$ and $P = -c$. Finally, in order to transition from the original to the perturbed problem, or vice versa, a correction step is mandatory (any predictor-corrector methods would be efficient since the perturbed solutions are very close to exact ones). The combination of different (positive and negative) values of the intensity of the perturbation and several correction steps allows one to explore the full bifurcation diagram of the slender elastic rod described by the equilibrium equations Eqs. (43)–(45).

So far, we have presented the ANM method in the context of thin elastic rods. We can now compute the bifurcation diagrams of our slender elastic rod under various mechanical and geometrical environments. However, the final crucial step of determining the stability of the solution branches is still missing, which is the focus of the following section.

### 4.3. Stability analysis

Determining the local stability of equilibrium branches is crucial for the physical understanding of the mechanical behavior of slender elastic rods, one of the main motivations being that locally unstable branches cannot be observed
<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) Schematic five continuation steps with the asymptotic numerical method. The step length of each portion of branch is determined a posteriori by analyzing the validity of the asymptotic solution. The highly nonlinear solution is known analytically for each portion and (b) branch switching through the perturbation method given in Eq. (62). Changing the sign of the intensity of perturbation, c, allows us to explore all the branches after a bifurcation point.</figcaption></figure>

experimentally, and must therefore be classified. Another advantage for gaining knowledge on the stability of a solution is that the loss of local stability is often associated with a bifurcation point. Assessing the stability is then useful to detect bifurcation points and navigate through the bifurcation diagram as illustrated in Fig. 5(b).

The equilibrium equations on their own are not sufficient to determine the local stability of the solutions; we also need to compute whether the solution is a local minimum or maximum of the system's energy. For practical purposes, we need to derive the second-order conditions of the geometrically constrained energy; a theoretical and numerical procedure that is well established (Luenberger, 1973; Manning et al., 1998; Luenberger and Ye, 2008).

In the previous sections, we have shown how to compute and follow the branches of solutions of the nonlinear algebraic equilibrium equations given in Eqs. (43)-(45). Since the resulting computed bifurcation diagrams are semi-analytical, we are therefore able to evaluate the solution $\mathbf{u}_e$ given in Eq. (48) at any finite value of the control parameter $\lambda_e$. Let $\mathbf{x}_e$, $\boldsymbol{\alpha}_e$ and $\boldsymbol{\mu}_e$ be the vectors of degrees of freedom and Lagrange multipliers respectively, associated with the vector solution $\mathbf{u}_e$. The vector $\mathbf{x}_e$ is a solution of the equilibrium equations Eqs. (43) and (44) for $\lambda = \lambda_e$ that satisfies the functional geometrical constraints given in Eq. (45). Rewriting Eqs. (43)-(45) in an energy minimization framework, $\mathbf{x}_e$ is the actual solution of the $n = 7N + 3$-dimensional constraint minimization problem

$$ \nabla(\mathcal{E}_e(\mathbf{x}_e) + W(\mathbf{x}_e)) + \sum_{i=1}^{N} \boldsymbol{\alpha}_e^i \nabla C_{\alpha}^i + \sum_{j=1}^{3N} \boldsymbol{\mu}_e^j \nabla C_{\mu}^j = \mathbf{0}, \qquad (63a) $$

subject to the $m = 4N$ functional constraints

$$ C_{\alpha}^{j}(\mathbf{x}_{e}) = \mathbf{q}^{j}\mathbf{q}^{j} - 1 = 0, \qquad (63b) $$

$$ C_{\mu x}^{j}(\mathbf{x}_{e}) = r_{x}^{j+1} - r_{x}^{j} - 2q_{1}^{j}q_{3}^{j} - 2q_{2}^{j}q_{4}^{j} = 0, \qquad (63c) $$

$$ C_{\mu y}^{j}(\mathbf{x}_{e}) = r_{y}^{j+1} - r_{y}^{j} - 2q_{2}^{j}q_{3}^{j} + 2q_{1}^{j}q_{4}^{j} = 0, \qquad (63d) $$

$$ C_{\mu z}^{j}(\mathbf{x}_{e}) = r_{z}^{j+1} - r_{z}^{j} + (q_{1}^{j})^{2} + (q_{2}^{j})^{2} - (q_{3}^{j})^{2} - (q_{4}^{j})^{2} = 0, \qquad (63e) $$

for the positional nodes $1 \le j \le N$. In Eq. (63), $\nabla$ is the gradient operator,<fnref n="1" /> $C_\alpha = [C_\alpha^1 C_\alpha^2 \dots C_\alpha^N]^T$ and $C_\mu = [C_\mu^1 C_\mu^2 \dots C_\mu^N]^T$.

According to the necessary and sufficient first order conditions of constrained minimization problems (Luenberger, 1973; Luenberger and Ye, 2008), the vector solution $\mathbf{x}_e$ is a local extremum (a minimum or maximum) of the total energy $\mathcal{E}(\mathbf{x}_e) = \mathcal{E}_e(\mathbf{x}_e) + W(\mathbf{x}_e)$ subject to the $m$ constraints in Eqs. (63b) and (63c). Supposing also that the $n \times n$ matrix

$$ \mathbf{L}(\mathbf{x}_e) = \nabla^2(\mathcal{E}_e(\mathbf{x}_e) + W(\mathbf{x}_e)) + \sum_{i=1}^{N} \boldsymbol{\alpha}_e^i \nabla^2 C_{\alpha}^i + \sum_{j=1}^{3N} \boldsymbol{\mu}_e^j \nabla^2 C_{\mu}^j, \qquad (64) $$

$$ \nabla f(\mathbf{x}) = \left[ \frac{\partial f(\mathbf{x})}{\partial \mathbf{x}_1} \quad \frac{\partial f(\mathbf{x})}{\partial \mathbf{x}_2} \quad \dots \quad \frac{\partial f(\mathbf{x})}{\partial \mathbf{x}_n} \right]^T. $$

<footnote n="1">¹ For a real-valued function $f \in C^1$ on $\mathbb{R}^n$ such that $f(\mathbf{x}) = f(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_n)$, we define the gradient of $f$ to be the $n$-dimensional vector,</footnote>
where $\nabla^2$ is the Hessian operator,<fnref n="2" /> is positive definite on the $m$-dimensional subspace $M = \{\mathbf{y} : \nabla\mathbf{h}(\mathbf{x}_e)\mathbf{y} = \mathbf{0}\}$ with $\mathbf{h}(\mathbf{x}_e) = [C_\alpha C_\mu]^T$, that is, for $\mathbf{y} \in M$ and $\mathbf{y} \neq \mathbf{0}$ that holds $\mathbf{y}^T\mathbf{L}(\mathbf{x}_e)\mathbf{y} > 0$, then, according to the second-order necessary and sufficient conditions, $\mathbf{x}_e$ is a strict local minimum of $\mathcal{E}(\mathbf{x}_e)$ subject to $C_\alpha = \mathbf{0}$ and $C_\mu = \mathbf{0}$.

The matrix $\mathbf{L}(\mathbf{x}_e)$ is the matrix of second partial derivatives, with respect to $\mathbf{x}$, of the discrete counterpart of the Lagrangian given in Eq. (24). When restricted to the subspace $M$ that is tangent to the constraint surface and which we denote by $\mathbf{L}_M$, $\mathbf{L}(\mathbf{x}_e)$ plays the role in second-order conditions directly analogous to that of the Hessian of the objective function in the unconstrained case (Luenberger, 1973; Luenberger and Ye, 2008). The eigenvalues, $\sigma_i$, and associated eigenvectors, $\mathbf{y}_i$, of $\mathbf{L}_M$, determine the local stability of the solutions of the constrained minimization problem. Mathematically, $\mathbf{L}_M$ is a $(n-m) \times (n-m)$ matrix defined, at each equilibrium point $\mathbf{x}_e$, as

$$ \mathbf{L}_M(\mathbf{x}_e) = (\ker(\nabla\mathbf{h}))^T \mathbf{L}(\ker(\nabla\mathbf{h})), \quad (65) $$

where $\ker(\bullet)$ denotes the kernel operator. Analyzing the $n-m$ eigenvalues $\sigma_i$ gives us information on the behavior of the associated perturbation $\mathbf{y}_i(t)$ in the neighborhood of the equilibrium $\mathbf{x}_e$. According to Lyapunov's theory (Bažant and Cedolin, 2010; Guckenheimer and Holmes, 1983):

* If $\sigma_i > 0$ for all $i \in [1...n-m]$, small perturbations remain small and the equilibrium is locally stable.

* If one index $i \in [1...n-m]$ exists, for which $\sigma_i < 0$, one perturbation diverges and the equilibrium is locally unstable.

* If $\sigma_i \ge 0$ for all $i \in [1...n-m]$ and if there exists one index $k$ such that $\sigma_k = 0$, the first order is insufficient to draw conclusions on the local stability of the equilibrium. In that case, a perturbation at higher order is necessary.

Applying this method to a sufficient number of fixed points $\mathbf{x}_e$ along the equilibrium branches, computed with the previous ANM method, allows us to determine the stability of the bifurcation diagram. The previous finite difference discretization presented in Section 4.1, together with the ANM algorithm described in Section 4.2 and the previous stability method, completes the semi-analytical continuation technique that we developed to compute and follow the equilibrium branches and the stability of an inextensible slender elastic rod undergoing extreme displacements and rotations. The combination of the conciseness and relative simplicity of our method offer the opportunity for it to be implemented in any programming language. In the following, we briefly present MANlab (Arquier, 2007; Karkar et al., 2010), an open-source bifurcation analysis software that provides a convenient framework to implement the previous numerical methods.

## 4.4. MANlab: an open-source bifurcation analysis software

MANlab is an interactive software package for the continuation and bifurcation analysis of algebraic systems, based on ANM continuation, and first released in 2009 (Arquier, 2007). Thanks to the implementation of most of the ANM equations in MATLAB using an object-oriented approach (Karkar et al., 2010), MANlab makes it simpler for the user to solve the system of Eqs. (43)–(45) and the stability of the solutions given by the second-order condition through Eq. (64). MANlab has a graphical user interface (GUI) with buttons, on-line inputs and graphical windows for generating, displaying and analyzing the bifurcation diagram and the solutions of the system. A unique identifying feature, when compared with other continuation codes, is that its computational efficiency, highlighted above, allows for interactive control of the continuation process. The full interactive and semi-automatic procedure consists of computation of a portion of a branch, choice of a new branch at a bifurcation point, reverse direction of continuation on the same branch, jump capability between solutions, visualization of user-defined quantities at a particular solution point, selection and deletion of a branch, or of one of its portion, possibility of correction step with a Newton-Raphson method and determination of the local stability of the solution.

To enter the system of equations, the user simply has to provide the three vector valued Matlab functions corresponding to the constant, linear and quadratic operators $\mathbf{L}_0$, $\mathbf{L}(\mathbf{u})$, and $\mathbf{Q}(\mathbf{u}, \mathbf{u})$ given in Eq. (49). To assess the local stability at each computed solution point in MANlab (Lazarus and Thomas, 2010), one can also provide the Hessian of the constrained Lagrangian restricted to $M$, $\mathbf{L}_M(\mathbf{x}_e)$ given in Eq. (65) and the package will automatically compute the eigenvalues of the linearized problem according to the previous section. Thanks to the flexibility offered by the MATLAB environment, users become rapidly familiar with MANlab. Calling of external routines such as finite elements codes is also possible.

In the following section, we validate our semi-analytical continuation method by using the MANlab package to simulate a precision model experiment; the quasi-static writhing of a double-clamped slender elastic rod, which equilibria are solutions of the discrete equilibrium equations given in Eqs. (43)–(45).

$$ \nabla^2 f(\mathbf{x}) = \left[ \frac{\partial^2 f(\mathbf{x})}{\partial \mathbf{x}_i \partial \mathbf{x}_j} \right]. $$

<footnote n="2">² We define the Hessian of $f$ at $\mathbf{x}$ ($f \in C^2$) to be the $n \times n$ symmetric matrix denoted $\nabla^2 f(\mathbf{x})$,</footnote>
5. Following the equilibria of an extremely twisted elastic rod

Having introduced the general theoretical and numerical framework to continue the equilibria and stability of slender
elastic rods, we proceed by implementing the specific problem of the writhing (extreme twisting) of a clamped heavy elastic
rod. Even though this fundamental problem appears seemingly simple, it can display an array of complex behavior with
intricate bifurcation diagrams and has received significant attention in the literature (Thompson and Champneys, 1996;
Furrer et al., 2000; Goriely and Tabor, 1998a; VanderHeijden and Thompson, 2000; Goyal et al., 2008; Liu and Center, 1975;
Yabuta, 1984; Goss et al., 2005; Goriely and Tabor, 1997b, 1997c, 1997a; Neukirch et al., 2002). The writhing of an elastic
rod is therefore an ideal scenario to challenge our theoretical and computational framework by contrasting the numerical
results with our own precision model experiments that were especially developed for the testing and validation of our
continuation method.

In this Section, we first present our apparatus and model experiments which consist of quasi-statically increasing the
rotation angle at one end of a slender elastic rod fixed between two concentrically aligned horizontal clamps. One of the
originalities of our experiments is that we fabricate our own elastic rods, enabling us to accurately target their material and
geometrical properties. In particular, we have full control in setting their intrinsic natural curvatures. After describing how
to account for the kinematic boundary conditions and control parameter specific to this writhing problem in the numerical
model described in Section 4, we compare a series of experimental and numerical results for two different elastic rods: a
straight rod with no natural curvature and a curved rod.

5.1. *Manufacturing of rods and experimental apparatus*

Our rods are cast by injecting vinylpolysiloxane (VPS), a two-part silicone-based elastomer, into a flexible PVC tube of inner and outer diameters $D_l=3.1$ mm and $D_o=5$ mm, respectively. The PVC mold is first wound around a cylinder of external radius $R_e$ and then injected with VPS, which eventually cross-links at room temperature (see inset of Fig. 6). After a setting period of 24 h, to ensure complete curing of the polymer, the outer flexible PVC pipe is cut to release the inner slender VPS elastic rod with a constant natural curvature $\hat{k}_1(s) = 1/(R_e + D_o/2)$ and a circular cross-section $R = D_l/2 = 1.55$ mm. The rod's second moments of area are $I_1 = I_2 = \pi R^4 / 4$ and $I_3 = J = \pi R^4 / 2$. We measure the Young's modulus of the elastomer to be $E = 1300 \pm 60$ KPa, a volumic mass $\rho = 1200$ kg/m$^3$ and a Poisson ratio of $\nu \approx 0.5$, so that its shear modulus is $G = E/2(1+\nu) = 433$ KPa.

The cast rod ($L=30$ cm long) is then attached between two horizontal concentric drill chucks of a lathe, separated by a
distance $d=22$ cm. A photograph of the side view of the experiment is presented in Fig. 6. The boundary conditions of the
rod are set to be rigidly clamped at both ends. For future representation of the rod configurations, we choose the origin of
the cartesian frame $(x, y, z)$ to be located at the clamp at the left extremity of the rod (see Fig. 6). The clamp located at the
origin, at the curvilinear coordinate $s=0$, is completely fixed but the other clamp, located at $s=L$, can be rotated with respect
to the $y$-axis, thereby imposing a rotation angle $\phi$ (see Fig. 6). Initially, for $\phi=0^\circ$, we ensure that the sign of the intrinsic
curvature $\hat{k}_1$ is such, that the rod naturally bends downwards, in the direction of gravity and that the difference between
twist angles, $\kappa_3(s)$, at both ends of the rod, is zero. In that configuration, the equilibrium shape of the clamped rod is close to
a planar inflectional elastica as theoretically described in VanderHeijden et al. (2003); the only difference arising from the
effects due to gravity which induces a catenary-like configuration. Our writhing experimental protocol then consists of
quasi-statically increasing the rotation angle, $\phi$, at $s=L$ and quantifying the evolution of equilibrium states as a function of
this control parameter, $\phi$. A variety of measurements on the configurations of the rod are performed by imaging the top of
the experiment (using a Nikon D90 SLR camera) and subsequent image processing.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** The writhing experiment. A $L=30$ cm long elastic rod with a circular cross section of radius 1.55 mm, Young's modulus $E=1300 \pm 60$ KPa, volumic mass $\rho = 1200$ kg/m$^3$ and a custom made constant natural curvature $\hat{k}_1$ is fixed at both ends between two concentrically aligned drill chucks separated by a distance $d=22$ cm. Our experiment consists of quasi-statically increasing the rotation $\phi$ at one end and investigating the evolution of equilibrium state with the control parameter $\phi$. Inset courtesy of Khalid Jawed: Fabrication process of an elastomeric rod. During working time, the PVC tubes containing the silicone-based elastomer are wound around cylinders of external radius $R_e$ and left in that position for 24 h. After demolding, it confers to the rod a constant natural curvature $\hat{k}_1(s) \approx 1/R_e$.</figcaption></figure>
## 5.2. Modeling of the boundary conditions for the writhing configuration

Before we can proceed with a direct comparison between experimental and numerical results, we first need to precise how to account for the specific kinematic boundary conditions and control parameter, $\Phi$, relevant to this specific writhing configuration, in our general numerical framework presented in Section 4. This specific implementation will serve as an example, which, following the series of procedures and rationale described below, can be extended to other kinematic conditions to solve a variety of other problems involving thin rods.

Representing the slender elastic rod of Fig. 6 by the discrete 3D Cosserat curve of Fig. 3 and applying the numerical method of Section 4, we can write its equilibrium equations in the quadratic form $f(\mathbf{u})$ defined in Appendix A. In the writhing experiment, gravity is the only external force applied to the rod. This gravitational force is represented by the weight of each element, reported at each node. In $f(\mathbf{u})$ given in Eq. (A.1), we can therefore write $\mathbf{P}^0 = \mathbf{P}^L = -1/2\rho g\pi R^2 L^2 / N^2 \mathbf{e}_z$ at both end nodes and $\mathbf{p}^i = -\rho g\pi R^2 L / N \mathbf{e}_z$ for all the internal nodes $1<i<N+1$, where $g=9.81 \text{ m/s}^2$ is the gravitational acceleration, $N$ is the number of segments and $\mathbf{e}_z$ is the unit vector in the z-direction. Since there are no external moments, we can also write $\mathbf{m}^i = \mathbf{0}$ for all the internal elements $1<j<N$ and $\mathbf{M}^0 = \mathbf{M}^L = \mathbf{0}$ at both ends.

In addition to the mechanical parameters which are general to thin rods, we also need to account for the kinematic boundary conditions and control parameter $\Phi$, specific to the writhing experiment, which are not included in the general formulation of $f(\mathbf{u})$ given in Eq. (A.1). At the left extremity of the rod ($s=0$), the first node must be fixed to the origin and the first element, which is naturally pointing in the z-direction so that $\mathbf{d}_3$ is parallel to $\mathbf{e}_z$, has to be re-oriented along the y-direction (see Fig. 7a). Mathematically, this translates into two functional constraints depending on the positional and rotational degrees of freedom, $\mathbf{r}^1$ and $\mathbf{q}_1$ alone,

$$C_r^0[\mathbf{r}^1] = \mathbf{r}^1 = \mathbf{0}, \qquad (66a)$$

$$C_q^0[\mathbf{q}^1] = \mathbf{q}^1 - \left[-\frac{\sqrt{2}}{2} \mathbf{0} \quad \mathbf{0} \quad \frac{\sqrt{2}}{2}\right]^T = \mathbf{0}. \qquad (66b)$$

Whereas the first condition (66a) can be physically interpreted as the fixed boundary conditions at $s=0$, $\mathbf{r}^1 = \mathbf{0}$, the second Eq. (66b) is more difficult to interpret due to the lack of direct physical significance of the quaternions. To determine the constraint set in Eq. (66b), we used the relation between Euler's principal geometric quantities and a set of unit quaternions, whose expression was given in Eq. (7). According to Euler's rotation theorem (Kuipers, 1999), the imposed rotation applied to the first segment (shown schematically in Fig. 7(a) can be represented by a rotation angle of $-\pi/2$ rad around the unit length vector $\mathbf{b}=[1 \ 0 \ 0]^T$. Following the conversion formula of Eq. (7), the equivalent representation in terms of unit quaternions is given by the rotation $\mathbf{q}^1 = [-\sin(\pi/4) \ 0 \ 0 \cos(\pi/4)]^T$ provided in Eq. (66b).

At the other end of the rod ($s=L$), where the rotation is being imposed, the last node is fixed at $y=d=22 \text{ cm}$. Furthermore, the last element has to be rotated by $-\pi/2$ rad around the x-axis to re-orient $\mathbf{d}_3$ along the y-direction, as explained above, but we also need to superimpose a rotation $\Phi$ with respect to the y-axis, to simulate writhing (see Fig. 7(b)). Mathematically, these conditions translate into two constraints depending on the positional and rotational degrees of freedom alone, $\mathbf{r}^{N+1}$ and $\mathbf{q}_N$

$$C_r^L[\mathbf{r}^{N+1}] = \mathbf{r}^{N+1} - [0 \ d \ 0]^T = \mathbf{0} \qquad (67a)$$

$$C_q^L[\mathbf{q}^N] = \mathbf{q}^N - \frac{\sqrt{2}}{2} [-\cos(\Phi/2) \sin(\Phi/2) \sin(\Phi/2) \cos(\Phi/2)]^T = \mathbf{0}. \qquad (67b)$$

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Schematics of the rotational boundary conditions of our writhing experiment. (a) At $s=0$, we need to orient the first element along the y-direction. To do so, we need to impose a rotation of $-\pi/2$ rad around the x-axis and (b) at $s=L$, we impose two rotations. We orient the Nth element in the y-direction as in (a) and superimpose a rotation $\Phi$, the control parameter, around the y-axis.</figcaption></figure>
Again, Eq. (67a) is a direct translation of the positional boundary conditions of our writhing experiment. Eq. (67b), however, is less intuitive; it is the quaternion representation of the composition of the two rotations shown in Fig. 7(b). The first rotation, A, has already been treated above and can be described by the quaternion vector $\mathbf{q}_A^N = [-\sin(\pi/4) 0 0 \cos(\pi/4)]^T$. The second re-orientation, B, can be represented by a rotation angle of $\Phi$ around the unit length vector $\mathbf{b} = [0 1 0]^T$ and translates as $\mathbf{q}_B^N = [0 \sin(\Phi/2) 0 \cos(\Phi/2)]^T$ in the quaternion basis. The total rotation imposed to the Nth element is therefore the composition of the two rotations A, then B. In terms of quaternions, this rotation $\mathbf{q}^N$ is represented by the multiplication of the two sets of quaternions $\mathbf{q}_A^N$ and $\mathbf{q}_B^N$ and reads $\mathbf{q}^N = \mathbf{q}_A^N \cdot \mathbf{q}_B^N$ which, following the multiplication rule in the quaternion basis (Kuipers, 1999), is given in Eq. (67b).

To properly account for the boundary conditions in our numerical model introduced in Section 4, in addition to the constraints set by Eqs. (66) and (67), we also need to include the quantities $\mu_r^0 \nabla C_r^0$, $\mu_r^L \nabla C_r^L$, $\mu_q^0 \nabla C_q^0$ and $\mu_q^L \nabla C_q^L$, into the quadratic vector of equilibrium equations $\mathbf{f}(\mathbf{u})$ given in Appendix A. Doing so involved the introduction of the Lagrange multipliers $\mu_r^0$, $\mu_r^L$, $\mu_q^0$ and $\mu_q^L$ in the vector of unknowns $\mathbf{u}$ of Eq. (48). Physically, $\mu_r^0$ and $\mu_r^L$ are the forces at each end, written in the Cartesian basis, required to impose the positional boundary conditions Eqs. (66a) and (67a) at equilibrium. Similarly, $\mu_q^0$ and $\mu_q^L$ are a set of quaternions representing the moments at each extremity, necessary to impose the rotational boundary conditions. Finally, these boundary conditions can be accounted for in the stability analysis by including the new Lagrange parameters and their associated constraints in the second-order condition given in Section 4.3 through Eq. (64).

Finally, before we are able to solve our modified nonlinear algebraic problem $\mathbf{f}(\mathbf{u}) = \mathbf{0}$, one last important step is needed. The condition given in Eq. (67b) is not quadratic in terms of the control parameter $\Phi$, and consequently, without further modification, the updated nonlinear vector valued function $\mathbf{f}(\mathbf{u})$ would not be adequate to the numerical framework posed in Section 4. Fortunately, there is an appropriate way in the ANM framework to quadratically recast Eq. (67b). The technique involves adding two new variables

$$c(a) = \cos(\Phi(a)/2) \qquad (68a)$$

$$s(a) = \sin(\Phi(a)/2), \qquad (68b)$$

into the vector of unknowns $\mathbf{u}(a)$. The procedure is based on the introduction of differential equations in terms of the path-parameter $a$ in $\mathbf{f}(\mathbf{u}(a))$. Differentiating Eq. (68) with respect to $a$, the unknowns $(\mathbf{q}^N, c, s)$ are now solutions of the quadratic algebraic system

$$\mathbf{0} = \mathbf{q}^N(a) - \frac{\sqrt{2}}{2} [-c(a) s(a) s(a) c(a)]^T, \qquad (69a)$$

$$0 = dc(a) + \frac{1}{2}s(a)d\Phi(a), \qquad (69b)$$

$$0 = ds(a) - \frac{1}{2}c(a)d\Phi(a). \qquad (69c)$$

To fully integrate Eq. (69) into the asymptotic numerical framework of Section 4.2, we now need to perform a minor modification to the identification technique of the power series explained in Eqs. (49)–(53). We recall that the fundamental idea behind the ANM is to express the vector of unknowns in a power series of $a$ such that, $\mathbf{u}(a) = \mathbf{u}_0 + \sum_{p=1}^{m} a^p \mathbf{u}_p$, whose differential version reads

$$d\mathbf{u}(a) = \mathbf{u}_1 + \sum_{p=2}^{m} p a^{p-1} \mathbf{u}_p. \qquad (70)$$

Substituting Eq. (70) into Eq. (69) and identifying the power of $a$ allow us to compute the contributions $\mathbf{u}_p$ of the semi-analytical vector of unknowns $\mathbf{u}(a)$, following the same procedure described in Section 4.2. We highlight the fact that this method of introducing differential equations in the ANM method is a convenient way to represent complex non-polynomial energy functions in our numerical model, e.g. to represent highly nonlinear phenomena such as contact forces (Karkar et al., 2013).

The kinematic boundary conditions, Eqs. (66) and (67), and rotational control parameter $\Phi$, Eqs. (50)–(68), for the writhing problem are all now correctly implemented into our algebraic equilibrium equations $\mathbf{f}(\mathbf{u}(a)) = \mathbf{0}$, where $\mathbf{f}(\mathbf{u})$ is given in Eq. (A.1). We can now compute the vector of unknowns, $\mathbf{u}(a)$, as an asymptotic expansion in terms of the path-parameter, $a$, using the method explained in Section 4, to perform the continuation of the solutions $\mathbf{u}(a)$ and assess their associated local stability.

## 5.3. Comparison between numerics and experiments

Having introduced, developed and described our theoretical and computational tools, we proceed by performing a direct comparison between numerics and experiments. In particular, we focus on quantifying the evolution of the equilibrium configurations and associated buckling instabilities, as a function of the control parameter, $\Phi$. We highlight that in this comparison, there are no fitting parameters; all material and geometric parameters of the experiments are independently measured and considered as input variables into the numerics.

In Figs. 8, 10 and 12b, we compare the top view of some representative experimental and numerical equilibrium shapes for a straight rod ($\hat{k}_1(s) = 0 \text{ m}^{-1}$) and a naturally curved rod ($\hat{k}_1(s) = 44.84 \text{ m}^{-1}$). From these images, we measure the
maximum transverse displacement of the rod, $X_{max}$, in the (x,y) plane (top view) as a function of the rotation angle, $\Phi$, which is treated as a control parameter. Experimentally, the quantity $X_{max}$ was measured from image analysis of the digital images taken by the camera located above the apparatus. Using these quantities, we then construct the bifurcation diagrams presented in Figs. 9(a), 11(a) and 12(a) (for experiments and numerics), for the straight and curved rods, respectively. We also analyze the stability of the equilibrium state by calculating the first eigenvalue of the stability problem as a function of the control parameter, $\Phi$, and the results are plotted in Figs. 9(a) and 11(a) (for numerics). In order to quantitatively validate our ANM continuation technique, the semi-analytical numerical curves (lines) are superposed onto the experimental results (data points, every $\Phi = 30^{\circ}$), for the same value of the control parameter. We highlight, once again, that there are no fitting parameters involved in this comparison; all quantities are measured in the experiments, independently from the numerics. The excellent quantitative agreement between experiments and numerics illustrates the sticking predictive power of our framework.

We now comment on the experimental and numerical results in more detail, focusing first on the case of the straight rod ($\dot{k}_1 = 0 \text{ m}^{-1}$), the results of which we plotted in Figs. 8 and 9. Initially, for $\Phi = 0^{\circ}$, the rod exhibits a planar equilibrium shape lying in the (y,z) plane due to the effect of gravity. This equilibrium configuration is calculated using a classic Newton-Raphson

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Top view of various equilibrium configurations of a straight elastic rod ($\dot{k}_1 = 0 \text{ m}^{-1}$) for increasing values of the rotation angle $\Phi$. The experimental pictures have a black background and the simulations have a white background. The simulation results are rendered to visualize twist by using bi-color rods. (a) Planar equilibrium shape at rest for $\Phi = 0^{\circ}$, (b) out-of-plane configuration for $\Phi = 720^{\circ}$, (c) out-of-plane configuration for $\Phi = 1470^{\circ}$ and (d) onset of formation of a plectoneme at the middle of the rod.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Bifurcation diagram of the straight elastic rod ($\dot{k}_1 = 0 \text{ m}^{-1}$). (a) Evolution of the maximum transverse displacement, $X_{max} = \max(\text{abs}(X))$ (defined in Fig. 8(a)), with the control parameter $\Phi$. Comparison between experimental results and the semi-analytical branches computed with MANlab. Solid/dashed lines represent stable/unstable branches, respectively and (b) evolution of the first eigenvalue of the stability problem with control parameter $\Phi$. The critical angle corresponding to the emergence of plectoneme is $\phi_{sim} = 2040^{\circ}$ for the simulations and $\phi_{exp} = 2025 \pm 15^{\circ}$ for the experiments.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Top view of various equilibrium configurations of a curvy elastic rod ($\hat{k}_1(s) = 44.84 \text{ m}^{-1}$) for increasing values of the rotating angle, $\phi$. The experimental pictures have a black background and the simulations have a white background. The simulation results are rendered to visualize twist by using bi-color rods. (a) Asymmetric out-of-plane equilibrium shape for $\phi = 150^\circ$, (b) one-twist-per-wave configuration for $\phi = 900^\circ$ with three wavelengths between the two clamps, (c) one-twist-per-wave configuration for $\phi = 1800^\circ$ with five wavelengths between the two clamps and (d) onset of formation of a plectoneme state at one extremity of the rod.</figcaption></figure>

algorithm (Riks, 1979) and taken as the initial fixed point $\mathbf{u}_0$, the solution of our equilibrium equations $\mathbf{f}(\mathbf{u}(a)) = \mathbf{0}$ (see Fig. 8(a)). When the rotation angle $\phi$ is increased, this initial planar shape evolves smoothly into an out-of-plane configuration, symmetric to the (y,z) plane, with an amplitude that grows due to an increasing internal twist (see Fig. 8(b), (c)]. At a critical value of the rotation angle, $\phi^c$, the out-of-plane shape loses stability and the rod buckles into a *plectoneme* state (Thompson and Champneys, 1996; VanderHeijden et al., 2003): a highly localized structure corresponding to a two-start right-handed helix with terminal loops. Beyond this point, our numerical model is no longer able to reproduce the rod's configurations since they involve self-contact which is not included in our description. To further quantify this process, in Fig. 9 we plot the maximum transverse displacement of the rod, $X_{\text{max}}$, as a function of the imposed rotation angle, $\phi$. Across the full range of $\phi$ explored, the experimental data plotted in Fig. 9 is in excellent quantitative agreement with the numerical prediction.

It is remarkable that the instability threshold for the formation of the plectoneme, $\phi^c$, is also well recovered by our local stability analysis showed in Fig. 9(b), where we plot the evolution of the first eigenvalue, $\sigma_1$, as a function of $\phi$. With no fitting parameters, the predicted critical threshold $\phi_{\text{sim}}^c = 2040^\circ$ is in excellent agreement (within 1%) with the experimental results $\phi_{\text{exp}}^c = 2025 \pm 15^\circ$. To plot the semi-analytical bifurcation diagram of Fig. 9, we computed 80 solution vectors, $\mathbf{u}(a)$, expressed in terms of power series expansions at the order $m=20$, as given in Eq. (50). Using a desktop computer with a standard processor (at the time of writing) of 2.71 Ghz and 2.75 Gb of RAM, the computation required 9 s to determine one asymptotic series. This represents a total running time of approximately 12 min to simulate the full problem, using MANlab. Note that this computational time is mostly due to the 80 inversions of the Jacobian matrix needed to solve the linear systems given in Eq. (53), which are of size $1415 \times 1415$ for $N=100$ elements. These computations could be made even more efficient by using a dedicated solver such as the ones offered by traditional finite difference codes but the time optimization of the ANM algorithm, which has been investigated (Cochelin et al., 1994), is beyond the scope of this paper.

Interestingly, for the case of the naturally curved rod ($\hat{k}_1(s) = 44.84 \text{ m}^{-1}$) the evolution of equilibrium configurations with the rotation angle $\phi$ is qualitatively different from the case of the naturally straight elastic rod, as shown in Figs. 10–12. Once again, the qualitative and quantitative agreement between the experimental and numerical results is remarkable. By introducing the natural curvature $\hat{k}_1(s)$, the previously symmetric out-of-plane solutions obtained for the case of straight rods, become asymmetric with respect to the (y,z) plane. For small rotation angles $\phi$, the initial planar shape exhibits an asymmetric out-of-plane configuration due to the competition between the imposed internal twist and the intrinsic twist naturally imposed by $\hat{k}_1(s)$ (see Fig. 10(a)). Above $\phi \approx 400^\circ$, our results confirm that the rod, jumps into a one-twist-per-wave mode due to the presence of natural curvature, as originally reported in Champneys et al. (1997). In this configuration, the number of waves is equal to the number of twists stored in the rod as shown in Fig. 10(b)-(c). For a critical rotation angle $\phi^c$, a plectoneme forms, superimposed onto the one-twist-per-wave equilibrium state which is no longer stable. It is interesting to note that, for the naturally curved rod, the plectoneme is located at one extremity of the rod rather than at its center (see Fig. 10(d)), as found above for the straight rod.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Bifurcation diagram of the curved elastic rod ($\dot{k}_1(s) = 44.84 \text{ m}^{-1}$). (a) Evolution of the maximum transverse displacement with the control parameter $\Phi$. Comparison between experimental results and the semi-analytical branches computed with MANlab. The presence of a constant natural curvature introduces a new wavy configuration with one-twist-per-wave. Solid/dashed lines represent stable/unstable branches, respectively and (b) evolution of the first eigenvalue of the stability problem with control parameter $\Phi$. The critical angle corresponding to the emergence of plectoneme is $\phi_{sim}^c = 2907^\circ$ for the simulations and $\phi_{exp}^c = 2945 \pm 15^\circ$ for the experiments.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Hysteretic behavior of the curved elastic rod ($\dot{k}_1(s) = 44.84 \text{ m}^{-1}$). (a) Evolution of the maximum transverse displacement with the control parameter $\Phi$. Comparison between experimental results obtained for increasing and decreasing $\Phi$ and the semi-analytical branches computed with MANlab. Solid/dashed lines represent stable/unstable branches, respectively and (b) the equilibrium configuration for $\Phi = 150^\circ$ (two helices with opposite handedness) is different from the configuration in Fig. 10(a) for the same control parameter $\Phi$, which is significant of multi-stability.</figcaption></figure>

For the naturally curved rod, our continuation method is able to robustly and efficiently follow the equilibrium branches across the full range of considered rotation angles, $\Phi$, exhibiting excellent agreement with the experimental data, as shown in Fig. 11. For this second test case with $\dot{k}_1(s) = 44.84 \text{ m}^{-1}$, we used $N=200$ elements and we computed 115 asymptotic expansions $u(a)$ to determine the full semi-analytical bifurcation diagram, leading to a computational time of approximately 35 minutes. Again, the instability threshold for the onset of a plectoneme is well recovered by our local stability analysis showed in Fig. 11(b). In this case of a more complicated bifurcation diagram, it is remarkable that the predicted $\phi_{sim}^c = 2907^\circ$ is within 2% of the experimentally measured value of $\phi_{exp}^c = 2945 \pm 15^\circ$. Counterintuitively, we find that imparting a constant natural curvature to our rods (essentially adding a finite imperfection to the stress-free configuration) results in postponing, by approximately 44%, the emergence of the plectoneme instability (often synonymous with failure in practical systems). To our knowledge, this interesting novel phenomenon has thus far been overlooked in the literature and deserves further investigation. A systematic study to quantify and rationalize the influence of natural curvature on the writhing of a slender elastic rod is beyond the scope of this paper, but is an aspect which we plan to investigate in future work.

Another interesting feature in the writhing of the naturally curved rod is the hysteretic behavior observed for small values of the rotation angle $\Phi$, before entering the one-twist-per-wave regime, as illustrated in Fig. 12. Initially, when we increased $\Phi$, the material twist originally stored due to the natural curvature $\dot{k}_1(s)$ is released until the rod jumps to the one-twist-per-wave mode at $\Phi \approx 400^\circ$. If we then decrease $\Phi$, the observed equilibrium configurations are not the asymmetric out-of-plane shapes we previously encountered as illustrated in Fig. 10(a). Instead, Fig. 12(b) shows an inversion of helix handedness, known as perversion (Goriely and Tabor, 1998b) and described as two helices with opposite handedness. If we decrease $\Phi$ even further to negative values, the rod jumps back on the previous stable equilibrium branches where the rod twists with the same handedness. In the regime of small rotations ($\Phi < 400^\circ$), our system is metastable; for the same control parameter $\Phi$, the rod can exhibit two different configurations depending on the loading path. Once again, our continuation method correctly predicts the different equilibrium states and stability threshold of this complex hysteretic behavior as shown
in Fig. 12. Moreover, this highly nonlinear feature emphasize the ability of our numerical technique to follow the equilibrium branches and stability of slender elastic rods independently of the complexity of the bifurcation diagram.

## 6. Conclusions and perspectives

We have presented an original theoretical and computational framework to follow the geometrically nonlinear 3D equilibrium shapes and stability of slender elastic rods. In our model, we account for the elastic energy due to changes of material curvatures and twist, as well as the work of external forces and moments, under the assumption that the rod is inextensible and unshearable. The main result of the paper lies in the combination of using quaternions with the Asymptotic Numerical Method (ANM) and stability analysis. Moreover, we have challenged and validated our computational framework by comparing numerical results against our own experiments on the writhing of a thin rod. Excellent quantitative agreement is found between the two. We were able to successfully and accurately calculate the geometrically-nonlinear configurations of the rods, as well as the critical thresholds for instability. We also note that our continuation algorithm is able to address regions of multi-stability and hysteresis, as in the regime of low rotation angles, when increasing or decreasing the control parameter.

A potential extension of this work would be to incorporate other additional mechanical ingredients into our model such as internal stretching, hydrostatic loading and contact forces arising, for instance, due to self-contact or when the rod interacts with external boundaries. One technical requirement in order to be able to introduce new energy terms in our description is that the resulting equilibrium equations have to be quadratic to match the present ANM framework. However, we have shown that even some cases of non-polynomial functions can easily be reduced to a quadratic form by introducing a limited number of new variables in the vector of unknowns in a process we call recasting. Otherwise, the introduction of these new features can be readily accomplished, as long as they derive from a potential energy since continuation methods only apply for conservative systems, where an equilibrium can be found.

Our novel method is relatively simple to implement, robust, accurate, flexible and computationally efficient. We hope that this technique will be invaluable in problems that demand the predictive understanding of complex mechanical phenomena intrinsic to the extreme deformation of slender elastic rods, whose timely revival is relevant in a variety of currently open problems in both nature and technology.

## Acknowledgments

We thank Basile Audoly for introducing us the breathtaking world of quaternions and Matt Metlitz for help with the writhing experiments. We are grateful to the support by the National Science Foundation (CMMI-1129894) and Schlumberger-Doll Research. Arnaud Lazarus acknowledges funding from a Battelle-MIT postdoctoral fellowship.

## Appendix A. Quadratic form of the vector of equilibrium equations

The quadratic form of the 14N-dimensional nonlinear vector valued function $f(u)$ given in Eq. (49) and representing the equilibrium equations (43)–(45) of our inextensible and unshearable slender elastic rod, can be written as

$$ f(u) = \begin{cases} \mathbf{P}^0 + \mu^i + \mathbf{0} & \text{for } i=1 \\ \frac{\mathbf{p}^i L}{N} + \mu_i - \mu^{i-1} + \mathbf{0} & \text{for } i=[1:N] \\ \mathbf{P}^L - \mu^{i-1} + \mathbf{0} & \text{for } i=N+1 \\ \mathbf{M}^0 + \mathbf{0} + 2 \sum_{k=1}^{3} \mathbf{B}_k \mathbf{q}^{i+1} G_k^i + (\mathbf{D}(\mathbf{q}^i)\mu^i - \alpha^i \mathbf{q}^i) & \text{for } j=1 \\ \mathbf{m}^i + \mathbf{0} + 2 \sum_{k=1}^{3} \mathbf{B}_k \mathbf{q}^{i-1} G_k^{i-1} - 2 \sum_{k=1}^{3} \mathbf{B}_k \mathbf{q}^{i+1} G_k^i + 2(\mathbf{D}(\mathbf{q}^i)\mu^i - \alpha^i \mathbf{q}^i) & \text{for } j=[1:N-1] \\ -\mathbf{M}^L + \mathbf{0} + 2 \sum_{k=1}^{3} \mathbf{B}_k \mathbf{q}^{i-1} G_k^{i-1} + 2(\mathbf{D}(\mathbf{q}^i)\mu^i - \alpha^i \mathbf{q}^i) & \text{for } j=N \\ 1 + \mathbf{0} - \mathbf{q}^i \mathbf{q}^i & \text{for } j=[1:N] \\ \mathbf{0} + \mathbf{r}^{i+1} - \mathbf{r}^i - L/N \mathbf{d}_3^i(\mathbf{q}^2) & \text{for } j=[1:N] \\ E_1 I_1 \bar{\kappa}_1^j + G_1^j L/N - E_1 I_1 \mathbf{B}_1(\mathbf{q}^i + \mathbf{q}^{i+1}).(\mathbf{q}^{i+1} - \mathbf{q}^i) & \text{for } j=[1:N-1] \\ E_2 I_2 \bar{\kappa}_2^j + G_2^j L/N - E_2 I_2 \mathbf{B}_2(\mathbf{q}^i + \mathbf{q}^{i+1}).(\mathbf{q}^{i+1} - \text{vect}\mathbf{q}^i) & \text{for } j=[1:N-1] \\ E_3 I_3 \bar{\kappa}_3^j + G_3^j L/N - E_3 I_3 \mathbf{B}_3(\mathbf{q}^i + \mathbf{q}^{i+1}).(\mathbf{q}^{i+1} - \mathbf{q}^i) & \text{for } j=[1:N-1] \end{cases} \quad (A.1) $$

where the expression of the director $d_3^i(\mathbf{q}^2)$, the projection matrix $\mathbf{D}(\mathbf{q}^i)$ and the unknown variables $u$ can be found in Eqs. (33), (42) and (48), respectively. Note that these algebraic equations are the general form of the equilibrium equations of the inextensible elastic rod under external forces and moments and do not account for particular kinematic boundary conditions or control parameter which can vary depending on the problem under study. An example of the additional
equations needed to impose some fixed boundary conditions and a rotation angle as control parameter for the case of writhing of a rod is given in Section 5.

## References

* Allgower, E.L., Georg, K., 2003. Introduction to Numerical Continuation Methods, vol. 45. Society for Industrial Mathematics.

* Altmann, S.L., 1986. Rotations, Quaternions and Double Groups, vol. 3. Clarendon Press, Oxford.

* Arquier, R., 2007. Une méthode de calcul des modes de vibrations non linéaires de structures. Ph.D. Thesis. Université de la Méditerranée (Aix-Marseille II).

* Arribas, M., Elope, A., Palacios, M., 2006. Quaternions and the rotation of a rigid body. Celestial Mech. Dynamical Astron. 96 (3), 239–251.

* Audoly, B., Pomeau, Y., 2010. Elasticity and Geometry: From Hair Curls to the Nonlinear Response of Shells. Oxford Press.

* Baguet, S., Cochelin, B., 2003. On the behaviour of the ANM continuation in the presence of bifurcations. Commun. Numer. Methods Eng. 19 (6), 459–471.

* Balaeff, A., Mahadevan, L., Schulten, K., 2006. Modeling DNA loops using the theory of elasticity. Phys. Rev. E 73 (3), 031919.

* Bažant, Z.P., Cedolin, L., 2010. Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories. World Scientific Publishing Company Incorporated.

* Bergou, M., Wardetzky, M., Robinson, S., Audoly, B., Grinspun, E., 2008. Discrete elastic rods. In: ACM Transactions on Graphics (TOG), vol. 27. ACM, pp. 63.

* Bertails, F., Audoly, B., Cani, M.P., Querleux, B., Leroy, F., Lévêque, J.L., 2006. Super-helices for predicting the dynamics of natural hair. In: ACM Transactions on Graphics (TOG), vol. 25. ACM, pp. 1180–1187.

* Champneys, A.R., VanderHeijden, G.H.M., Thompson, J.M.T., 1997. Spatially complex localization after one-twist-per-wave equilibria in twisted circular rods with initial curvature. Philos. Trans. R. Soc. London Ser. A: Math. Phys. Eng. Sci. 355 (1732), 2151–2174.

* Chou, J.C.K., Kamel, M., 1991. Finding the position and orientation of a sensor on a robot manipulator using quaternions. Int. J. Robotics Res. 10 (3), 240–254.

* Cochelin, B., 1994. A path-following technique via an asymptotic-numerical method. Comput. Struct. 53 (5), 1181–1192.

* Cochelin, B., Damil, N., Potier-Ferry, M., 1994. Asymptotic-numerical methods and Padé approximants for non-linear elastic structures. Int. J. Numer. Methods Eng. 37 (7), 1187–1213.

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

* Crisfield, M.A., 1991. Nonlinear Finite Element Analysis of Solids and Structures. Essentials, vol. 1. Wiley, New York, NY (United States).

* Dam, E.B., Koch, M., Lillholm, M., 1998. Quaternions, Interpolation and Animation. Datalogisk Institut, Københavns Universitet.

* Damil, N., Potier-Ferry, M., 1990. A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures. Int. J. Eng. Sci. 28 (9), 943–957.

* Davies, M.A., Moon, F.C., 1993. 3D spatial chaos in the elastica and the spinning top: Kirchhof analogy. Chaos 3, 93.

* Dichmann, D.J., Li, Y., Maddocks, J.H., 1996. Hamiltonian formulations and symmetries in rod mechanics. IMA Volumes Math. Appl. 82, 71–114.

* Dill, E.H., 1992. Kirchhoff's theory of rods. Arch. Hist. Exact Sci. 44 (1), 1–23.

* Doedel, E.J., 1981. Auto: a program for the automatic bifurcation analysis of autonomous systems. Congressus Numerantium 30, 265–284.

* Furrer, P.B., Manning, R.S., Maddocks, J.H., 2000. DNA rings with multiple energy minima. Biophys. J. 79 (1), 116–136.

* Goriely, A., Tabor, M., 1997a. Nonlinear dynamics of filaments. III. Instabilities of helical rods. Proc. R. Soc. London Ser. A: Math. Phys. Eng. Sci. 453 (1967), 2583–2601.

* Goriely, A., Tabor, M., 1997b. Nonlinear dynamics of filaments I. Dynamical instabilities. Phys. D: Nonlinear Phenom. 105 (1), 20–44.

* Goriely, A., Tabor, M., 1997c. Nonlinear dynamics of filaments II. Nonlinear analysis. Phys. D: Nonlinear Phenom. 105 (1), 45–61.

* Goriely, A., Tabor, M., 1998a. Nonlinear dynamics of filaments. IV Spontaneous looping of twisted elastic rods. Proc. R. Soc. London Ser. A: Math. Phys. Eng. Sci. 454 (1980), 3183–3202.

* Goriely, A., Tabor, M., 1998b. Spontaneous helix hand reversal and tendril perversion in climbing plants. Phys. Rev. Lett. 80 (7), 1564–1567.

* Goss, V.G.A., vanderHeijden, G.H.M., Thompson, J.M.T., Neukirch, S., 2005. Experiments on snap buckling, hysteresis and loop formation in twisted rods. Exp. Mech. 45 (2), 101–111.

* Goyal, S., Perkins, N.C., Lee, C.L., 2008. Non-linear dynamic intertwining of rods with self-contact. Int. J. Non-Linear Mech. 43 (1), 65–73.

* Guckenheimer, J., Holmes, P., 1983. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, vol. 42. Springer-Verlag, New York.

* Hamilton, W.R., 1847. On quaternions. In: Proceedings of the Royal Irish Academy, vol. 3, pp. 1–16.

* Hamilton, S.W.R., 1853. Lectures on Quaternions. Hodges and Smith.

* Hanson, A.J., 2005. Visualizing Quaternions. Morgan Kaufmann Pub.

* Healey, T.J., Mehta, P.G., 2005. Straightforward computation of spatial equilibria of geometrically exact cosserat rods. Int. J. Bifurcation Chaos Appl. Sci. Eng. 15 (3), 949–966.

* Horn, B.K.P., 1987. Closed-form solution of absolute orientation using unit quaternions. J. Opt. Soc. Am. A 4 (4), 629–642.

* Karkar, S., Arquier, R., Lazarus, A., Thomas, O., Vergez, C., Cochelin, B., 2010. Manlab: an interactive path-following and bifurcation analysis software. http://manlab.lma.cnrs-mrs.fr/, November 2010.

* Karkar, S., Cochelin, B., Vergez, C., 2013. A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: the case of non-polynomial nonlinearities. J. Sound Vib. 332 (4), 968–977.

* Kehrbaum, S., Maddocks, J.H., 1997. Elastic rods, rigid bodies, quaternions and the last quadrature. Philos. Trans: Math. Phys. Eng. Sci. 2117–2136.

* Koiter, W.T., 1970. The stability of elastic equilibrium. Technical Report, DTIC Document, 1970.

* Kuipers, J.B., 1999. Quaternions and Rotation Sequences. Princeton University Press, Princeton, NJ, USA.

* Lazarus, A., Thomas, O., 2010. A harmonic-based method for computing the stability of periodic solutions of dynamical systems. Comptes Rendus Mécanique 338, 510–517.

* Lazarus, A., Barois, T., Perisanu, S., Poncharal, P., Manneville, P., deLangre, E., Purcell, S.T., Vincent, P., Ayari, A., 2010a. Simple modeling of self-oscillation in NEMS. Appl. Phys. Lett. 96 (193114).

* Lazarus, A., de Langre, E., Manneville, P., Vincent, P., Perisanu, S., Ayari, A., Purcell, S., 2010b. Statics and dynamics of a nanowire in field emission. Int. J. Mech. Sci. 52 (11), 1396–1406.

* Levien, R., 2008. The Elastica: A Mathematical History. Electrical Engineering and Computer Sciences University of California at Berkeley.

* Liu, F.C., Center, N.C.B., 1975. Kink Formation and Rotational Response of Single and Multistrand Electromechanical Cables. Defense Technical Information Center.

* Luenberger, D.G., 1973. Introduction to Linear and Nonlinear Programming. Addison-Wesley Publishing Company.

* Luenberger, D.G., Ye, Y., 2008. Linear Nonlinear Programming, vol. 116. Springer.

* Manning, R.S., Rogers, K.A., Maddocks, J.H., 1998. Isoperimetric conjugate points with application to the stability of DNA minicircles. Proc. R. Soc. London Ser. A: Math. Phys. Eng. Sci. 454 (1980), 3047–3074.

* Marko, J.F., Neukirch, S., 2012. Competition between curls and plectonemes near the buckling transition of stretched supercoiled DNA. Phys. Rev. E 85 (1), 011908.

* Neukirch, S., VanderHeijden, G.H.M., Thompson, J.M.T., 2002. Writhing instabilities of twisted rods: from infinite to finite length. J. Mech. Phys. Solids 50 (6), 1175–1191.

* Riks, E., 1979. An incremental approach to the solution of snapping and buckling problems. Int. J. Solids Struct. 15 (7), 529–551.

* Spillmann, J., Harders, M., 2010. Inextensible elastic rods with torsional friction based on lagrange multipliers. Comput. Animation Virtual Worlds 21 (6), 561–572.

* Spillmann, J., Teschner, M., 2009. Cosserat nets. IEEE Trans. Visualization Comput. Graphics 15 (2), 325–338.
Sun, Y., Choi, W.M., Jiang, H., Huang, Y.Y., Rogers, J.A., 2006. Controlled buckling of semiconductor nanoribbons for stretchable electronics. Nat. Nanotechnol. 1 (3), 201–207.

Thompson, J.M.T., Champneys, A.R., 1996. From helix to localized writhing in the torsional post-buckling of elastic rods. Proc. R. Soc. London Ser. A: Math. Phys. Eng. Sci. 452 (1944), 117–138.

Tweed, D., Cadera, W., Vilis, T., 1990. Computing three-dimensional eye position quaternions and eye velocity from search coil signals. Vision Res. 30 (1), 97–110.

VanderHeijden, G.H.M., Thompson, J.M.T., 2000. Helical and localised buckling in twisted rods: a unified analysis of the symmetric case. Nonlinear Dyn. 21 (1), 71–99.

VanderHeijden, G.H.M., Neukirch, S., Goss, V.G.A., Thompson, J.M.T., 2003. Instability and self-contact phenomena in the writhing of clamped rods. Int. J. Mech. Sci. 45 (1), 161–196.

Vannucci, P., Cochelin, B., Damil, N., Potier-Ferry, M., 1998. An asymptotic-numerical method to compute bifurcating branches. Int. J. Numer. Methods Eng. 41 (8), 1365–1389.

Waldvogel, J., 2008. Quaternions for regularizing celestial mechanics: the right way. Celestial Mech. Dyn. Astron. 102 (1), 149–162.

Wang, M., Alexander-Katz, A., Olsen, B.D., 2012. Diffusion of entangled rod-coil block copolymers. ACS Macro Lett. 1 (6), 676–680.

Wicks, N., Wardle, B.L., Pafitis, D., 2008. Horizontal cylinder-in-cylinder buckling under compression and torsion: review and application to composite drill pipe. Int. J. Mech. Sci. 50 (3), 538–549.

Yabuta, T., 1984. Submarine cable kink analysis. Bull. JSME 27 (231), 1821–1828.

Zahrouni, H., Cochelin, B., Potier-Ferry, M., 1999. Computing finite rotations of shells by an asymptotic-numerical method. Comput. Methods Appl. Mech. Eng. 175 (1), 71–85.

Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z., 2005. The Finite Element Method: Its Basis and Fundamentals, vol. 1. Butterworth-Heinemann.