Numerical solution of a bending-torsion model for elastic rods

Aiming at simulating elastic rods, we discretize a rod model based on a general theory of hyperelasticity for inextensible and unshearable rods. After reviewing this model and discussing topological effects of periodic rods, we prove convergence of the discretized functionals and stability of a corresponding discrete flow. Our experiments numerically confirm thresholds e.g. for Michell's instability and indicate a complex energy landscape, in particular in the presence of impermeability.


Introduction
Long slender objects-such as springy wires made of plastic or metalcan be approximated by curves. In many cases, equilibrium shapes are characterized in terms of the bending energy, i.e., (half of) the total squared curvature. The latter has a long history, dating back to Bernoulli, and can be seen as the starting points of elasticity theory.
The bending energy depends just on the centerline of an object and does not incorporate other physical effects such as twisting, friction, or shear.
For instance, only relying on the bending energy one cannot explain why a telephone cable tends to curl. It also does not preclude self-penetration.
In this paper we extend the study of inextensible elastic curves by the first author [3] to inextensible and unshearable elastic rods. To this end we discretize the minimization problem and devise a numerical scheme in order to simulate a suitable H 2 -like gradient flow.
Here c b , c t > 0 are bending and torsion rigidities that are determined by the Lamé coefficients of the material and geometrical properties of the rod. Furthermore, L rod bc : H 2 × H 1 → Y encodes the boundary data rod bc in some finite-dimensional linear space Y . We assume that it only involves linear combinations of boundary points of y, y , and b. Therefore L rod bc is continuous with respect to weak convergence in H 2 × H 1 . In particular L rod bc can be used to incorporate periodicity in case of a closed curve y.
For ease of readability, we will rescale I rod by 1/c t from now on and abbreviate κ = c b /c t .
We will always assume A to be nonempty which is guaranteed if the boundary data rod bc is compatible with the frame condition and implies that the distance between the endpoints is strictly less than L. Any boundary data on a frame F can be matched by adjusting a reference frame F 0 using a suitable cumulative angle function ϕ (see Section 2.3 below).
Elastic rods. Based on the work of Mora and Müller [49] for general rods, the minimization problem (P rod ) can be rigorously derived from a general three-dimensional hyperelastic model, see [4] for a short formal derivation. In the situation of rods with circular cross-section, made of some isotropic and homogeneous material, we find that c b ≥ 2c t . According to Coleman and Swigon [15, p. 195] there is some indication that values less than κ = 3 2 are appropriate for modeling DNA.
The study of elastic rods has a long history. It is closely related to elasticae, i.e., stationary points of the bending energy, see Levien [41] and references therein. A comprehensive presentation on the subject from the perspective of elasticity theory is provided by Antman [1].
We find applications in different fields such as the modeling of coiling and kinking of submarine cables (Zajac [71]; Goyal et al. [32,31]), cell filaments (Manhart et al. [44]), and computer graphics (Bergou et al. [7]; Spillmann and Teschner [60]). Modeling in molecular biology has stimulated a lot of activity in this field as well.
A prototypical model for DNA supercoiling which has received considerable attention is the twisted elastic ring investigated by Maddocks in various collaborations [22,37,37,43,45,46]. The solution of the corresponding minimization problem leads to an intrinsically straight uniform rod with equal bending stiffness. The analysis bases on Hamiltonian formulations of rod mechanics. The general idea is to impose a twist rate β on a unit loop.
For small values of β the round circle (with a uniform twist) remains an equilibrium. This phenomenon is known as Michell's instability [47], see Goriely [28] and references therein. Larger perturbations lead to instability and bifurcation phenomena, cf. Goriely and Tabor [29,30].
Ivey and Singer [36] reconsidered the problem from a variational point of view, obtaining a complete description of the space of closed and quasiperiodic minimizers. Recently, a reformulation in terms of symplectic geometry has been given by Needham [50]. Regarding the discretization of elasticae we refer to Scholtes et al. [58].
Gradient flows. We aim at numerically detecting configurations of framed curves with low bending and twisting energy. For this we consider the gradient flow of the energy functional in (P rod ), a weighted sum of an elastic bending energy term and a functional that tracks the twisting of the frame about its centerline.
We implement a constraint ensuring that the curves stay close to arclength parametrization if the initial curve is arclength parametrized. Moreover the bending energy can be replaced by the squared L 2 norm of the second derivative of the curve which is a crucial point in the analysis of the discretization.
Impermeability. Based on earlier work aiming at modeling DNA plasmids [14,17,66,68], Coleman and Swigon [15] take self-contact phenomena into account and discuss the interaction between certain topological quantities such as writhe, excess link, and the number of self-contact points. In [16] they include the case of (two-bridge) torus knots. In contrast to a related approach by Starostin [61] they impose a (small) positive thickness.
The corresponding case of open curves with appropriate boundary conditions has also been studied by several authors. Van der Heijden et al. [69] provide a comprehensive study of jump phenomena in clamped rods with and without self-contact. A more detailed classification of the respective equilibrium configurations is given by Neukirch and Henderson [51]. Clauvelin, Audoly, and Neukirch [13] modeled the situation of a small loosely knotted arc with open end-points and studied the shape of the set of self-contact and the influence of twist applied to the end-points as well. Starostin and van der Heijden [62] model the situation of so-called two-braids, i.e., structures formed by two elastic rods winding around each other, which also covers the case of (2, b)-torus knots [63,64]. The dynamic evolution of intertwined clamped loops subject to varying loads has been addressed by Goyal et al. [31,32]. Some of the above-mentioned models involve initial assumptions on the geometry, especially regarding the contact situation, focussing on explicit constructions for modeling and simulation. Our approach of treating (P rod ) does not rely on any precondition.
We will redefine (P rod ) in Section 6.6 to incorporate impermeability. To this end, we rely on the tangent-point energies whose impact on the evolution of (unframed) curves has been discussed in [5]. We thereby extend a regularization ansatz due to von der Mosel [70] with O'Hara's energies [52] in place of the tangent-point functional. In fact, one might conjecture that any self-avoiding functional will qualitatively produce the same results.
Computationally, this case is particularly challenging since strong forces related to bending effects have to by compensated by repulsive forces related to the tangent-point functional to avoid self-intersections. Regularization approaches guaranteeing global injectivity have been successfully implemented in different fields, see Krömer and Valdman [38] for an example in the context of elasticity.
The existence of curves minimizing (P rod ) follows via the direct method of the calculus of variations or, equally, from the Gamma-convergence (Proposition 4.1) together with the coercivity of the functionals. In the presence of uniform thickness bounds, however, we cannot rely on these reasoning, see Gonzalez et al. [27]. While (P rod ) covers the "uniform symmetric case" of the Kirchhoff rod, which constitutes maybe the simplest model that involves both bending and twisting, the setting discussed in [27] offers more flexibility and especially also covers the cases of extensible shearable rods. Schuricht and von der Mosel [59] derived the Euler-Lagrange equations for elastic rods with self-contact. A similar approach has been followed by Hoffman and Seidman [34,33].
The evolution of impermeable rods preserves isotopy classes, so topology aspects come into play. Here we will encounter a more involved picture compared to the analysis of the twist-free setting, see Langer and Singer [39] and Gerlach et al. [25]. A complete characterization of minimizers is wide open.
Outline. We review the geometry of elastic rods in Section 2. In Section 3 we derive an approximation result (Lemma 3.1) that is used in Section 4 in order to prove Gamma-convergence of the discrete problem to the continuous one (Proposition 4.1). We prove stability of the numerical scheme in Section 5 (Proposition 5.2). Several experiments discussed in Section 6 indicate a complex energy landscape.

