A quaternion-based continuation method to follow the equilibria and stability of slender elastic rods

# A quaternion-based continuation method to follow the equilibria and stability of slender elastic rods

A. Lazarus J. T. Miller P. M. Reis Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge 02139, USA Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge 02139, USA
###### Abstract

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 simple 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 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.

###### 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 [1] and Kirchhoff’s kinetic analogy in [2] to the burgeoning of numerical approaches such as finite element-based methods in the late century [3], and the more recent algorithms based on discrete differential geometry [4]. 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 range from the supercoiling of DNA [5], self-assembly of rod-coil block copolymers [6], design of nano-electromechanical resonators [7, 8], development of stretchable electronics [9], computed animation of hairs [10] and coiled tubing operations in the oil-gas industries [11].

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 [12]. 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 [13] or standard minimization techniques [14] are often inappropriate since, depending on the initial guess, they may not converge towards the desired solution, or any solution. Addressing these computational difficulties calls for alternative numerical techniques, such as well-known continuation methods [15, 16]. Continuation techniques are based on coupling nonlinear algorithms (e.g. predictor-corrector [15] or perturbation methods [16]) 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 [13]. 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 [17]. Quasi-static deformations of slender elastic rods have been intensively studied using this software [18, 19], mostly due to the analogy between the rod’s equilibrium equations with the spinning top’s dynamic equations [20]. 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 1990’s [21, 22]. 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 [23, 24, 25].

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 nonlinear equilibrium equations, by minimizing its geometrically-constrained mechanical energy including internal bending and twisting energy, as well as the work of external forces and moments. Expressing the flexural and torsional internal moments in a quaternion basis yields differential equilibrium equations that are simply quadratic in terms of the unknowns. In Section 4, we proceed by presenting the numerical method developed to compute the equilibrium solutions. Using a finite element approach, 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 [26, 27]. In Section 5, we develop our own precision model experiment for the fundamental problem of the writhing of a clamped elastic rod [18, 28, 29, 30] and challenge our numerical model against experimental results. The 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 and the other length scale by . At large scales, the rod can be regarded as an adapted material curve: its centerline. If 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 in the reference configuration,

 r(s)=[rx(s)ry(s)rz(s)]T=[x(s)y(s)z(s)]T. (1)

We consider unstretchable rods whose centerline remains inextensible upon deformation. As explained in detail in [12], this assumption is physically justified for a wide range of loading conditions, provided that the aspect ratio of the rod, , is small. Under this assumption, the variable 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 in the deformed configuration. At each particular location , we associate an orthonormal basis , attached to the centerline. The centerline, together with this set of material frames, form what is called a Cosserat curve [31]. We choose the orientation of these material frames in a way such that the directors and lie in the plane of the cross-section, while the third director is always parallel to the tangent of the curve [see Fig. 1(a)]. Considering the case of small strains, the triad 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 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, which we now express using the framework of differential geometry of curves in 3D space [12]. The previous condition of orthonormality (Euler-Bernoulli assumption) yields the relations,

 d′i(s).di(s)=0andd′i(s).dj(s)=−d′j(s).di(s), (2)

for all indices and varying from to (there is no implicit sum over in the first equation) and where denotes differentiation with respect to . The most general set of first-order linear equations conserving the orthonormal character of the material frame represented in Eq. (2) can be expressed as,

 d′1(s)=κ3(s)d2(s)−κ2(s)d3(s), (3a) d′2(s)=−κ3(s)d1(s)+κ1(s)d3(s), (3b) d′3(s)=κ2(s)d1(s)−κ1(s)d2(s). (3c)

where , and are scalar functions interpreted below. These equations describe the rigid-body rotation of the frame. Using the notation for the cross product of two vectors, Eqs. (3) can be rewritten as,

 d′1(s)=Ω(s)×d1(s),d′2(s)=Ω(s)×d2(s),d′3(s)=Ω(s)×d3(s), (4)

where we have introduced the Darboux vector ,

 Ω(s)=κ1(s)d1(s)+κ2(s)d2(s)+κ3(s)d3(s). (5)

