Chaire d'Analyse Appliquée

Computational and theoretical developments in rod mechanics

Table of Contents

The problem: static equilibria of elastic rods

An elastic rod can be described by (r(s),d1(s),d2(s),d3(s)), where:
s is arclength along the rod,
the centerline r(s) in R3 runs through the centers of mass of the rod cross-sections,
the frames (d1(s),d2(s),d3(s)) in SO(3) are the principal axes of inertia of the rod cross-sections.
If the rod is inextensible, then d3(s) = r'(s), so the rod is completely specified by r(s) (shown below as a green tube) and a normal vector field d1(s) tracking the rod's twist (depicted below as a blue ribbon).

The energy of a rod configuration is a functional of the form:

I have been studying the possible equilibrium configurations of an inextensible rod subject to various boundary constraints; for example, the twisted ring, in which the centerline closes smoothly and the ribbon forms an angle  at the point of closure (see below). Mathematically, this is a calculus of variations problem: we seek critical points of the functional E subject to prescribed boundary conditions on (r,di). These critical points are found via the classic Euler-Lagrange equations for E, which yield a two-point boundary value problem for a system of ODEs. Since we want to solve this 2-point BVP as various parameters are varied (either  or various parameters appearing in W), our numerical computations use the parameter-continuation package AUTO .

The energy functional

The function W incorporates all material-dependent properties (e.g., intrinsic curvature, bending stiffnesses, twisting stiffness). For example, a common energy model used today is:

where the strains ui are defined by:

and the expression  denotes the strains defined from the rod's intrinsic shape. Thus,  and  describe intrinsic curvatures while  describes intrinsic twist. The coefficients Ki are called the stiffnesses of the rod (bend stiffnesses for i=1,2 and twist stiffness for i=3). The computations described below can be easily extended to more elaborate energy models, such as ones with higher-than-quadratic terms or twist-bend coupling. 

Solution set for "perfect" rod (cf. Li & Maddocks, 1998 )

For the simplest case in which Ki are independent of s, K1=K2, and the  are 0, which we call the perfect rod (physically: uniform, circular cross-section, and intrinsically straight), the set of solutions to the twisted ring problem is highly symmetric. As  varies the solutions sweep out the bifurcation diagram shown below. The energy of each solution is plotted against the torque m3 at s=0 (both scaled to be unitless). Coloring indicates stability of the solutions, with green denoting local minima of E, red denoting saddles with one negative direction, etc. (see section on stability below for further information).

Degenerate orbits of solutions to the perfect problem (cf. Manning & Maddocks, 1998 )

Due to symmetries of the perfect rod, each point on the above bifurcation diagram actually represents an entire manifold (or orbit) of twisted ring solutions (all with the same energy and torque). For example, see below. Here, we take a given rod equilibrium, and in Step 1 do a rigid-body rotation by an angle  (about 45 degrees in this case), and in Step 2 rotate the rod at every point about its tangent vector by . In the end, we will still be at an equilibrium of the energy (because of the isotropy of the rod) and the boundary conditions on r and di are the same as in the original.
  
On the parabolic branches in the perfect diagram, the orbit of degenerate solutions is homeomorphic to a circle (generated by the operation shown above for 0    ). On the other branches, the orbit is homeomorphic to a torus (generated by the operation shown above for 0    , plus another operation relating to translation of the rod along its arclength). 

Splitting of bifurcation diagram for imperfect problems (cf. Manning & Maddocks, 1998 )

If we remove the symmetries of the perfect rod, say by introducing intrinsic curvature into the rod, then the orbits of degenerate solutions for the perfect rod yield in general a finite set of solutions, each with different energy, to the imperfect problem. In fact, through a perturbation computation, we have determined exactly which points on the orbit yield solutions to the imperfect problem. Most commonly, the parabolic branches in the perfect diagram (whose orbits are circles) yield 2 branches in the imperfect diagram, while the remaining branches (whose orbits are tori) yield 4 branches. See for example, the diagram below.

We can make even stronger statements using the perturbation expansion. For example, for infinitesimal intrinsic curvatures  , a solution on a parabolic branch in the perfect diagram always yields two solutions in the imperfect diagram, unless the following two conditions hold:

For non-infinitesimal perturbations, and especially those nearly satisfying the above conditions, the nice "generic" splitting shown above need not occur. Thus, the perturbation expansion suggests which intrinsic shapes will lead to more unusual imperfect diagrams. See for example, the diagram below. To better illustrate the unusual splitting, just the bottom portion of the diagram is shown, and a combination of energy and bending moment m1 is plotted on the y axis.

In the category of half-science, half-entertainment, here are a few movies demonstrating the transition from perfect to imperfect diagrams. Be forewarned that they're somewhat large (400 KB to 1 MB).

Movie 1: a superposition of the perfect diagram (in red) and an imperfect diagram (in purple and yellow) as the imperfection (which in this case intrinsic curvature of a 157-base-pair DNA) is gradually turned on. As above, energy is plotted against m3.

Movie 2: same as Movie 1 (but without the perfect diagram superimposed) for a larger imperfection representing protein-bound DNA. The dots mark solutions for which d1(0)=d1(1).

Movie 3: same as Movie 2, except that m3 is plotted versus .

Stability of DNA equilibria (cf. Manning, Rogers & Maddocks, 1998 )

The above computations were all made by solving the Euler-Lagrange equilibrium equations for the elastic rod numerically. Therefore, the computed equilibria are critical points of the strain energy functional but not necessarily stable local minima. We have thus designed a conjugate point test which determines a necessary condition for a given equilibrium to be a local minimum. This test is a generalization of the standard Jacobi/Morse conjugate point theory in the unconstrained calculus of variations to the isoperimetrically constrained calculus of variations problem which describes the twisted ring. The number of conjugate points is called the index, and should give the number of downward-turning directions on the potential energy surface in the neighborhood of the critical point (assuming the appropriate sufficiency conditions hold in addition to the necessary condition we have just discussed)

In the diagrams above, the index is denoted by color: green = 0 (local minimum), red = 1 (saddle with one downward-turning direction), light blue = 2, orange = 3, purple = 4, and dark blue = 5.

In the perfect diagram, each critical point on a parabolic branch also contains a single flat direction (in addition to the number of downward-turning directions indicated by the color); this flat direction arises exactly because of the existence of a degenerate orbit of solutions. When a perturbation is added, this flat direction can either turn upward or downward. Thus, in the generic imperfect diagram, one of the two children of the perfect parabolic branches always has the same index as the corresponding perfect branch (in this case, the flat direction perturbed to an upward-turning direction), whereas the other child has index one higher (in this case, the flat direction perturbed to a downward-turning direction).

Similarly, critical points on the nonparabolic branches contain two flat directions. In the generic imperfect diagram, one of the four children has the same index as in the perfect diagram (both flat directions perturb up), two of the children have index one higher (one flat direction perturbs up, one down), and the remaining child has index two higher (both flat directions perturb down).

Self-contact

An ongoing project is to extend the above equilibrium and stability computations to incorporate rod self-contact. The mathematics and computation are significantly more challenging since integro-differential equations are involved. With self-contact incorporated, a rich variety of equilibrium structures will be possible, including stable knotted configurations and self-interwound plectonemes.
 

Return to the homepage of "Chaire d'Analyse Appliquée"
Send your email or comments about this page to Patrick Furrer
last update : PF, july 8, 1998