Elastic rods
Here we provide a short presentation of the geometry of elastic rods which is inspired by Langer, Singer and Ivey [40,36]. It is not essential for the analysis of the numerical scheme in the subsequent sections but sheds some light on the interpretation of the experiments in the last section.
2.1. Framed curves. A rod is modeled by a curve y : [0, L] → R 3 which corresponds to its centerline and an orthonormal frame F : where × denotes the vector cross product. In the following we consider y ∈ H 2 and F ∈ H 1 .
We will assume that the first column of F coincides with the unit tangent The idea is that the directors b and d track the twisting of the material about the centerline.
The assumed energy regime for bending stiffness in (P rod ) imposes inextensibility as a physical property. Therefore we can prescribe arclength parametrization which leads to the unit tangent vector t = y and the curvature k = |t | = |y |.
Our analysis also covers the case of closed rods where [0, L] is understood to be the periodic interval R/LZ. We will realize the latter by imposing suitable (periodic) boundary conditions at 0 and L. In general a twist-free frame (see Section 2.3 below) of a closed curve will not close up, i.e., there can be a discontinuity at one point of R/LZ. One has to take care of this fact when defining boundary conditions. An important frame that will always be well-defined and continuous for sufficiently smooth (both open and closed) curves y with nonvanishing curvature is the Frénet frame where b F = t / |t |.
A rod is assumed to have some small diameter which can be considered infinitesimal; however, self-penetrations are not excluded at this stage (see Section 6.6 below for a discussion on modeling impermeability).

