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

Authors: Ademir L. Xavier Jr.
Journal: Journal of Computational Methods in Physics
 (2014)

---

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


---

## Body

## 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 (Section2) 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 (Figure1) 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 A 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.Figure 1
Illustration of a rod segment with cross-section areaA actuated by external and contact forces.With this in mind, we considered a segment of a rod or thin filament with lengthd
s as shown in Figure 1. A given segment can be assessed by the arclength parameter s which runs from 0 to L, the total rod length. The rod has cross-section A which is symmetrical in relation to the centroid position defined by the vector C
(
s
,
t
) written in relation to a fixed frame with origin O in the laboratory.Each segment is subject to internal and external forces which, in general, are functions ofs. Internal forces p act upon the segment differential area d
A and external forces f
ext on the differential volume element d
A
d
s. Newton's laws for the segment can then be written as
(1)
∫
‍
∫
A
‍
p
n
d
A
+
∫
‍
∫
‍
∫
V
‍
f
ext
d
v
=
∫
‍
∫
‍
∫
V
‍
ρ
∂
2
R
∂
t
2
d
v
,
(2)
∫
‍
∫
A
‍
R
×
p
n
d
A
+
∫
‍
∫
‍
∫
V
‍
R
×
f
ext
d
v
=
∫
‍
∫
‍
∫
V
‍
ρ
R
×
∂
2
R
∂
t
2
d
v
,
where ρ is the material density of the rod and R
(
s
,
t
) is the position of the element in the external reference system O. 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
d
A of the segment and the external force act upon a volume element d
v.To each segment we attach a local reference system of unitary vectorsd
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 s and t as the rod evolves in time. Since
(3)
∂
∂
s
(
d
i
·
d
j
)
=
0
⟶
∂
d
i
∂
s
=
k
×
d
i
,
∂
∂
t
(
d
i
·
d
j
)
=
0
⟶
∂
d
i
∂
t
=
w
×
d
i
,the last relations are twist and spin equations of the rod. Writingk in the d-basis, k
1 and k
2 are the curvature vectors and k
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 positionr
(
s
,
t
) in relation to the local reference frame given by the d-basis and the position of the segment central axis in relation to the external reference system:
(4)
R
(
s
,
t
)
=
C
(
s
,
t
)
+
r
(
s
,
t
)with(5)
r
(
s
,
t
)
=
r
1
d
1
+
r
2
d
2
.As a consequence of the chosen basis, the centroid position is given by(6)
C
(
s
,
t
)
=
∫
0
s
‍
d
3
(
s
′
,
t
)
d
s
′
,and the curve can be followed by the time dependence of thed
3 vector.Defining the total contact force integrated on the rod section as(7)
F
≡
∫
‍
∫
A
‍
p
n
d
aand the segment moment about the central axis as(8)
M
≡
∫
‍
∫
A
‍
r
×
p
n
d
a
,the next step is to differentiate (1) and (2) with respect to s what, using (7) and (8), gives
(9)
∂
F
∂
s
+
∫
‍
∫
A
‍
f
ext
d
a
=
∫
‍
∫
A
‍
ρ
∂
2
R
∂
t
2
d
a
,
∂
M
∂
s
+
∂
C
∂
s
×
F
+
∫
‍
∫
A
‍
r
×
f
ext
d
a
=
∫
‍
∫
‍
ρ
(
r
×
∂
2
R
∂
t
2
)
d
a
.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 A. The result is given in terms of the local reference system, the d-basis:
(10)
∂
F
∂
s
+
∫
‍
∫
A
‍
f
ext
d
a
=
ρ
A
∫
0
s
‍
∂
2
d
3
∂
t
2
d
s
′
,
(11)
∂
M
∂
s
+
d
3
×
F
+
∫
‍
∫
A
‍
r
×
f
ext
d
a
=
ρ
(
I
1
d
1
×
∂
d
1
∂
t
+
I
2
d
2
×
∂
d
2
∂
t
)
,withI
1 and I
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]
(12)
M
=
E
I
1
k
1
d
1
+
E
I
2
k
2
d
2
+
G
J
k
3
d
3
=
(
E
I
1
k
1
E
I
2
k
2
G
J
k
3
)
,withE and G being Young's and the shear moduli of the material, respectively, and J the axial moment of inertia. Equation (12) was written in such a way that if k
1
=
0 and k
2
=
0, the filament is naturally straight.In (10) and (11), the derivatives of F and M with respect to s can be written as
(13)
∂
F
∂
s
=
∂
F
i
∂
s
d
i
+
F
i
∂
d
i
∂
s
,
∂
M
∂
s
=
∂
M
i
∂
s
d
i
+
M
i
∂
d
i
∂
s
,
where we have used the summation convention. With the help of (3), we find
(14)
∂
F
∂
s
=
(
∂
F
1
∂
s
-
k
3
F
2
+
k
2
F
3
∂
F
2
∂
s
+
k
3
F
1
-
k
1
F
3
∂
F
3
∂
s
-
k
2
F
1
+
k
1
F
2
)
,
∂
M
∂
s
=
(
E
I
1
∂
k
1
∂
s
-
(
E
I
2
-
G
J
)
k
2
k
3
E
I
2
∂
k
2
∂
s
+
(
E
I
1
-
G
J
)
k
1
k
3
G
J
∂
k
3
∂
s
-
E
(
I
1
-
I
2
)
k
1
k
2
)
,so that the equations for the static filaments (removing the terms proportional to the derivative with respect to time in (10) and  (11)) read:
(15)
d
F
1
d
s
-
k
3
F
2
+
k
2
F
3
+
F
1
ext
=
0
,
d
F
2
d
s
+
k
3
F
1
-
k
1
F
3
+
F
2
ext
=
0
,
d
F
3
d
s
-
k
2
F
1
+
k
1
F
2
+
F
3
ext
=
0
,
E
I
1
d
k
1
d
s
-
(
E
I
2
-
G
J
)
k
2
k
3
-
F
2
+
M
1
ext
=
0
,
E
I
2
d
k
2
d
s
+
(
E
I
1
-
G
J
)
k
1
k
3
+
F
1
+
M
2
ext
=
0
,
G
J
d
k
3
d
s
-
E
(
I
1
-
I
2
)
k
1
k
2
+
M
3
ext
=
0
,whereF
i
ext and M
i
ext, i
=
1,2
,
3, are the components, in d-basis, of the vectors F
ext and M
ext, respectively, defined by
(16)
F
ext
=
∫
‍
∫
A
‍
f
ext
d
a
,
M
ext
=
∫
‍
∫
A
‍
r
×
f
ext
d
a
.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 d-basis vectors. However, solutions are best understood if written in the external and fixed system represented by the origin O 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
(17)
d
(
s
,
t
)
=
S
(
s
,
t
)
e
0
,
where e
0 is a vector in Cartesian basis (e
^
1
,
e
^
2
,
e
^
3) of the fixed system and S
(
s
,
t
) is a 3
×
3 real matrix representing a sequence of pure rotations, so that det
[
S
(
s
,
t
)
]
=
1. Since
(18)
∂
d
(
s
,
t
)
∂
s
=
∂
S
(
s
,
t
)
∂
s
e
0
,
(19)
∂
d
(
s
,
t
)
∂
t
=
∂
S
(
s
,
t
)
∂
t
e
0
,and (3) are valid for the d-basis, we have
(20)
∂
S
∂
s
=
K
S
,
(21)
∂
S
∂
t
=
W
S
,
(22)
∂
S
∂
t
=
W
K
-
1
∂
S
∂
s
,with(23)
K
=
(
0
k
3
-
k
2
-
k
3
0
k
1
k
2
-
k
1
0
)
,
(24)
W
=
(
0
w
3
-
w
2
-
w
3
0
w
1
w
2
-
w
1
0
)
.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 fork by (15), then S is obtained by (20) using (23), d
(
s
) is directly given by (17), and, finally, the centroid curve can be generated by the integration of
(25)
d
C
d
s
=
(
S
31
(
s
)
S
32
(
s
)
S
33
(
s
)
)
,
where (6) was used and now the vector C is given in the fixed (e
^
1
,
e
^
2
,
e
^
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
(26)
C
(
0
)
=
(
C
1
(
0
)
C
2
(
0
)
C
3
(
0
)
)
,the orientation of the centerline of the rod at the tip with respect to the fixed system given by(27)
d
(
0
)
=
S
(
0
)
e
0
,and thes-“dynamics” given by (20) in accordance to the K 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 W
=
0 which correspond to equilibrium solutions as described in [28]. These form a system of first order differential equations with initial value
(28)
F
(
0
)
=
(
F
1
(
0
)
F
2
(
0
)
F
3
(
0
)
)
,
k
(
0
)
=
(
k
1
(
0
)
k
2
(
0
)
k
3
(
0
)
)
.From the point of view of theinitial value problem, the set of conditions (28) generates solutions for (15). However, these solutions are “unreferenced” in as much as the s-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 Q
i, i
=
[
1
,
…
,
6
], as
(29)
Q
t
=
(
F
1
F
2
F
3
k
1
k
2
k
3
)
,
(30)
T
t
=
(
S
11
S
12
S
13
S
21
S
22
S
23
S
31
S
32
S
33
)
,then (20) becomes
(31)
(
d
T
d
s
)
=
(
Q
6
(
s
)
T
4
-
Q
5
(
s
)
T
7
Q
6
(
s
)
T
5
-
Q
5
(
s
)
T
8
Q
6
(
s
)
T
6
-
Q
5
(
s
)
T
9
Q
4
(
s
)
T
7
-
Q
6
(
s
)
T
1
Q
4
(
s
)
T
8
-
Q
6
(
s
)
T
2
Q
4
(
s
)
T
9
-
Q
6
(
s
)
T
3
Q
5
(
s
)
T
1
-
Q
4
(
s
)
T
4
Q
5
(
s
)
T
2
-
Q
4
(
s
)
T
5
Q
5
(
s
)
T
3
-
Q
4
(
s
)
T
6
)
,which is the reference system problem to be solved subjected to the initial condition(32)
T
(
0
)
t
=
(
S
11
(
0
)
S
12
(
0
)
S
13
(
0
)
S
21
(
0
)
S
22
(
0
)
S
23
(
0
)
S
31
(
0
)
S
32
(
0
)
S
33
(
0
)
)
.Rewriting the set of first order differential equations (15), Q
(
s
) should be obtained from(33)
(
d
Q
d
s
)
=
(
Q
6
Q
2
-
Q
5
Q
3
Q
4
Q
3
-
Q
6
Q
1
Q
5
Q
1
-
Q
4
Q
2
[
(
E
I
2
-
G
J
)
Q
5
Q
6
+
Q
2
]
E
I
1
-
[
(
E
I
1
-
G
J
)
Q
4
Q
6
+
Q
1
]
E
I
2
(
I
1
-
I
2
)
Q
4
Q
5
G
J
)
.A solution of (31) is such that the matrix S
(
s
) is orthogonal. To show this, we note that (23) implies that K
t
+
K
=
0, or S
t
K
t
S
+
S
t
K
S
=
0. Using   (20) we find that
(34)
d
d
s
(
S
t
S
)
=
0
.Since we chooseS
t
(
0
)
S
(
0
)
=
1, then
(35)
S
t
S
=
1for alls with 1 the unit matrix. We also remark that, for the t-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),
(36)
d
C
d
s
=
(
T
7
(
s
)
T
8
(
s
)
T
9
(
s
)
)
,should be solved subjected to the position of the initial filament tipC
(
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
(37)
X
=
(
Q
T
C
)
,with 18 components. The derivative ofX with respect to s will then be given by
(38)
d
X
d
s
=
(
d
Q
d
s
d
T
d
s
d
C
d
s
)subjected to the initial condition(39)
X
(
0
)
=
(
Q
(
0
)
T
(
0
)
C
(
0
)
)
.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 tipC
(
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 S
(
s
) is orthogonal and represents a reference basis, one can parametrize S 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
(40)
d
3
(
0
)
=
(
x
0
y
0
z
0
)
,
x
0
2
+
y
0
2
+
z
0
2
=
1Figure 2
Orientation of thed-basis vectors and the fixed reference frame. Vector d
3 is always tangent to the filament curve.withx
0, y
0, and z
0 three real numbers. The initial vector X may be written as
(41)
X
(
0
)
=
(
F
(
0
)
k
(
0
)
U
1
U
2
d
3
(
0
)
0
)with(42)
U
1
t
=
(
γ
x
0
z
0
γ
y
0
z
0
-
γ
)
,
U
2
t
=
(
-
γ
y
0
γ
x
0
0
)
,whereγ
=
1
/
x
0
2
+
y
0
2. 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 S matrix as
(43)
S
=
(
D
11
(
θ
,
ϕ
)
D
12
(
θ
,
ϕ
)
D
13
(
θ
,
ϕ
)
D
21
(
θ
,
ϕ
)
D
22
(
θ
,
ϕ
)
D
23
(
θ
,
ϕ
)
D
31
(
θ
,
ϕ
)
D
32
(
θ
,
ϕ
)
D
33
(
θ
)
)with(44)
D
11
(
θ
,
ϕ
)
=
cos
θ
cos
φ
cos
ψ
-
sin
φ
sin
ψ
,
D
12
(
θ
,
ϕ
)
=
cos
θ
cos
φ
sin
ψ
,
D
13
(
θ
,
ϕ
)
=
-
cos
φ
sin
θ
,
D
21
(
θ
,
ϕ
)
=
-
cos
θ
sin
φ
cos
ψ
-
cos
φ
sin
ψ
,
D
22
(
θ
,
ϕ
)
=
-
cos
θ
sin
φ
sin
ψ
+
cos
φ
cos
ψ
,
D
23
(
θ
,
ϕ
)
=
sin
φ
sin
θ
,
D
31
(
θ
,
ϕ
)
=
sin
θ
cos
φ
,
D
32
(
θ
,
ϕ
)
=
sin
θ
sin
ψ
,
D
33
(
θ
)
=
cos
θ
.The total number of degrees of freedom in the specification of a solution is therefore 12.


## 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 d-basis vectors. However, solutions are best understood if written in the external and fixed system represented by the origin O 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
(17)
d
(
s
,
t
)
=
S
(
s
,
t
)
e
0
,
where e
0 is a vector in Cartesian basis (e
^
1
,
e
^
2
,
e
^
3) of the fixed system and S
(
s
,
t
) is a 3
×
3 real matrix representing a sequence of pure rotations, so that det
[
S
(
s
,
t
)
]
=
1. Since
(18)
∂
d
(
s
,
t
)
∂
s
=
∂
S
(
s
,
t
)
∂
s
e
0
,
(19)
∂
d
(
s
,
t
)
∂
t
=
∂
S
(
s
,
t
)
∂
t
e
0
,and (3) are valid for the d-basis, we have
(20)
∂
S
∂
s
=
K
S
,
(21)
∂
S
∂
t
=
W
S
,
(22)
∂
S
∂
t
=
W
K
-
1
∂
S
∂
s
,with(23)
K
=
(
0
k
3
-
k
2
-
k
3
0
k
1
k
2
-
k
1
0
)
,
(24)
W
=
(
0
w
3
-
w
2
-
w
3
0
w
1
w
2
-
w
1
0
)
.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 fork by (15), then S is obtained by (20) using (23), d
(
s
) is directly given by (17), and, finally, the centroid curve can be generated by the integration of
(25)
d
C
d
s
=
(
S
31
(
s
)
S
32
(
s
)
S
33
(
s
)
)
,
where (6) was used and now the vector C is given in the fixed (e
^
1
,
e
^
2
,
e
^
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
(26)
C
(
0
)
=
(
C
1
(
0
)
C
2
(
0
)
C
3
(
0
)
)
,the orientation of the centerline of the rod at the tip with respect to the fixed system given by(27)
d
(
0
)
=
S
(
0
)
e
0
,and thes-“dynamics” given by (20) in accordance to the K 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 W
=
0 which correspond to equilibrium solutions as described in [28]. These form a system of first order differential equations with initial value
(28)
F
(
0
)
=
(
F
1
(
0
)
F
2
(
0
)
F
3
(
0
)
)
,
k
(
0
)
=
(
k
1
(
0
)
k
2
(
0
)
k
3
(
0
)
)
.From the point of view of theinitial value problem, the set of conditions (28) generates solutions for (15). However, these solutions are “unreferenced” in as much as the s-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 Q
i, i
=
[
1
,
…
,
6
], as
(29)
Q
t
=
(
F
1
F
2
F
3
k
1
k
2
k
3
)
,
(30)
T
t
=
(
S
11
S
12
S
13
S
21
S
22
S
23
S
31
S
32
S
33
)
,then (20) becomes
(31)
(
d
T
d
s
)
=
(
Q
6
(
s
)
T
4
-
Q
5
(
s
)
T
7
Q
6
(
s
)
T
5
-
Q
5
(
s
)
T
8
Q
6
(
s
)
T
6
-
Q
5
(
s
)
T
9
Q
4
(
s
)
T
7
-
Q
6
(
s
)
T
1
Q
4
(
s
)
T
8
-
Q
6
(
s
)
T
2
Q
4
(
s
)
T
9
-
Q
6
(
s
)
T
3
Q
5
(
s
)
T
1
-
Q
4
(
s
)
T
4
Q
5
(
s
)
T
2
-
Q
4
(
s
)
T
5
Q
5
(
s
)
T
3
-
Q
4
(
s
)
T
6
)
,which is the reference system problem to be solved subjected to the initial condition(32)
T
(
0
)
t
=
(
S
11
(
0
)
S
12
(
0
)
S
13
(
0
)
S
21
(
0
)
S
22
(
0
)
S
23
(
0
)
S
31
(
0
)
S
32
(
0
)
S
33
(
0
)
)
.Rewriting the set of first order differential equations (15), Q
(
s
) should be obtained from(33)
(
d
Q
d
s
)
=
(
Q
6
Q
2
-
Q
5
Q
3
Q
4
Q
3
-
Q
6
Q
1
Q
5
Q
1
-
Q
4
Q
2
[
(
E
I
2
-
G
J
)
Q
5
Q
6
+
Q
2
]
E
I
1
-
[
(
E
I
1
-
G
J
)
Q
4
Q
6
+
Q
1
]
E
I
2
(
I
1
-
I
2
)
Q
4
Q
5
G
J
)
.A solution of (31) is such that the matrix S
(
s
) is orthogonal. To show this, we note that (23) implies that K
t
+
K
=
0, or S
t
K
t
S
+
S
t
K
S
=
0. Using   (20) we find that
(34)
d
d
s
(
S
t
S
)
=
0
.Since we chooseS
t
(
0
)
S
(
0
)
=
1, then
(35)
S
t
S
=
1for alls with 1 the unit matrix. We also remark that, for the t-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),
(36)
d
C
d
s
=
(
T
7
(
s
)
T
8
(
s
)
T
9
(
s
)
)
,should be solved subjected to the position of the initial filament tipC
(
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
(37)
X
=
(
Q
T
C
)
,with 18 components. The derivative ofX with respect to s will then be given by
(38)
d
X
d
s
=
(
d
Q
d
s
d
T
d
s
d
C
d
s
)subjected to the initial condition(39)
X
(
0
)
=
(
Q
(
0
)
T
(
0
)
C
(
0
)
)
.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 tipC
(
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 S
(
s
) is orthogonal and represents a reference basis, one can parametrize S 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
(40)
d
3
(
0
)
=
(
x
0
y
0
z
0
)
,
x
0
2
+
y
0
2
+
z
0
2
=
1Figure 2
Orientation of thed-basis vectors and the fixed reference frame. Vector d
3 is always tangent to the filament curve.withx
0, y
0, and z
0 three real numbers. The initial vector X may be written as
(41)
X
(
0
)
=
(
F
(
0
)
k
(
0
)
U
1
U
2
d
3
(
0
)
0
)with(42)
U
1
t
=
(
γ
x
0
z
0
γ
y
0
z
0
-
γ
)
,
U
2
t
=
(
-
γ
y
0
γ
x
0
0
)
,whereγ
=
1
/
x
0
2
+
y
0
2. 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 S matrix as
(43)
S
=
(
D
11
(
θ
,
ϕ
)
D
12
(
θ
,
ϕ
)
D
13
(
θ
,
ϕ
)
D
21
(
θ
,
ϕ
)
D
22
(
θ
,
ϕ
)
D
23
(
θ
,
ϕ
)
D
31
(
θ
,
ϕ
)
D
32
(
θ
,
ϕ
)
D
33
(
θ
)
)with(44)
D
11
(
θ
,
ϕ
)
=
cos
θ
cos
φ
cos
ψ
-
sin
φ
sin
ψ
,
D
12
(
θ
,
ϕ
)
=
cos
θ
cos
φ
sin
ψ
,
D
13
(
θ
,
ϕ
)
=
-
cos
φ
sin
θ
,
D
21
(
θ
,
ϕ
)
=
-
cos
θ
sin
φ
cos
ψ
-
cos
φ
sin
ψ
,
D
22
(
θ
,
ϕ
)
=
-
cos
θ
sin
φ
sin
ψ
+
cos
φ
cos
ψ
,
D
23
(
θ
,
ϕ
)
=
sin
φ
sin
θ
,
D
31
(
θ
,
ϕ
)
=
sin
θ
cos
φ
,
D
32
(
θ
,
ϕ
)
=
sin
θ
sin
ψ
,
D
33
(
θ
)
=
cos
θ
.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 nth 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 s, which runs from 0 to L (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]:
(45)
Q
(
0
)
=
(
0
κ
τ
0
κ
τ
(
G
J
-
E
I
1
)
E
I
1
τ
2
(
G
J
-
E
I
1
)
E
I
1
)
,andI
1
=
I
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, I
1
=
I
2
=
2.0, G
=
1.0, E
=
0.5, κ
=
1.0, L
=
60.0, k
1
(
0
)
=
0, k
2
(
0
)
=
1, k
3
(
0
)
=
0.03, F
1
(
0
)
=
0, F
2
(
0
)
=
0.09, F
3
(
0
)
=
2.7
×
1
0
-
3, and the integration was made with 3000 points. As the value of τ is increased, the Frenet helix is distended indicating an increase in the force along the z-axis (vertical). The orientation of the initial local basis is
(46)
(
d
1
(
0
)
d
2
(
0
)
d
3
(
0
)
)
=
(
-
1
0
0
0
0
1
0
1
0
)
(
e
1
e
2
e
3
)
,Figure 3
Frenet helices for circular cross section filaments and increasing values oft. As usual, z is the vertical axis of the box.and the coil filament starts from the origin.  Space periodic solutions,C
(
0
)
=
C
(
L
), can also be generated. Figure 4 contains filaments with 3000 steps with circular cross-section, I
1
=
I
2 and the same initial orientation. The parameters used to generate each curve are as follows: (a) E
=
0.2, G
=
0.5, k
1
(
0
)
=
0.0, k
2
(
0
)
=
0.1, k
3
(
0
)
=
0.1, F
1
(
0
)
=
0.0, F
2
(
0
)
=
0.01, F
3
(
0
)
=
0.0, I
1
=
1.0, L
=
63.0, (b) E
=
1.0, G
=
2.5, k
1
(
0
)
=
0.0, k
2
(
0
)
=
1.3, k
3
(
0
)
=
0.01, F
1
(
0
)
=
0.0, F
2
(
0
)
=
0.9885, F
3
(
0
)
=
1.0, I
1
=
2.0, L
=
71.0, (c) E
=
0.1, G
=
2.5, k
1
(
0
)
=
0.0, k
2
(
0
)
=
1.3, k
3
(
0
)
=
0.01, F
1
(
0
)
=
0.0, F
2
(
0
)
=
0.9885, F
3
(
0
)
=
1.0, I
1
=
2.0, L
=
13.0, and (d) E
=
0.2, G
=
2.5, k
1
(
0
)
=
0.0, k
2
(
0
)
=
1.3, k
3
(
0
)
=
0.01, F
1
(
0
)
=
0.0, F
2
(
0
)
=
0.99, F
3
(
0
)
=
2.0, I
1
=
2.0, L
=
35.0. (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).A variety of solutions for a circular cross section filament and different values of initial conditions. (a) and (b) are closed solutions. See text for the parameters of each figure;z is the vertical direction.
(a)
(b)
(c)
(d)Since in the last examples,I
1
=
I
2, then k
3
(
s
)
=
k
3
(
0
) for all s according to the last element of  (38). Young's modulus E multiplies I
1 and I
2 as shown by  (15), so that the cases for which I
1
≠
I
2 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 I
1
=
3.0 and I
1
=
2.0 and decreasing values of G, the shear modulus, as indicated in the figure. For this case we have E
=
1.0, k
1
(
0
)
=
0.0, k
2
(
0
)
=
1.0, k
3
(
0
)
=
0.07, F
1
(
0
)
=
0.0, F
2
(
0
)
=
0.21, F
3
(
0
)
=
0.0147, L
=
63.0 and 3000 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(47)
f
ext
=
f
1
d
1
+
f
2
d
2
+
f
3
d
3
,with(48)
f
1
=
α
10
+
α
11
s
+
α
12
s
2
+
⋯
,
f
2
=
α
20
+
α
21
s
+
α
22
s
2
+
⋯
,
f
3
=
α
30
+
α
31
s
+
α
32
s
2
+
⋯
,
where the α's are real coefficients describing the dependence of the force as a function of parameter s 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 α
10, α
20, and α
30. 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
(49)
F
ext
=
-
ρ
g
A
e
^
3withg the gravity constant. This term must be inserted in   (15), but first we need to write e
^
3 in terms of the d-basis. Using   (17) we write
(50)
e
^
0
=
S
-
1
(
s
)
d
(
s
)
,
S
-
1
(
s
)
=
(
S
23
S
32
-
S
22
S
33
S
12
S
33
-
S
13
S
32
S
13
S
22
-
S
12
S
23
S
33
S
21
-
S
23
S
31
S
13
S
31
-
S
33
S
11
S
23
S
11
-
S
13
S
21
S
22
S
31
-
S
21
S
32
S
32
S
11
-
S
31
S
12
S
12
S
21
-
S
22
S
11
)
.Therefore, using definition (30), we write (49) as
(51)
F
ext
=
-
ρ
g
A
[
(
T
5
T
4
-
T
8
T
4
)
d
1
+
(
T
8
T
1
-
T
7
T
2
)
d
2
F
ext
=
-
ρ
g
A
+
(
T
2
T
3
-
T
5
T
1
)
d
3
]
.Equation (33) then becomes
(52)
(
d
Q
d
s
)
=
(
Q
6
Q
2
-
Q
5
Q
3
+
ρ
g
A
(
T
5
T
4
-
T
8
T
4
)
Q
4
Q
3
-
Q
6
Q
1
+
ρ
g
A
(
T
8
T
1
-
T
7
T
2
)
Q
5
Q
1
-
Q
4
Q
2
+
ρ
g
A
(
T
2
T
3
-
T
5
T
1
)
[
(
E
I
2
-
G
J
)
Q
5
Q
6
+
Q
2
]
E
I
1
[
(
-
E
I
1
+
G
J
)
Q
4
Q
6
-
Q
1
]
E
I
2
E
(
I
1
-
I
2
)
Q
4
Q
5
G
J
)
.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 d-basis through   (50).In order to illustrate gravity effects, we consider two Frenet helices with distinct orientations. In the first case (Figure5), the gravity zero helix is vertical. The parameters of Figure 5 are as follows: E
=
1.0, G
=
0.5, k
1
(
0
)
=
0.0, k
2
(
0
)
=
2.0, k
3
(
0
)
=
0.1, F
1
(
0
)
=
0.0, F
2
(
0
)
=
0.0, F
3
(
0
)
=
0.0, I
1
=
1.0, I
2
=
1.0, L
=
30.0, r
=
1.0, A
=
1.0 and g is varied as indicated in the figure. These simulations were obtained with 3000 steps. We note that, as g 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.Figure 5
Effect of anisotropy in a filament with decreasing shear modulus. The originalI
1
=
I
2 case is a classic Frenet helix with g
=
1.0. For all cases above I
1
=
3.0 and I
2
=
2.0. As g decreases the filament becomes less rigid against torsion. The parameters used to generate each curve are in the text; z is the vertical direction.Another possible solution comes from horizontal helices under the action of gravity. In Figure7, a sequence of Frenet helix is shown with the same initial conditions of Figure 6 but with L
=
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 -
z as intuitively expected by simple experiments with coils and ribbons.Figure 6
“Curly-hair” solutions or Frenet helices under the action of gravity for increasing values ofg, stretching the filament vertically. Gravity points downward in the -
z direction. The parameters used to generate each curve are in the text.Figure 7
Horizontal Frenet helix under the action of gravity. The parameters are the same as in Figure6, with exception of d
(
0
) which is orthogonal to the one used in previous examples.


## 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 N first order differential equations subjected to N initial conditions can be transformed to a boundary value problem with n (n
<
N) initial conditions and N
-
n  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
(53)
(
C
′
,
T
′
,
F
′
)
⟶
(
C
′′
)
,
(54)
(
C
′
,
T
′
,
k
′
)
⟶
(
C
′′
)
,
(55)
(
C
′
,
T
′
)
⟶
(
C
′′
,
F
′′
)
,
(56)
(
C
′
,
T
′
)
⟶
(
C
′′
,
k
′′
)
,
(57)
(
C
′
,
F
′
,
k
′
)
⟶
(
C
′′
)
,
(58)
(
C
′
,
F
′
)
⟶
(
C
′′
,
F
′′
)
,
…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
(59)
B
(
k
′
)
|
C
′
,
T
′
,
F
′
=
C
(
k
′
,
s
=
L
)
-
C
′′
.The search is operated iteratively for the(
i
+
1
)-step by the globally convergent Newton-Raphson method [29] and by giving an initial k
i
′
(60)
k
i
+
1
'
⟵
k
i
'
+
δ
k
i
'
,such that
(61)
J
i
·
δ
k
i
'
=
-
B
i
,
where B
i is the value of   (59) evaluated at the i-step and J is the Jacobian matrix
(62)
J
i
j
=
∂
B
i
∂
k
jevaluated numerically for the current step. In order to provide a numerical example, we again go back to Figure4(a) and consider the problem of “closing” a single loop. We fix L
=
62.835 and let I
1
=
I
2
=
1.0 and E
=
0.2, G
=
0.5, k
1
(
0
)
=
0.0, k
2
(
0
)
=
0.1, k
3
(
0
)
=
0.1, F
1
(
0
)
=
0.0, F
2
(
0
)
=
0.01, F
3
(
0
)
=
0.0. Instead of calculating ∂
B
i
/
∂
k
j, we determine the surface generated by |
B
(
k
′
)
| and inspect the topology of ∂
B
i
/
∂
k
j on the (
k
1
′
,
k
2
′
) plane since k
3 is constant. The following set of points on this plane
(63)
k
¯
1
(
i
)
=
k
1
(
0
)
+
Δ
k
1
(
2
i
-
m
)
2
m
,
i
=
0
,
…
,
m
,
k
¯
2
(
i
)
=
k
2
(
0
)
+
Δ
k
2
(
2
j
-
m
)
2
m
,
j
=
0
,
…
,
m
,
with m
=
20, Δ
k
1
=
Δ
k
2
=
0.001, will result in Figure 8(a) while a magnified view is shown with Δ
k
1
=
Δ
k
2
=
0.00001 in Figure 8(b). These plots exhibit a clear minimum corresponding to the (
k
1
′
,
k
2
′
) pair that minimizes |
B
(
k
′
)
|.Plots oflog
|
B
(
k
¯
1
(
i
)
,
k
¯
2
(
j
)
)
| for 20
×
20 solutions of Kirchoff loops around k
1
(
i
=
10
)
=
0.0, k
2
(
j
=
10
)
=
0.1 with i and j given by   (63). In (a), Δ
k
1
=
Δ
k
2
=
0.001 and in (b) Δ
k
1
=
Δ
k
2
=
0.00001.
(a)
(b)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 S as stated earlier. Let us admit that α, β, and γ are the set p of three parameters. These can be the x
0, y
0, and z
0 of   (40) or the Euler angles. For problem (58) one can write the discrepancy vector as
(64)
B
(
T
′
,
k
′
)
=
(
C
(
s
=
L
)
-
C
′′
F
(
s
=
L
)
-
F
′′
)
,which can be expanded as a functionT
′ and k
′
(65)
B
(
T
′
+
δ
T
′
,
k
′
+
δ
k
′
)
=
B
(
T
′
,
k
′
)
B
(
T
′
+
δ
T
′
,
k
′
+
δ
k
′
)
=
+
∑
i
=
1
9
‍
∂
B
∂
T
i
′
∂
T
i
′
+
∑
i
=
1
3
‍
∂
B
∂
k
i
′
∂
k
i
′
.However, one can also write(66)
∂
T
i
′
=
∑
j
=
1
3
‍
∂
T
i
′
∂
p
j
δ
p
j
,so that(67)
B
(
T
′
+
δ
T
′
,
k
′
+
δ
k
′
)
=
B
(
T
′
,
k
′
)
+
∑
j
=
1
3
‍
∑
i
=
1
9
‍
∂
B
∂
T
i
′
∂
T
i
′
∂
p
j
δ
p
j
+
∑
i
=
1
3
‍
∂
B
∂
k
i
′
δ
k
i
′
.Imposing the conditionB
(
T
′
+
δ
T
′
,
k
′
+
δ
k
′
)
=
0 to get a solution, we find the equation system
(68)
(
∑
i
=
1
9
‍
∂
B
∂
T
i
′
∂
T
i
′
∂
p
0
0
∂
B
∂
k
′
)
(
δ
p
δ
k
′
)
=
-
B
(
T
′
,
k
′
)
,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 theS 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 asE and G. In summary, several classes of static solutions can be stud by using the approach here presented.



---

*Source: 650365-2014-01-09.xml*