Simple deformation measures for Discrete elastic rods and ribbons

The Discrete elastic rod method (Bergou et al., 2008) is a numerical method for simulating slender elastic bodies. It works by representing the center-line as a polygonal chain, attaching two perpendicular directors to each segment, and defining discrete stretching, bending and twisting deformation measures and a discrete strain energy. Here, we investigate an alternative formulation of this model based on a simpler definition of the discrete deformation measures. Both formulations are equally consistent with the continuous rod model. Simple formulas for the first and second gradients of the discrete deformation measures are derived, making it easy to calculate the Hessian of the discrete strain energy. A few numerical illustrations are given. The approach is also extended to inextensible ribbons described by the Wunderlich model, and both the developability constraint and the dependence of the energy of the strain gradients are handled naturally.


Introduction
The geometric non-linearity of thin elastic rods gives rise to a rich range of phenomena even when the strains are small, see e.g. [9,31] for recent examples. So, the non-linear theory of rods has traditionally combined geometrically non-linearity with linear constitutive laws [1,6]. However, recent interest has expanded beyond the linearly elastic regime, including viscous threads [14,33], plastic and visco-plastic bars [17,3,4], viscoelastic rods [26], capillary elastic beams made of very soft materials [25]. Thin elastic ribbons may also be viewed in this class with a non-linear constitutive law that captures the complex deformation of the cross-sections [34,41,36,38,18,5].
The study of instabilities, especially in the presence of complex constitutive relations, requires an accurate but efficient numerical method. Here, we build on the work of Bergou et al. [12] to propose a numerical method applicable to slender elastic structures in general. To keep the presentation focused, we limit our presentation to elastic rods: both linearly elastic and non-linear elastic constitutive laws are covered. Our main contribution consists in providing a discrete geometric description of slender rods. This kinematic building block is independent of the elastic constitutive law in our formulation, making the extension to inelastic constitutive laws relatively straightforward, as discussed in Section 4.
We follow the classical kinematic approach, and use the arc-length s in the undeformed configuration as a Lagrangian coordinate. We denote the center-line of the rod in the current configuration as x(s) (boldface symbols denote vectors). We introduce an orthonormal set of vectors (d I (s)) 1 I 3 , called the directors, to describe the orientation of the cross-section. We impose the adaptation condition that the director d 3 matches the unit tangent t to the center-line: This kinematic description is common to all variants of the rod model. It is complemented by constitutive equations specifying either the stored energy density (in the case of a hyperelastic theory) or the reaction forces and moments as functions of the four deformation measures or their histories. The formulation is completed by imposing either equilibrium or balance of momenta. The resulting equations for linear elastic constitutive relations are known as the Kirchhoff equations for rods, and they can be derived variationally, see [39,6]; we will not discuss them further. we will not discuss them further.
Various strategies have been proposed to simulate the equations for thin rods numerically. In approaches based on the finite-element methods, it is challenging to represent the kinematic constraint of adaptation (1.1) between the unknown center-line x(s) and the unknown rotation representing the orthonormal directors d I (s). Another approach is based on super-helices or super-clothoids: in these high-order approaches, the bending and twisting strain measures κ (I) (s) are discretized into constant or piecewise linear functions. The result is a highly accurate method which has been successfully applied to several challenging problems [13,15,16]. The price to pay is that the reconstruction of the center-line in terms of the degrees of freedom is non-trivial and non-local. Additionally, some common boundary conditions, such as clamped-clamped conditions, must be treated using non-linear constraints.
A new approach called the Discrete elastic rods method was introduced by Bergou et al. [12]; see [24] for a recent primer. The Discrete elastic rod method is a low-order method, which starts out by discretizing the center-line into a polygonal chain with nodes (x 0 , . . . , x N ). The tangents and material frames d i I are defined on the segments, see Figure 1.1. The adaptation condition (1.1) is used to parameterize the material frames (d i I ) 1 i 3 in terms of the positions (x i−1 , x i ) of the adjacent nodes and of a single twisting angle ϕ i , as described in Section 22.4. A discrete rotation gradient is obtained by comparing the orthonormal directors from adjacent segments: this yields a differential rotation at a vertex between the segments. This must now be projected onto a material frame to yield the bending and twisting strain measures according to equation (1.3). The material frame, however, lives on segments. The original Discrete elastic rod formulation worked around this difficulty by introducing an additional director frame living on the nodes, obtained by averaging the director frames from the adjacent segments [12,24]. In the present work, a different definition of the discrete bending and twisting strain measures is used, see Equations (2.11) and (2.13). This small change simplifies the formulation of model considerably. We note that a similar measure was introduced independently in a recent work on shearable rod models [21].
Overall, the proposed formulation offers the following advantages: • As in the original Discrete rod model, the proposed formulation eliminates two out of the three degrees of freedom associated with the directors at each node using of the adaptation condition (1.1); this leads to a constraint-free formulation that uses degrees of freedom sparingly.
• The formulation of the model is concise: in particular the gradient and Hessian of the discrete elastic energy are given by the simple, closed form formulas listed in Section 3.
• The proposed deformation measures have a clear geometric interpretation: in the context of inextensible ribbons, for example, a discrete developability condition can easily be formulated in terms of the new set of discrete strains, see Section 24.2.
• The kinematic description can easily be combined with various constitutive models to produce discrete models for elastic rods, inextensible ribbons, viscous or visco-elastic rods, etc., as discussed in Section 4.
2 Discrete bending and twisting deformation measures