The physical interpretation of Eqs. (4) is that the material frame rotates with a rotation velocity, , when following the centerline at unit speed. The quantities and in Eq. (5), called the material curvatures, illustrated in Fig. 1(b1)-(b2), represent the extent of rotation of the material frame, with respect to the directions and of the cross-section. The quantity quantifies the rotation of the material frame with respect to the tangent , 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, , and the rotation matrix of the Euclidean 3D space , this rotated Darboux vector is,

 κ(s)=RT(s)Ω(s), (6)

or, in terms of the directors ,

 κk(s)=dk(s).Ω(s), (7)

since the directors constitute the columns of the rotation matrix .

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 , and . 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 [12]. 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 or the directors, , and the strain rate vector, , 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 [32]. They were first described by Hamilton in 1843 [33, 34] and were extensively used in many physics and geometry problems before loosing prominence in the late century following the development of numerical analysis. Quaternions were then revived in the late 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 [35], computer graphics [36, 37], optics [38, 39], robotics [40] and orbital mechanics [41, 42]. It is beyond the scope of this article to discuss a detailed evaluation of the advantages and disadvantages of using quaternions over other rotation parametrizations. 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 , is,

 i2=j2=k2=ijk=−1, (8)

where , , and are the basis elements of . A quaternion number in is written in the form where the imaginary part is an element of the vector space and the real part is a scalar. Using the basis , , , of makes it possible to write a quaternion as a set of quadruples, usually expressed as a vector in ,

 q=[q1q2q3q4]T. (9)

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=[bxsin(Φ/2)bysin(Φ/2)bzsin(Φ/2)cos(Φ/2)]T, (10)

where is the Euler principal angle and is the unit length principal vector such that [see Fig. 2]. Given that four Euler parameters are needed to define a three-dimensional rotation, a natural constraint equation prescribing that is indeed a unit quaternion follows from Eq. (10),

 q21+q22+q23+q24=1. (11)

The orthogonal matrix representation corresponding to a rotation by the quaternion with is,

 R(q)=⎡⎢ ⎢⎣q21−q22−q23+q242(q1q2−q3q4)2(q1q3+q2q4)2(q1q2+q3q4)−q21+q22−q23+q242(q2q3−q1q4)2(q1q3−q2q4)2(q2q3+q1q4)−q21−q22+q23+q24⎤⎥ ⎥⎦. (12)

Returning to the context of a thin elastic rod discussed above, its material frame remains orthonormal upon deformations and its rigid-body re-orientation can be expressed by the rotation matrix given in Eq. (12) such that,

 R(q(s))=[d1(q(s))d2(q(s))d3(q(s))]. (13)

The local frame is now parametrized in terms of the curvilinear unit quaternion coordinates vector along the slender rod.

We proceed by relating the strain rate vector of Eq. (6) to the Euler parameters . Multiplying each of the three geometric relations given in Eqs. (3) by the relevant director yields expressions for the material curvatures and twist in terms of the directors alone,

 κ1(s)=−d2(s).d′3(s),κ2(s)=−d3(s).d′1(s),κ3(s)=−d1(s).d′2(s). (14)

To compute in terms of quaternions, we note that is a function of , which is itself a function of the curvilinear coordinate . Upon employing the chain rule of partial differentiation, we obtain

 d′k(s)=d′k(q(s))=∂dk(q(s))∂s=∂dk∂q∂q∂s=Jk(q(s))q′(s) (15)

for the three directions , where is the Jacobi matrix . Replacing and by their respective expressions (12)-(13) and (15) in Eq. (14) allows us to express the material curvatures, and , and the twist, , solely in terms of the unit quaternions,

 κk(s)=2Bkq(s)q′(s)for k=1,2,3, (16)

 B1=⎡⎢ ⎢ ⎢⎣000100100−100−1000⎤⎥ ⎥ ⎥⎦,B2=⎡⎢ ⎢ ⎢⎣00−10000110000−100⎤⎥ ⎥ ⎥⎦,B3=⎡⎢ ⎢ ⎢⎣0100−1000000100−10⎤⎥ ⎥ ⎥⎦. (17)

