Static Kirchhoff Rods under the Action of External Forces: Integration via Runge-Kutta Method

This paper shows how to apply a simple Runge-Kutta algorithm to get solutions of Kirchhoff equations for static filaments subjected to arbitrary external and static forces. This is done by suitably integrating at once Kirchhoff and filament reference system equations under appropriate initial boundary conditions. To show the application of the method, we display several numerical solutions for filaments including cases showing the effect of gravity.


Introduction
Filamentary bodies such as rods, ribbons, and wires are present in many scales from microscopic to astrophysical bodies. Marine cables [1], hair beams [2], rubber bands, microfibers [3], nanowires [4], and DNA molecules [5] share common dynamics, which is testified by the onset of rich elastic phenomena such as buckling, coiling, and looping [6]. The writhing instability for example [7,8] is a looping process in which the spatial configuration of an elastic filament changes remarkably as a response to an applied external torque, the change itself as a minimization process of the filament elastic energy.
The dynamics of filamentary bodies has been the subject of extensive research in continuum mechanics with the aim of establishing either exact or approximate solutions for thin filaments. The theory itself started with Kirchhoff [9,10] in 1892 and Clebsch and Love [11] who considered small deformation of thin rods within the elastic range. The theory was developed under special assumptions that later became the approximations adopted by the so-called simple bending and torsion theory [12]. Although the deformations of infinitesimal portions of a rod in the theory may be small, the overall deformation is quite large and allows the coiling and looping phenomena above mentioned.
Exact solutions of Kirchhoff equations for static rods are well known. Shi and Hearst found the most general analytical solution to the static Kirchhoff rod, which they called helix on a linear helix [13]. Nizzete and Goriely [14] used the equivalence between the static Kirchhoff equations for rods with circular cross-sections and the Euler equations for spinning tops to provide a classification of all different equilibrium solutions for a filament. Chouaieb et al. [15] completed the study of existence and stability of helical structures within the Kirchhoff rod model. Goriely and Tabor [7,16] developed a dynamical method to study the stability of equilibrium solutions of the Kirchhoff rod model, based on a time-dependent perturbation scheme, thus providing a framework for classification of new solutions. This method has been used, for example, to study DNA rings [17,18] and nanowires [19] under large stresses.
In biology and nanoscience, for example, most of the interesting filaments are the ones subject to diverse boundary conditions. Tobias et al. [20] derived analytic expressions for the equilibrium configurations of filaments representing long segments of a DNA double helix under different end conditions. da Fonseca et al. [19,21] recently used the Kirchhoff model to propose schemes for the study of the elastic properties of nanowires and nanosprings. McMillen and Goriely [22] obtained equilibrium solutions for the tendril perversion structure within the framework of the Kirchhoff rod model, which are heteroclinic orbits joining asymptotically two helices with opposite torsion. Gottlieb and Perkins [23] investigated spatially complex forms in 2 Journal of Computational Methods in Physics a BVP governing the equilibrium of slender cables subjected to thrust, torsion and gravity. In a very comprehensive work, Powers [24] treats the dynamics of filaments in viscous fluids and solves some interesting cases.
Despite the success of the above examples of equilibrium solutions of Kirchhoff model subject to some particular types of boundary conditions, boundary value problems are still a great challenge in the study of the static and dynamics of filaments. The main approach to find out equilibrium solutions for filaments subject to arbitrary boundary conditions is to rely on numerical computation. Examples are the use of finite element methods to obtain the equilibrium solutions of supercoiled DNAs [25], numerical package tools to solve differential equations based on a Hamiltonian formalism of the rod model [26], and a method of linearization of Kirchhoff equations, to find equilibrium solutions of open rods with fixed ends [27].
In this paper, we explore the use of straightforward methods, as the use of implicit or explicit iterative methods such as Runge-Kutta of a given order subjected to boundary conditions, to solve static filaments under the action of (static) external forces. This is done by a numerical procedure in which the Kirchhoff equations and the reference frame of the filament centerline are combined to solve the boundary value problem. The organization of the paper is as follows: we first make an introduction to the general Kirchhoff equations of rods (Section 2) and show how to use the external force to obtain a transformation matrix that couples the system of Kirchhoff equations to the filament reference frame (Sections 2.1 and 2.2). In Section 2.2 we introduce the general formulation of boundary value problems. In Section 3 we present numerical examples of closed and open rods subjected to boundary conditions, including illustrations of the effect of gravity in several kinds of filamentary solutions obtained in the absence of external forces. Section 4 is in fact a note on how to transform the filament problem into several types of boundary value problems. Finally, Section 5 (conclusion) summarizes our main results.