A compendium on quaternions
Rod models make use of rotations in the three-dimensional space. These rotations are conveniently represented using quaternions. Here, we provide a brief summary of quaternions and their main properties. A complete and elementary introduction to quaternions can be found in [28]. A quaternion q ∈ Q can be seen as a pair made up of a scalar s ∈ R and a vector v ∈ R 3 , q = (s, v). Identifying the scalar s and the vector v with the quaternions (s, 0) and (0, v) respectively, one has the quaternion decomposition q = s + v.
The product of two quaternions q 1 = (s 1 , v 1 ) and q 2 = (s 2 , v 2 ) is defined as The product is non-commutative. A unit quaternion r = s + v is a quaternion such that s 2 + |v| 2 = 1. Unit quaternions represent rotations in the three-dimensional Euclidean space, in the following sense. Definer = s−v as the quaternion conjugate to r. Define the action of the unit quaternion r on an arbitrary vector w as where the left-hand side defines a linear map on the set of vectors w, and the right-hand side is a double product of quaternions. It can be shown that (i) the quaternion r * w is a pure vector, (ii ) the mapping w → r * w is a rotation in Euclidean space, (iii ) the quaternion r can be written as r = ±r n (θ) where r n (θ) = cos θ 2 + n sin θ 2 = exp n θ 2 , (2.2) Figure 2.1: A node x i , its adjacent segments, and the adjacent nodes x i±1 in reference (gray background) and current (white background) configurations. Director frames, shown in purple, are represented by a unit quaternion, whose action on the Cartesian frame e I yields the director frame.
θ is the angle of the rotation, and n is a unit vector subtending the axis of the rotation. Note that both unit quaternions +r n (θ) and −r n (θ) represent the same rotation.
Given two unit quaternions r 1 and r 2 , consider the product r 2 r 1 : for any vector w, the equality (r 2 r 1 ) * w = r 2 r 1 w r 2 r 1 = r 2 r 1 w r 1 r 2 = r 2 * (r 1 * w) shows that the unit quaternion r 2 r 1 represents the composition of the rotations associated with r 1 applied first, and r 2 applied last. The multiplication of unit quaternions is therefore equivalent to the composition of rotations. In view of this, we will identify rotations with unit quaternions. The inverse of the rotation r will accordingly be identified with the conjugate r.

Parallel transport
Parallel transport plays a key role in the Discrete elastic rods model, by allowing one to define twistless configurations of the material frame in an intrinsic way. For two unit vectors a and b such that b = −a, the parallel transport from a to b is the rotation mapping a to b, whose axis is along the binormal a×b. Parallel transport can be interpreted geometrically as the rotation mapping a to b and tracing out the shortest path on the unit sphere [12].
An explicit expression of the parallel transport from a to b in terms of unit quaternions is [27] The proof is as follows. First it can be verified that p b a is a unit quaternion, as can be shown by using the identity |a×b| 2 1+a·b = 1−(a·b) 2 1+a·b = 1 − a · b. Second, the rotation p b a indeed maps a to as can be checked by explicit calculation. Finally, the axis of p b a is indeed about the binormal a × b: equation (2.2) shows that the vector part of the unit quaternion is aligned with the rotation axis and equation (2.3) shows that the vector part of p b a is aligned with a × b. For two units vectors a and b such that a = −b, the parallel transport p b a is ill-defined.

Reference and current configurations
A configuration of the discrete rod is defined by a set of nodes x i indexed by an integer i, 0 i N . We consider an open rod having unconstrained endpoints x 0 and x N for the moment; alternate boundary conditions such as periodic or clamped boundary conditions are discussed later. For simplicity, we limit attention to the case where the nodes are equally spaced in the undeformed configuration, i.e., the undeformed length j is independent of the segment index j: it is denoted as In addition to the undeformed configuration, the simulation deals with two configurations shown in Figure 2.1: • Reference configuration (shown with a gray background in the figure). The only role of the reference configuration is to allow a parameterization of the current configuration. It does not bear any physical meaning and its choice does not affect the results of the simulations. It is chosen for convenience.
In the reference configuration, the position of node i is denoted by x i . The orthonormal frame of directors on segment i connecting nodes x i and x i+1 is denoted as (d i I ) I∈{1,2,3} . The adaptation condition from equation (1.1) requires that the third director d j 3 coincides with the unit tangent T j to the segment in reference configuration, • Current configuration (shown with a white background). The current configuration is the physical configuration of the rod and is the unknown in a simulation. It is parameterized by the degrees of freedoms (see Section 22.7).
In the current configuration, the center-line of the rod is defined by the node positions x i . On segment i connecting the nodes x i and x i+1 , the directors are denoted as (d i I ) I∈{1,2,3} . The adaptation condition from equation (1.1) requires As shown in the figure, the orthonormal director frames (d j I ) 1 I 3 and (d j I ) 1 I 3 are represented by unit quaternions D j and d j , respectively, that yield the directors when applied to the Cartesian basis e I : (2.7) The quaternions d j and d j therefore represent the rotations 3 I=1 d j I ⊗ e I and 3 I=1 d j I ⊗ e I , respectively. They fully describe their respective frames.
The reference and current configurations are not assumed to be close to one another. However, our parameterization introduces a weak restriction: the reference configuration must be chosen such that the angle of the rotation d j D j mapping d j I to d j I does not come close to π, in any of the segments j. This condition is fulfilled by resetting periodically the reference configuration to the current configuration: • in dynamic simulations, this reset is typically done at the end of any time step; • in equilibrium problems, it is typically done whenever an equilibrium has been found and the load is incremented.
In principle, it is even possible to reset the reference configuration in the middle of the Newton-Raphson iteration used to update a time step (in the dynamic case) or the non-linear equilibrium (in the static case), but special care is required as this amounts to changing the parameterization of the unknown during iteration. All the applications shown at the end of this paper deal with the static case, i.e., they involve the calculation of equilibria for a series of load values: our simulations are initialized with the reference configuration x i , d j I representing a simple starting point which is typically a straight or circular equilibrium configuration without any load (see the example description for further details). The reference configuration is reset each time an equilibrium is found.