With the new expression of in Eq. (16), we have been able to write the total strains in terms of the locally perturbable variables , 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 , , and are not geometrically independent. First, to represent a three-dimensional rotation with four coordinates, the unit quaternion assumption given in Eq. (11), must be verified. Secondly, whereas thus far we have treated the centerline position and the orientations as separate entities, the positions and the orientations cannot be considered independently. Indeed, the material frames parametrized by and are coupled by the constraint that the third director is always parallel to the tangent ,

 r′(s)=d3(q(s)), (18)

where is the unit tangent vector to the Cosserat curve and along the centerline since we assumed inextensibility. The three constraints set by Eq. (18) 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. (16), which links the strain rates to the local orientation of the material frame, together with the four geometrical constraints given in Eq. (11) and Eq. (18). For the remainder of this article, the three positions and four quaternion coordinates constitute the seven degrees of freedom of our slender elastic rod [31, 43]. 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 [12]. 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 and a constant cross-section reads,

 Ee=EI12∫L0(κ1(s)−^κ1(s))2ds+EI22∫L0(κ2(s)−^κ2(s))2ds+GJ2∫L0(κ3(s)−^κ3(s))2ds, (19)

where we used the previously defined rotational strain rate vector , and where the quantities , and are the intrinsic natural curvature and twist of the rod along the directors , and , respectively. In this expression, is the Young’s modulus of the material and is the shear modulus of the material with Poisson’s ratio . The constants and are the moments of inertia along the principal directions of curvature in the plane of the cross-section and and is the moment of twist which, similarly to and for the bending energy, depends only of the geometry of the cross-section. Replacing the material curvatures , and twist by their expression given in Eq. (16) allows us to write the elastic energy in a more compact form, in term of the rotational degrees of freedom alone,

 Ee(q(s))=3∑k=1EkIk2∫L0(2Bkq(s)q′(s)−^κk(s))2ds, (20)

where , and .

### 3.2 Variation of the energy

We now follow a variational approach for the elastic energy in Eq. (20), and consider an infinitesimal perturbation from an arbitrary configuration of the rod. The perturbed quantities are preceded by . Carrying out the first variation of Eq. (20), the corresponding variation of the energy is,

 δEe=3∑k=1EkIk∫L0(2Bkqq′−^κk)(2Bkδqq′+2Bkqδq′)ds, (21)

where is the vector of the arbitrary perturbations of the rotational degrees of freedom . Upon integration by parts, we transform Eq. (21) into an integral that depends on alone to arrive at,

 (22)

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 parametrized from to . 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 where,

 T(s)=G1(s)2B1q+G2(s)2B2q+G3(s)2B3q. (23)

The vector 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 , and given by,

 Gk(s)=EkIk(κk(s)−^κk(s))=EkIk(2Bkq(s),q′(s)−^κk(s)) (24)

are respectively the two flexural and torsional moments, defined as the components of in the local material frame. The second term in Eq. (22) 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 where as a four-dimensional vector written in the quaternion basis that reads,

 τ(s)=3∑k=1(Gk(s)4Bkq′(s)+G′k(s)2Bkq(s)), (25)

where is the differential of with respect to . The quantity is the net moment applied on an infinitesimal element of the rod located between the cross-sections at and .

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 and torques that are applied at the two ends and , and distributed forces and torques that are applied along the length of the rod, with linear densities and , respectively. The density of forces, , can represent, for instance, the weight of the rod, and the density of moments, , 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,

 δW=P(0)δr(0)+M(0)δq(0)+P(L)δr(L)+M(L)δq(L)+∫L0(p(s)δr(s)+m(s)δq(s))ds, (26)

