Christian Hafner
Bernd Bickel
Abstract
The Kirchhoff rod model describes the bending and twisting of slender elastic rods in three dimensions and has been widely studied to enable the prediction of how a rod will deform, given its geometry and boundary conditions. In this work, we study a number of inverse problems with the goal of computing the geometry of a straight rod that will automatically deform to match a curved target shape after attaching its endpoints to a support structure. Our solution lets us finely control the static equilibrium state of a rod by varying the cross-sectional profiles along its length.
We also show that the set of physically realizable equilibrium states admits a concise geometric description in terms of linear line complexes, which leads to very efficient computational design algorithms. Implemented in an interactive software tool, they allow us to convert three-dimensional hand-drawn spline curves to elastic rods and give feedback about the feasibility and practicality of a design in real time. We demonstrate the efficacy of our method by designing and manufacturing several physical prototypes with applications to interior design and soft robotics.
1 INTRODUCTION
Rodlike elements are ubiquitous in natural and human-made environments and form the building blocks of the world around us from beams used in the construction of buildings, machines, and tools to electrical cables, clothing fibers, tree branches, human hair, muscle tissue, and DNA strands. This omnipresence has led to a long and rich scientific history concerned with developing the tools to predict how an elastic rod of a given length and shape will bend and twist under different loads and constraints.
The history of investigating the inverse problem—to determine the initial shape of a rod whose equilibrium state is known—is considerably younger but no less ambitious: It lays the theoretical foundation for designing rods with a programmable deformed shape that is hard coded into its geometry.
Motivation. Most rods, beams, and ribbons used in construction and design are straight in their undeformed state and acquire curvature only through the process of bending and twisting. Straight rods are preferred, because they can be manufactured more cheaply and with less material waste and packed more easily for transport. To enable the inverse design of straight rods with a programmable curved equilibrium shape, it is useful to know exactly which equilibrium shapes are physically realizable in the first place.
This problem was recently solved by Hafner and Bickel [2021] for the special case of planar deformations, leading to a concise description of all plane curves that appear as static equilibria of straight elastic rods with spatially varying thickness. This result also gives rise to a computational design algorithm, which solves for the thickness distribution necessary to achieve a given curved shape, and has applications in architecture and design.
A limitation of the aforementioned work is the restriction to planar shapes and the inability to reproduce shapes that are curved in three dimensions (3D). In this work, we show that a geometric characterization can also be attained for the case of fully three-dimensional bending and twisting of straight elastic rods with spatially varying cross sections. This allows the computational design of structures such as the ones shown in Figure 1, in which the curvature of the rods is induced purely by the elastic response of the material and the cross sections computed by our algorithm.
Fig. 1. Photographs of physical prototypes. Our computational design algorithm takes as input a curve in \(\mathbb {R}^3\) and computes the geometry of a straight elastic rod that will deform to match the input curve once installed in a support structure. Left: Six silicone rods that match six input curves subject to elastic forces and gravity. Right: Free-form light sculpture made from three straight silicone rods and electroluminescent wire, attached to a 3D-printed fixture.
Problem Statement. We study the design problem associated with rods of vanishing natural curvature from a geometric and from an algorithmic perspective. To formalize this problem, we use the large-displacement small-strain model developed by Kirchhoff and Clebsch, called the Kirchhoff rod model for short.
The deformed state of a Kirchhoff rod can be described as a framed curve in \(\mathbb {R}^3\), so a natural question to ask is which framed curves occur as deformed equilibrium states of rods—supposing we can freely choose the cross-sectional profile at every point. We study three flavors of this question, which differ by how much of the deformed state is prescribed, i.e., whether the twist of the rod is constrained to vanish, constrained to be a prescribed function, or not constrained at all.
Design Algorithm. Our main theoretical contribution is to show that there is a close connection between these sets of equilibrium curves and classical projective line geometry. This connection leads to a characterization of these curves that is both concise and computationally convenient. Furthermore, it directly translates into an algorithm that checks whether a given curve has the equilibrium property and lets us construct the geometry of a rod that will realize this curve as one of its equilibrium states.
This algorithm takes as input a design consisting of a number of curves in three-dimensional space, which we would like to turn into a physical model where elastic rods take the place of these curves. For each curve, the algorithm outputs the geometry of a straight rod with a spatially varying cross-sectional profile. This rod will adopt the shape of the input curve at static equilibrium once its endpoints are mounted in a support structure at the correct locations, orientations, and twist, as illustrated in Figure 2.
Fig. 2. Algorithm Overview. Left: The input to our computational design algorithm is a curve in \(\mathbb {R}^3\) , which we want to convert to a deformed elastic rod. Center-left: We show that this problem has a solution if and only if there exists a helical motion (green) that is nowhere orthogonal to the binormal lines (yellow) of the curve. Center-right: Based on the parameters of the helical motion, we compute a moving frame along the curve and the geometry of a rod that solve the Kirchhoff rod equilibrium equation. Right: Photograph of a silicone rod that deforms according to the computed solution in the real world.
The basic version of our algorithm, which applies when the gravitational effect on the deformed shape is negligible, comes with the rigorous guarantee that it terminates successfully if and only if the input is feasible. This is in contrast to existing iterative form-finding algorithms, which cannot guarantee that existing solutions will always be found, as they may encounter local minima in the complex energy landspace of Kirchhoff rods or terminate early due to an excessively shallow gradient close to a minimum. Our feasibility test, form-finding, and geometry generation take 10 ms all together, so they are fast enough to be part of an interactive design pipeline.
Extensions. In preparation of our main results, we formulate and show two interesting properties of the design space of Kirchhoff rods. First, we prove convexity of the set of stiffness matrices that describe the resistance to bending and twisting at a point of the rod. This indicates that many design problems involving Kirchhoff rods can be posed as convex optimization problems. Second, we investigate how using cross sections with different shapes affects the equilibrium curves that can be achieved. We prove that elliptical cross sections suffice to reach every possible deformed state, and other cross sections do not add any more design freedom.
Finally, we show how external forces and the dead load of a rod can be modeled as part of the inverse problem, thereby opening up applications in design on a larger scale. To take these forces into account, we present an iterative algorithm that post-processes the result obtained with the basic, load-free algorithm. Convergence of this algorithm is at least linear on all of our examples, as we show by empirical means.
We apply our design algorithm to several fields of application, such as fixture design, interior design, and soft robotics. In this context, we discuss our fabrication pipeline for producing physical copies of rods using 3D-printed molds and silicone casting. The efficacy of our approach is validated by comparing photographs of our manufactured examples to renderings of the target designs.
2 RELATED WORK
Our work fits into a body of research on fabrication-aware design, a field that aims to provide computational tools for the design of physical artifacts that are subject to special constraints derived from its material, function, or fabrication pipeline. Elastic rods are popular building blocks in computational fabrication, and they have been studied for their ability to take on complex shapes through bending and to approximate surfaces through the formation of networks. We give an overview of these lines of research, followed by an account of works that study mathematical properties and physical limits of the design space of rods, which is the primary concern of this text.
Fabrication-aware Design. The inverse design of free-form objects from elastic materials has received wide attention since the early 2010s and covers a range of applications, such as inverse shape design under gravity [Chen et al. 2014], design of objects with prescribed natural frequencies [Bharaj et al. 2015], and design of curved surfaces with elastic membranes [Pérez et al. 2017; Guseinov et al. 2017]. Other strategies for surface fabrication include modeling the deformation of wire mesh with Chebyshev nets [Sageman-Furnas et al. 2019], approximation with Origami [Dudte et al. 2016], and inflation of air channel networks [Panetta et al. 2021]. A comprehensive overview of fabrication-aware design can be found in recent surveys [Bermano et al. 2017; Bickel et al. 2018].
Design of Rod Networks. Grids of elastic rods have been used in an architectural context for the economical construction of curved façades and pavilions. One successful concept is that of a gridshell, an elastic lattice that is assembled in a planar configuration and curved during deployment [Lienhard et al. 2013]. Gridshells have been built based on asymptotic nets [Schling et al. 2018], special geodesic nets and sliding joints [Pillwein et al. 2020; Pillwein and Musialski 2021], and using deployment simulation with inverse design optimization [Panetta et al. 2019]. Another approach for constructing curved surfaces from flat elements at scale is the use of spiraling microstructures [Malomo et al. 2018].
Outside of architecture, rod networks have also been used to approximate deformable surfaces with multiple target poses at a low fabrication cost [Pérez et al. 2015] and to cover a surface in a structurally sound, decorative network of user-defined patterns [Zehnder et al. 2016]. A recent application of ribbon networks concerns the design of tri-axially woven free-form surfaces, based on near-geodesic families of curves [Vekhter et al. 2019] or using naturally curved ribbons [Ren et al. 2021].
Design of (and with) Rods. Solutions to design problems are often based on forward simulation methods, which have been studied extensively. The earliest rod simulation method in graphics literature treats the static clamped-free case [Pai 2002], while the first dynamic rod simulation method is based on a piecewise-helical discretization [Bertails et al. 2006]. The popular discrete elastic rod model supports arbitrary boundary conditions, constraints, and anisotropic cross sections [Bergou et al. 2008; , 2010].
Rod simulation is frequently used to compute the dynamics of hair in computer animation. Modeling hair styles that look lifelike under gravity is among the earliest inverse design problems studied in computer graphics [Hadap 2006], with later works taking into account collisions and frictional contact [Derouet-Jourdan et al. 2013]. In computational fabrication, rod models have been employed to predict the elastic response of wires for instance. Applications include the design of structurally stable wire sculptures [Miguel et al. 2016] and cable-actuated wire characters with multiple target poses [Xu et al. 2019b], as well as optimizing the motor trajectories of robotic wire characters to minimize mechanical oscillations [Hoshyari et al. 2019]. Recently, Duenser et al. [2020] used a rod model to robotically control a hot-wire foam cutter with a flexible elastic cutting wire. That the manifold of all equilibrium configurations attained by such a wire can be parametrized by a single coordinate chart was shown by Bretl and McCarthy [2014].
Rods and ribbons also play a significant role in self-formation processes, in which mechanical properties are controlled and exploited to realize a design. Examples are spatially varying thermal properties to induce controlled curvature in rods with a straight initial state [Wang et al. 2019], ribbons with spatially varying width to approximate a target shape given double-clamped boundary conditions [Liu et al. 2020] and ribbons whose natural curvature is modified to control their buckled shape [Xu et al. 2019a].
Properties of Design Spaces. Some works go beyond solving inverse problems numerically and also provide formal conditions for the feasibility of design problems. Pioneering this type of contribution in graphics literature, Derouet-Jourdan et al. [2010] study the problem of determining the natural curvature of an isotropic clamped-free elastic rod in the plane whose deformed shape under gravity is known. In addition to providing a numerical inversion algorithm, the authors give a sufficient stability condition on the stiffness and density of the rod. Bertails-Descoubes et al. [2018] study the same inversion problem for Kirchhoff rods in three dimensions and show that the natural configuration is uniquely determined by the equilibrium curve up to framing. Closest to our work, Hafner and Bickel [2021] give a geometric characterization of all plane curves that appear as equilibria of clamped-clamped elastic rods with spatially varying stiffness. We provide a direct generalization of this result to the setting of Kirchhoff rods in three dimensions with spatially varying anisotropy.
Beyond the design of elastic rods, Konaković Luković et al. [2018] study pneumatically actuated deployment of auxetic structures and show that this design space is described exactly by surfaces with positive mean curvature almost everywhere. Finally, Wang and Solomon [2021] represent skinning weights on meshes as solutions to a parametrized family of elliptic partial differential equations, which are shown to possess the same properties as high-quality weights. We use a similar characterization, based on solutions to a parametrized differential equation, to describe twist-free equilibria of elastic rods.
3 OVERVIEW
The technical sections of this article are organized in three parts. Our contributions are presented in Sections 5 and 7–9, while Sections 4 and 6 serve as technical introductions.
1. Kirchhoff Rods. We summarize the classical Kirchhoff rod model, with an emphasis on how to compute the stiffness matrix for a given cross section, and introduce the equilibrium equation for clamped-clamped boundary conditions (Section 4).
These tools are used to characterize the set of all stiffness matrices achievable within Kirchhoff rod theory and to show that this set is convex. Furthermore, we prove that a much smaller set of stiffness matrices, namely that induced by elliptical cross sections, suffices to reach every possible equilibrium state (Section 5).
2. Equilibrium Curves. We give a summary of Plücker coordinates and linear line complexes (Section 6), our main tools to study the design space of Kirchhoff rods. Then we present a geometric characterization of twist-free equilibrium states, which directly generalizes a previous result in the plane setting [Hafner and Bickel 2021]. We explore interactive design of these curves and discuss how to heuristically avoid stability issues (Section 7).
Next, we drop the twist-free constraint and present a geometric characterization of all Frénet curves that appear as deformed centerlines of Kirchhoff rods at equilibrium. This leads to a computational design algorithm that decides for a given (framed) curve whether it has the equilibrium property and, in case it does, computes the geometry of the corresponding Kirchhoff rod (Section 8).
3. Design under Load. We discuss a generic load model that supports point loads and line loads, such as the dead load of a rod, and the resulting modification to the equilibrium equation. The problem of finding solutions to this equation is then posed as a fixed-point problem, which can be solved via iteration (Section 9).
Finally, in Section 10, we show several physical examples that have been designed with our approach and manufactured using silicone casting, before we draw conclusions in Section 11.
All algorithms presented in this article rely only on linear programs and initial value problems, both of which can be solved numerically within a few milliseconds. The result is a set of techniques fast enough to produce results within a fraction of a second, so they can be used to give immediate feedback during an interactive editing session, or as part of more complex, iterative algorithms.
3.1 Main Result
Much of this article builds toward an algorithm that lets us interpret a Frénet curve \(\gamma :(0,\ell)\rightarrow \mathbb {R}^3\) as the static equilibrium of a Kirchhoff rod. We show in Proposition 8 that this is possible for any given curve if and only if there exist \(c,\bar{c}\in \mathbb {R}^3\) such that \(\begin{equation*} \langle \gamma ^{\prime } \times \gamma ^{\prime \prime }, c\times \gamma +\bar{c} \rangle \gt 0 \end{equation*}\) holds along the entire curve. This inequality is linear in the unknowns \(c\) and \(\bar{c}\), so it can be checked by a linear program. Furthermore, we will see that choosing concrete values for \(c,\bar{c}\) allows us to construct the three-dimensional geometry of an elastic rod that is guaranteed to deform in such a way that it will match \(\gamma\) at static equilibrium under suitable boundary conditions.
We can attach a geometric meaning to this inequality: The vector \(\gamma ^{\prime }\times \gamma ^{\prime \prime }\) gives the direction of the binormal line at a point of the curve, which is the normal line of the osculating plane. Vector fields of the form \(x\mapsto c\times x+\bar{c}\), for some fixed \(c\) and \(\bar{c}\), are so-called helical vector fields and arise as the velocity fields of a simultaneous translation along and rotation around a fixed axis.
The inequality states that \(\gamma\) can be converted to the equilibrium state of a Kirchhoff rod if and only if one can find a helical vector field that is not orthogonal to the binormal line in any point of the curve.1 Figure 2 illustrates an example in which the trajectories of a helical vector field (green) and the binormal lines of the curve (yellow) are very well aligned (center-left). This leads to a rod geometry that can reproduce the equilibrium state with a low amount of twist and that has a cross-sectional distribution with a high level of uniformity, which simplifies fabrication (center-right, right). The mechanical arguments for this will be made precise in Section 8.
4 KIRCHHOFF RODS: TECHNICAL PRELIMINARIES
The goal of this section is to introduce the mathematical tools to describe a popular model for the deformation of thin elastic rods, developed by Kirchhoff and Clebsch. Readers familiar with this subject may skip ahead to Equation (6) and consult the supplemental cheat sheet for an overview of our notation.
A rod is a slender three-dimensional body whose extent in one direction is much greater than that in the orthogonal directions, such as beams, cables, yarn, and hair. It can be shown that deformations of a rod under moderate loads admit a number of kinematic simplifications, such as near-inextensibility of the centerline and cross sections remaining nearly planar in bending. These assumptions lead to a reduction of the full three-dimensional displacement field to a one-dimensional description via curvatures and twist, which is captured mathematically by the concept of framed curves.
Fig. 3. Kirchhoff rod. Left: Initial, undeformed state. Right: Deformed state with bending and twist. The cross sections (gray) are anisotropic and vary across the length. Their centroids form the centerline \(\gamma\) (black), to which the material normals \(n_1\) and \(n_2\) (black) are attached. The purple-green stripe pattern on the surface visualizes the twist in the deformed state.
4.1 Framed Curves
We will focus our investigation on rods with a center line that is straight in the initial state. To keep track of bending and twisting modes, we rigidly attach a standard orthonormal frame \((e_1,e_2,e_3)\) to every point of the center line, assuming that the direction of the rod coincides with \(e_3\), as shown in Figure 3 (left).
Once forces are applied to the rod, the center line deforms isometrically into an arc-length parametrized curve \(\gamma :(0,\ell)\rightarrow \mathbb {R}^3\), with \(\ell\) the length of the rod. Likewise, the frames rotate such that \(e_1\) and \(e_2\) map onto the material normals \(n_1\) and \(n_2\), and \(e_3\) onto the curve tangent \(\gamma ^{\prime }\). These three orthonormal vectors form the columns of the moving frame \(F=(n_1,n_2,\gamma ^{\prime }):(0,\ell)\rightarrow SO(3)\). The pair \((\gamma ,F)\) is called a framed curve. Figure 3 shows an anisotropic Kirchhoff rod in its undeformed and deformed state.
Darboux Vector and Curvature. The derivative of \(F\) wrt the arc-length parameter \(s\) is described geometrically by the Darboux vector \(\omega :(0,\ell)\rightarrow \mathbb {R}^3\), which satisfies \(F^{\prime } = [\omega ]_\times F\). Here \([\omega ]_\times\) denotes the skew-symmetric matrix such that \([\omega ]_\times v = \omega \times v\) for all \(v\in \mathbb {R}^3\). The coordinates of \(\omega\) with respect to \(F\) are called the curvature vector \(k:(0,\ell)\rightarrow \mathbb {R}^3\), so \(\omega = Fk\). The components \(k=(\kappa _1,\kappa _2,\tau)^T\) are known as the material curvatures \(\kappa _i\), for \(i=1,2\), and the twist \(\tau\).
The normal component of \(\omega\), denoted by \(\omega _n\), is the same for all frames adapted to \(\gamma\), and can be computed as \(\begin{equation*} \omega - \langle \omega , \gamma ^{\prime } \rangle \gamma ^{\prime } = \omega _n = \gamma ^{\prime } \times \gamma ^{\prime \prime }. \end{equation*}\) Likewise, the (geometric) curvature \(\kappa =\sqrt {\kappa _1^2 +\kappa _2^2}\ge 0\) is determined by \(\gamma\) alone and can be computed as \(\kappa = \Vert \gamma ^{\prime \prime }\Vert = \Vert \omega _n\Vert\). Analogously to \(\omega _n\), we also introduce notation for the first two columns of \(F\) and the first two entries of \(k\): \(\begin{align*} F_n:=(n_1,n_2)=FS, \quad k_n:=(\kappa _1,\kappa _2)^T=S^T k, \quad \text{with} \quad S=\left(\begin{array}{ll}1 & 0 \\ 0 & 1\\ 0 & 0 \end{array}\right). \end{align*}\) This lets us write \(\omega _n = F_n k_n\) for any frame.
It will be useful to express \(k\) purely in terms of \(F\). To do this, we can use the following identity: For any \(Q\in \mathrm{SO}(3)\) and \(v\in \mathbb {R}^3\), it holds that \(Q[v]_\times Q^T = [Qv]_\times\).2 Then we compute \(F^{\prime }F^T=[\omega ]_\times =[Fk]_\times =F[k]_\times F^T\), which implies that \([k]_\times =F^T F^{\prime }\).
Special Frames. A frame with \(\tau \equiv 0\) is called parallel and describes a rod deformation without twist. For any given curve, a parallel frame \(F\) is uniquely determined by the image under \(F\) of a single point, for example \(F(0)\). Two parallel frames differ only by a constant rotation in the normal plane. The Darboux vector of a parallel frame is contained in the normal plane, so \(\omega =\omega _n\).
Away from inflection points, a curve has a unique Serret–Frénet frame defined by the principal normal \(n_1=\gamma ^{\prime \prime }/\Vert \gamma ^{\prime \prime }\Vert\) and characterized by \(\kappa _1\equiv 0\). The curve \({binormal}\) \(n_2=(\gamma ^{\prime }\times \gamma ^{\prime \prime })/\Vert \gamma ^{\prime \prime }\Vert\) is parallel to \(\omega _n\). A curve that has a Serret–Frénet frame everywhere is called a Frénet curve and is characterized by \(\kappa \gt 0\).
Relating Adapted Frames. Two frames \(F\) and \(F_\beta\) adapted to the same curve \(\gamma\) share the same third basis vector at every point, i.e., \(Fe_3=\gamma ^{\prime }=F_\beta e_3\). Thus, they are related through a rotation around the third basis vector by an angle \(\beta :(0,\ell)\rightarrow \mathbb {R}\) that may vary as a function of \(s\). Likewise, we can relate their curvature vectors:
(1) \(\begin{equation} \begin{gathered}F_{\beta ,n} = F_n Q_\beta , \quad \text{with} \quad Q_\beta =\left(\begin{array}{ll} \cos \beta & -\sin \beta \\ \sin \beta & \cos \beta \end{array}\right), \\ k_{\beta ,n} = Q_\beta ^T k_n, \quad \tau _\beta = \tau + \beta ^{\prime }. \end{gathered} \end{equation}\) Figure 4 shows a curve with a parallel frame and its Serret–Frénet frame, as well the rotation between the two frames.
Fig. 4. Rotation between a parallel and the Serret–Frénet frame. Any two frames adapted to the same curve are related through a rotation with some angle \(\beta :(0,\ell)\rightarrow \mathbb {R}\) around the curve tangent. In this example, \(\beta\) relates a parallel frame (purple) and the Frénet frame (green) of the curve.
4.2 Elastic Energy
The elastic energy of a deformed Kirchhoff rod is defined based on the assumption that bending and twisting modes of deformation can be decoupled and contribute separately to the energy. The contribution of each mode is derived by assuming a state of uniform bending or twisting and computing the respective energy from three-dimensional elasticity. We will only discuss aspects of the resulting formulation that are relevant to subsequent sections—for a detailed derivation, see Audoly and Pomeau (2010)3.3–3.5].
The elastic energy of a Kirchhoff rod is defined as (2) \(\begin{equation} \!\!W=\frac{1}{2}\!\int _0^\ell \!\! \langle k, K k\rangle , \text{ with } K= \begin{array}{ll} EI & 0_{2\times 1}\\ 0_{1\times 2} & \mu J \end{array}, \, I=\begin{array}{ll} I_{xx} & I_{xy} \\ I_{xy} & I_{yy} \end{array}. \end{equation}\) The Young’s modulus \(E\gt 0\) and shear modulus \(\mu \gt 0\) only depend on the base material of the rod and are assumed to be constant across the length. The stiffness matrix \(K(s)\) contains the bending stiffness \(EI(s)\) and twisting stiffness \(\mu J(s)\), with \(I(s)\) and \(J(s)\) dependent on the cross section of the rod at \(s\in (0,\ell)\).
In particular, the bending rigidity \(I\) is the area moment of inertia tensor, whose coordinates are given by (3) \(\begin{equation} I(s)=\int _{\mathcal {D}(s)} \begin{array}{ll} y^2 & -xy\\ -xy & x^2 \end{array}\, dA(x,y). \end{equation}\) Here \(\mathcal {D}(s)\subset \mathbb {R}^2\) is the cross section of the rod at \(s\), with the centerline passing through its centroid. The integrand can be written as the outer product of \((-y,x)^T\) with itself, so \(I(s)\in S^2_{++}\), i.e., \(I(s)\) is symmetric positive-definite whenever \(\mathcal {D}(s)\) has positive area measure. The bending energy is given by \(\frac{1}{2}E\langle k_n, I k_n\rangle\) and is due to normal stresses away from the centerline, as seen in Figure 5 (center).
If \(\mathcal {D}(s)\) is simply connected, then the torsional rigidity \(J\) is given by the Dirichlet energy of the solution to a Poisson equation,
(4) \(\begin{equation} J(s)=4\int _{\mathcal {D}(s)} \Vert \nabla \chi \Vert ^2, \quad \text{with} \quad \begin{array}{ll} \Delta \chi = -1 & \text{in } \; \mathcal {D}(s), \\ \chi = 0 & \text{on } \; \partial \mathcal {D}(s). \end{array} \end{equation}\) In general, neither \(I\) nor \(J\) have closed-form solutions, but they do for special cases such as circular and elliptical disks. The twisting energy is given by \(\frac{1}{2} \mu J \tau ^2\) and is due to shear stresses away from the centerline. Twisting causes cross sections to deform out of plane to reach the minimum-energy state, as shown in Figure 5 (right).
Fig. 5. Deformation modes. Left: Undeformed rod. Center: Uniform bending. Right: Uniform twist, with cross sections warping out of plane.
4.3 Equilibrium Equations
We assume kinematic, or double-clamped, boundary conditions: Both \(\gamma\) and \(F\) are fixed at \(s=0\) and \(s=\ell\). These boundary conditions are often encountered in architectural and interior design applications, and they make for the richest design space of equilibrium states.
To set up the variational problem, we choose \(F\) as the primary variable, so fixing \(F\) at both ends imposes Dirichlet boundary conditions. Assuming that \(\gamma (0)\) coincides with the origin, we can express \(\gamma\) as a function of \(F\) via \(\gamma (s)=\int _0^s \gamma ^{\prime }=\int _0^s Fe_3\). Thus, the endpoint constraint \(\gamma (\ell)=\gamma _\ell\) takes the form of an integral constraint \(\int _0^\ell Fe_3 = \gamma _\ell\). Constrained extremals of the Kirchhoff energy are characterized by extremals of the Lagrangian, (5) \(\begin{equation} \mathcal {L}= \int _0^\ell \left(\frac{1}{2}\langle k, K k\rangle - \langle c, Fe_3\rangle \right), \end{equation}\) with Lagrange multiplier \(c\in \mathbb {R}^3\).
Section 1 of the supplemental document shows how to derive the constrained Euler–Lagrange equation of this variational problem by investigating variations of \(F\) that respect its orthogonal structure. This results in the following statement: A framed curve \((\gamma ,F)\) is a static equilibrium state of a Kirchhoff rod with stiffness \(K\) and kinematic boundary conditions if and only if there exist \(c,\bar{c}\in \mathbb {R}^3\) such that (6) \(\begin{equation} FKk=c\times \gamma + \bar{c}, \end{equation}\) called the moment equilibrium equation.
Our main contribution is to give geometric characterizations of (framed) curves having this property at three levels of generality, which are captured by the following definitions:
Let \((\gamma ,F)\) be an arc-length parametrized framed curve of length \(\ell\) with curvature vector \(k\). Assume there exist \(c,\bar{c}\in \mathbb {R}^3\) and domains \(\mathcal {D}(s)\subset \mathbb {R}^2\) centered at the origin for all \(s\in (0,\ell)\) such that \(FKk=c\times \gamma +\bar{c}\) holds, where \(K\) denotes the stiffness matrix induced by \(\mathcal {D}\). Then \((\gamma ,F)\) is called a framed equilibrium curve.
Let \(\gamma\) be an arc-length parametrized curve. If there exists a frame \(F:(0,\ell)\rightarrow \text{SO}(3)\) adapted to \(\gamma\) such that \((\gamma ,F)\) is a framed equilibrium curve, then \(\gamma\) is called an equilibrium curve. If, additionally, \(F\) can be chosen to be a parallel frame, then \(\gamma\) is called a parallel equilibrium curve.
Our geometric characterizations of parallel equilibrium curves, equilibrium curves, and framed equilibrium curves are given in Sections 7.1, 8.1, and 8.5, respectively.
5 KIRCHHOFF RODS: THE CONSTITUTIVE RELATION
To uniquely define the geometry of a straight elastic rod, we have to choose a cross section \(\mathcal {D}(s)\subset \mathbb {R}^2\) at every \(s\in (0,\ell)\), such that the center line passes through the centroid of \(\mathcal {D}(s)\). This choice determines \(K(s)\) and thus the mechanical behavior of the rod at this point. Ultimately, we want to study the set of equilibrium states, but to do this, we first have to characterize the set of all stiffness matrices \(K\) that are induced by admissible cross sections.
We carry out this characterization in Section 5.1 and then infer two results that are relevant in the context of design. First, we show in Section 5.2 that the set of admissible stiffness matrices is convex, thus allowing design problems with convex objectives to be cast as convex optimization problems. Second, we prove in Section 5.3 that the subset of admissible stiffness matrices corresponding to elliptical cross sections is sufficient to span the entire design space of framed equilibrium curves. Thus, we do not lose any design freedom upon neglecting exotic cross sections that would be hard to fabricate or fail to satisfy the Kirchhoff assumptions.
5.1 Admissible Stiffness Matrices
We have already seen in Equations (2)–(4) that every stiffness matrix \(K\) consists of a 2-by-2 block \(EI\in S^2_{++}\), the bending stiffness, and \(\mu J\gt 0\), the torsional stiffness. But is every matrix that obeys these constraints induced by an admissible cross section?
Regarding the bending stiffness, it is easy to see that every \(I\in S^2_{++}\) can be attained, for example using elliptical cross sections. In particular, we can choose the radii \(a,b\) of an ellipse to determine the eigenvalues of \(I\) and the orientation \(\varphi\) to determine the eigenvectors, as shown in the inset. The resulting area moment matrix is (7) \(\begin{equation} I=\frac{\pi }{4}Q\begin{array}{ll} ab^3 & 0 \\ 0 & a^3 b \end{array}Q^T, \quad \text{with} \quad Q=\begin{array}{ll} \cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{array}. \end{equation}\)
For the torsional stiffness, we see from Equation (4) that any \(J\gt 0\) can be attained in principle, for example, by circular disks of different radii. However, if we restrict our attention to cross sections with fixed \(I\), then the range of attainable \(J\) is limited by (8) \(\begin{equation} J\le 4\psi (I), \quad \text{with} \quad \psi :S^2_{++}\rightarrow \mathbb {R}:X\mapsto \frac{\det X}{\operatorname{tr}X}, \end{equation}\) as shown by Diaz and Weinstein (1948)p. 5]. This upper bound on \(J\) is tight, and the unique maximizer for any given \(I\) is an ellipse. We can see from Equation (7) that in this case, \(J=\frac{\pi a^3b^3}{a^2+b^2}\).
However, we prove in Section 2 of the supplemental document that there is no positive lower bound:
Given \(I\in S^2_{++}\) and \(J\gt 0\), there is a bounded domain \(\mathcal {D}\subset \mathbb {R}^2\) with bending rigidity \(I\) and torsional rigidity \(J\) if and only if \((I,J)\) is an element of (9) \(\begin{equation} \mathcal {K} = \lbrace (I,J) \in S^2_{++}\times \mathbb {R} : 0 \lt J \le 4 \psi (I) \rbrace . \end{equation}\)
A heuristic argument that we can decrease \(J\) without changing \(I\) is illustrated in Figure 6: We start with a cross section realizing \(I\) and progressively add cuts from the boundary to the interior of the domain. This does not change \(I\), because we only remove zero-measure sets from \(\mathcal {D}\). However, \(J\) can be made arbitrarily small as the boundary points of the domain move closer and closer together.
Fig. 6. \(J\) has no positive lower bound. The green gradient shows the solution \(\chi\) to Equation (4) for circular disks with an increasing number and length of cuts. Progressively adding cuts makes \(J\) arbitrarily small while maintaining \(I\) . Greatest incircles are drawn in yellow.
We assume \(E\) and \(\mu\) to be fixed, so a pair \((I,J)\in \mathcal {K}\) uniquely defines a stiffness matrix \(K=\mathrm{diag}(EI,\mu J)\) and vice versa. For convenience, we will sometimes write \(K=(I,J)\) and \(K\in \mathcal {K}\) by abuse of notation.
5.2 Convexity of 𝒦
Designing a Kirchhoff rod with a prescribed equilibrium state \((\gamma ,F)\) amounts to finding \(c,\bar{c}\in \mathbb {R}^3\) and \(K:(0,\ell)\rightarrow \mathcal {K}\) that solve Equation (6). Because Equation (6) is linear in \(c,\bar{c}\), and \(K\), we need only convexity of \(\mathcal {K}\) to show that this is a convex problem for arbitrary convex objectives. Our proof is based on:
The function \(\psi :X\mapsto \frac{\det X}{\operatorname{tr}X}\) is concave on \(S^2_{++}\).
This can be shown either with any computer-algebra system, or by hand, as we do in Section 3 of the supplemental document. The lemma implies the following:
Let \((I_0,J_0),(I_1,J_1)\in \mathcal {K}\), and \(t\in (0,1)\). Because \(I_0,I_1\in S^2_{++}\), and \(S^2_{++}\) is convex, we have \((1-t)I_0+tI_1\in S^2_{++}\). From concavity of \(\psi\), it follows that \(\begin{equation*} 0 \lt (1-t)J_0+tJ_1 \le 4((1-t)\psi (I_0) + t\psi (I_1)) \le 4\psi ((1-t)I_0+tI_1), \end{equation*}\) which shows that \(((1-t)I_0+tI_1, (1-t)J_0+tJ_1)\in \mathcal {K}\).□
This result guarantees that we can numerically find the global optimizer of optimization problems over \(\mathcal {K}\) constrained by Equation (6), as long as the objective function is convex. In particular, we can solve the constraint satisfaction problem (CSP) of determining whether a given framed curve is a framed equilibrium curve. In Section 8, we improve this result and show that the CSP can even be solved by a linear program.
Fig. 7. Elliptification. Left: Rod with cross-shaped cross sections (light blue) and rod with elliptical cross sections (yellow) having the same equilibrium state. Right: Comparison between cross-shaped and elliptical cross sections at marked locations, with \(k_n\) visualized as arrow.
5.3 Elliptical Cross Sections
While \(\mathcal {K}\) is convenient for numerical optimization, it contains stiffness matrices that are only induced by cross sections like the one in Figure 6 (right), which are impractical to manufacture. Moreover, cross sections like this will buckle even under moderate loads and thus break the assumptions of Kirchhoff rods. Therefore, it would be best to avoid using cross sections with \(J \ll 4\psi (I)\) in design.
This raises two questions: How is the design space of Kirchhoff rods reduced if we restrict ourselves to a proper subset of \(\mathcal {K}\) that has impractical cross sections removed? And what is a good subset to use? Surprisingly, a very restrictive choice, namely, (10) \(\begin{equation} \mathcal {K}^\ast := \partial \mathcal {K}\cap \mathcal {K} = \lbrace (I,J)\in S^2_{++} \times \mathbb {R}: J=4\psi (I)\rbrace , \end{equation}\) is still optimal. This set contains exactly the stiffness matrices induced by elliptical cross sections and realizes all framed equilibrium curves, with no reduction of the design space, as we show here:
Let \((\gamma ,F)\) be a framed equilibrium curve such that \(FKk=c\times \gamma +\bar{c}\) holds with \(c,\bar{c}\in \mathbb {R}^3\) and \(K:(0,\ell)\rightarrow \mathcal {K}\). Then there exists \(K^*:(0,\ell)\rightarrow \mathcal {K}^\ast\) such that \(FK^*k=c\times \gamma +\bar{c}\).
We rewrite the equilibrium equation by splitting it up into its normal and tangential components, which gives (11) \(\begin{equation} EIk_n=F_n^T(c\times \gamma +\bar{c}), \quad \mu J\tau = \langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle . \end{equation}\) To show the statement, we compute \(I^*\in S^2_{++}\) for each \(s\in (0,\ell)\) such that \((I^*,J)\in \mathcal {K}^*\) and \(I^*k_n=Ik_n\).
Let \(Q\in SO(2)\) such that \(Qk_n = \kappa e_1\), and define \(\tilde{I}=QIQ^T\). Let \(\begin{equation*} \tilde{I}^* = \begin{array}{ll}\tilde{I}_{xx} & \tilde{I}_{xy} \\ \tilde{I}_{xy} & \tilde{I}_{yy}^* \end{array}, \quad \text{with} \quad \tilde{I}_{yy}^*=\frac{J\tilde{I}_{xx}+4\tilde{I}_{xy}^2}{4\tilde{I}_{xx} - J}, \end{equation*}\) which only has its bottom-right entry different from \(\tilde{I}\). Then \(I^*=Q^T\tilde{I}^*Q\) satisfies all requirements, as we will now show.
The equilibrium equation is satisfied, because \(\begin{equation*} I^* k_n = Q^T \tilde{I}^* Qk_n = \kappa Q^T \tilde{I}^* e_1 = \kappa Q^T\tilde{I} e_1 = Q^T\tilde{I}Qk_n = Ik_n. \end{equation*}\) To verify that \(I^*\in S^2_{++}\), it suffices to show that \(\tilde{I}^*\in S^2_{++}\). We have \(\tilde{I}_{xx}\gt 0\), because \(I\) and thus \(\tilde{I}\) is in \(S^2_{++}\). Furthermore, \(\det \tilde{I}^*=\frac{J(\tilde{I}_{xx}^2+\tilde{I}_{xy}^2)}{4\tilde{I}_{xx} - J} \gt 0\) if \(J \lt 4\tilde{I}_{xx}\). But this always holds, because \(\begin{equation*} J \le 4\psi (I)=4\psi (\tilde{I}) = \frac{4(\tilde{I}_{xx}\tilde{I}_{yy}-\tilde{I}_{xy}^2)}{\tilde{I}_{xx}+\tilde{I}_{yy}} \le 4\tilde{I}_{xx} \frac{\tilde{I}_{yy}}{\tilde{I}_{xx}+\tilde{I}_{yy}} \lt 4\tilde{I}_{xx}, \end{equation*}\) where we have used that \(\psi\) only depends on the tensor invariants of its argument. Finally, \(J=4\psi (I^*) = 4\psi (\tilde{I}^*)\), and thus \((I^*,J)\in \mathcal {K}^*\) is shown by direct computation, from the definition of \(\tilde{I}^*\).□
This result shows that any Kirchhoff rod equilibrium state is also attainable by a Kirchhoff rod consisting exclusively of elliptical cross sections. Figure 7 shows an example in which a rod with spatially varying cross-shaped cross sections is converted to a rod with elliptical cross sections. Even though the resulting elliptical cross sections do not have the same bending rigidity, it is guaranteed that the equilibrium state is maintained.
Using \(\mathcal {K}\) in optimization and in proofs, and \(\mathcal {K}^*\) for design and fabrication gives us the best of both worlds: The convexity of \(\mathcal {K}\) gives optimality guarantees in optimization and simplifies mathematical analysis. Meanwhile, any solution obtained using \(\mathcal {K}\) can be converted to a solution in \(\mathcal {K}^*\), which avoids cross-sectional buckling issues, and is guaranteed to yield moldable geometries, because elliptical cross sections are convex.
6 EQUILIBRIUM CURVES: TECHNICAL PRELIMINARIES
Our next goal is to give a thorough analysis of the types of equilibrium curves introduced in Definitions 1 and 2. In particular, we provide characterizations that relate the shape of these curves to existing concepts from projective line geometry and that pave the way to very efficient computational design algorithms. We will show that some fundamental inverse design problems related to Kirchhoff rods can be solved exactly, in the sense that we can decide computationally whether there is a solution for a given input, and find a solution near-instantaneously if it exists.
Looking back at the previous section, we found that choosing domains \(\mathcal {D}(s)\subset \mathbb {R}^2\) in Definition 1 is equivalent to choosing a stiffness function \(K:(0,\ell)\rightarrow \mathcal {K}\). Proposition 4 further guarantees that we can restrict \(\mathcal {K}\) to \(\mathcal {K}^*\), thus using only elliptical cross sections, without changing the design space of Kirchhoff rod equilibrium states. To characterize this space, we ask the following: For which (framed) curves is it possible to choose \(c,\bar{c}\in \mathbb {R}^3\) and a stiffness function so Equation (6) is satisfied? Our answers, which we give in Sections 7 and 8, rely on concepts from the geometry of lines in \(\mathbb {R}^3\), and in particular on linear line complexes, which we introduce below. For a more extensive treatment, we point the reader to Pottmann and Wallner (2001)Ch. 2.1, 3.1]PLXBIB0033.
6.1 Line Geometry
Let \(L\) be a line in \(\mathbb {R}^3\), with direction \(v\in \mathbb {R}^3\) and passing through a point \(x\in \mathbb {R}^3\). We call \((l,\bar{l}):=(v,x\times v)\in \mathbb {R}^6\) the Plücker coordinates of \(L\) and write \(L=\lambda (l,\bar{l})\) to denote the line defined by \((l,\bar{l})\) as a subset of \(\mathbb {R}^3\). Plücker coordinates satisfy \(\langle l,\bar{l}\rangle = 0\) and are homogeneous, i.e., every non-zero scalar multiple refers to the same line. The projective extension of this set is known as the Klein quadric, \(\begin{equation*} \Lambda _\text{kl}:=\lbrace (l,\bar{l})\in \mathbb {P}^5 : \langle l,\bar{l}\rangle = 0\rbrace . \end{equation*}\) Plücker coordinates of the form \((0,\bar{l})\) represent an ideal line, which contains exactly the ideal points (“points at infinity”) with directions orthogonal to \(\bar{l}\). One useful application of Plücker coordinates is the computation of intersections: two lines \(\lambda (l,\bar{l})\) and \(\lambda (c,\bar{c})\) intersect if and only if \(\langle l,\bar{c}\rangle +\langle \bar{l}, c \rangle =0\).
The set of all lines whose Plücker coordinates satisfy a homogeneous linear equation with (fixed) coefficients \(c,\bar{c}\in \mathbb {R}^3\), \(\begin{equation*} \mathcal {C}:=\lbrace (l,\bar{l})\in \Lambda _\text{kl} : \langle l,\bar{c}\rangle +\langle \bar{l}, c \rangle =0 \rbrace , \end{equation*}\) is known as a linear line complex. If \((c,\bar{c})\in \Lambda _\text{kl}\), i.e., if \((c,\bar{c})\) are themselves Plücker coordinates of a line, then the complex is said to be singular and consists of all lines intersecting \(\lambda (c,\bar{c})\), as shown in Figure 8 (left). In case \((c,\bar{c})\) is ideal, so \(c=0\), the complex \(\mathcal {C}\) consists of all lines with directions orthogonal to \(\bar{c}\), see Figure 8 (right).
We can also interpret a linear complex geometrically if \((c,\bar{c})\notin \Lambda _\text{kl}\), in which case \(\mathcal {C}\) is called regular. To do this, consider the vector field \(h:x\mapsto c\times x+\bar{c}\) in \(\mathbb {R}^3\). A vector field of this form is called helical, because every field line, i.e., every curve tangent to \(h\) in every point, is a helix. The family of all field lines of \(h\) is called a helical motion and consists of all helices with a fixed axis and a fixed pitch, see Figure 9 (left). Next, we consider at every \(x\in \mathbb {R}^3\) the pencil of lines through \(x\) that are normal to the helix passing through this point, see Figure 9 (right). The union of all such line pencils gives exactly the lines contained in \(\mathcal {C}\). In summary: A regular linear complex is the set of path normals of a helical motion.
Fig. 8. Singular linear complexes. Left: A Euclidean singular linear complex ( \(c \ne 0\) ) contains all lines (black) incident to a fixed Euclidean line (purple) with Plücker coordinates \((c,\bar{c})\) . Right: An ideal singular linear complex ( \(c=0\) ) contains all lines (black) orthogonal to a fixed direction \(\bar{c}\) (purple).
7 PARALLEL EQUILIBRIUM CURVES
Our first application of linear line complexes is a geometric characterization of parallel equilibrium curves—curves that appear as twist-free equilibria of Kirchhoff rods. This is a natural starting point, because the result directly generalizes the characterization of plane elastic curves, which we summarize here for convenience:
Let \(\gamma :(0,\ell)\rightarrow \mathbb {R}^2\) with signed curvature \(\kappa\) and a finite number of inflection points, and let \(a\in \mathbb {R}^2\), \(b\in \mathbb {R}\) not all zero. Then there exists \(K:(0,\ell)\rightarrow \mathbb {R}\) with \(0\lt \inf K\le \sup K\lt \infty\) such that \(K\kappa =\langle a,\gamma \rangle + b\) if and only if
In Proposition 6, this theorem corresponds exactly to the special case in which the constants from the right-hand side of the equilibrium equation \(FKk=c\times \gamma +\bar{c}\) define a singular complex \((c,\bar{c})\). Then the line \(\lambda (c,\bar{c})\) plays the role of \(L\) in Theorem 5. In contrast, regular complexes \((c,\bar{c})\) correspond exactly to parallel equilibrium curves that are not plane.
Fig. 9. Regular linear complexes. Left: A helical motion is a family of helices (green) with fixed axis (purple) and fixed pitch. Right: For every helix (green) in the helical motion and its axis (purple), take the pencil of normal lines (black) at every point. The union of all such line pencils forms a regular linear complex.
Before stating the characterization theorem, we outline the remainder of this section: The class of parallel equilibrium curves proves to be rigid in the sense that only one curvature can be chosen independently at every point, unlike general curves in \(\mathbb {R}^3\), which have two independent curvatures. To facilitate the design of curves that lie in this class by construction, we derive a second characterization based on a family of ordinary differential equations and parametrized directly by curvature. In this context, we also show that all parallel equilibrium curves are essentially warped superhelices, which gives an intuitive understanding of the shapes achievable within this class. Finally, we discuss how to avoid unstable solutions by judiciously choosing cross sections.
7.1 Geometric Characterization
The central objects in the characterization are the tangent lines and binormal lines of a curve, i.e., the lines passing through \(\gamma (s)\) in directions \(\gamma ^{\prime }(s)\) and \(\omega _n(s)=\gamma ^{\prime }(s)\times \gamma ^{\prime \prime }(s)\), respectively, as shown in Figure 10. Binormal lines are not defined at inflection points, and this is reflected in the fact that non-plane parallel equilibrium curves are always Frénet curves, as we will show. In this sense, planar solutions are exceptional: They are the only parallel equilibrium curves that may contain inflections.
Let \(\gamma\) be an arc-length parametrized curve in \(\mathbb {R}^3\) with a finite number of inflection points. Then \(\gamma\) is a parallel equilibrium curve if and only if one of the following holds:
(1) | The curve \(\gamma\) is not plane, and there is a linear complex \(\mathcal {C}\) that contains all tangent lines of \(\gamma\) but none of its binormal lines. | ||||
(2) | The curve \(\gamma\) is plane and satisfies the assumptions of Theorem 5. | ||||
Solutions of type (1) and (2) correspond to regular and singular values of \((c,\bar{c})\) in Equation (6), respectively. Solutions of type (1) are Frénet curves.
The proof of this statement is given in Appendix A and shows that the geometric conditions on tangent and binormal lines are equivalent to (12a) \(\begin{align} \langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle &= 0, \end{align}\) (12b) \(\begin{align} \langle \gamma ^{\prime }\times \gamma ^{\prime \prime },c\times \gamma +\bar{c}\rangle &\gt 0, \end{align}\)
for some \(c,\bar{c}\in \mathbb {R}^3\). The parallel frame \(F\) does not appear in these conditions, so we need not explicitly compute it to check if they are satisfied. Furthermore, the conditions are linear in the unknowns \(c\) and \(\bar{c}\), so they can be checked with a linear program in principle.
However, Equation (12a) is a pointwise equality constraint on \(\gamma\), which means that the number of independent curvatures of \(\gamma\) is reduced from two to one. Non-plane curves satisfying this equation generally have no spline-based representation, because finitely many control points do no provide enough degrees of freedom to satisfy infinitely many constraints. This rules out most curve-drawing tools as a viable option for interactive design. Below, we propose a different way of exploring the design space of parallel equilibrium curves, but we will find another use for Equation (12a) in Section 8.2 for the design of equilibrium curves that need not be parallel.
Fig. 10. Tangent and binormal lines. Left: Curve (green) and its family of tangent lines (black). Right: Curve (green) and its family of osculating planes (green) and binormal lines (black).
7.2 ODE Characterization
In Section 4 of the supplemental document, we derive two more characterizations of non-plane parallel equilibrium curves:
Let \(\gamma\) be an arc-length parametrized curve in \(\mathbb {R}^3\). Then the following are equivalent:
(1) | There exist \(c,\bar{c}\in \mathbb {R}^3\) such that \(\gamma\) satisfies Equation (12). | ||||
(2) | There exist \(c,\bar{c}\in \mathbb {R}^3\) such that the ordered set \(\lbrace \gamma ^{\prime },\gamma ^{\prime \prime },c\times \gamma +\bar{c}\rbrace\) is a right-handed orthogonal basis at every point. | ||||
(3) | There exist \(c,\bar{c}\in \mathbb {R}^3\) and \(m:(0,\ell)\rightarrow \mathbb {R}_{\gt 0}\) such that (13) \(\begin{equation} \gamma ^{\prime \prime }(s)=m(s)\cdot (c\times \gamma (s)+\bar{c})\times \gamma ^{\prime }(s) \end{equation}\) and \(\langle \gamma ^{\prime }(0),c\times \gamma (0)+\bar{c}\rangle = 0\). It holds that \(\kappa =m\cdot \Vert c\times \gamma +\bar{c}\Vert\). | ||||
Equation (13) is well suited for the interactive exploration of parallel equilibrium curves, because the user can directly modify \(\kappa\), which results in a predictable change of the curve shape. A useful way of thinking about the space of feasible designs is that parallel equilibrium curves are warped superhelices, the shapes obtained by coiling one helix around another helix. Constant-curvature solutions of Equation (13) constitute a family of superhelices (with an exact discrete helical symmetry, as discussed in Section 5 of the supplemental document) and contain regular helices as a special case, as shown in Figure 11. Circles and lines are obtained in the limit as \(\langle c,\bar{c}\rangle \rightarrow 0\) and \(\kappa \rightarrow 0\). By modifying \(\kappa\) to be non-constant, the user can warp portions of the superhelix to explore the shape space.
Apart from \(\kappa\), we can also choose \(c\), \(\bar{c}\), and initial conditions \(\gamma (0)\) and \(\gamma ^{\prime }(0)\) satisfying \(\langle \gamma ^{\prime }(0),c\times \gamma (0)+\bar{c}\rangle = 0\). After eliminating rigid motions, these reduce to three scalars \(p,r,\alpha \in \mathbb {R}\) such that
(14) \(\begin{equation} c=e_3, \quad \bar{c}=pe_3, \quad \gamma (0)=re_1, \quad \gamma ^{\prime }(0) = Q_\alpha e_1, \end{equation}\) where \(Q_\alpha \in SO(3)\) is the rotation by angle \(\alpha\) around the axis in direction \(re_2+pe_3\). Modifying these scalars controls the ratio between the outer and inner radius of the superhelix, as well as the tightness of the windings.
Fig. 11. Constant-curvature parallel equilibrium curves. The one-parametric family of solutions to Equations (13) and (14) for \(r=\frac{1}{2}=p\) and \(\kappa \equiv 1\) , obtained by sweeping \(\alpha \in (0,2\pi ]\) . The family interpolates smoothly between superhelices and a special helical solution at \(\alpha =\pi /2\) .
7.3 Cross Sections and Stability
The last step in turning a parallel equilibrium curve into fabricable geometry is to choose cross sections satisfying the equilibrium equation. Because \(\tau \equiv 0\), the torsional rigidity \(J\) drops out of the equation, so we need only choose \(I\in S^2_{++}\) at every point to satisfy \(EIk_n = F_n^T(c\times \gamma +\bar{c})\). We can then convert \(I\) into a cross section \(\mathcal {D}\subset \mathbb {R}^2\), for example an elliptical one. Naturally, we can scale all cross sections uniformly along the rod to control the overall thickness. This is because a rescaled cross section \(t\mathcal {D}\) has bending rigidity \(t^4 I\), which solves the equilibrium equation with \((t^4 c, t^4 \bar{c})\).
To characterize all \(I\in S^2_{++}\) that solve \(EIk_n = F_n^T(c\times \gamma +\bar{c})\), note that \(k_n\) and \(F_n^T(c\times \gamma +\bar{c})\) are parallel: From Proposition 7(2), we have that \(\omega _n=\gamma ^{\prime }\times \gamma ^{\prime \prime }\) is parallel to \(c\times \gamma +\bar{c}\), and we know that \(k_n = F_n^T\omega _n\). The factor of proportionality is given by \(m\), as seen from \(\Vert k_n\Vert =\kappa =m\cdot \Vert c\times \gamma +\bar{c}\Vert\). Thus, we can parametrize all admissible \(I\) by the spectral decomposition \(\begin{gather*} I=\lambda _1 \, v_1\otimes v_1 + \lambda _2 \, v_2\otimes v_2, \quad \text{with}\\ \lambda _1 = 1/(Eq), \quad v_1 = k_n/\kappa , \quad v_2 = (-\kappa _2,\kappa _1)^T / \kappa , \end{gather*}\) and free parameter \(\lambda _2 \gt 0\). If we consider an elliptical cross section in the coordinate system spanned by \(n_1\) and \(n_2\), then the eigenvectors \(v_1\) and \(v_2\) give the semiaxes of the ellipse. The bending axis has direction \(k_n\) and coincides with \(v_1\).
Even though all choices for \(\lambda _2 \gt 0\) will give an equilibrium state, they differ greatly in their stability properties. This is because a rod will bend easily around its weak axis, but trying to bend it around its strong axis will usually result in loss of stability by buckling. It is thus essential to pick cross sections in such a way that \(v_1\) coincides with the major semiaxis of the ellipse, as shown in Figure 12. This is achieved by choosing \(a\gt b\) with the notation of Equation (7), which is equivalent to \(\lambda _2 \gt \lambda _1\). Violating this rule will almost surely result in an unstable equilibrium state.
We expose the ratio \(u=a/b\) to the user and restrict its domain to \(u\ge 1\). In computation, the ratio is achieved by setting \(\lambda _2 = u^2\lambda _1\). Even though this rule excludes a common source of buckling, it is a heuristic and not a formal stability guarantee. Indeed, for curves with many helical windings, a stable choice of \(I\) will likely not exist. If it is not intuitively clear whether a design is stable or not, then it may be necessary to check for stability issues numerically prior to fabrication, or to verify the formal stability conditions for the solution [Manning et al. 1998].
Fig. 12. Cross sections and stability. From left to right: undeformed rod; stable equilibrium attained by aligning the major axis of the cross section with the bending axis; unstable equilibrium attained by aligning the minor axis with the bending axis; and the resulting buckled solution.
8 GENERAL EQUILIBRIUM CURVES
To take advantage of the full design space offered by Kirchhoff rods, we need to bring twist into the equation. This goal requires a fundamental design decision: Does the burden to specify the twist fall on the designer, or is it done computationally? Both choices have a drawback: The designer may have a good idea of the curve shape they want but may lack the intuition to know what twist will work to realize it. However, relying on a computational solution makes it more difficult to ensure that the twist falls into physically plausible bounds, as we will see.
For most of this section, we explore the latter choice, in which the designer only prescribes the deformed center line of the rod and leaves the computation of twist and cross sections to the computer. This leads to the characterization of general equilibrium curves—curves in \(\mathbb {R}^3\) that can be framed to yield a solution to Equation (6), the equilibrium equation for Kirchhoff rods. Like in the parallel case, the characterization yields conditions that are linear in the unknowns and can thus be checked with a linear program. However, general equilibrium curves prove to be free of pointwise equality constraints, so the linear program is a useful tool for computational design and will work on manually drawn input curves, for example splines.
The challenge with this workflow is to find a solution that not only satisfies the equilibrium equation but is also practical for fabrication. Two relevant criteria are that cross sections should not vary too much in size across the length and that the twist should stay within reasonable bounds. We will show how these goals can be achieved without jeopardizing linearity.
At the end of this section, we also present a theoretical result about framed equilibrium curves, in which the curve and the moving frame are both prescribed. While we do not use this result for design applications, it provides a general insight into the structure of equilibrium states attainable within the theory of Kirchhoff rods.
8.1 Geometric Characterization
As a proper superset of parallel equilibrium curves, the conditions for general equilibrium curves must be strictly weaker than those in Proposition 6. Indeed, for Frénet curves, the characterization emerges simply by removing the tangent line condition, so only the binormal line condition remains.
We can view curves with continuous geometric curvature as elements of \(C^2((0,1);\mathbb {R}^3)\) by relaxing the requirement of an arc-length parametrization. With the standard \(C^2\)-norm, Frénet curves are dense within the set of all regular curves. This implies that inflection points of curves in \(\mathbb {R}^3\) are non-essential features that can be removed by a perturbation of arbitrarily small magnitude—this is in contrast to plane curves, where inflections are not removable.
We opt to prove the characterization only for Frénet curves, because the resulting theory suffices for the application explored in this work: to turn user-drawn spline curves into Kirchhoff rod equilibrium states.
Let \(\gamma :(0,\ell)\rightarrow \mathbb {R}^3\) be an arc-length parametrized Frénet curve. Then \(\gamma\) is an equilibrium curve if and only if there is a linear complex that contains none of the binormal lines of \(\gamma\).
“\(\Rightarrow\)”: Assume there exist \(c,\bar{c}\in \mathbb {R}^3\) and \(K:(0,\ell)\rightarrow \mathcal {K}\) such that \(FKk=c\times \gamma +\bar{c}\). We can write the normal component of this equation as \(EIk_n=F_n^T(c\times \gamma +\bar{c})\). From \(I\in S^2_{++}\) and the Frénet assumption \(\kappa = \Vert k_n\Vert \gt 0\), it follows that (15) \(\begin{equation} \begin{aligned}0 &\lt \langle k_n, F_n^T(c\times \gamma +\bar{c}) \rangle = \langle F_nk_n, F_nF_n^T(c\times \gamma +\bar{c}) \rangle \\ &= \langle \omega _n, c\times \gamma +\bar{c}\rangle = \langle \omega _n,\bar{c}\rangle + \langle \gamma \times \omega _n,c\rangle . \end{aligned} \end{equation}\) This shows that \((c,\bar{c})\) is a linear complex that does not contain any of the binormal lines, which have Plücker coordinates \((\omega _n, \gamma \times \omega _n)\).
“\(\Leftarrow\)”: Assume there exists a linear complex \(c,\bar{c}\in \mathbb {R}^3\) that contains no binormal line of \(\gamma\). By continuity, we have \(0 \lt \langle \omega _n,c\times \gamma +\bar{c}\rangle\), possibly after flipping the signs of \(c\) and \(\bar{c}\). Choosing a parallel frame \(F\) adapted to \(\gamma\) yields \(0\lt \langle k_n,F_n^T(c\times \gamma +\bar{c})\rangle\), with \(k\) the curvature vector of \(F\). Thus, there exists \(I\in S^2_{++}\) such that \(EIk_n=F_n^T(c\times \gamma +\bar{c})\).
Next, choose \(J\) such that \(0\lt J\le 4\psi (I)\), and define (16) \(\begin{equation} \tau _\beta =\frac{1}{\mu J}\langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle , \quad \text{and} \quad \beta (s)=\int _0^s \! \tau _\beta . \end{equation}\) Then the frame \(F_\beta\) defined by Equation (1) has material curvatures \(k_{\beta ,n}=Q_\beta ^Tk_n\) and twist \(\tau _\beta\). Assembling the tangential and normal parts of \(c\times \gamma +\bar{c}\) gives \(\begin{equation*} EF_{\beta ,n}Q^T_\beta I Q_\beta k_{\beta ,n} + \mu J \tau _\beta \gamma ^{\prime } = c\times \gamma +\bar{c}, \end{equation*}\) so \(F_\beta\) solves the equilibrium equation with \((Q^T_\beta I Q_\beta , J)\in \mathcal {K}\).□
In summary, every choice of \(c,\bar{c}\in \mathbb {R}^3\) that satisfies (17) \(\begin{equation} \langle \gamma ^{\prime }\times \gamma ^{\prime \prime },c\times \gamma +\bar{c}\rangle \gt 0 \end{equation}\) enables us to pick cross sections that solve the equilibrium equation. Even if we limit the choice of \((I,J)\) to elliptical cross sections, we have a one-dimensional family of ellipses to choose from at every point. This choice will influence both the twist and the ratio between the ellipse radii \(a\) and \(b\). In the next section, we discuss how to compute \(c\) and \(\bar{c}\) in a way that favors low twist and a reasonably small gap between \(a\) and \(b\).
8.2 Linear Program
To choose \(c\) and \(\bar{c}\), we recall a result from Section 7 about parallel equilibrium curves: A non-plane Frénet curve has a parallel equilibrium frame if and only if we can find \(c\) and \(\bar{c}\) that, in addition to Equation (17), also satisfy \(\langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle = 0\). As discussed in Section 7.3, this implies that \(k_n\) and \(F_n^T(c\times \gamma +\bar{c})\) are parallel at every point, with \(m=\kappa /\Vert c \times \gamma +\bar{c}\Vert\) the factor of proportionality. Thus, one solution to \(EIk_n=F_n^T(c\times \gamma +\bar{c})\) is given by \(I=\mathcal {I}_2 / (Eq)\), with \(\mathcal {I}_2\) the 2-by-2 identity matrix. This solution corresponds to a circular cross section, so \(a=b\). This shows that parallel equilibrium curves satisfy both criteria for fabricability to perfection: They have zero twist and allow us to choose \(a/b\) as close to unity as we like.
Unfortunately, the pointwise constraint \(\langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle = 0\) cannot be satisfied exactly for most curves, but we can attempt to satisfy it approximately. To keep the problem linear, a natural choice is to minimize \(\sup |\langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle |\). Using Equation (16) (left), we can also interpret this objective mechanically as the minimization of the twist couple \(\mu J\tau _\beta\). Next, we recall from Proposition 7 that parallel equilibrium curves also satisfy \(\langle \gamma ^{\prime \prime }, c\times \gamma +\bar{c}\rangle =0\). Likewise, we can interpret this expression mechanically by computing \(\begin{equation*} \langle \gamma ^{\prime \prime },c\times \gamma +\bar{c}\rangle = \langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle ^{\prime } = \mu (J\tau _\beta)^{\prime }, \end{equation*}\) so low values of \(|\langle \gamma ^{\prime \prime },c\times \gamma +\bar{c}\rangle |\) correspond to a twist couple function with high \(C^1\)-regularity. We combine both objectives by minimizing \(\sup |B^T(c\times \gamma +\bar{c})|_*\), with \(B=(\gamma ^{\prime },\gamma ^{\prime \prime }/\kappa)\) the orthonormal matrix having the tangent and principal normal as its columns, and \(\begin{equation*} |(x,y)^T|_* = \max \lbrace |x|, w_\text{reg}|y| \rbrace , \end{equation*}\) where \(w_\text{reg}\ge 0\) is a user-controlled regularization weight. Finally, our linear program for determining \(c\) and \(\bar{c}\) reads (18a) \(\begin{align} \text{minimize } &R, \nonumber \nonumber\\ \text{subject to } & 1 \le \langle \gamma ^{\prime }\times \gamma ^{\prime \prime },c\times \gamma +\bar{c} \rangle , \nonumber \nonumber\\ & {-R} \le \langle \gamma ^{\prime }, c\times \gamma + \bar{c} \rangle \le R, \end{align}\) (18b) \(\begin{align} & {-R} \le w_\text{reg}\langle \gamma ^{\prime \prime }/\kappa , c\times \gamma + \bar{c} \rangle \le R, \end{align}\)
where the constraints are enforced at a dense set of samples along the curve. This linear program in the variables \((R,c,\bar{c})\) can be solved near-instantaneously, so the user can interactively browse the family of solutions generated by varying the only parameter, \(w_\text{reg}\).
8.3 Geometry Generation
Given \(c\) and \(\bar{c}\), the next step is to choose a matrix \(I\in S^2_{++}\) at every point that solves \(EIk_n=F_n^T(c\times \gamma +\bar{c})\), where \(k_n\) and \(F_n^T(c\times \gamma +\bar{c})\) are not necessarily parallel. The equation constrains two out of three independent entries of \(S^2_{++}\), which leaves a one-parametric family of solutions. We can compute an explicit representation of the form \(I=I_1+t I_2\), where \(I_1\) and \(I_2\) are symmetric positive semi-definite rank-1 matrices, and \(t\gt 0\) is the free parameter.
Figure 13 shows the family of solutions for different angles between \(k_n\) and \(F_n^T(c\times \gamma +\bar{c})\). A canonical choice that exists in every family and benefits fabricability is the “most circular” ellipse, for which the ratio \(a/b\) is closest to unity. A symbolic computation shows that this solution can be found easily by choosing \(t\) as \(t^*=\operatorname{tr}I_1 / \,\operatorname{tr}I_2\).
For \(t\rightarrow 0\), the ellipse elongates orthogonally to \(F_n^T(c\times \gamma +\bar{c})\), while for \(t\rightarrow \infty\), it elongates along \(k_n\). The latter property is useful for curves that are similar to parallel equilibrium curves, so the family of solutions looks like the one in Figure 13 (left) at most points of the curve. Then choosing \(t\gt t^*\) instead of \(t=t^*\) avoids stability issues for the reasons discussed in Section 7.3. After choosing \(I\), we compute \(J=4\psi (I)\), \(\tau _\beta\) and \(\beta\) via Equation (16). Finally, we can compute \(Q_\beta ^T I Q_\beta\) at every point, which gives the elliptical cross section with the correct rotation.
Fig. 13. Elliptical families. Family of elliptical cross sections that solve the equilibrium equation for fixed \(k_n\) (black) and \(F_n^T(c\times \gamma +\bar{c})\) (gray). The angle between the two vectors increases from left to right. The most circular ellipse in each family ( \(t=t^*\) ) is marked purple.
8.4 Interactive Design
Figure 14 shows the resulting geometry for a horseshoe-shaped input curve that bends out of plane. Even though it is not intuitively obvious what geometry and forcing mechanism will result in a curve like this, our algorithm finds a solution in which the out-of-plane deformation is induced by torque applied to the endpoints. The user interface provides control over \(w_\text{reg}\) to explore the tradeoff between the objective in Equation (18a) and the regularization term in Equation (18b). Prioritizing the former yields the solution on the left, with smaller overall twist, while increasing the regularization weight reduces the total variation, shown on the right.
To validate the solution, we perform a forward simulation of the Kirchhoff rod with the cross sections computed by the computational design algorithm. The boundary conditions are applied by first bending the rod and then twisting its endpoints to cause the out-of-plane deformation. Figure 15 shows frames from the resulting animation, verifying that the desired equilibrium state can indeed by reached by continuous deployment.
8.5 Framed Equilibrium Curves
We conclude this section with a characterization of the set of all framed equilibrium curves. Here we consider the moving frame as part of the “input” instead of deriving it from the input curve. This results in a description of all equilibrium states that can be attained within the theory of Kirchhoff rods with zero natural curvature.
Let \((\gamma ,F)\) be an arc-length parametrized framed Frénet curve. Then the following are equivalent:
Compared to Equation (17), this characterization introduces a sharper lower bound for \(\langle \gamma ^{\prime }\times \gamma ^{\prime \prime },c\times \gamma +\bar{c}\rangle\). This bound is a consequence of the non-linear constraint \(J \le \psi (I)\), which was not used in the characterization of general equilibrium curves. Surprisingly, the new condition is linear in the unknowns \(c\) and \(\bar{c}\) despite this non-linearity and can thus be checked by a linear program.
Fig. 14. Hook curve. Result of the geometry generation algorithm for general equilibrium curves used on a spline curve (top center) without (left) and with regularization (right). Top: Deformed configuration. Center: Twist as a function of arc length. Bottom: Undeformed configuration.
To make the statement of this result rigorous, one needs to control the singularities in the fractional term that may appear at zeroes of \(\tau\). This is done in Appendix B, where we restate the proposition with the level of technical detail necessary to provide a proof.
9 KIRCHHOFF RODS UNDER LOAD
The geometric characterizations in the preceding sections are subject to the limitation that we consider only the elastic energy of Kirchhoff rods and neglect other factors such as the dead load and external loads. However, our computational design algorithms can be adapted to scenarios in which external loads have a small to moderate effect on the equilibrium state.
We will preface this section with the caveat that our proposed algorithm is, in contrast to the ones presented previously, a heuristic and comes without a formal guarantee for finding solutions to all feasible inputs. However, it terminates successfully for all reasonable inputs we have tried, including the examples shown in Section 10. We include the algorithm here it for its practical relevance and as a starting point for future work on this subject.
9.1 Load Model
To model line loads (such as the dead load) and point loads (such as a weight hanging from a specific point) together, we use a load distribution \(q(s)=p(s)+\sum _{i=1}^n\! \delta (s-s_i) f_i\). Here \(p:(0,\ell)\rightarrow \mathbb {R}^3\) models line loads (force per length) applied to \(\gamma\), and \(f_i\in \mathbb {R}^3,i=1,\ldots ,n,\) are concentrated forces applied at \(\gamma (s_i)\), which are modeled as delta distributions. Integrating \(q\) on an interval \(I\subset (0,\ell)\) gives the accumulated force applied to the rod on this interval.
The effect of this load distribution is captured by adding the potential \(-\int _0^\ell \langle q,\gamma \rangle\) to the Lagrangian from Equation (5). This modification leads to the equilibrium equation
(19) \(\begin{equation} (FKk)(s)=(c+Q(s))\times \gamma (s) + (\bar{c}+M(s)), \end{equation}\) where \(Q(s):=\int _0^s \!q\) is the accumulated force and \(M(s):=\int _0^s \!\gamma \times q\) the accumulated moment.
Fig. 15. Hook curve validation. Forward simulation of a Kirchhoff rod generated by the computational design algorithm to verify that the equilibrium state can be reached. The boundary conditions are applied gradually by twisting the endpoints.
9.2 Design Algorithm
The feasibility condition from Equation (17) can easily adapted to a given load case by substituting \(c\) and \(\bar{c}\) with the expressions appearing in Equation (19), which yields \(\begin{equation*} \langle \gamma ^{\prime }\times \gamma ^{\prime \prime },(c+Q)\times \gamma + (\bar{c}+M)\rangle \gt 0. \end{equation*}\) If the load is thought to be fixed, then we can check for the existence of \(c,\bar{c}\in \mathbb {R}^3\) satisfying this condition using a linear program as usual.
However, if we consider the dead load of a rod, then the rod geometry and the load are coupled, because the cross sections determine the weight per unit arc-length. More precisely, the choice of \((c,\bar{c})\) determines the cross sections, which in turn determine the load, and thus \(M\) and \(Q\), which appear in the equilibrium equation. This dependence of \(M\) and \(Q\) on \((c,\bar{c})\) is non-linear, so we can no longer check for feasibility with a linear program.
Geometry and Dead Load. In preparation of our computational design algorithm under dead load, we discuss the relationship between rod geometry and dead load in more detail. The geometry is encoded in the bending stiffness \(I(s)\) via Equation (7), and the dead load is encoded as a line load \(p(s)=A(s)\rho g\), where \(A(s)\) is the cross-sectional area at \(s\in (0,\ell)\), \(\rho \gt 0\) is the material density, and \(g\in \mathbb {R}^3\) is the gravitational acceleration.
Next, we define two maps: one from \(I\) onto \(p\) and one from \(p\) onto \(I\). The first map, \(p(I)\), computes the dead load for a (fixed) geometry via \(p(s)=2\sqrt {\pi }(\det I(s))^{1/4}\rho g\). This holds, because \(A(s)=2\sqrt {\pi }(\det I(s))^{1/4}=\pi a(s) b(s)\) is the area of an ellipse with radii \(a(s),b(s)\). The second map, \(I(p)\), computes the geometry that equilibrates a (fixed) load \(p\). This map is defined by the steps in Section 8.3, with the difference that \(c\) is replaced by \(c+Q\), and \(\bar{c}\) by \(\bar{c}+M\).
Note that the maps \(p(I)\) and \(I(p)\) are not generally inverses of each other: We only have \(p^* = p(I(p^*))\) for some \(p^*:(0,\ell)\rightarrow \mathbb {R}^3\) if \(\gamma\) is an equilibrium curve for the geometry \(I(p^*)\) under its own dead load. Thus, solving the inverse design problem under dead load is equivalent to finding a fixed point of \(p\circ I\).
Algorithm. We propose to first fix \((c,\bar{c})\) heuristically (Step 1), and then apply a fixed-point iteration procedure (Step 2):
Step 1. We compute the load-free solution to Equation (18), which yields constants \(c_0\) and \(\bar{c}_0\) and a first guess of the rod geometry, encoded by the bending stiffness \(I_0(s)\). We compute the corresponding dead load \(p_0=p(I_0)\) and the accumulated force \(Q_0\) and moment \(M_0\). Our heuristic choice for the constants is \(c:=c_0-\frac{1}{2}(\inf Q_0 + \sup Q_0)\) and \(\bar{c}:=\bar{c}_0-\frac{1}{2}(\inf M_0 + \sup M_0)\). The rationale behind this is that the functions \(c+M_0(s)\) and \(\bar{c}+Q_0(s)\)—which replace \(c_0\) and \(\bar{c}_0\)—will be as close as possible to \(c_0\) and \(\bar{c}_0\), in the uniform norm.
Step 2. To find a fixed point of \(p\circ I\), we iterate \(p_{i+1}=p(I(p_i))\), until \(\sup |p_{i+1}-p_i|\lt \varepsilon\), where we set \(\varepsilon =10^{-10}\) in our examples. We cannot offer a formal proof of convergence but see experimentally that convergence is at least linear, and our examples terminate after fewer than 10 iterations, which we show in Section 10.2.
10 RESULTS
We implemented the computational design algorithms for parallel equilibrium curves (Section 7) and general equilibrium curves (Section 8) in a software tool that allows users to interactively design Kirchhoff rods. Edits made by the user trigger the computational design algorithms to re-run and display the resulting geometry near-instantaneously to enable fast prototyping.
Once a design has been finalized, the resulting geometry can be exported for fabrication as parametric CAD geometry. This allows our tool to be integrated seamlessly into a CAD workflow, for example to create molds for fabrication or design a support structure on which the Kirchhoff rods can be mounted.
Below we show objects designed with our algorithms and fabricated by casting silicone in 3D-printed two-part molds. The first example demonstrates the fabrication method and the design space of parallel equilibrium curves; subsequent examples use the algorithm for general equilibrium curves (with and without dead load and external loads) and show applications in soft robotics and design with active bending.
10.1 Parallel Equilibrium Curve and Fabrication
Rod Design. To design a Kirchhoff rod within the constrained space of parallel equilibrium curves (Section 7), the user is given control over the geometric curvature function \(\kappa\) and the three scalar parameters discussed in Section 7.2. Together, these quantities uniquely determine a parallel equilibrium curve up to rigid motion.
Initially, \(\kappa\) is set to be constant, which generates a segment of a double helix, but the user can modify \(\kappa\) by adding and dragging sample points to bend or straighten the curve in certain locations. Figure 16 shows three snapshots from a design session, in which the user progressively edits a parallel equilibrium curve. The preview of the rod geometry that realizes this curve is updated in real time as the control points are being dragged.
Fig. 16. Parallel curve design. Top: In-app preview of the deformed rod geometry. Bottom: Editable geometric curvature. From left to right: Starting from a constant-curvature rod, the user adds more samples to the curvature graph to locally straighten or tighten the windings of the double helix. The curves from the helical motion (gray) serve as a visual guide.
Mold Design. Once a design is finished, it can be exported as a FeatureScript for the CAD system Onshape to generate solid parts of the undeformed and deformed rod geometries. The undeformed rod serves as a starting point for designing a 3D-printable mold, which is used for silicone casting during fabrication. To simplify this, our app also exports a parting surface that splits the mold into two parts, each guaranteed to be a height field.
The deformed rod is used to design a support structure, with sockets that enforce the kinematic boundary conditions at both endpoints of the rod. The automatically generated CAD geometry and the manually designed mold and support structure are shown in Figure 18.
Fabrication. We 3D print the two-part mold from PLA on an Ultimaker S5, close it, and seal the seams with gaffer tape. For casting, the mold is placed vertically, which is important to prevent the formation of air bubbles. Next, we use a syringe to inject liquid silicone (SmoothSil 945) through the injection hole located near the bottom of the mold. The silicone rises through the cavity until it reaches the air vent at the top, at which point we seal the injection hole with a rubber plug.
After a 16-hour curing period, the rod is ready to be unmolded and mounted in the 3D-printed support structure, as shown in the photograph in Figure 22. The helical windings, small thickness of the rod, and material chosen for this example all contribute to a low overall stiffness, which makes the rod sag visibly under gravity, compared to the target design. This illustrates the relevance of accounting for gravity in inverse design even on this scale—which we do in the following examples.
10.2 Loop Array
The addition of twist opens up the design space of Kirchhoff rods and gives more creative freedom to the designer by enabling direct control via spline editing tools or specifying curves analytically. These curves serve as input to the general equilibrium curve algorithm discussed in Section 8 and can be post-processed by the algorithm in Section 9 to account for the dead load.
In this example, we design and fabricate an array of six curves from a smooth family, shown in Figure 17. The rods in purple are the result of solving the linear program in Equation (18), which neglects gravity, and then forward-simulating the geometry with gravity. This makes the rods sag noticeably under their own weight and prevents a faithful reproduction of the input curve.
Dead Load. To improve reproduction of the input curve under gravity, we apply the dead load optimization algorithm from Section 9.2. The inset shows the fixed-point error \(\sup |p_{i+1}-p_{i}|\) as a function of the iteration count \(i\) on a log scale, which provides numerical evidence that convergence is at least linear. Every graph in green corresponds to a rod from this example, while the ones in yellow and purple correspond to the rods from Sections 10.3 and 10.5, respectively.
The rods shown in green in Figure 17 are the result of our optimization and reproduce the input curves perfectly (up to numerical error). The undeformed rod geometries before and after dead load optimization, as shown in Figure 19, are visually very similar. This is because the geometry shown in purple serves as the initial value for the fixed-point iteration, which naturally converges to an attractive fixed point (green) close to it. Nevertheless, this small change suffices to counteract the effect of gravity. Renderings and photographs of the result are shown in Figures 1 and 22 for direct comparison.
Performance. Solving the linear program from Equation (18) for one of the input curves (with 400 sample points) takes about 5 ms and computing the rod geometry according to Section 8.3 an additional 6 ms. The dead load fixed-point iteration takes between six and eight iterations to converge, and each iteration takes about 6 ms. Therefore, the total computation time per curve is about 60 ms, fast enough to show geometry changes due to user edits in real time.
10.3 Fixture Design
A popular application of bending-active materials is the design of structures that take their final shape only under the effect of a weight, such as a lampshade hanging from it, or a person sitting on it. Our system supports the design of objects like this by using the load distribution described in Section 9, which models external line loads and point loads in addition to the dead load.
We demonstrate this feature by designing a fixture that is in equilibrium under its own weight plus a weight hanging from a specific point. The input curve is designed manually by manipulating the control points of a quartic B-spline curve. As usual, the user sees immediately after each edit how the geometry of the rod and the twist of the equilibrium state were affected by the change.
Fig. 18. Loop array deformed. Front view (top) and side view (bottom) of rods from the loop array example. The final result, taking into account gravity for inverse design, is shown in green; the result of neglecting gravity during design, but forward-simulating with gravity, in purple.
Load Optimization. The external weight hanging from the rod is modeled as a point load \(f_1\in \mathbb {R}^3\) at a curve point \(\gamma (s_1)\) as described in Section 9.1. During the fixed-point procedure, the load distribution \(q(s)=p(s)+\delta (s-s_1) f_1\) takes into account both the dead load, which is updated in every iteration, and the external weight, which remains constant.
Figure 20 shows the evolution of the dead load \(i\) during the fixed-point procedure. After a single iteration, the solution is converged enough for the graph of \(p\) to remain visually unchanged afterwards, and after seven iterations, the pointwise change is below \(\varepsilon =10^{-10}\). As seen in the top-right part of Figure 20, the algorithm adds twist near the bottom endpoint. This extra twist has the effect of lifting up the hook enough to counter the force introduced by the external weight.
The inset shows the final geometry, forward-simulated without the external weight (purple), in which case there is a large deviation from the target curve, and with the external weight (green), in which case the target curve is matched precisely. Figure 22 offers a side-by-side comparison between photographs of the physical model and renderings of the target shape, showing a good agreement.
Fig. 17. Parallel curve mold. Top and bottom right: Undeformed rod geometry (green) with parting surface (purple) and two-part mold (gray). Bottom left: Rendering of deformed rod (green) with support structure (gray).
10.4 Soft Robotics
Elastomers like silicone are used in soft-robotics applications for the design of soft grippers or tools for minimally invasive surgery. Most of the soft-robotic mechanisms currently in use are actuated pneumatically or by cables [Runciman et al. 2019]. Both actuation systems add considerable complexity to the compliant part of the mechanism, through a sequence of air chambers or a network of cables around or inside the part.
With our algorithm, we can design simple compliant mechanisms that are actuated by twisting the endpoints of a rod. This does not add any complexity to the compliant part itself, because one only needs to add twisting joints to the fixture. The actuation itself can be performed by motors in the fixture, to follow a prespecified trajectory, or manually if human fine-motor skills are required. Either way, this shifts the complexity away from the compliant part of the mechanism to an outside controller.
Fig. 19. Loop array undeformed. Undeformed geometry of rods from the loop array before (purple) and after (green) applying the dead load fixed-point iteration.
Lifting Tool. We show a macro-scale version of a tool that can be used to lift objects out of an inaccessible location. The deformed configuration of the tool is given by the geometry shown in Figure 14 (right), where the out-of-plane deformation of the horse shoe is caused by twist applied to the endpoints. We 3D print a mechanical fixture that connects the endpoints to rigid bars that a human user can turn to change the twisting angle from a distance, and to control the amount of out-of-plane bending.
The sequence of photographs in Figure 23 shows a usage scenario in which the tool is inserted into a tunnel with a sudden drop at the end. Using the twisting actuation, the flexible part of the tool bends downwards to grab a box and slide it up along the wall of the protrusion, so it can be pulled back through the tunnel.
10.5 Free-form Light Sculpture
We show an application to interior design, inspired by the Freeform Light Sculpture series of New York artist John Procario.3 Our design mimics the organic wooden shapes of the original sculptures with black rubber beams that are made to follow a three-dimensional curved target shape by taking advantage of elastic and gravitational forces.
For this example, we forgo spline curves in favor of specifying the target design with a mathematical expression for a closed curve, as shown in the inset: \(\begin{equation*} \gamma (t)=\begin{pmatrix} r_1\cos t+r_2\cos (t/2+p) \\ r_1\sin t+r_2\sin (t/2+p) \\ a \cos (3t/2) \end{pmatrix}, \end{equation*}\) with \(t\in (0,4\pi)\), and \(r_1=1\), \(r_2=1/4\), \(p=9/20\), and \(a=2/5\). We split the curve into three segments, \((0,4\pi /3)\), \((4\pi /3,8\pi /3)\), and \((8\pi /3,4\pi)\), and compute separately for each the geometry of a Kirchhoff rod, taking into account the dead load as per Section 9.
Figure 21 illustrates that gravity takes a central role in shaping the final deformed shape of this model. If we neglect the dead load while computing the rod geometry (top-left), then the resulting deformed shape under gravity is saggy and very far from the target. If the dead load is accounted for during design, then we have perfect agreement (top-right). The figure also shows the iterations of the dead-load optimization for each segment: After two iterations, the solutions are converged almost completely.
We manufactured the model using SmoothSil 960 silicone, with black pigments added for coloring. The design features a small indentation running along the inside of the rod, in which we place an electroluminescent wire. Figures 1 and 22 show a rendering of the target design, and photos of the physical model with and without external lighting from a similar perspective, to allow for a visual comparison.
Fig. 20. Load optimization. Top: Deformed geometry before (left) and after (right) accounting for external load and dead load. Bottom: Vertical component of the dead load before and after fixed-point iteration.
11 DISCUSSION
In this work, we characterize the design space of Kirchhoff rods with spatially varying cross sections and vanishing natural curvature. This geometric characterization gives rise to computational design algorithms that translate a curve in three-dimensional space to a Kirchhoff rod that attains this curve at equilibrium, given appropriate boundary conditions. We also discuss an extension that takes into account gravity, to enable applications on a larger scale.
A current limitation of our algorithm is that stability of the target equilibrium state is only enforced heuristically and verified after the design stage through numerical means. However, stability cannot be rigorously guaranteed by our computational design algorithm. One could approach this issue by combining the adjoint method for stability optimization [Hafner and Bickel 2021] with the constrained Jacobi equation for Kirchhoff rods [Manning et al. 1998]. Likewise, if the user specifies an infeasible input curve, or one that leads to geometry that is too difficult to fabricate, then there is currently no fully automatic system in place to improve the design computationally. Doing this with iterative optimization methods is an avenue for future work.
Fig. 21. Light sculpture design. Top: Deformed configuration under gravity if dead load is neglected (left) or accounted for (right) by the design algorithm. Center: Undeformed configuration of final design. Bottom: Iterations of the dead-load optimization for all three segments.
Fig. 22. Results. Comparison between renderings of the target designs and photos of our elastic physical models. (1) The loop array from the front (a), cf. Figure 1 (left), and from the side (b). (2) The fixture from the side (a) and from the top (b), conforming to the target shape when equipped with a weight of 100 g. (3) The free-form light sculpture, cf. Figure 1 (right). (4) The parallel equilibrium curve example, which sags due to a lack of gravity correction.
Fig. 23. Soft robotics. The lifting tool is pushed into the cavity from the left (a) and bends downwards out of plane once the user applies twist to the endpoints of the mechanism (b). This allows the tool to grab and lift the box (c) and pull it out of the cavity (d). Sequence shown in supplemental video.
Allowing arbitrary elliptical cross sections for rods has the disadvantage that one needs to use molding or three-axis CNC milling for fabrication. It would therefore be of practical relevance to study the design of rods that can be produced with simpler means, e.g., laser cutting or water jetting rods from a sheet material or lathing to produce circular cross sections. Geometries that can be manufactured using these methods give rise to a variety of design subspaces whose exploration could be of great practical relevance and may enable computational design algorithms with a wider range of applications. In the same vein, some fabrication methods can produce rods with non-vanishing natural curvature, for example, by cutting a strip with a curved centerline from a sheet material. Treating this natural curvature as an additional design variable may enlarge the design space of equilibrium states, and lead to a greatly simplified fabrication pipeline.
The authors consider Proposition 10, the characterization of framed equilibrium curves, the most surprising result in this article. This is because it constitutes a linear representation of a design space that is subject to non-linear constraints with no apparent special structure except convexity. Is this a lucky coincidence specific to Kirchhoff rods or could it hold more generally for inverse problems associated with other mechanical systems? Answering this question for two- or three-dimensional systems, for example to understand the design space of general elasticity, may be a rewarding endeavor with immediate applications to the design of compliant mechanisms.
APPENDICES
A PROOF OF PROPOSITION 6
We can rewrite Equation (6) as \(EF_nIk_n+\mu J \tau \gamma ^{\prime } = c\times \gamma +\bar{c}\).
“\(\Rightarrow\)”: By assumption, there exist a parallel frame \(F\) adapted to \(\gamma\), constants \(c,\bar{c}\in \mathbb {R}^3\), and \(I:(0,\ell)\rightarrow S^2_{++}\) such that \(EF_n Ik_n=c\times \gamma +\bar{c}\), because \(\tau \equiv 0\). Taking the inner product with \(\gamma ^{\prime }\) implies \(\begin{equation*} 0 = \langle \gamma ^{\prime }, c\times \gamma +\bar{c}\rangle = \langle \gamma ^{\prime }, \bar{c}\rangle + \langle \gamma \times \gamma ^{\prime }, c\rangle , \end{equation*}\) so the linear complex \(\mathcal {C}\) associated with \((c,\bar{c})\) contains all tangent lines of \(\gamma\), with Plücker coordinates \((\gamma ^{\prime }, \gamma \times \gamma ^{\prime })\). The normal part of the equilibrium equation gives \(EIk_n = F_n^T(c\times \gamma +\bar{c})\), so (20) \(\begin{equation} \begin{aligned}0 &\le \langle k_n, F_n^T(c\times \gamma +\bar{c}) \rangle = \langle F_nk_n, F_nF_n^T(c\times \gamma +\bar{c}) \rangle \\ &= \langle \omega _n, c\times \gamma +\bar{c}\rangle = \langle \omega _n,\bar{c}\rangle + \langle \gamma \times \omega _n,c\rangle , \end{aligned} \end{equation}\) where the inequality follows from \(I\in S^2_{++}\). Because \(I\) has full rank, we have \(k_n=0\) exactly where \(F^T_n(c\,\times \,\gamma +\bar{c})=0\), which is equivalent to \(c\times \gamma +\bar{c}=0\), because \(c\times \gamma +\bar{c}\) is in the column space of \(F_n\).
(1) Assume \(\mathcal {C}\) is regular, so \(\langle c,\bar{c}\rangle \ne 0\). The column space of \([c]_\times\) is exactly the orthogonal complement of \(c\), so there are no curve points satisfying \(c\times \gamma +\bar{c}=0\), and in consequence no points with \(k_n=0\). Thus, \(\gamma\) is a Frénet curve. This also shows \(0\lt \langle \omega _n,\bar{c}\rangle + \langle \gamma \times \omega _n,c\rangle\), so no binormal line is contained in \(\mathcal {C}\).
For the sake of contradiction, assume that \(\gamma\) was plane. Then all tangent lines of \(\gamma\) are coplanar, and we can choose three that are not concurrent in a point. However, no such set of three lines is contained in a regular linear complex, so \(\gamma\) must not be plane.
(2) Assume \(\mathcal {C}\) is singular, so \(\langle c,\bar{c}\rangle = 0\). If \(c=0\), then all tangent lines of \(\gamma\) are orthogonal to \(\bar{c}\), and \(\gamma\) must be plane. From \(c\times \gamma +\bar{c}=\bar{c}\ne 0\), we see that \(k_n\ne 0\), so \(\gamma\) has no inflection points.
We can show that the case \(c\ne 0\) reduces exactly to Theorem 5. We know by assumption that all tangent lines of \(\gamma\) intersect the Euclidean line \(\lambda (c,\bar{c})\), as shown in Figure 9 (left). This implies that every connected component of \(\gamma ((0,\ell)) \setminus \lambda (c,\bar{c})\) is contained in a plane \(P\) with \(\lambda (c,\bar{c})\subset P\). The points \(\gamma (s_0)\in \lambda (c,\bar{c})\) are exactly the inflection points of \(\gamma\), so there are finitely many connected components and planes.
Consider an interval \((s_0-\varepsilon , s_0+\varepsilon)\) small enough such that \(\gamma ((s_0-\varepsilon ,s_0))\) is contained in a plane \(P\) and \(\gamma ((s_0,s_0+\varepsilon))\) in a plane \(\bar{P}\). Assume for the sake of contradiction that \(P \ne \bar{P}\), so \(P\cap \bar{P} = \lambda (c,\bar{c})\). Then the tangent line of \(\gamma\) at \(s_0\) must coincide with \(\lambda (c,\bar{c})\).
Assume wlog that \(\gamma (s_0)\) coincides with the origin, \((c,\bar{c})=(e_2,0)\), and \(P\) is normal to \(e_3\). Next, consider the parallel frame \(\bar{F}=(\bar{n}_1, \bar{n}_2, \gamma ^{\prime })\) such that \(\bar{n}_1\equiv -e_3\) on the interval \((s_0-\varepsilon , s_0)\). This frame is related to \(F\) by a constant rotation \(Q\in SO(2)\) and satisfies the equilibrium equation with \(\bar{I}=Q^T I Q\): \(\begin{equation*} E\bar{F}_n \bar{I} \bar{k}_n = E (F_n Q) (Q^T I Q) (Q^T k_n) = E F_n I k_n = c\times \gamma +\bar{c}, \end{equation*}\) according to Equation (1). Noting that \(\bar{\kappa }_2\equiv 0\), \(\langle \bar{n}_2,e_1\rangle = -\langle \gamma ^{\prime }, e_2\rangle\), and \(\langle \bar{n}_2,e_2\rangle = \langle \gamma , e_1\rangle\) on \((s_0-\varepsilon , s_0)\), we can write out the equilibrium equation in coordinates, which yields \(\begin{equation*} \bar{I}_{xy}\bar{\kappa }_1\langle \gamma ^{\prime },e_1\rangle = 0 = \bar{I}_{xy}\bar{\kappa }_1\langle \gamma ^{\prime },e_2\rangle , \quad E \bar{I}_{xx}\bar{\kappa }_1 = \langle \gamma , e_1\rangle . \end{equation*}\) The two equations on the left imply \(\bar{I}_{xy}\bar{\kappa }_1\equiv 0\) and thus \(\bar{I}_{xy}\equiv 0\), because \(\bar{\kappa }_1(s)=0\) only at \(s=s_0\). The equation on the right is exactly the equilibrium equation from Theorem 5, with \(\lambda (e_2,0)\) equivalent to \(L\). That \(\lambda (e_2,0)\) is tangent to \(\gamma\) contradicts Theorem 5(1), and shows \(P=\bar{P}\). By repeating the argument, we see that every connected component of \(\gamma ((0,\ell)) \setminus \lambda (c,\bar{c})\) is contained in the same plane.
“\(\Leftarrow\)”: Assume either (1) or (2) holds, and let \(F\) be a parallel frame adapted to \(\gamma\). For (2), there is nothing to prove, because the result follows directly from Theorem 5. If (1) holds for some linear complex defined by \(c,\bar{c}\in \mathbb {R}^3\), then we have \(0=\langle \gamma ^{\prime }, c\times \gamma +\bar{c}\rangle\). We also have \(0\lt \langle \omega _n, c\times \gamma +\bar{c}\rangle\), after possibly replacing \((c,\bar{c})\) with \((-c,-\bar{c})\), because \(\omega _n(s)\) and \(c\times \gamma (s)+\bar{c}\) are continuous. This implies \(0 \lt \langle k_n, F_n^T(c\times \gamma +\bar{c})\rangle\), so we can find \(I:(0,\ell)\rightarrow S^2_{++}\) such that \(EIk_n = F_n^T(c\times \gamma +\bar{c})\). This implies \(EF_nIk_n=c\times \gamma +\bar{c}\).
B FRAMED EQUILIBRIUM CURVES
Let \(\gamma \in C^2((0,\ell);\mathbb {R}^3)\) be an arc-length parametrized Frénet curve and \(F:(0,\ell)\rightarrow SO(3)\) a moving frame adapted to \(\gamma\) such that all zeros of \(\tau\) are isolated. Formally define \(\begin{equation*} v_1=\frac{\gamma ^{\prime }}{\tau }, \quad v_2=\frac{4\omega _n}{E\kappa ^2}-\frac{\gamma ^{\prime }}{\mu \tau }. \end{equation*}\) Then the following are equivalent:
(1) | The framed curve \((\gamma ,F)\) is an equilibrium curve with a constitutive relation \((I,J):(0,\ell)\rightarrow \mathcal {K}\) that is coercive and bounded.4 | ||||
(2) | There exist \(c,\bar{c}\in \mathbb {R}^3\) such that \(\langle v_i,c\times \gamma +\bar{c}\rangle \gt 0\) for \(i=1,2\) in the sense that the inequalities hold away from the zeros of \(\tau\), and the upper and lower limits are finite and positive as \(\tau \rightarrow 0\). | ||||
\((1)\Rightarrow (2)\): By (1), we assume that there exists \((I,J):(0,\ell)\rightarrow \mathcal {K}\) satisfying the equilibrium equation \(EF_nIk_n+\mu J\tau \gamma ^{\prime }=c\times \gamma +\bar{c}\), such that \(J\) and the eigenvalues of \(I\) are bounded from below and above by positive constants. Furthermore, \(\kappa =\Vert k_n\Vert\) is also bounded from below and above by positive constants, because it is positive and continuous on \([0,\ell ]\).
Taking the inner product between the equilibrium equation and \(\gamma ^{\prime }\) yields \(\mu J\tau =\langle \gamma ^{\prime },c\times \gamma +\bar{c}\rangle\), which shows that \(\langle \gamma ^{\prime }/\tau , c\times \gamma +\bar{c}\rangle\) has the required properties by the coercivity and boundedness of \(J\).
To show the second inequality, take the inner product between the equilibrium equation and \(\omega _n/(E\kappa ^2)\), which yields \(\begin{equation*} \lambda _1 \le \left\langle \frac{k_n}{\kappa }, I\frac{k_n}{\kappa } \right\rangle = \left\langle \frac{\omega _n}{E\kappa ^2},c\times \gamma +\bar{c} \right\rangle \le \lambda _2, \end{equation*}\) with \(\lambda _1\le \lambda _2\) the eigenvalues of \(I\). Furthermore, from \(\begin{equation*} J\le 4\psi (I)\le 4\lambda _1\frac{\lambda _2}{\lambda _1+\lambda _2} \end{equation*}\) and the coercivity and boundedness of \(I\), we find that there exists \(\varepsilon \gt 0\) such that \(J+\varepsilon \le 4\lambda _1\). Combining these inequalities with the expression for \(\mu J\tau\) from above, we get \(\begin{equation*} \varepsilon = (J+\varepsilon)-J\le \left\langle \frac{4\omega _n}{E\kappa ^2}-\frac{\gamma ^{\prime }}{\mu \tau },c\times \gamma +\bar{c} \right\rangle \le 4 \sup \lambda _2, \end{equation*}\) which shows the required bounds.
\((2)\Rightarrow (1)\): Define \(J=\langle \frac{\gamma ^{\prime }}{\mu \tau },c\times \gamma +\bar{c} \rangle\) to satisfy the tangential component of the equilibrium equation. From the limit properties of \(\langle \gamma ^{\prime }/\tau , c\,\times \,\gamma +\bar{c}\rangle\) at zeros of \(\tau\), we see that \(J\) is coercive and bounded, and we can extend it with arbitrary positive values at these zeros.
Satisfying the normal component equation \(EIk_n=F_n^T(c\times \gamma +\bar{c})\) with an appropriate choice of \(I\) is equivalent to satisfying \(E \tilde{I} Qk_n = QF_n^T (c\times \gamma +\bar{c})\) for some \(Q\in SO(2)\) with an appropriate choice of \(\tilde{I}\), because we can transform \(\tilde{I} = QIQ^T\). Furthermore, \(I\) inherits coerciveness and boundedness from \(\tilde{I}\), so it suffices to show these properties for the latter.
We choose \(Q\in SO(2)\) such that \(Qk_n=\kappa e_1\), so it suffices to find \(\tilde{I}\) such that \(E\kappa \tilde{I} e_1 = QF_n^T(c\times \gamma +\bar{c})\). Choosing \(\tilde{I}_{xx}=\frac{1}{E\kappa ^2}\langle \omega _n,c\times \gamma +\bar{c}\rangle\) satisfies the first component of this equation, and we can uniquely determine \(\tilde{I}_{xy}\) from the second. Next, we pick \(\tilde{I}_{yy}=\frac{J\tilde{I}_{xx}+4\tilde{I}_{xy}^2}{4 \tilde{I}_{xx}-J}\) to satisfy \(J=4\psi (\tilde{I})\), which is checked by direct computation.
To verify that \(\tilde{I}\) is coercive and bounded, it suffices to show that \(\operatorname{tr}\tilde{I}\) and \(\det \tilde{I}=\frac{J(\tilde{I}_{xx}^2+\tilde{I}_{xy}^2)}{4\tilde{I}_{xx}-J}\) are bounded from below and above by positive constants. This can be seen from the formula \(2\lambda _{1,2}=\operatorname{tr}\tilde{I} \pm \sqrt {\operatorname{tr}^2 \tilde{I} - 4\det \tilde{I}}\). For \(\tilde{I}_{xx}\), boundedness is clear by continuity and from \(4\tilde{I}_{xx} = \langle v_1/\mu \, + \, v_2, c \, \times \, \gamma +\bar{c}\rangle \gt 0\). For \(\tilde{I}_{yy}\) and \(\det \tilde{I}\), boundedness of the numerator is clear from the boundedness of \(\tilde{I}_{xx}\) and \(J\), and boundedness of the denominator can be seen from \(\begin{equation*} 4\tilde{I}_{xx}-J = \left\langle \frac{4\omega _n}{E\kappa ^2}-\frac{\gamma ^{\prime }}{\mu \tau },c\times \gamma +\bar{c} \right\rangle , \end{equation*}\) and the limit properties of the expression on the right-hand side, which follow from (2). We conclude that \(\lambda _1\) and \(\lambda _2\) are also bounded from below and above by positive constants, which gives the required coercivity and boundedness of \(I\).□
Footnotes
1 The expression \(\langle \gamma ^{\prime } \times \gamma ^{\prime \prime }, c\times \gamma +\bar{c} \rangle\) is continuous in the curve parameter. Thus, non-orthogonality is equivalent to a constant sign of this inner product. Should the sign be negative, we can replace \((c,\bar{c})\) with \((-c,-\bar{c})\) to produce an inner product that is strictly positive everywhere.
Footnote2 Proof: Let \(w\in \mathbb {R}^3\). Then \(Q[v]_\times Q^T w=Q(v\times (Q^T w)) = (Qv)\times w=[Qv]_\times w\).
Footnote3 http://www.johnprocario.com/.
Footnote4 There exist \(c,C\gt 0\) such that for all \(s\in (0,\ell)\), it holds that \(c\le J(s)\le C\) and \(c\le \lambda _1(s),\lambda _2(s) \le C\) with \(\lambda _1\le \lambda _2\) the eigenvalues of \(I\).
Supplemental Material
Available for Download
Supplementary materials
- . 2010. Elasticity and Geometry: From Hair Curls to the Non-linear Response of Shells. Oxford University Press, Oxford, UK.
2010017278 Google Scholar
- . 2010. Discrete viscous threads. ACM Trans. Graph. 29, 4 Article
116 (07 2010), 10 pages.DOI: Google ScholarDigital Library
- . 2008. Discrete elastic rods. In ACM SIGGRAPH 2008 Papers (Los Angeles, California) (SIGGRAPH’08). Association for Computing Machinery, New York, NY, Article 63, 12 pages. Google ScholarDigital Library
- . 2017. State of the art in methods and representations for fabrication-aware design. Comput. Graph. Forum 36, 2 (
May 2017), 509–535.DOI: Google ScholarDigital Library - . 2006. Super-helices for predicting the dynamics of natural hair. ACM Trans. Graph. 25, 3 (
July 2006), 1180–1187.DOI: Google ScholarDigital Library - . 2018. Inverse design of an isotropic suspended Kirchhoff rod: theoretical and numerical results on the uniqueness of the natural shape. Proc. Roy. Soc. A: Math. Phys. Eng. Sci. 474, 2212 (
April 2018), 20 pages.DOI: Google ScholarCross Ref
- . 2015. Computational design of metallophone contact sounds. ACM Trans. Graph. 34, 6, Article
223 (October 2015), 13 pages.DOI: Google ScholarDigital Library - . 2018. State of the art on stylized fabrication. Comput. Graph. Forum 37, 6 (
February 2018), 325–342.DOI: Google ScholarCross Ref - . 2014. Quasi-static manipulation of a Kirchhoff elastic rod based on a geometric analysis of equilibrium configurations. Int. J. Robot. Res. 33, 1 (2014), 48–68.Google ScholarDigital Library
- . 2014. An asymptotic numerical method for inverse elastic shape design. ACM Trans. Graph. 33, 4, Article
95 (July 2014), 11 pages.DOI: Google ScholarDigital Library - . 2013. Inverse dynamic hair modeling with frictional contact. ACM Trans. Graph. 32, 6 (
November 2013), 1–10.DOI: Google ScholarDigital Library - . 2010. Stable inverse dynamic curves. ACM Trans. Graph. 29, 6, Article
137 (December 2010), 9 pages.DOI: Google ScholarDigital Library - . 1948. The torsional rigidity and variational methods. Am. J. Math. 70, 1 (1948), 107–116. http://www.jstor.org/stable/2371935Google ScholarCross Ref
- . 2016. Programming curvature using origami tessellations. Nat. Mater. 15, 5 (2016), 583–588.Google ScholarCross Ref
- . 2020. RoboCut: Hot-wire cutting with robot-controlled flexible rods. ACM Trans. Graph. 39, 4, Article
98 (July 2020), 15 pages.DOI: Google ScholarDigital Library - . 2017. CurveUps: Shaping objects from flat plates with tension-actuated curvature. ACM Trans. Graph. 36, 4 (
July 2017), 1–12.DOI: Google ScholarDigital Library - . 2006. Oriented strands: Dynamics of stiff multi-body system. In Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA’06). Eurographics Association, Goslar, DEU, 91–100. Google Scholar
- . 2021. The design space of plane elastic curves. ACM Trans. Graph. 40, 4, Article
126 (July 2021), 20 pages.DOI: Google ScholarDigital Library - . 2019. Vibration-minimizing motion retargeting for robotic characters. ACM Trans. Graph. 38, 4, Article
102 (July 2019), 14 pages.DOI: Google ScholarDigital Library - . 2018. Rapid deployment of curved surfaces via programmable auxetics. ACM Trans. Graph. 37, 4 (
August 2018), 1–13.DOI: Google ScholarDigital Library - . 2013. Active bending, a review on structures where bending is used as a self-formation process. Int. J. Space Struct. 28, 3-4 (
September 2013), 187–196.DOI: Google ScholarCross Ref - . 2020. Tapered elasticæ as a route for axisymmetric morphing structures. Soft Matter 16, 33 (2020), 7739–7750.
DOI: Google ScholarCross Ref - . 2018. FlexMaps: Computational design of flat flexible shells for shaping 3D objects. ACM Trans. Graph. 37, 6, Article
241 (December 2018), 14 pages.DOI: Google ScholarDigital Library - . 1998. Isoperimetric conjugate points with application to the stability of DNA minicircles. Proc. Math. Phys. Eng. Sci. 454, 1980 (1998), 3047–3074. http://www.jstor.org/stable/53424Google ScholarCross Ref
- . 2016. Computational design of stable planar-rod structures. ACM Trans. Graph. 35, 4 (
July 2016), 1–11.DOI: Google ScholarDigital Library - . 2002. STRANDS: Interactive simulation of thin solids using cosserat models. Comput. Graph. Forum 21, 3 (
September 2002), 347–352.DOI: Google ScholarCross Ref - . 2021. Computational inverse design of surface-based inflatables. ACM Trans. Graph. 40, 4, Article
40 (July 2021), 14 pages. Google ScholarDigital Library - . 2019. X-shells: A new class of deployable beam structures. ACM Trans. Graph. 38, 4 (
July 2019), 1–15.DOI: Google ScholarDigital Library - . 2017. Computational design and automated fabrication of kirchhoff-plateau surfaces. ACM Trans. Graph. 36, 4 (
July 2017), 1–12.DOI: Google ScholarDigital Library - . 2015. Design and fabrication of flexible rod meshes. ACM Trans. Graph. 34, 4, Article
138 (July 2015), 12 pages.DOI: Google ScholarDigital Library - . 2020. On elastic geodesic grids and their planar to spatial deployment. ACM Trans. Graph. 39, 4, Article
125 (July 2020), 12 pages.DOI: Google ScholarDigital Library - . 2021. Generalized deployable elastic geodesic grids. ACM Trans. Graph. 40, 6, Article
271 (December 2021), 15 pages.DOI: Google ScholarDigital Library - . 2001. Computational Line Geometry. Springer-Verlag, Berlin and Heidelberg.Google ScholarCross Ref
- . 2021. 3D weaving with curved ribbons. ACM Trans. Graph. 40, 4, Article
127 (July 2021), 15 pages.DOI: Google ScholarDigital Library - . 2019. Soft robotics in minimally invasive surgery. Soft Robot. 6, 4 (2019), 423–443.
DOI: PMID: 30920355. Google ScholarCross Ref - . 2019. Chebyshev nets from commuting PolyVector fields. ACM Trans. Graph. 38, 6, Article
172 (November 2019), 16 pages.DOI: Google ScholarDigital Library - . 2018. Design and construction of curved support structures with repetitive parameters. In Advances in Architectural Geometry 2018. Klein Publishing, Vienna, Austria, 140–165.Google Scholar
- . 2019. Weaving geodesic foliations. ACM Trans. Graph. 38, 4, Article
34 (July 2019), 22 pages.DOI: Google ScholarDigital Library - . 2019. A-Line: 4D printing morphing linear composite structures. In Proceedings of the 2019 CHI Conference on Human Factors in Computing Systems (Glasgow, Scotland Uk) (CHI’19). Association for Computing Machinery, New York, NY, 1–12. Google ScholarDigital Library
- . 2021. Fast quasi-harmonic weights for geometric data interpolation. ACM Trans. Graph. 40, 4, Article
73 (July 2021), 15 pages.DOI: Google ScholarDigital Library - . 2019b. Bend-it: Design and fabrication of kinetic wire characters. ACM Trans. Graph. 37, 6 (
January 2019), 1–15.DOI: Google ScholarDigital Library - . 2019a. Optimization-based approach for the inverse design of ribbon-shaped three-dimensional structures assembled through compressive buckling. Phys. Rev. Appl. 11, 5 (
May 2019), 14.DOI: Google ScholarCross Ref - . 2016. Designing structurally-sound ornamental curve networks. ACM Trans. Graph. 35, 4, Article
99 (July 2016), 10 pages.DOI: Google ScholarDigital Library
Index Terms
- The Design Space of Kirchhoff Rods
Recommendations
The design space of plane elastic curves
Elastic bending of initially flat slender elements allows the realization and economic fabrication of intriguing curved shapes. In this work, we derive an intuitive but rigorous geometric characterization of the design space of plane elastic rods with ...
On the Geometric Flow of Kirchhoff Elastic Rods
Recently, rod theory has been applied to the mathematical modeling of bacterial fibers and biopolymers (e.g., DNA) to study their mechanical properties and shapes (e.g., supercoiling). In static rod theory, an elastic rod in equilibrium is the critical ...
Optimization of cross-section of actively bent grid shells with strength and geometric compatibility constraints
A method for cross-section optimization of actively bent grid shells is introduced.A co-rotational formulation allows both Form-finding and Structural analyses.The configuration of static equilibrium is sought by Dynamic Relaxation scheme.Strength and ...
Comments