Centerline-twist representation
In this section, we introduce a parameterization that provides a concise representation of the current configuration that is at the heart of the Discrete elastic rod method. All quantities from the reference configuration, such as the node positions x i , unit tangents T j , material frames d j 3 and associated rotations D j , are known. We proceed to analyze the current configuration. A key observation is that equation (2.6) yields the tangent director d j 3 as a function of the node positions x i : if the nodes are prescribed, the full frame of directors d j I can only twist about this tangent. The three directors (d j I ) 1 I 3 on segment j, as well as the associated unit quaternion d j by equation (2.7), can therefore be parameterized in terms of • the adjacent nodes positions x j and x j+1 , • a scalar twist angle ϕ j .
The parameterization used by the Discrete elastic rod method may be written as [12,11,2] where x j and x j+1 are the positions of the adjacent nodes, ϕ j is the twisting angle, is the parallel transport from the reference unit tangent T i to the current unit tangent t i (x j , x j+1 ) given as a function of the node positions by equation (2.6), r T j (ϕ j ) = cos ϕ j 2 + T j sin ϕ j 2 is the rotation about T j with angle ϕ j (see equation (2.2)), and D j is the unit quaternion associated with the reference configuration of the directors (see equation (2.7)).
This yields a parameterization of the rod in terms of the degrees of freedom vector X = (x 0 , ϕ 0 , x 1 , ϕ 1 , x 2 , · · · , x n−1 , ϕ n−1 , x n ), (2.10) where the nodes positions x i are read off directly from X and the directors are reconstructed using equations (2.7) and (2.8). It is called the centerline-twist representation. As observed in Section 22.2, the parallel transport in equation (2.9) is singular if t i (x j , x j+1 ) = −T i , i.e., if any one of the tangents flips by an angle π between the reference and current configuration. The periodic reset of the reference configuration described earlier in Section 22.3 prevents this from happening.
Note that in the original paper of [12], parallel transport was used to move the directors from one segment to an adjacent segment (spatial parallel transport). This makes the directors dependent on the degrees of freedom associated with all the nodes and segments located on one side of the directors. Here, like in subsequent work by the same authors [11,2], we use parallel transport 'in time': in equation (2.8), p j (x j , x j+1 ) serves to parameterize the directors in current configuration in terms of the same set of directors in reference configuration. With this approach, the directors are a function of the local degrees of freedom, as implied by the notation d j (x j , ϕ j , x j+1 ) in equation (2.8).

Lagrangian rotation gradient
This rotation is an Eulerian quantity: like its continuous counterpart κ(S), it is not invariant when the rod rotates rigidly. The following, however, is a Lagrangian version q i of the rotation gradient that is invariant by rigid-body rotations, Here, we depart from earlier work on Discrete elastic rods [12] who used q i := q avg i is some average of the adjacent frames d i−1 and d i . A definition of the rotation gradient similar to (2.11) has been used in the context of shearable rods [21] and in a purely geometric analysis of discrete rods [27].
We now explain why this definition represents a Lagrangian rotation gradient. One way to define a Lagrangian rotation gradient, is to pull back the Eulerian rotation gradient d i d i−1 to the reference configuration. However, the discreteness of our representation raises a difficulty: the frames are defined on the segment while the Eulerian rotation gradient d i d i−1 is defined on the nodes. So, we could use the frame associated with the segment on the left of the node for the pull back by defining q left but this biases the choice on the left. Or, we could use the right counter-part, q right this biases the choice to the right. However, these biases are apparent only: elementary calculations shows that these are in fact identical thereby justifying our definition. The unit quaternion q i introduced in equation (2.11) is the discrete analogue of the pull-back (e I ⊗ d I (s)) · κ(s) of the rotation gradient κ(s) used in the continuous rod theory, whose components κ J (s) = e J · [(e I ⊗ d I (s)) · κ(s)] = d J (s) · κ(s) define the bending and twisting measures. In the following section, bending and twisting are similarly extracted from the unit quaternion q i .

Bending and twisting deformation measures
The discrete bending and twisting deformation measures are defined as the components of the pure vector, This κ i is twice the vector part I(q i ) = qi−q i 2 of the quaternion q i , which shows that it is indeed a vector. Let κ i,I denote its components in the Cartesian basis, such that κ i = 3 I=1 κ i,I e I . The first two components κ i,1 and κ i,2 can be interpreted as measures of bending about the transverse directors d j 1 and d j 2 , while the third component κ i,3 is a discrete measure of twisting. Like q i , these are integrated versions of their smooth counterparts, that are proportional to the discretization length ; this will be taken into account when setting up a discrete strain energy.

Summary
The current configuration is reconstructed in terms of the degrees of freedom X from equation (2.10) as follows: • the node positions x i are directly extracted from X, see equation (2.10), • the unit tangents t j (x j , x j+1 ) are obtained from equation (2.6), • parallel transport p j (x j , x j+1 ) is obtained by combining equations (2.9) and (2.3), • the director frames d j (x j , ϕ j , x j+1 ) are obtained from equation (2.8), • the bending and twisting deformation vector Finally, a possible definition of the discrete stretching measure on segment j joining nodes x j and x j+1 is see for instance [26]. Here, denotes the undeformed length of the segments, which is different from the length |x j+1 − x j | in reference configuration. This discrete stretching measure is an integrated version of the continuous strain ε(S), like the discrete bending and twisting deformation measures κ i,I . The particular definition of the stretching measure ε j in equation (2.14) requires the evaluation of the squared norm and not of the norm itself, which simplifies the calculation of the gradient significantly.

Interpretation of the discrete deformation measures
We now show that the discrete deformation measures (up to a minor rescaling) may be interpreted as the rotation that transports the director frame from one segment to the next. Consider the function ψ 15) and note that ψ(t) ≈ 1 for t 1 (See supplementary information for a plot of this function). Define the adjusted deformation measure to be Even for moderate values of |κ i |, the original and adjusted deformations measures are not very different, ω i,J ≈ κ i · e J , as the variations of the function ψ are bounded by 1 ψ(t) π/2.
The adjusted deformation measure has a simple geometric interpretation. We start from the decomposition (2.2) of the rotation gradient q i = r ni (θ i ) = cos θi 2 + n i sin θi 2 = exp ni θi 2 , where n i is a unit vector aligned with the axis of the rotation q i , and θ i is the angle of this rotation, 0 θ π. In view of equation (2.13), κ i = q i − q i = 2 sin θi 2 n i . In particular, |κ i | = 2 sin θi 2 and so ψ(|κ i |) = θi/2 sin(θi/2) from equation (2.15). The adjusted strain is then ω i,J e J = ψ(|κ i |) κ i = θi/2 sin(θi/2) 2 sin θi 2 n i = θ i n i : in effect, the adjustment factor ψ(|κ i |) transforms κ i = 2 I(q i ) (twice the vector part of q i ) into ω i,J e J = θ i n i = 2 log q i (twice its logarithm).
: as is well known, the conjugate rotation d i d i−1 has the same angle θ i as the original rotation q i and its axis is obtained by applying the rotation d i−1 to the original axis. This can be rewritten as where is a (finite) rotation vector. Similar relations have been derived in the work of [27]. Repeating the same argument with i , one can show that the vector Ω has the same decomposition in the other directors frame, Ω i = ω i,J d i J : The benefit is that ω i,J have an even simpler interpretation, see equations (2.17-2.18). The drawback is that the function ψ gets involved in the calculation of the strain, resulting in cumbersome formulas for the strain gradients (Section 3). Therefore, we continue to use the original deformation measures.