where is the vector of the small arbitrary perturbations of the translational degrees of freedom . According to Eq. (26), the external forces and are defined in terms of the global directions , , whereas the external moments and are expressed in the quaternion basis , , , of , defined by Eq. (8). 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 can be chosen freely. This is, however, not the case since our rod is subject to the kinematical constraints introduced previously in Eqs. (11) and (18). 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 and external loads . In this Lagrangian formalism, the enforcement of the unicity of quaternion in Eq. (11), translates as the continuous functional constraint,

 Cα[q(s)]=q(s)q(s)−1=0, (27)

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

 Cμ[r(s),q(s)]=r′(s)−d3(q(s))=0. (28)

With the expressions for the energy of an arbitrary configuration of the rod in Eqs. (22) and (26) 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. (11) and (18). This is equivalent to requiring that the first order variation of the functionals and , combined linearly with the variation of the constraints and over the interval from to (i.e. the Lagrangian) vanish,

 δELag=δEe−δW+∫L0α(s)δCα(s)ds+∫L0μ(s)δCμ(s)ds=0. (29)

In this equation, the variation of the constraint given in Eq. (27) takes the form,

 ∫L0α(s)δCα(s)ds=∫L02αqδqds, (30)

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

 ∫L0μ(s)δCμ(s)ds=[μδr]L0−∫L0μ′δr+2D(q)μδqds, (31)

where the terms of the vector valued function are the Lagrange multipliers ensuring the condition of inextensibility of the slender elastic rods and the operator reads,

 D(q(s))=⎡⎢ ⎢ ⎢ ⎢⎣q3(s)−q4(s)−q1(s)q4(s)q3(s)−q2(s)q1(s)q2(s)−q3(s)q2(s)−q1(s)−q4(s)⎤⎥ ⎥ ⎥ ⎥⎦. (32)

Now, substituting Eqs. (22), (26), (30) and (31) into the Lagrangian of Eq. (29), we arrive at the first variation of the geometrically constraint elastic energy of the slender elastic rod,

 δELag=[T(s)δq(s)+μ(s)δr(s)]L0−M(0)δq(0)−M(L)δq(L)−P(0)δr(0)−P(L)δr(L)−∫L0(p(s)+μ′(s))δr(s)ds−∫L0(τ(s)+m(s)−2α(s)q(s)+2D(q(s))μ(s))δq(s)ds. (33)

The condition that the variation in Eq. (33) must vanish for an arbitrary perturbations and yields the strong form of the equilibrium equations for our elastic rod as second-order differential equations,

 0 =p(s)+μ′(s) (34a) 0 =τ(s)−m(s)+2α(s)q(s)−2D(q(s))μ(s). (34b)

When projected along the three directions of the global cartesian frame , the vector equation Eq. (34a) yields a set of three differential equations that can be interpreted as the balance of forces. The vector of Lagrange multiplier 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 , we find that the element is submitted to the contact forces and from the neighboring elements, and to the external force . At equilibrium, the total forces is zero as described by Eq. (34a).

When projected along the four elements of the quaternion basis , the vector equation Eq. (34b) 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 .

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

 0 =q21(s)+q22(s)+q23(s)+q24(s)−1 (35a) 0 =r′x(s)−2q1(s)q3(s)−2q2(s)q4(s) (35b) 0 =r′y(s)−2q2(s)q3(s)+2q1(s)q4(s) (35c) 0 =r′z(s)+q21(s)+q22(s)−q23(s)−q24(s). (35d)

In the seven differential equilibrium equations of Eqs. (34) plus the four differential equations in Eqs. (35), the eleven unknowns are the four Lagrange multipliers , , and , the four rotational degrees of freedom , , and and the three translational degrees of freedom , and . 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 given in Eq. (23), which is cubic in . 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. (33), 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 and yield,

 0 =(T(0)+M(0))δq(0), (36a) 0 =(T(L)−M(L))δq(L). (36b) The remaining boundary terms associated with displacements δr(0) and δr(L) of the ends s=0 and s=L, respectively, yield, 0 =(μ(0)+P(0))δr(0), (36c) 0 =(μ(L)−P(L))δr(L). (36d)