2.2.
Twist rate. Using the orthonormality of the frame, we may express the variation of the director b by The first term tracks the change of b that is imposed by the spatial behavior of the curve. It is just a component of the curvature vector as y = (y · b)b+ (y · d)d. Only the second one actually provides information about the twisting of the frame about the centerline. Therefore we will call b · d the twist rate of the frame.
We may also characterize a frame by a 9 × 9 linear system, namely for scalar coefficient functions k b , k d , and β where 0, 1 ∈ R 3×3 denote the zero and identity matrices. Here k b = y ·b and k d = y ·d are the components of the curvature k of y and β = b · d is the twist rate. For instance, the Frénet frame is characterized by k d ≡ 0.
A more detailed discussion of the impact of the twist rate is given in Section 2.4 below.
2.3. Reference frame. For any curve y, a point ξ ∈ [0, L], and b ∈ S 2 , b ⊥ y (ξ), we obtain by integration a unique frame ] whose twist rate is constantly zero. We call it synonymously a Bishop frame, natural frame, reference frame, or twistfree frame for y as it is a frame in rest position subject to a fixed curve. Therefore, up to a rotation of the initial vector b (which corresponds to an element of S 1 ) there is a unique twist-free frame for any given curve.
A twist-free frame F 0 provides a useful reference configuration. Denoting the (cumulative) angle between the director b of any other frame and b 0 by ϕ, we arrive at b = (cos ϕ)b 0 + (sin ϕ)d 0 and d = −(sin ϕ)b 0 + (cos ϕ)d 0 .
Consequently, the rate of change of ϕ is just the twist rate Two frames that just differ by a constant angle ϕ may be considered equivalent, in particular when modeling rods with a circular diameter where there is no natural choice of a director. This is of course different for small ribbons with lateral extension in a particular direction.
Note that even for closed curves with frames that close up, ϕ(L) − ϕ(0) does not need to be an integer multiple of 2π, unless the twist-free reference frame closes up. The latter applies in particular to rods with planar centerline where the vector being perpendicular to the respective plane provides a "canonical" twist-free frame.
In general, there is no direct correlation between ϕ(L) − ϕ(0) and the angle enclosed by b(0) and b(L). One can think of a revolute joint that controls the latter angle.
2.4. Total twist. An important quantity, in the literature often simply referred to as "twist", is the total twist (more precisely, total twist rate) where s is an arclength parameter and the last expression is parametrization invariant.
As the first identity suggests, the total twist can be interpreted as the number of rotations the director b (or, equivalently, d) performs about the curve, i.e., the centerline of the rod.
The total twist takes integer values on any closed curve for which both frame and twist-free reference frame close up. In particular, this holds for any planar closed curve with a closed frame.
2.5. Cȃlugȃreanu's identity. A given (sufficiently smooth) embedded closed curve y together with a closed frame b (i.e., there are no discontinuities of the frame) defines a link consisting of y and y + εb for some small ε > 0. For embedded closed curves, Cȃlugȃreanu's identity [11,12] Lk = Tw + Wr provides a decomposition of the Gauss linking number Lk (an integer topological invariant, a special case of the mapping degree) into the sum of two geometric terms (that can take arbitrary real values each), namely the total twist Tw and the writhe functional Wr. See Moffatt and Ricca [48] for an account on the history of this result.
The Gauss linking number amounts to half of the sum of all signed crossings of y and y + εb with respect to a regular projection direction. The latter guarantees that we can identify any crossing either as an over-or undercrossing, respectively. To this end, the intersection of (the images of) y and y + εb must be empty. It was Gauss' seminal discovery that this quantity can also be expressed by a double integral over [0, L], see, e.g., the nice review by Ricca and Nipoti [56]. For two curves y, y forming a link, such as y = y + εb, we have which is independent under reparametrization.
Any crossing of y(x) and y(x ) + εb(x ) where the preimages x, x ∈ [0, L] are close will be called local and global otherwise. Given any regular projection direction, we can distinguish between local and global contributions to the linking number, see Figure 1. Note that, in contrast to the linking number, its local and global contributions are not invariant for any regular projection direction in general. By Sard's theorem, almost every direction is regular, so we can just compute the average of local and global contributions over the sphere of all projection directions. The first one agrees with the total twist, the second one is the writhe functional. For details we refer to Dennis and Hannay [21] and references therein.
Writhe can equally be computed by summing up the signs of self -crossings of any projection of the curve and averaging over all projection directions. In particular, it only depends on the curve, not on the frame, but in contrast to linking number and total twist it requires y to be embedded. It can also be represented by the Gauss double integral via Wr(y) = Lk(y, y).
Writhe measures non-planarity (and non-sphericity) of a curve; in particular, it vanishes on planar (and spherical) embedded curves.
Consequently, the total twist will be close to the Gauss linking number for embedded curves that are nearly planar (which occur at several stages of the Of course, there is no "writhe-free" frame, but, as suggested by Maddocks, we can geometrically construct a "writhe frame" whose linking number is zero, see [21].
Compared to the other two functionals in Cȃlugȃreanu's identity, which can be evaluated using double integrals, the total twist is much easier to compute. Furthermore it does not require embeddedness of the curve.
2.6. Energies. We assume that the behavior of the rod is driven by a linear combination of the bending energy (half of the total squared curvature) and the twisting energy (half of the total squared twist rate), cf. Mora and Müller [49]. More precisely we consider the functional where s is an arclength parameter, k(s) denotes the curvature of y at y(s), and c b , c t > 0 are material constants. Minimizers are called elastic rods.
As mentioned in the introduction, we rescale the energy functional by 1/c t and define κ = c b /c t .
At the end of this section, we will briefly discuss two related minimization issues.

2.7.
Optimal frames. For a given curve y, we may consider the problem to find a director b minimizing I rod [y, ·] subject to the boundary condition L rod bc [y, b] = rod bc . In first place, if b is a stationary point of I rod [y, ·] for some fixed y then is constant due to du Bois-Reymond's lemma. According to the Cauchy-Schwarz inequality, the twisting energy is bounded below by 2π 2 L Tw(y, b) 2 . This minimum is attained if and only if (3) holds such that the twisting energy then amounts to 2π 2 L Tw(y, b) 2 = L 2 β 2 . In particular, we can check whether a given rod has a uniform twist rate by computing the quotient of total squared twist rate over squared total twist rate.
Note that for any global minimizer (y, b) of I rod , the director b is a global minimizer of I rod [y, ·] as well.
In case L rod bc [y, b] does not affect both points 0 and L, minimizing I rod [y, ·] is equivalent to constructing a twist-free frame.
Otherwise we face a clamped problem, i.e., b has to satisfy b (If the boundary condition just forces the frame to close up, i.e., b(0) = b(L), we may just let b − = b + for an arbitrary vector perpendicular to y (0) and y (L).) In this case there is a global minimizer b min with constant twist rate (3). Using a twist-free reference frame L . If ϕ(L) = π there are precisely two minimizers with twist rate ϕ ≡ ± π L . Topological restrictions can enforce arbitrary angles ϕ(L) ∈ R, however, this does not apply to (P rod ) which does not preserve this sort of condition throughout the evolution. Keeping track of topology enforces modeling impermeability-quite a natural feature which we will address in Section 6.6.
2.8. Releasing total twist. In light of Section 2.7 we must have |Tw(y, b)| ≤ 1 2 for any global minimizer (y, b) of I rod . In general, we have however, the absolute value of the total twist does not have to be decreasing throughout the evolution.
At the final stage of an evolution of a closed curve (the frame does not have to close up), all we can hope for, however, is |Tw(y, b)| ≤ 1. We briefly explain how this bound can be realized.
One can change the Gauss linking number of a given (embedded) rod by ±2 by locally forming a small loop, performing a suitable self-penetration and moving the curve back to the original position. The value of the writhe functional is not affected as it does not depend on the frame. So we have changed the total twist by ±2 as well according to the Cȃlugȃreanu identity (cf. Section 2.4).
A self-penetration of the curve will in general lead to a change in topology resulting in a discontinuity of the linking number. While the total twist is continuous throughout the evolution, the writhe functional is not welldefined on non-embedded curves and thereby compensates the change of the linking number.
An evolution does not necessarily realize the bound |Tw| ≤ 1. First of all, it is in general unclear whether it will in fact converge to a (local) minimizer at all. Another obstruction is discussed in the next section.
2.9. Michells instability. Among all closed curves, the round circle framed by its normal vector is the unique global minimizer of I rod (up to a constant rotation of the frame).
It is a remarkable fact that the round circle remains a minimizer (at least a local one, cf. Section 2.8) when we add some twist by increasing the (constant) twist rate β (which results in a discontinuity of the frame at one point). This phenomenon which is referred to as Michell's instability has been discovered 130 years ago [47] and then been rediscovered several times, see Goriely [28] for more details.
Zajac [71] has found the threshold β * = 2π √ 3κ/L that separates the stable and unstable regime. As before, L denotes the length of the curve. More precisely, the circular rod is stable as long as |β| < β * and unstable if |β| > β * . The dependency on κ = c b /c t is quite intuitive: If κ is very small, the bending energy dominates which always prefers the circle. A proof of Zajac's result adapted to our setting can be found in Ivey and Singer [36,Sect. 6].
Starting an evolution with β ini ∈ β * , 2π L , we have |Tw| = Tw < 1. Therefore the rod cannot reduce twist by self-penetration, so we will merely face some buckling of the rod-which is difficult to detect numerically. In Experiment 6.2 we chose κ = 3 2 , for which we measure a drastic change of the twisting energy by self-penetration of the curve.

Density
We can smoothly approximate any framed curve in A, i.e., A∩(C ∞ × C ∞ ) is dense in A with respect to the H 2 ×H 1 -topology, preserving given boundary conditions.
Proof. Our strategy is as follows. We first construct smooth approximizers (y δ , b δ ) ∈ C ∞ × C ∞ . In a second step we correct the boundary values by adding (smooth) functions v δ and c δ . The new curve y δ + v δ will not have length L. We balance the length by adding another smooth function w δ compactly supported in (0, L) \ supp v δ . Now we reparametrize the curve y δ +v δ +w δ to arclength and apply the same reparametrization to the vector field b δ + c δ . Renormalizing it by the usual Gram-Schmidt scheme produces the required director.
To begin with, we perpare the length correction. If |y(L) − y(0)| = L, the curve u just parametrizes the segment from y(0) to y(L), so y ∈ C ∞ and we only have to treat the director b (which can be done similarly as outlined below). If |y(L) − y(0)| < L we infer from H 2 ⊂ C 1 that y ∈ A cannot only move on a straight line. So we may assume that there is some constant unit vector v ∈ S 2 , v ⊥ (y(L) − y(0)) (this condition is empty for closed curves), such that y · v > 0 on some interval I + ⊂ [0, L]. Due to the fact that´L 0 y (x) · v dx = (y(L) − y(0)) · v = 0 there has to be another interval I − ⊂ [0, L] on which y · v < 0. Diminishing I ± if necessary, we may assume that they are both contained in (µ, L − µ) for some µ ∈ (0, L/2). Moreover we can assume that there is some λ ∈ (0, 1 2 ] such that ±y · v ≥ 2λ on I ± . We choose δ ∈ (0, δ 0 ] for some δ 0 ∈ (0, 1] which will be fixed later on only depending on (y, b) and ε. Using a standard mollifier, we obtain (y We may assume that ±y δ · v ≥ λ on I ± for all δ ∈ (0, δ 0 ].
In order to match the boundary conditions, we subtract suitable functions. More precisely, we let where C µ only depends on µ and C > 0 on the embedding H 1 → C 0 .
To make this more precise, let w = ωφ v for some ω ∈ R with |ω| ≤ λ φ C 0 .
Recalling that v δ and w have disjoint support, we infer if ω ≤ 0, if ω ≥ 0, which allows for the desired length correction depending on the sign of ω. More precisely, we can change the length of y δ + v δ by at least ±α where α = . Diminishing δ 0 if necessary, we can ensure that The embedding H 1 → C 0 guarantees that the curveȳ δ is immersed and min b δ ≥ 1 2 if δ 0 is small enough. So we may apply the reparametrization operator from Lemma 3.2 and let y δ =ȳ δ • ψ −1 and b δ − b H 1 tend to zero as δ 0. Using the continuity of the reparametrization and |y | ≡ 1 we find that y δ − y H 2 tends to zero as well. Choosing δ 0 sufficiently small, we may assume that y δ − y H 2 ≤ ε and b if δ 0 is sufficiently small. Indeed, using Leibniz rule vw H 1 ≤ v H 1 w H 1 and the fact that both y δ − y H 2 and b δ − b H 1 get arbitrarily small provided δ 0 is chosen accordingly, the same applies to Arguing as above, we find that the term 1 b δ H 1 is uniformly bounded.
sufficiently small. This allows for bounding 1 b δ as well. In the same way we prove boundedness of 1  that reparametrizes an immersed curve to constant speed is continuous with respect to the H 2 -norm. Moreover, A proof can be found in [55,Appendix]; the argument applies without rescaling and the additional requirement of embeddedness. It applies to non-periodic intervals [0, L] as well. The last statement can be derived from the formula.

Discretization
Our discretization is based on cubic and linear finite element spaces. We consider a partition of [0, L] by a set of nodes N h that contains the endpoints 0 and L. We define nodal bases (ϕ z ) z∈N h and (ψ z,j ) z∈N h , j = 0, 1 with the following properties. If z ± ∈ N h are neighboring nodes of z ∈ N h then ϕ z and ψ z,j are supported in [z − , z + ]. On the intervals [z − , z] and [z, z + ] the functions ϕ z are (affine) linear with ϕ z (z) = 1 while ψ z,j are cubic polynomials satisfying ψ z,j (z) = δ j,0 and ψ z,j (z) = δ j,1 , j = 0, 1. We define nodal interpolation operators on C 0 ([0, L]) and C 1 ([0, L]) respectively by letting Furthermore, we will employ the averaging operator Q h which is piecewise defined on any element [z, z ] (i.e., z, z ∈ N h are neighboring) by Both for the discrete approximation result and the stability of our numerical scheme presented in Section 5 below it is crucial to exploit the structure of the dimensionally reduced functionals.
To identify convex and concave terms, we observe that the orthonormality of the frame F = [t, b, d] implies b · b = 0 and b · t = −b · t . Therefore the integrand of the twisting functional becomes To obtain a coercivity property (under the restriction b h L ∞ ≤ 1) we set which ensures κ ≥ 2θ. This will allow for controlling (part of) the second term of I rod [y, b] by the first one even if κ < 2. Now we decompose By V h rod ⊂ H 2 × H 1 we denote the cross product of piecewise cubic and piecewise linear functions subject to N h . With the product finite element space V h rod and the operator Q h we consider the following discretization of the minimization problem (P rod ) in which the pointwise orthogonality relation y · b = 0 is approximated via a penalty term: The operators Q h and I 1,0 h are included in a way that leads to a simple assembly of the corresponding matrices and avoids quadrature.
Note that the constraints are imposed on particular degrees of freedom which makes the method practical. In order to ensure A h = ∅, the distance between the endpoints of y h must be less than L and the spatial mesh has to be chosen sufficiently fine.
For establishing the Gamma-convergence result, it is useful to write where P b and P Q h b h denote the square roots of the positive semidefinite The only difference to I rod is the penalization term and the fact that b is replaced by Q h b h in the third term and once in the fourth term.
Our first task is to show that minimizers of I h,ε rod within A h approximate I rod -minimizers within A. In the following statement, we assume that I rod and I h,ε rod attain the value +∞ outside of A and A h , respectively. Proposition 4.1. As ε, h 0, the functional I h,ε rod Gamma-converges to I rod with respect to the weak H 2 × H 1 -topology.
Note that there is no restriction on the ratio of ε and h.
Proof. To establish the Lim-inf inequality, we consider an arbitrary sequence ((y h , b h )) h>0 ⊂ A h and (y, b) ∈ H 2 × H 1 with y h y in H 2 and b h b in H 1 as h 0. In particular, we have y h → y in C 1 and b h → b in C 0 .
We have to show that if lim inf (h,ε) 0 I h,ε rod [y h , b h ] < ∞ then the limit point (y, b) belongs to A (such that the formula for I rod is applicable) and the lim inf-inequality holds.
As the mesh size tends to zero as h 0, the condition |y | = |b| = 1 is satisfied everywhere due to the uniform convergence y h → y , b h → b. Furthermore L rod bc is continuous with respect to weak convergence. It remains to verify that b is perpendicular to y . By the interpolation estimate we have As all terms of I h,ε rod are non-negative, the term ε −1´L 0 I 1,0 h (y h · b h ) 2 dx is uniformly bounded which implies´L 0 (y · b) 2 dx = 0. This gives y ⊥ b, so (y, b) ∈ A.
The second term on the right-hand side tends to zero as h 0 due to the boundedness of weakly converging sequences and This estimate also yields Now the lim inf-inequality follows from the lower semicontinuity of the L 2norm and the fact that the penalization term is non-negative (since it is the linear interpolation of a non-negative term).
We turn to the Lim-sup inequality. Let (y, b) ∈ A and δ ∈ (0, 1]. Of course, I rod [y, b] < ∞. We aim at constructing a recovery sequence. We apply Owing to the smoothness of the regularized frame, we have above by an expression that tends to zero. For the first three terms of For the last two terms in the previous line, we infer and, recalling (5), We treat the fourth term of I h,ε rod [y h , b h ] similarly as above. More precisely, we infer Finally we deal with the penalty term. Due to y δ ⊥ b δ we find that y h (z) · b h (z) = 0 for all z ∈ N h . Therefore I 1,0 h (y h · b h ) 2 (x) = 0 for all x ∈ I which implies that the penalty term vanishes on the entire sequence.

Iterative minimization
We linearize the pointwise constraints in our iterative algorithm for computing minimizers of I h,ε rod . For a vector field y h ∈ S 3, The functionals L rod bc,y and L rod bc,b are the components of L rod bc corresponding to the variables y and b, respectively, assuming that the boundary conditions can be appropriately separated.
In what follows we abbreviate the partial torsion term weighted by the factor (1 − θ) by the nonlinear operator The nonlinear term N h does not occur if κ ≥ 2. In this case the negative contribution in the energy functional can be treated using its separate concavity properties. Otherwise, an inductive argument is used in the stability analysis which requires a different treatment. We therefore set We abbreviate We generate a sequence (y k h , b k h ) k=0,1,... that approximates a stationary configuration for I h,ε rod with the following algorithm. It approximates a gradient flow that is determined by the metrics (·, ·) and (·, ·) † . These can be chosen quite general, however, our stability result relies on certain embeddings, see (6) below.
Algorithm 5.1 (Gradient descent for elastic rods). Choose an initial pair (y 0 h , b 0 h ) ∈ A h and a step size τ > 0, set k = 1.
Note that the N h -terms on the right-hand sides vanish in case θ = 1 which corresponds to κ ≥ 2.
It is useful to view d t y k h and d t b k h as the unknowns in Steps (1) and (2) in- The algorithm exploits the fact that the penalty term is separately convex while the nonquadratic contribution to the torsion term is separately concave. Therefore, the decoupled semi-implicit treatment of these terms is natural and unconditionally energy stable if θ = 1, i.e., in the bending-dominated case.

Proposition 5.2 (Convergent iteration). Assume that we have
If c 0 τ ≤ 1/2 then the iteration is energy decreasing, convergent, and the unit-length violation is controlled via where e 0,h = I h,ε rod [y 0 h , b 0 h ] < ∞ and c , † > 0 only depends on the metrics. If θ = 1 then we may choose c 0 = 0, i.e., stability and convergence hold unconditionally.
Remark 5.3. In case θ = 1, i.e., in case of low torsion rigidity, we actually only require the second line of (6). Moreover, if we even replace it by w h ≤ c w h and r h ≤ c † r h † the estimates for the constraint violation still hold in L 1 instead of L ∞ .
In general, the condition (6) can be satisfied if · and · † are H 2 -and H 1seminorms and if L rod bc imposes suitable Dirichlet conditions on one endpoint of the interval. An L 2 -gradient flow, however, requires stronger assumptions on the step size.
Proof of Proposition 5.2. (a) We first consider the case θ = 1 for which k = k and the nonlinearities related to the operator N h disappear and the asserted estimate holds with c 0 = 0. For this we note that the functional is separately convex, i.e., convex in y h and in b h . Therefore, we have that , which by summation leads to the inequality is separately convex and we have (1) and (2) of Algorithm 5.1 and find that . Since for θ = 1 we have that

Repeated application leads to
Using (6) and the previously established bound for the discrete time derivatives proves the estimate for the constraint violation if θ = 1.

Experiments
Here we report on some numerical experiments. The parameters that have been used are shown in Table 1. The length of the curve is denoted by L, the number of nodes by N , the maximum step size h max is normally close to L/N where L is the length of the curve. Except for Experiment 6.1, the initial curve has constant twist rate β ini . While the director b will always be clamped on both ends of the rod, the boundary conditions (bc) for the curve are either periodic (p) or clamped on both ends (c). Unless otherwise stated the penalization parameter has been chosen to be ε = 2π N . We always use the time step size τ = 1 8 h max . In some cases we also added some small perturbation to the initial curve which is specified in the text.
We briefly comment on the legend of the corresponding energy plots where the horizontal axis shows the iteration steps of the evolution. Of course, "bending" refers to the scaled bending energy κ 2 y h 2 L 2 , "twisting" to the scaled twisting energy θ while "total twist" means the functional Tw defined in Section 2.4, "potential" to the tangent-point energy (see Section 6.6), "penalizing" to P h,ε , and "total" is I h,ε rod . In general the total twist is scaled differently, i.e., there is a second axis on the right margin of the energy plots while all other values refer to the axis on the left margin of the coordinate system. A coloring of a curve represents curvature values.
For typical discretizations with 400 nodes the overall runtime of our implementation in Matlab on a server (3.1 GHz) took about 3 minutes per 10,000 iteration steps. An impermeable rod with 800 nodes in Experiment 6.6 (b), with an assembly of the self-avoidance potential in C, requires about 15 minutes per 10,000 iteration steps.
We would like to stress that we observe energy monotonicity for all our experiments, confirming the stability features discussed in Proposition 5.2. 6.1. Uniform twist rate. According to Section 2.7, stationary frames have constant twist rate. We expect that a non-uniform twist configuration will become constant within the evolution.
We start with a straight line curve of length 2π with uniform twist rate 4 on [0, π 2 ] and 0 on [ π 2 , 2π]. No perturbation of the initial rod was used. Initial stage and some intermediate steps from the evolution are visualized in Figure 2. After 4,000 steps the twist rate is nearly constant and amounts to 1.
All energy values are neglectable except for the twisting energy which virtually coincides with the total energy. Furthermore, at initial and final step the twisting energy data matchs quite closely the expected values of 1 2´2 π 0 β(s) 2 ds which amount to 4π and π. Looking at the following experiments whose initial configurations all have a uniform twist rate, we find this property being violated in the course of the iteration. In first place, this is due to the spatial behavior of the curve which does not seems to allow for an simultaneous reaction by the director. Eventually, uniform twist rate is restored, at least when the evolution has reached a stationary configuration, cf. Section 2.7.
We can test for uniform twist rate by computing the quotient of 2π 2 L Tw 2 over the twisting energy. As to Experiments 6.2 to 6.4, throughout the evolutions this number stays close to 1 for most of the time where we detect a relative deviation below 1 200 . The twist rate for Experiment 6.3 is plotted in Figure 6 (right). In Experiment 6.4 we ignore the initial steps where the twist is zero. In the other cases, we also see that the twist rate eventually dissipates but it can last a relatively long time.  Left: In order to experimentally confirm Zajacs threshold, we repeat the evolution depicted in Figure 3 above for different twist rates β ini > β * . The evolutions turn out to be very similar; essentially they only differ in speed. The (logarithmically scaled) twisting profiles reveal a drastic reduction of the twisting energy (due to the self-penetration of the curve). The smallest iteration step at which the twisting energy is below 27 4 π serves as a threshold for the speed of the evolution. Right: The energy plot for the evolution from Figure 3 where β ini = 4.2. Here the twisting energy decay occurs around iteration step 114,200. 6.2. Michells instability. We aim at experimentally confirming Zajac's threshold β * = 2π √ 3κ/L, see Section 2.9.
We consider the initial curve y(s) = (cos s, sin s, 0) with the frame b(s) = cos(β ini s) In order to break symmetry which seems to prevent rod configurations from leaving even energetically unfavorable states, a slight perturbation has been applied to the initial curve, namely (9) s → 1 1000 sin(7s) perpendicularly to the plane (i.e., in z-direction); the frame is corrected accordingly.
In order to quantitatively verify Zajac's threshold, we had to choose a rather high penalty coefficient, namely ε = 10 −5 . For κ = 3/2, we obtain In Figure 4 (left) we plot the twisting energy of several evolutions using logarithmic scales on both axes. From top to bottom, the profiles correspond to initial values of β ini = β * + 2 10 for = −1, 0, 1, 2, 3, 4. More precisely, we plot the twisting energy of the evolutions corresponding to different values of β ini minus the twisting energy of the configuration at β * , i.e., 1 2´2 π 0 β(s) 2 ds− 27 4 π, which seems to be stationary. Of course, values less than 27 4 π ≈ 21.2058 are ignored.
Experimentally we find that evolutions for different initial values of β ini seem to be very similar and essentially only vary in speed. The region of iteration steps where the twisting energy drastically falls is a good indication for the latter. There is one caveat-probably due to symmetry it turned out that the evolution corresponding to the initial values β ini = 3 is remarkably slower than expected. Therefore we chose β ini = 2.98 instead.
The plot in Figure 4 (left) indicates a reciprocal dependence of the evolution speed on β ini − β * . For values β ini ≤ β * the initial configuration remained unchanged (at least until step k = 5 · 10 6 ). This confirms Zajac's threshold as desired.
We show a typical full energy profile in Figure 4 (right) for the initial value β ini = 4.2. Some iteration steps are visualized in Figure 3. The corresponding plots for the other initial values of β ini essentially differ by the scaling of the horizontal axis and the simulations looks very similar. Initial and final stage are round circles which correspond to the (unique) global bending energy minimum of L = 2π among all closed curves.
Note that the total twist is reduced by precisely 2 from β ini = 4.2 to β = 2.2. In light of Section 2.8 a second reduction would be possible as well. However, in contrast to Experiment 6.3 below, the gradient of the energy does not seem to be steep enough to invest the amount of additional bending required for  Figure 5. Evolution of a round circle twisted by five full rotations in Experiment 6.3. Twist is reduced due to selfpenetrations of the rod around iteration steps k = 46,400 and k = 48,800. Initial and final curves are round circles which appear to lie in the same plane. Plots of energy profile and twist rate are shown in Figure 6 below. another self-penetration. As 2.2 < β * the evolution becomes stationary due to Michell's instability.
6.3. Reducing twist by self-penetration. We repeat the experiment from Section 6.2 with κ = 2 and β ini = 5, i.e., for a continuous frame.
Here we have β * ≈ 3.4641. The same slight perturbation has been added to the initial curve as before in (9).
In this case we observe twist reduction by two consecutive self-penetrations although the evolution could possibly stop after the first one since the twist value is then already below the threshold β * . Obviously, the gradient of the energy of the noncircular configuration around iteration step k = 48,000 is so steep that Michell's instability does not play any role here.
The evolution is visualized in Figure 5 and the energy values are plotted in Figure 6 (left). As the initial and final configurations are circular and the frame closes up (because of β ini ∈ Z), the total twist is integer at the beginning and end of the evolution (cf. Section 2.5). Right: Plot of the twist rate β for the iteration steps visualized in Figure 5. Apart from the initial configuration, the twist rate is non-uniform throughout the evolution. It eventually dissipates around iteration step 100,000. . Evolution of a twist-free planar elastic figure eight from Experiment 6.4 with κ = 0.7. The initial curve is (almost) planar and evolves to a circle in a plane which seems to be perpendicular to the initial configuration. The corresponding frame performs one full rotation. Obviously the bending forces dominate in this case. The energy plot is shown in Figure 8 (right).
The twist rate, however, does not stay uniform throughout the evolution as can be seen from Figure 6 (right). Eventually the twist will be balanced similarly to Experiment 6.1.
6.4. Planar figure eight. Any closed planar elastica (i.e., a critical point of the bending energy) is either a circle or a planar figure-eight curve, possibly several times covered, cf. Sachkov [57].
Explicit formulae for elastica based on special functions have been computed in the 19th century, see the references in Levien [41]. Here we make use of an arclength parametrization given by Dall'Acqua and Pluda [19] which relies Here E denotes the incomplete elliptic integral of the second kind and K the complete elliptic integral of the first kind while am is the Jacobi amplitude function and cn the elliptic cosine function, cf. [19]. The (signed) curvature amounts to s → −2 √ m cn(s, m). The figure-eight curve corresponds to m ≈ 0.82611 which is the uniquely defined number in (0, 1) with 2E( π 2 , m) = K(m) ≈ 2.321.
In order to break the symmetry which may result in an unstably stationary configuration, a slight perturbation similarly to (9) has been added to the initial curve, namely s → 1 1000 sin 7 · 2π 4K(m) s , perpendicularly to the plane (i.e., in z-direction); the frame is corrected accordingly.
According to Ivey and Singer [36,Sect. 7] the twist-free planar figure-eight is stable if κ = c b /c t < 1 2 and unstable if κ > 1 2 . We performed several evolutions whose energy plots are depicted in Figure 8 (left). As in Experiment 6.2, the evolutions are very similar and essentially only differ in speed. In each case the total twist is raised from zero to one which is still in accordance with the bound |Tw| ≤ 1 in Section 2.8. Mind the semi-logarithmic scaling of the horizontal axis. Negative values of the twisting energy are due to discretization errors and tend to zero when choosing a larger number of nodes.
Snapshots of the evolution for κ = 0.7 can be found in Figure 7. We observe an evolution to a round circle with one full twist. For the same reason as in Experiment 6.3 we face integer values of Tw at beginning and end of the evolution. The full energy plot is depicted in Figure 8 (right).
The parameters κ of the profiles shown in Figure 8 (left) amount to κ = 1 2 + 2 10 for = −2, −1, 0, 1, 2, 3, 4. The red solid line corresponding to κ = 0.52 is just about to lift at the right margin. The speed of the evolution seems to reciprocally depend on κ − 1 2 , suggesting that the evolutions for κ ≤ 1 2 will be stationary. This confirms the threshold by Ivey and Singer as desired. This curve is not parametrized by arclength, with the notation of elliptic integrals introduced in Experiment 6.4 we have L = 4 √ 2E( π 2 , −2) ≈ 12.357. No perturbation of the initial rod was used.
Choosing β ini = 4 gives an initial total twist of 8 according to Section 2.4.
The evolution is depicted in Figure 9, the corresponding energy plot can be found in Figure 10. It seems to become stationary after 30,000 steps although the total twist could be further reduced, see Section 2.8. 6.6. Implementing impermeability. In order to preclude rods from selfpenetrations we consider the modified total energy I rod + TP where ≥ 0 and TP denote the tangent-point functional (10) u → 1 2 q qˆL 0ˆL 0 ds ds r(y(s), y(s)) q , q > 2.
Here s,s denotes arclength parameters, and r(y(s), y(s)) is the radius of the circle that is tangent to y at the point y(s) and that intersects with y in y(s).
As many so-called knot energies [53], the tangent-point energies provide a monotonic uniform bound on the bi-Lipschitz constant. This implies in particular that the energy values of a sequence of embedded curves converging to a curve with a self-intersection will necessarily blow up.
The tangent-point energies have been proposed by Gonzalez and Maddocks [26]; the scale invariant case q = 2 has already been introduced by Buck and Orloff [10]. They are defined on (smooth) embedded curves y : [0, L] → R n and take values in [0, +∞], see Strzelecki and von der Mosel [65] and references therein. Blatt [8] has characterized the energy spaces in terms of Sobolev-Slobodeckiȋ spaces; regularity aspects are discussed in [9].
More information on the discretization of the tangent-point functional can be found in [5,6]. We cut out a strip of radius 2h max off the diagonal in [0, L] 2 . Figure 9.
Evolution of an open clamped rod from Experiment 6.5. The total twist amounts to 8 initially and is then reduced to about 2 by two self-penetrations which occur around iteration steps 6,000 and 16,500. The energy plot is shown in Figure 10 below.
We repeat Experiments 6.3 and 6.5 in the presence of self-avoidance. No perturbation was added to the initial curves. The parameters are chosen similarly, see Table 1.
Note that for closed curves with a closed frame the linking number is preserved throughout the evolution. Changes of the total twist will be entirely compensated by the writhe functional.   Figure 11. In the first part of Experiment 6.6 we repeat Experiment 6.3 in the presence of impermeability. As selfpenetrations are excluded, we observe the formation of coilings. The energy plot is shown in Figure 13 (left). In the second part of Experiment 6.6 we repeat Experiment 6.5 in the presence of impermeability. As self-penetrations are excluded, we observe the formation of a plectoneme. The viewer's position has been rotated by 90 degrees. The energy plot is shown in Figure 13 (right).