Variations of the discrete deformation measures
In this section, we present explicit formulae for the first and second derivatives of the deformation measures κ i (summarized in Section 22.7) with respect to X. The first gradient is required for determination of the internal forces, which are the first gradient of the strain energy. The availability of the second gradient in analytical form makes it possible to use implicit time-stepping methods (in dynamic problems) or to evaluate the Hessian for second order methods (in static problems).
Our notation for variations is first introduced based on a simple example. For a function y = f (x) taking a vector argument x and returning a vector y, the first variation is the linear mapping δx → δy = f (x) · δx, where δx is a perturbation to x and f (x) is the gradient matrix. To compute the second variation, we start from δy = f (x) · δx, perturb the argument x of f as x + δx and linearize the result as Here, the second variation is defined as the second order term δ 2 y := f (x) : (δx ⊗ δx), where f (x) is the Hessian. By construction, δ 2 y is a quadratic form of δx.
In this section, the reference configuration is fixed and the degrees of freedom are perturbed by δX = (· · · , δx i , δϕ i , · · · ). We simply present the final results; the detailed calculations are cumbersome but straightforward, and provided as supplementary material.
where I is the identity matrix, τ i is the third-order tensor τ i = (I − t i ⊗ t i ) ⊗ t i , the colon denotes the double contraction of the last two indices of the rank-three tensor on the left-hand side. For any permutation (n 1 , n 2 , n 3 ) of (1, 2, 3), T (n 1 , n 2 , n 3 ) denotes the generalized transpose of a rank-three tensor µ such that µ • parallel transport p i = p t i T i from equations (2.9) and (2.3), where for any vector a, a × is the linear operator and k i is the binormal defined by (3.5) • rotation gradient q i from equation (2.11), • discrete bending and twisting strain measure vector κ i from equation (2.13), where I(q) = q−q 2 denotes the vector part of a quaternion q. • stretching measure ε i from equation (2.14), In these formula, the first and second variations of the rotations p i , d i and q i are not captured by quaternions but by regular vectors, bearing a hat, such as δp i , δ 2pi , δd i , etc. Equations (3.1-3.8) involve standard calculations from Euclidean geometry: the more advanced quaternion calculus is only required in the proof given in the supplementary materials.
Equations (3.1-3.8) suffice to calculate the strain gradients. They can be implemented easily and efficiently using standard libraries for vector and matrix algebra. These formulas for the first and second gradient of strain are considerably simpler than those applicable to the discrete strain measures used in earlier work on Discrete elastic rods [12,2,31,26].
In equations (3.1-3.8), the perturbations to the degrees of freedom such as δx i and δϕ i are dummy variables. The first-order variations such as δt i , δp i , must be represented numerically as linear forms, by storing their coefficients as vectors. Similarly, the second-order variations such as δ 2 t i , δ 2pi , etc. are represented as quadratic forms, whose coefficients are stored as sparse symmetric matrices; the reader is referred to [26] for further details on this aspect of implementation. All these coefficients depend on the current configuration and must be updated whenever the degrees of freedom X or the reference configuration change.
These vectors and symmetric matrices should be stored at an appropriate place in the data structure representing the Discrete elastic rod. The tensors representing δt i , δp i , δ 2pi and δ 2d i depend on the perturbations δx i and δx i+1 to the nodes adjacent to a given segment, and therefore best stored in the data structure representing segments, which have access naturally to the degrees of freedom of the adjacent nodes. The quantities δd i and δ 2d i make use of the twisting angle δϕ i in addition to the adjacent nodes δx i and δx i+1 , and should be stored in the data structure representing the material frame attached to particular segment. The quantities δq i , δκ i , δ 2q i and δ 2 κ i are best stored in a data structure representing an elastic hinge at a node, that depends on the material frames at the adjacent segments.

Constitutive models
The discrete kinematics from Sections 2 and 3 can be combined with a variety of constitutive laws to produce discrete numerical models for rods that are elastic, viscous, visco-elastic, etc.: the procedure has been documented in previous work, and it is similar to the general approach used in finite-element analysis. Elastic problems are treated by introducing a strain energy function U (X), whose gradient with respect to X yields the negative of the discrete elastic forces [12,26]; while viscous problems are treated by introducing a discrete Rayleigh potential U (X,Ẋ) ,whose gradient with respect to velocitiesẊ yields discrete viscous forces [11,14,2]. More advanced constitutive models such as visco-elastic laws can be treated by variational constitutive updates of a discrete potential that makes use of the same discrete deformation measures [26]. In [26], it is emphasized that these different constitutive models can be implemented independently of the geometric definition of discrete deformation measure. Using this decoupled approach, it is straightforward to combine the kinematic element proposed in the present work with constitutive element from previous work. We illustrate this with the classical, linearly elastic rod in Section 44.1 (Kirchhoff rod model), and a discrete inextensible ribbon model in Section 44.2 (Wunderlich model). The latter is a novel application of the Discrete elastic rod method.