To provide a physical interpretation of the behavior at the boundary conditions, we first consider Eq. (36a). If the endpoint is free to rotate, the vector is arbitrary and one is led to the boundary condition . This is the total torque applied on the section , which is the sum of the internal moments transmitted by the downstream part of the rod, , and of the moment applied by the operator. At equilibrium, the total torque should vanish when the end is free to rotate. If the endpoint is fixed, the perturbations that are consistent with the kinematics are such that 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. (36b) near the opposite end, , although the total torque is now , since, in this case the internal moment is applied by the downstream part of the rod, .

The two other boundary conditions written in Eqs. (36c)-(36d) 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 near the end , by a similar reasoning as above. However, the total force is near the opposite end, , given that the internal forces are now applied by the downstream part of the rod, . This remark validates our previous interpretation as for the physical interpretation of the Lagrange multipliers , and ; 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. (34)-(36) 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 [13, 17]. 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 [21, 22]. Taking advantage of the particular cubic form of the geometrically-exact equilibrium Eqs (34)-(36), 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. (34)-(36) 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 (34)-(36), we first explain how to discretize the main function unknowns such as strain rate vector , material frame , positional and rotational degrees of freedom and or Lagrange multipliers and .

The position of the rod is represented by discretizing its centerline into elements separated by spatial control points, in , located by the discrete curvilinear coordinate 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,

 r′(si)≈r(si+dsi)−r(si)dsi=ri+1−ridsi (37)

where is the length of the element and is the total length of the rod. To ensure the inextensibility condition of our rods, is constant upon deformation and the stretch along the centerline is forced to verify .

The orientations of the centerline elements are represented by employing material frames in where is the set of quaternions associated with each element . According to Eq. (12), the directors of the element are vectors in represented at the midpoints on the centerline segments [see Fig. 3(a)] such that,

 dj1=⎡⎢ ⎢ ⎢⎣qj1qj1−qj2qj2−qj3qj3+qj4qj42(qj1qj2+qj3qj4)2(qj1qj3−qj2qj4)⎤⎥ ⎥ ⎥⎦,dj2=⎡⎢ ⎢ ⎢⎣2(qj1qj2−qj3qj4)−qj1qj1+qj2qj2−qj3qj3+qj4qj42(qj2qj3+qj1qj4)⎤⎥ ⎥ ⎥⎦, dj3=⎡⎢ ⎢ ⎢⎣2(qj1qj3+qj2qj4)2(qj2qj3−qj1qj4)−qj1qj1−qj2qj2+qj3qj3+qj4qj4⎤⎥ ⎥ ⎥⎦. (38)

Replacing the quaternions functions by their discrete counterparts in the expression of strain rates given in Eq. (16), we can write the discrete material curvatures , and the twist expressing the extent of rotation around the directors, , and , between two successive elements [see Fig. 3(b)] in the form

 κjk=2Bk¯qjq′j% for k=1,2,3. (39)

In Eq. (39), we introduced the average and the spatial derivative of the rotational degrees of freedom of the element as,

 ¯qj=qj+1+qj2andq′j=qj+1−qjdsj, (40)

where , taking into consideration the rod’s inextensibility condition.

In a similar fashion, replacing the continuous function and its derivative by their discretized counterparts, and , respectively, given in Eq. (40), the first variation of elastic energy previously given in Eq. (21) can be approximated by a Riemann sum over the elements from to ,

 δEe≈3∑k=1N−1∑j=12Gjk(Bkqjδqj+1−Bkqj+1δqj)dsj. (41)

In this equation, the vectors are the discrete version of the perturbed rotational degrees of freedom and are associated with each element . The constants are an approximation of the flexural and torsional internal torques given in Eq. (24), which read

 Gjk=EkIk(Bk(qj+qj+1)1dsj(qj+1−qj)−^κjk), (42)

and are defined between two successive elements for .

Replacing and by their discrete counterparts at each element and each node , respectively, the variation of the work done by the external forces and torques previously given in Eq. (26), can be approximate by its discrete version,

 δW≈P0δr1+M0δq1+PLδrN+1+MLδqN+N∑i=2piδridsi+N−1∑j=2mjδqjdsj, (43)

where