Kirchhoff Model for Thin Rods
In order to write the static equation for rods, we first make a short summary of the main dynamical equations of rods ( Figure 1) that are valid under the following definitions and assumptions [10]: (i) a rod or filament is a 3D body in which two of the dimensions are much smaller than a third one which defines the axis or centerline of the rod; (ii) a rod is viewed as a sequence of short segments of area which are held together by contact forces, p, between adjacent segments sharing the same elastic properties; (iii) the cross-section of a segment remains undistorted and perpendicular to the axis defined by the third dimension. The axis may be taken as the centroid of the section; (iv) displacements are always small and the breakdown of the elastic limit is never achieved. In this way, stresses and torques are always proportional to strains in full accordance with Hook's law.
With this in mind, we considered a segment of a rod or thin filament with length as shown in Figure 1. A given  can be assessed by the arclength parameter which runs from 0 to , the total rod length. The rod has crosssection which is symmetrical in relation to the centroid position defined by the vector C( , ) written in relation to a fixed frame with origin in the laboratory.
Each segment is subject to internal and external forces which, in general, are functions of . Internal forces p act upon the segment differential area and external forces f ext on the differential volume element . Newton's laws for the segment can then be written as where is the material density of the rod and R( , ) is the position of the element in the external reference system . Equations (1) and (2) are, respectively, the conservation of linear momentum and angular momentum of the rod in which all internal forces were substituted by contact stresses on the normal n of the segment and the external force act upon a volume element V.
To each segment we attach a local reference system of unitary vectors d 1 , d 2 , and d 3 , so that d 3 is parallel to the segment axis and d 1 and d 2 are orthogonal to d 3 . These are often called Frenet-Serret basis vectors [8] and are all dependent on and as the rod evolves in time. Since the last relations are twist and spin equations of the rod.
Writing k in the -basis, 1 and 2 are the curvature vectors and 3 is the twist density. Vector w is the spin vector and represents the angular momentum of the section in the local basis. It is also common to separate the element position in two terms: the position r( , ) in relation to the local reference Journal of Computational Methods in Physics 3 frame given by the -basis and the position of the segment central axis in relation to the external reference system: with As a consequence of the chosen basis, the centroid position is given by and the curve can be followed by the time dependence of the d 3 vector. Defining the total contact force integrated on the rod section as and the segment moment about the central axis as the next step is to differentiate (1) and (2) with respect to what, using (7) and (8), gives Note that (9) are valid for constant density rods which are the cases to be treated here, and they can be further simplified by using definitions (4), (5), and (6) and integrating on the section area . The result is given in terms of the local reference system, the -basis: with 1 and 2 the moment of inertia of the section in relation to the axis formed by the d 1 and d 2 vectors. To proceed further, we need a constitutive relation for the rod moment M. As usual, we write [6, 10] with and being Young's and the shear moduli of the material, respectively, and the axial moment of inertia. Equation (12) was written in such a way that if 1 = 0 and 2 = 0, the filament is naturally straight. In (10) and (11), the derivatives of F and M with respect to can be written as where we have used the summation convention. With the help of (3), we find so that the equations for the static filaments (removing the terms proportional to the derivative with respect to time in (10) and (11)) read: where ext and ext , = 1, 2, 3, are the components, inbasis, of the vectors F ext and M ext , respectively, defined by Note that the static condition does not imply that the force will be constant throughout the filament since there is an interaction among different components, arising from terms of the rotated coordinate system that is attached to each filament segment.

Integrals for the External Reference
System. If we manage to integrate (10) and (11), we obtain a solution which is still written in the segment reference system given by the -basis vectors. However, solutions are best understood if written in the external and fixed system represented by the origin as in Figure 1. In particular, the integral of the third vector d 3 gives the parametric curve of the centroid as stated by (6). It is possible to find a transformation that takes vectors from the fixed system to the local one and the corresponding inverse operation [14]. In general we can write where e 0 is a vector in Cartesian basis (̂1,̂2,̂3) of the fixed system and ( , ) is a 3 × 3 real matrix representing a sequence of pure rotations, so that det[ ( , )] = 1. Since and (3) are valid for the -basis, we have After integration of (10) and (11), Equations (20) and (21) must be integrated in order to obtain a description of the rod in the fixed system. If we restrict ourselves to static problems, first we obtain the solution for k by (15), then is obtained by (20) using (23), d( ) is directly given by (17), and, finally, the centroid curve can be generated by the integration of where (6) was used and now the vector C is given in the fixed (̂1,̂2,̂3) basis. The description of the filament is now complete: we have as initial condition the coordinate of one of the filament tips given by the vector the orientation of the centerline of the rod at the tip with respect to the fixed system given by and the -"dynamics" given by (20) in accordance to the matrix.