Elastic rods (Kirchhoff model)
The classical, continuous theory of elastic rods uses a strain energy functional U where κ (I) (s) = κ(s) · d I (s) are the components of the rotation gradient in the frame of directors, see equation (1.3). For an inextensible, linearly elastic rod made of a Hookean material with natural curvature κ (0) , for instance, the strain energy density is where Y and µ are the Young modulus and the shear modulus of the material, I 1 and I 2 are the geometric moments of inertia of the cross-section, and J is the torsional constant.
In the discrete setting, we introduce a strain energy i E i (κ i ) where the sum runs over all interior nodes i. The strain energy assigned to an interior node i is defined in terms of the strain energy density as (no implicit sum over i), where is the undeformed length of the segments for a uniform mesh. The factor in the argument of E takes care of the fact that κ i is an integrated quantity, i.e., it is κi · e J and not just κ i · e J that converges to the continuous strain κ (J) (s); for a non-uniform grid, this would need to be replaced with the Voronoi length associated with the interior vertex i in undeformed configuration. The factor in factor of E in equation (4.2) ensures that the discrete sum i E i = i E converges to the integral L 0 E ds = U [12]. Consider for instance an equilibrium problem with dead forces F i on the nodes: it is governed by the total potential energy Φ(X) defined in terms of X = (x 0 , ϕ 0 , . . . , ϕ N −1 , x N ) as This energy is minimized subject to the inextensibility constraints In equations (4.3-4.4), the elastic deformation measures κ i and ε j is reconstructed in terms of the unknown X by the method described in Section 2, as expressed by the notation In the case of dead forces, the first and second variations of the total potential energy is derived as see for instance [26]. Here, ∂Ei ∂κi and ∂ 2 Ei ∂κ 2 i are the internal stress and tangent elastic stiffness produced by the elastic constitutive model E i (κ i ). The two terms appearing in the parentheses in the right-hand side of δ 2 Φ are known as the elastic and geometric stiffness, respectively. The first and second variations of the strain, δκ i and δ 2 κ i , are available from Section 3: the equilibrium can be solved using numerical methods that require evaluations of the Hessian of the energy. Note that the Hessian can be represented as a sparse matrix thanks to the local nature of the energy contributions In the applications presented in the forthcoming sections, we find equilibrium configurations by minimizing Φ(X) in equation (4.3) using the sequential quadratic programming method (SQP) described by [30]; it is an extension of the Newton method for non-linear optimization problems which can handle the non-linear constraints in equation (4.4). It requires the evaluation of the first and second gradient of the energy Φ, see equation (4.5), and of the first gradient of the constraints that are available from equation (3.8). We used an in-house implementation of the SQP method in the C++ language, with matrix inversion done using the SimplicialLDLT method available from the Eigen library [23].

Inextensible elastic ribbons (Wunderlich model)
Ribbons made up of material that are sensitive to light [43,22] or temperature change [8] have been used to design lightweight structures that can be actuated. They are easy to fabricate, typically by cutting out a thin sheet of material, and their thin geometry can turn the small strains produced by actuation into large-amplitude motion. For this reason, there has been a surge of interest towards mechanical models for elastic ribbons recently. When the width-to-thickness ratio of a ribbon cross-section is sufficiently large, its mid-surface is effectively inextensible. Sadowsky has proposed a one-dimensional mechanical model for inextensible ribbons [34]. Sadowsky model is one-dimensional but differs from classical rod models in two aspects: one of the two bending modes is inhibited due to the large width-to-thickness aspect-ratio, and the two remaining twisting and bending modes are governed by an non-quadratic strain energy potential that effectively captures the inextensible deformations of the ribbon mid-surface. Sadowsky's strain energy is non-convex which can lead to the formation of non-smooth solution representing a micro-structure [20,32]; to avoid these difficulties, we use the higher-order model of Wunderlich that accounts for the dependence of the energy on the longitudinal gradient of bending and twisting strain [41].
The Wunderlich model has been solved numerically by a continuation method, see for instance the work of [38]. The continuation method is an extension of the shooting method that can efficiently track solutions depending on a parameter [19]. It requires the full boundary-value problem of equilibrium to be specified spelled out, which is quite impractical in the case of Wunderlich ribbons. A recent and promising alternative is the high-order method of [16] that starts from linear and quadratic interpolations of the bending and twisting strains, and treats the center-line position and the directors as secondary (reconstructed) quantities. In the present work, we explore an alternative approach, and show that simulations of the Wunderlich model are possible with limited additional work on top of the generic Discrete elastic rod framework.
We build on the work of [18] who have shown that the Wunderlich model can be viewed as a special type of a non-linear elastic rod, see also [37]. Accordingly, simulations of the Wunderlich model can be achieved using a simple extension of the Discrete elastic rod model, which we describe now. We first introduce a geometric model for a discrete inextensible ribbon, in which the inextensibility of the mid-surface is fully taken into account. We start from a rectangular strip lying in the plane spanned by (e 1 , e 3 ), as shown in Figure 4.1a. Through every node (shown as black dots in the figure), we pick a folding direction within the plane of the strip (brown dotted line in the figure); we denote by π/2 − γ i the angle of the fold line relative to the centerline. Next, we fold along each one of these lines by an angle θ i , as shown in Figure 4.1b. We call the resulting surface a discrete inextensible ribbon. By construction, it is isometric to the original strip.
Let us now introduce the director frames d i I following rigidly each one of the faces: the planar faces are spanned by the directors d i 1 and d i 3 . By construction the vector Ω i for the rotation that maps one frame, d i−1 I , to the next, d i I , see equation (2.17), is aligned with the fold line. We observe that the unit tangent along the fold direction is e 3 sin γ i + e 1 cos γ i in the flat configuration of the strip; it is therefore mapped to In view of this, we conclude Comparing with equation (2.18), we obtain the discrete deformation measure in the developable ribbon as ω i,1 = θ i cos γ i (bending mode), ω i,2 = 0 (inhibited bending mode) and ω 3,i = 0 (twisting mode).
Eliminating θ i , we find ω i,2 = 0 and ωi,3 ωi,1 = tan γ i , which can be rewritten in terms of the original discrete strain κ i = (κ i,1 , κ i,2 , κ i,3 ) with the help of equation (2.15) as The continuous version of the developability conditions is κ 2 (s) = 0 and κ 3 (s) = η(s) κ 1 (s), where η(s) = tan γ(s) and π/2 − γ(s) is the angle between the generatrix and the tangent, see for instance [18]. It is remarkable that the discrete developability conditions (4.6) are identically satisfied. This is a consequence of the simple geometric interpretation for the discrete deformation measures introduced in Section 2.
To simulate inextensible ribbons, we introduce the unknown η i as an additional degree of freedom at each one of the interior nodes, and we use in equation (4.3) a strain energy density directly inspired from that of Wunderlich [18,38] In equation (4.7), D = Y h 3 12 (1−ν 2 ) is the bending modulus from plate theory, h is the thickness, w is the width and is the discretization length. The quantity η i is calculated by a central-difference approximation of the gradient of η, where is the mesh size. The constraint (4.6) 2 is imposed at each node using the SQP method. Introducing the nodal degrees of freedom η i together with the constraint (4.6) 2 allows us to work around calculating η i = κ i,3 /κ i,1 , which is a division with a potentially small denominator; in addition, this approach warrants that κ i,3 = 0 whenever κ i,1 = 0, which is necessary for the Wunderlich energy to remain finite. It is a feature of the Wunderlich model that η can take on arbitrary values in intervals where κ 1 vanishes identically. To work around this, we have introduced an artificial drag on the η i 's between iterations of the solve. When convergence is reached, the drag force is identically zero.
The discrete potential energy Φ(X) is minimized by the same numerical method as described in Section 44.1, taking into account the kinematic constraints (4.6) and the centerline inextensibility (4.4).

