Abstract

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.

1. 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 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.

2. 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, , 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 segment can be assessed by the arclength parameter which runs from 0 to , the total rod length. The rod has cross-section which is symmetrical in relation to the centroid position defined by the vector 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 act upon the segment differential area and external forces 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 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 of the segment and the external force act upon a volume element .

To each segment we attach a local reference system of unitary vectors , , and , so that is parallel to the segment axis and and are orthogonal to . 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 in the -basis, and are the curvature vectors and is the twist density. Vector 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 in relation to the local reference 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 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 and the moment of inertia of the section in relation to the axis formed by the and vectors. To proceed further, we need a constitutive relation for the rod moment . 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 and , the filament is naturally straight.

In (10) and (11), the derivatives of and 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 and , , are the components, in -basis, of the vectors and , 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.

2.1. 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 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 is a vector in Cartesian basis () of the fixed system and is a real matrix representing a sequence of pure rotations, so that . Since

and (3) are valid for the -basis, we have

with

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 by (15), then is obtained by (20) using (23), 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 is given in the fixed () 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.

2.2. 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 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 must still be obtained by integrating (20). In what follows, for simplicity, we will initially take all external forces as zero, and . If we define the vector with components , , 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), should be obtained from

A solution of (31) is such that the matrix is orthogonal. To show this, we note that (23) implies that , or . Using (20) we find that

Since we choose , then

for all with 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 . 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 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 , the direction of the three reference vectors that compose the basis, , and the initial force and curvature vectors 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 as

with , , and three real numbers. The initial vector may be written as

with

where . In (41) the third component of was set to zero as shown in Figure 2. Note that there is some freedom in the choice of and . 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

with

The total number of degrees of freedom in the specification of a solution is therefore 12.

3. 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 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 . 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 integration was made with points. As the value of is increased, the Frenet helix is distended indicating an increase in the force along the -axis (vertical). The orientation of the initial local basis is

and the coil filament starts from the origin. Space periodic solutions, , can also be generated. Figure 4 contains filaments with steps with circular cross-section, and the same initial orientation. The parameters used to generate each curve are as follows: (a) , , , , , , , , , , (b) , , , , , , , , , , (c) , , , , , , , , , , and (d) , , , , , , , , , . (a) and (b) are almost closed solutions. 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).

Since in the last examples, , then for all according to the last element of (38). Young's modulus multiplies and as shown by (15), so that the cases for which either represent homogeneous filaments with noncircular cross-sections or circular cross-section filaments with anisotropic bending stiffness. To illustrate this case, we show in Figure 4 a helix with and and decreasing values of , the shear modulus, as indicated in the figure. For this case we have , , , , , , , and points. As the shear modulus decreases, the filament stretches as a response to less torsion rigidity.

General forms of external forces may be written as

with 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 is zero, since gravity can not make the segment turn about its symmetric axis. For gravity all 's are zero except , , and . We will restrict ourselves to gravity in the numerical examples that follow. For the segment, the evaluation of the first integral in (16) leads to

with the gravity constant. This term must be inserted in (15), but first we need to write in terms of the -basis. Using (17) we write

Therefore, using definition (30), we write (49) as

Equation (33) then becomes

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 5 are as follows: , , , , , , , , , , , , and is varied as indicated in the 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 . The orientation of 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.

4. Note about Static Filament Solutions as a Boundary Value Problem

The initial value problem defined by (38) and (39) can be cast into a variety of boundary value problems. Due to space limitations, it is not in the scope of this paper to 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 in which we explicitly represent (29) by its constituents and (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 space for a solution which minimizes the discrepancy vector

The search is operated iteratively for the -step by the globally convergent Newton-Raphson method [29] and by giving an initial such that where is the value of (59) evaluated at the -step and is the Jacobian matrix

evaluated numerically for the current step. In order to provide a numerical example, we again go back to Figure 4(a) and consider the problem of “closing” a single loop. We fix and let and , , , , , , , . Instead of calculating , we determine the surface generated by and inspect the topology of on the plane since is constant. The following set of points on this plane with , , will result in Figure 8(a) while a magnified view is shown with in Figure 8(b). These plots exhibit a clear minimum corresponding to the pair that minimizes .

Problems (54)–(58) can be solved by the same reasoning. However, for problems (57) and (58) the actual number of search variable is six because of the parametric dependency of matrix as stated earlier. Let us admit that , , and are the set of three parameters. These can be the , , and of (40) or the Euler angles. For problem (58) one can write the discrepancy vector as

which can be expanded as a function and

However, one can also write

so that

Imposing the condition to get a solution, we find the equation system

which is essentially similar to (61) with an extended Jacobian.

5. 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.

Acknowledgments

The author would like to thank Christine Fernandes for the preparation of the LaTeX version of this paper, Marcus A. M. de Aguiar and Alexandre F. da Fonseca (IFGW/UNICAMP) for useful discussions.