Coupling the Kirchhoff Equations to the Reference Frame
System. From now on we treat static rods only, whose solutions are given by solving (15) and = 0 which correspond to equilibrium solutions as described in [28]. These form a system of first order differential equations with initial value ) .
From the point of view of the initial value problem, the set of conditions (28) generates solutions for (15). However, these solutions are "unreferenced" in as much as the -evolution of the basis (d 1 , d 2 , d 3 ) must still be obtained by integrating (20). In what follows, for simplicity, we will initially take all external forces as zero, F ext = 0 and M ext = 0. If we define the vector Q with components , = [1, . . . , 6], as then (20) becomes which is the reference system problem to be solved subjected to the initial condition Rewriting the set of first order differential equations (15), Q( ) should be obtained from ) ) ) ) ) ) ) ) ) ) .
A solution of (31) is such that the matrix ( ) is orthogonal.
To show this, we note that (23) implies that + = 0, or + = 0. Using (20) we find that Since we choose (0) (0) = 1, then for all with 1 the unit matrix. We also remark that, for the -evolution in the dynamical case, the same result holds given (21) and (24). Finally, in order to get the filament curve, the following equation, derived from (6) and (26), should be solved subjected to the position of the initial filament tip C(0). Only after solving the sequence of these three sets of first order differential equations a complete solution of Kirchhoff 's problem is obtained. Such a procedure can be implemented numerically, but, instead of running each set of equations sequentially we compose the great vector with 18 components. The derivative of X with respect to will then be given by subjected to the initial condition ) .
Therefore, a single numerical procedure can be used to integrate the three sets of problems at once. We again emphasize the physical meaning of such approach. In order to get a complete filament solution, we must give the coordinate of the initial tip C(0), the direction of the three reference vectors that compose the basis, T(0), and the initial force and curvature vectors Q(0) in such a predefined basis. The Kirchhoff problem of rods as an initial value problem therefore spans a space with 18 dimensions. However, this does not mean that there are 18 independent values. In fact, since ( ) is orthogonal and represents a reference basis, one can parametrize by only 3 parameters. For example, the initial orientation of a filament tip (Figure 2) at the origin may be set by d 3 (0) as with 0 , 0 , and 0 three real numbers. The initial vector X may be written as where = 1/√ 2 0 + 2 0 . In (41) the third component of d 2 (0) was set to zero as shown in Figure 2. Note that there is some freedom in the choice of d 1 (0) and d 2 (0). These may be chosen so as to coincide with the main axis of the filament cross-section which leads to an easy visualization of the segment orientation if the cross-section is not circular. One can also use Euler angles [18], , , and , and write the initial matrix as The total number of degrees of freedom in the specification of a solution is therefore 12.

Solving the Initial Value Problem with and without External Force
Equation (38) can be straightforwardly integrated by Runge-Kutta algorithm of th order, [29] subjected to the initial condition (39). The method is a well-known technique of integration which symmetrically advances the solution step bystep using information of the current derivative such as the one given by (38). In the case of static filaments, the independent variable is the parameter , which runs from 0 to (the filament length) and plays exactly the same role of time in dynamical problems. The interaction of crossed terms, which is relevant in external forces cases, is completely taken into account through the scheme of (37). To illustrate some solutions, we use a 5th order Runge-Kutta integrator with the following initial condition [6]: and 1 = 2 . This produces Frenet helices with constant curvature and force and the solution for the centroid is shown in Figure 3 for the indicated values of . In this case, and the coil filament starts from the origin. Space periodic solutions, C(0) = C( ), can also be generated. Figure 4 contains filaments with 3000 steps with circular cross-section, In order to produce them, specially tuned boundary conditions were chosen as described in [28]. In this sense, the method is unable to produce closed solutions without further tuning of boundary values and other parameters. By changing the value of Young's modulus, the curve opens itself and becomes a complex helix as shown in Figure 4(d).
General forms of external forces may be written as with 1 = 10 + 11 + 12 2 + ⋅ ⋅ ⋅ , 2 = 20 + 21 + 22 2 + ⋅ ⋅ ⋅ , where the 's are real coefficients describing the dependence of the force as a function of parameter along the rod [10]. Gravity is an example of force for which the external torque with the gravity constant. This term must be inserted in (15), but first we need to writê3 in terms of the -basis. Using (17) we writê0 ) .
(50) Therefore, using definition (30), we write (49) as ) ) ) ) ) ) ) ) ) ) ) ) ) ) . (52) Equation (52) shows the appearance of crossed terms in the integration scheme of (38). For other external forces pointing to arbitrary directions, or under the action of external torques, the same procedure should be applied by writing whatever fixed components in terms of the -basis through (50). In order to illustrate gravity effects, we consider two Frenet helices with distinct orientations. In the first case ( Figure 5), the gravity zero helix is vertical. The parameters of Figure  These simulations were obtained with 3000 steps. We note that, as is increased, the helix distorts itself resembling a "curly-hair" filament. The helix pitch at the lower turns is smaller than the upper ones, because these must support most of the filament weight and are therefore more stretched.
Another possible solution comes from horizontal helices under the action of gravity. In Figure 7, a sequence of Frenet helix is shown with the same initial conditions of Figure 6 but with = 50.0. The orientation of d(0) is now orthogonal to the previous case so that gravity is orthogonal to the helix main axis. The natural spring is deformed toward − as intuitively expected by simple experiments with coils and ribbons.