Illustrations
In this section, the Discrete elastic rod model is used to simulate • a linearly elastic model for an isotropic beam, Section 55.1, • a linearly elastic model for an anisotropic beam with natural curvature, Section 55.2, • Sano and Wada's extensible ribbon model, Section 55.3, These examples serve to illustrate the capabilities of the model. In addition, comparison with reference solutions available from the literature provide a verification of its predictions.

Euler buckling
We consider Euler buckling for a planar, inextensible elastic rod that is clamped at one endpoint. We consider two types of loading: either a point-like force f p at the endpoint opposite to the clamp, or a force f d distributed along the length of the rod. In both cases, the force is applied along the initial axis of the rod, is invariable (dead loading), and is counted positive when compressive. A sketch is provided in Figure 5 The typical simulation time is about 1/10s for each equilibrium on a personal computer, and the results are shown in Figure 5.1, and compared to that obtained by solving (5.1) using the bvp4c solver from Matlab. A good agreement on the position of the endpoint of the rod is found in the entire post-bifurcation regime. In addition, the onset of bifurcation agrees accurately with the prediction (5.3) from the linear stability analysis.

Folding of an over-curved ring
A circular elastic ring with length L can buckle out of plane if its natural natural curvature κ (0) does not match the curvature 2π/L of the circle with length L. In the case of an over-curved ring, such that κ (0) > 2 π/L, a buckled shape featuring two symmetric lobes has been reported [29,8,7]. Here, we simulate the buckling of over-curved rings using the Discrete elastic rod model and compare the results to the experimental shapes reported by [29].
In the experiments of [29], a commercial slinky spring with a width w = 5 mm, thickness t = 2 mm and length L = 314 mm is used; Poisson's ratio has been measured as ν = 0.41. Note that the aspect-ratio t/w = 0.4 is not small. In our simulations, we use a discrete version of the quadratic strain energy for a linearly elastic rod having an anisotropic cross-section (I 1 = I 2 ), see equations (4.1-4.3). We use the elastic moduli reported in the supplement of [29]: The value 0.256 in the numerator was obtained by [29] from the book of [40], and applies to the particular commercial Slinky used in their experiments. In the absence of applied loading, the value of the Young modulus is irrelevant and we set Y = 1 in the simulations. The equilibria of the Discrete elastic rod are calculated numerically for different values of the dimensionless loading parameter O = 2 π κ (0) /L, with O > 1 corresponding to the over-curved case. We use N = 400 nodes. We start from a circular configuration having curvature κ (0) = 2 π/L. The Discrete elastic rod model is closed into a ring as follows: the first two nodes and the last two nodes are prescribed to x 0 = x N −1 = 0 and x 1 = x N = e x , respectively; the first and last frames are also fixed, such that d 0 1 = d N −1 1 = e y . Next, the over-curvature κ (0) is varied incrementally. For each value of κ (0) , an equilibrium configuration is sought, and we extract the minimal distance D between pairs of opposite points on the ring. In Figure 5

Buckling of a bent and twisted ribbon
We now turn to an effective rod model applicable to thin ribbons. Sano and Wada [35] have proposed an effective beam model that accounts for the stretchability of the ribbon having moderate width, thereby Here, is the uniform segment length in undeformed configuration, A 1 = Y h w 3 /12 and A 2 = Y h 3 w/12 are the initial bending moduli, A 3 = Y h 3 w/[6 (1 + ν)] is the initial twisting modulus and The parameter ξ is the typical length-scale where the stretchability of the mid-surface starts to play a role.
The potential E i from equation (5.5) is non-quadratic, meaning that the equivalent rod has non-linear elastic constitutive laws. The elastic model (5.5) of Sano and Wada is applicable to thin ribbons, for w h. It is based on kinematic approximations. A refined version of their model has been obtained very recently by [5], by asymptotic expansion starting from shell theory; in the latter work, a detailed discussion of the validity of the various models for thin ribbons can also be found. We do not expect any difficulty in applying the present numerical model to the ribbon model in [5]. Both the models of Sano and Wada, and of Audoly and Neukirch improve on Wunderlich model by addressing the stretchability of the ribbon; unlike the Wunderlich model, however, they ignore the dependence of the energy on η , and therefore account less accurately for the 'conical' singularities often observed in ribbons [42] as η varies quickly there.
Following [35], we consider the buckling of a ribbon with length L = π R bent into half a circle, whose ends are twisted in an opposite senses by an angle φ, see Figure 5.3. Specifically, they identified a snapping instability which occurs for moderately wide ribbons, when the width w < w * is below a threshold w * ≈ 1.24 √ L h, but not for wider ribbons, when w > w * ; they showed that their equivalent rod model can reproduce this instability, as well as its disappearance for larger widths. In Figure 5.3, we compare the predictions of a Discrete elastic rod model using (5.5) with the original experiments and simulations from [35]. Our simulations use N = 350 vertices each. Our simulation results are in close agreement with both their experimental and numerical results. In particular, we recover the instability when w < w * only.

The elastic Möbius band
An extension of the Discrete elastic rod model that simulates the inextensible ribbon model of Wunderlich has been described in Section 44.2, see equation (4.7). With the aim to illustrate and verify this discrete model, we simulate the equilibrium of a Möbius ribbon, and compare the results to those reported in the seminal paper of Starostin and van der Heijden [36]. In our simulations, the inextensible strip is first bent into a circle, and the endpoints are turned progressively twisted by an angle of 180 • to provide the correct topology. The final equilibrium shapes are then recorded for all possible values of the aspect-ratio w/L. For   = −e y . The equilibrium shape for a particular aspect-ratio w/L = 1/ (2 π) is shown in Figure 5.4a, with arclength L = 1, width w = 1/ (2 π) and N = 150 simulation nodes. A detailed comparison with the results of [36] is provided in Figure 5.4b, where the scaled bending and twisting strains κ i,1 / and κ i,3 / from the discrete model with N = 250 vertices are compared to the strains κ 1 (s) and κ 3 (s) obtained by [36] using numerical shooting, for different values of the width w.

Conclusion
We have presented a new formulation of the Discrete elastic rod model. The formulation is concise and uses only the minimally necessary degrees of freedom: the position of the nodes and the angle of twist of the segments between the nodes. It naturally incorporates the adaptation condition without the need for any constraint, penalty or Lagrange multiplier. We use bending and twisting deformation measures that are different from those used in earlier work on Discrete elastic rods, are equally consistent with their continuum counterparts, and have a simple physical interpretation in the discrete setting. Consequently, the formulation is versatile in the sense that it can be combined with a variety of linear and nonlinear as well as elastic and inelastic constitutive relations. In fact, ribbons can be incorporated as generalized rods with a nonlinear constitutive model. Similarly, the formulation can be used both for static and dynamic simulations.
We have presented explicit formulae for the first and second derivatives of the deformation measures that eases implementation. We have demonstrated our method with four examples, and verified our results against prior experimental and theoretical findings in the literature.
The source code used for the numerical simulation is available through CaltechDATA at https://data. caltech.edu/records/2147.
All three authors conceived of the work and the formulation. KK conducted the theoretical and numerical calculations with advice from BA and KB. KK and BA took the lead in writing the manuscript and all three authors finalized it.
The authors declare that there are no competing interests.

Appendices A Plot of function ψ
The function ψ(t) from equation (2.15) is plotted in figure A.1.

B Detailed derivation of the strain gradients
In this appendix, we provide a detailed derivation of the first and second gradients of the strain appearing in section 3.
To derive the first gradient, we continue to use the conventions of section 3: we use a perturbation δX of the degrees of freedom, and we denote by δy = f (x) · δx the first variation of a generic quantity y = f (x) entering in the reconstruction of the discrete strain, where x depends indirectly on the degrees of freedom X.
For the second variation, however, we work here in a slightly more general setting than in the main text, as we consider two independent perturbations δ 1 X and δ 2 X of the degrees of freedom. We denote by δ 1 x and δ 2 x the corresponding perturbations to the variable x, and by δ 1 y and δ 2 y the first-order variations of the functions: δ 1 y = f (x) · δ 1 x and δ 2 y = f (x) · δ 2 x are simply obtained by replacing the generic increment δx appearing in the first order variation δy with δ 1 x and δ 2 x, respectively. To obtain the second variation, we perturb the argument x appearing in δ 1 y = f (x) · δ 1 x as x + δ 2 x, leaving δ 1 x untouched, and we expand the result to first order in δ 2 x. This yields a quantity denoted as δ 12 y, which we can write formally as δ 12 y = f (x) : (δ 1 x ⊗ δ 2 x), where f (x) is the Hessian. By a classical result in the calculus of variations, the quantity δ 12 y is bilinear and symmetric with respect to δ 1 x and δ 2 x. The second variation δ 2 y given in the main text is the quadratic form obtained by ultimately condensing the variations δ 1 x and δ 2 x appearing in δ 12 y into a single perturbation δx = δ 1 x = δ 2 x.

B.1 Infinitesimal rotation vectors
As an important preliminary result, we show that the first variation of a rotation represented by a unit quaternion s can be characterized by means of first-order vector-valued increment δŝ ∈ R 3 , and that the second variation of s can be represented by means of a second-order vector-valued increment δ 12ŝ ∈ R 3 . These vectors will be referred as the infinitesimal rotation vectors. They are connected to the variations δs and δ 12 s of the quaternion by δs = 1 2 δŝ s δ 12 s = The increment δŝ is linear with respect to the variation δX of the degrees of freedom, and the increment δ 12ŝ is bilinear with respect to the independent variations δ 1 X and δ 2 X of the degrees of freedom. As usual in our notation, δ 1ŝ and δ 2ŝ denote the first-order variation δŝ, evaluated on the increment δ 1 X and δ 2 X, respectively. This representation of the first and second variations of a parameterized quaternion is equivalent to that proposed by [10]. The proof is as follows. By taking the first variation of the condition 2 (s s − 1) = 0 that s is a unit quaternion, we have 0 = 2 δs s + 2 s δs = 2 δs s + 2 δs s. This shows that the quaternion 2 δs s is a pure vector: this the vector δŝ introduced in equation (B.1) above. Now, by inserting the increment δ 1 X in the relation just derived, we have 2 δ 1 s s ∈ R 3 ; perturbing this expression as s ← s + δ 2 s, one shows that the following quaternion is a pure vector: 2 δ 12 s s + 2 δ 1 s δ 2 s = 2 δ 12 s s + 1 2 (δ 1ŝ s) (δ 2ŝ s) = 2 δ 12 s s − 1 2 δ 1ŝ δ 2ŝ = 2 δ 12 s s + 1 2 δ 1ŝ · δ 2ŝ − 1 2 δ 1ŝ × δ 2ŝ ; here, the quaternion product δ 1ŝ δ 2ŝ has been evaluated using the definition (2.1). Adding the vector quantity 1 2 δ 1ŝ × δ 2ŝ , the quantity 2 δ 12 s s + 1 2 δ 1ŝ · δ 2ŝ appears to be another pure vector: this is the vector δ 12ŝ introduced in equation (B.1).
The second-order infinitesimal rotation vector δ 12ŝ can be calculated directly from the first-order one δŝ as Here, δ 1 (δ 2ŝ ) denotes the first-order variation of δ 2ŝ when s is perturbed into s + δ 1 s; this quantity is not symmetric with respect to the perturbations δ 1 s and δ 2 s. Similarly, δ 2 (δ 1ŝ ) denotes the first-order variation of δ 1ŝ when s is perturbed into s + δ 2 s. The proof of equation (B.2) is as follows. Take the second variation of δ 1 s = 1 2 δ 1ŝ s from equation (B.1) as The left-hand side is symmetric with respect to the perturbations δ 1 s and δ 2 s, by definition of the second variation. Symmetrizing the righthand side, we obtain δ 12 s = δ1(δ2ŝ)+δ2(δ1ŝ) In the following sections, the first and second variations of the rotations that enter into the Discrete elastic rod model, such as the parallel transport p i and the director rotation d i , will be systematically represented using the corresponding infinitesimal rotation vectors, such as δp i , δ 12p i , δd i and δ 12d i .