Note about Static Filament Solutions as a Boundary Value Problem
The initial value problem defined by (38)  fully describe the variety of Cauchy's problems arising from filaments in the Kirchhoff approximation, even in its purely static cases. However, in order to give some directions on how solutions could be built, we briefly make here some very straightforward comments. As is known [29], a problem consisting of a set of first order differential equations subjected to initial conditions can be transformed to a boundary value problem with ( < ) initial conditions and − final conditions. Let us first define the initial value problem by the sequence of conditions (C , T , F , k ) in which we explicitly represent Q (29) by its constituents F and k (note that now subscripts do not represent time differentiation). Then, the following sequence of boundary value problems examples can be written as In order to solve Problem (53) for example, one can use the shooting method by searching the {k } space for a solution which minimizes the discrepancy vector The search is operated iteratively for the ( + 1)-step by the globally convergent Newton-Raphson method [29] and by giving an initial k such that where B is the value of (59) evaluated at the -step and J is the Jacobian matrix evaluated numerically for the current step. In order to provide a numerical example, we again go back to Figure 4 with = 20, Δ 1 = Δ 2 = 0.001, will result in Figure 8(a) while a magnified view is shown with Δ 1 = Δ 2 = 0.00001 in Figure 8 of matrix as stated earlier. Let us admit that , , and are the set p of three parameters. These can be the 0 , 0 , and 0 of (40) or the Euler angles. For problem (58) one can write the discrepancy vector as which can be expanded as a function T and k B (T + T , k + k ) = B (T , k ) Journal of Computational Methods in Physics

11
However, one can also write so that Imposing the condition B(T + T , k + k ) = 0 to get a solution, we find the equation system which is essentially similar to (61) with an extended Jacobian.

Conclusions
In this paper we have obtained static solutions for Kirchhoff 's equations by numerical integration, using a scheme in which both Kirchhoff and reference system equations are integrated at once. By incorporating additional dimensions, the final vector that describes the filament curve is also obtained. The scheme allows a quick and easy determination of solutions under the action of external (static) forces. These appear as interaction terms in both force and moment equations, (10), provided the components of force and torque are expressed in the local reference system. This was the case of gravity and is also the case of many other external forces such as the centrifugal acceleration for a filament in a rotating system, the electric force for a charged filament [30], or a magnetized filament in the presence of magnetic fields [31]. In order to obtain solutions, we have treated the integration of the static equations as initial and boundary value problems. As an initial value problem, the filaments can be generated as solutions of a set of first order differential equations subjected to 18 initial values. Some of these values are linked to each other (due to the parametric dependency of the matrix), so that the actual number of degrees of freedom for the initial values is 12. This sets the maximum number of dimensions in the constraints for the boundary value version of the filament problem. As a boundary value problem, several types of initial and final constraints can be written.
In deriving Kirchhoff 's equations for static filaments, we have kept the filament density constant. The effect of nonhomogeneities in rods can be easily incorporated into Runge-Kutta system as an additional term in both force and moment laws. The same is valid for wires with nonuniform elastic constants such as and . In summary, several classes of static solutions can be stud by using the approach here presented.

Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.