B.2 Variation of parallel transport
We start by deriving the variations of the parallel transport p b a from the unit vector a to the unit vector b defined in equation (2.3), assuming b = −a. As a represents the fixed unit tangent T i in reference configuration, it remains unperturbed, δa = 0 δ 12 a = 0.
Since b remains a unit vector during the perturbation, we have 1 2 (|b| 2 − 1) = 0. Taking the first and second variation of this constraint, we have

B.2.1 First variation of parallel transport
As a preliminary step, we consider the case of parallel transport from b to its perturbation b + δb. Using b · δb = 0, we find from equation (2.3), We now return to the calculation of p b+δb a . Following the work of [12], as well as equations [3.7] and [A.2] from [26], one can use a holonomy reasoning to shows that, to first order in δb, We rewrite this as where δθ = − k 2 · δb and k is the scaled binormal that characterizes the holonomy (see [12]), The infinitesimal rotation r a (δθ) from equation (B.3) can be found from equation ( In view of this, the first order variation of parallel transport writes as Identifying with equation (B.1), we find that it is captured by the infinitesimal rotation vector

B.2.2 Second variation of parallel transport
From equation (B.6), we have Using equation (B.4), the variation of the binormal is found as Inserting into equation (B.7) and reordering the terms, we find In view of equation (B.2), we can obtain the second-order infinitesimal rotation vector by symmetrizing this with respect to the increments δ 1 b and δ 2 b: (B.8)

B.2.3 Application to a Discrete elastic rod
In a Discrete elastic rod, the transport is from the undeformed tangent a = T i to the deformed tangent b = t i , see equation (

B.3 Variation of unit tangents
With E i = x i+1 −x i as the segment vector, the variation of the unit tangent Next, the second variation is calculated as Here, we have used δ 12 E i = 0 since E i = x i+1 − x i depends linearly on the degrees of freedom. Inserting the expression of the first variations from equation (3.1), the second variation δ 12 t i can be rewritten as where the third-order tensor τ i = (I − t i ⊗ t i ) ⊗ t i and its generalized transpose are defined below equation (3.1). The expression of δ 2 t i announced in equation (3.1) is obtained by condensing δ 1 x i = δ 2 x i = δx i and identifying δ 2 t i = δ 12 t i .

B.4 Variation of directors rotation
In view of equation (B.1), the infinitesimal rotation vector δd i associated with the directors rotation d i is Differentiating the expression of d i from equation (2.8), we have δd i = δ p i r T i (ϕ i ) D i = δp i r T i (ϕ i ) D i + p i δ(r T i (ϕ i )) D i . Equation (2.2) shows that, with a fixed unit vector T i , δ(r T i (ϕ i )) = 1 2 δϕ i T i r T i (ϕ i )here, the vector in square bracket is an infinitesimal rotation vector, see equation (B.1). This yields δd i = δp i r T i (ϕ i ) D i + 1 2 p i δϕ i T i r T i (ϕ i ) D i . Inserting into the equation above, and using d i = D i r T i (−ϕ i ) p i from equation (2.8), we find as announced in equation (3.5).
The second-order infinitesimal rotation vector is then obtained from equation (B.2) as Here, we have used δ 12 ϕ i = 0 as ϕ i is a degree of freedom and the variations δ 1 ϕ i and δ 2 ϕ i are independent. Upon condensation of the two variations, the equation leads to the expression of δ 2d i announced in equation (3.5). Symmetrizing with respect to the independent variations δ 1 and δ 2 and using equation (B.2), we obtain the second infinitesimal vector as

B.5 Rotation gradient
Upon condensation of the two variations, the equation leads to the expression of δ 2q i announced in equation (3.6).

B.6 Strain vector
Equation (2.13) can be rewritten as κ i = 2 I(q i ), where I(q) = q−q 2 denotes the vector part of a quaternion. The operator I being linear, we have δκ i = 2 I(δq i ) = I (δq i q i ) as well as δ 12 κ i = 2 I(δ 12 q i ) as announced in equation (3.7). In the equation above, the second variation of the unit quaternion δ 12 q i has been expressed using equation (B.1).

B.7 Numerical verification
We verify the gradient and Hessian of the elastic energy, by considering a Kirchhoff rod having 80 nodes. Starting from a straight rod, we increment the magnitude of the natural curvature, magnitude of gravity, and a point load applied at the ends over 100 iterations. At each iteration we compute the equilibrium, disabling the update of the reference configuration discussed in Section 2.3. This allows us to verify the gradient in the generic setting where the reference and current configurations differ significantly from each other. The computed equilibrium solution is denoted by the vector X. We introduce a second configuration vectorX by adding a random perturbation to X where each perturbation is chosen randomly between (−0.1, 0.1). This magnitude of perturbation ensures that the configurationX is sufficienlty far from an equilibrium. By starting with the different equilibrium solutions X, we ensure that the variations are taken at different locations in the configuration space. The gradient of the discrete strain energy E = N −1 i=1 E i is evaluated at the pointX either as ∇E a computed based on the analytical formula given in the main text, or as ∇E f d using finite differences as (∇E f d ) i = (E(X + he i ) − E(X + he i ))/2h where h = 10 −7 and e i is a unit vector where the i th component is 1.
We then calculate the relative gradient error as: Similarly for the Hessian, we calculate the hessian ∇ 2 E of the strain energy gradient at the pointX, either analytically (∇ 2 E a ) using the methods described in the manuscript or using finite differences (∇ 2 E f d ).
We calculate ∇ 2 E f d using finite differences on the analytical form of the gradient. The relative hessian error is calculated as: At every iteration, we calculate a different random perturbation and calculate the errors ∇E err and ∇ 2 E err at that point. The results are shown in Figure (