Variational principles for nonlinear Kirchhoff rods

TL;DR

This paper develops and clarifies variational principles for nonlinear Kirchhoff rods, highlighting the role of differential geometry and proposing simplified models for static and dynamic problems.

Contribution

It systematically presents new and existing variational principles for nonlinear Kirchhoff rods, emphasizing geometric aspects and approximation methods.

Findings

New variational principles for nonlinear Kirchhoff rods

Clarification of geometric roles in rod modeling

Identification of approximations for simplified formulations

Abstract

The present article studies variational principles for the formulation of static and dynamic problems involving Kirchhoff rods in a fully nonlinear setting. These results, some of them new, others scattered in the literature, are presented in a systematic way, helping to clarify certain aspects that have remained obscure. In particular, the study of transversely isotropic models reveals the delicate role that differential geometry plays in their formulation and unveils consequently some approximations that can be made to obtain simplified formulations.

Equations

$$S^2:=\left\{\mathbf{d}\in {\mathbb{R}}^3\mid \mathbf{d}\cdot \mathbf{d}=1\right\},$$

$$T{S}^2:=\left\{(\mathbf{d},\mathbf{c}),\mathbf{d}\in S^2,\mathbf{c}\in {\mathbb{R}}^3,\mathbf{d}\cdot \mathbf{c}=0\right\}.$$

$$\mathbf{c}=\mathbf{w}\times \mathbf{d},\mathrm{with}\mathbf{w}\cdot \mathbf{d}=0,$$

$${\nabla }_{\mathbf{w}}\mathbf{v}:=(\mathbf{I}-\mathbf{d}\otimes \mathbf{d})D\mathbf{v}\cdot \mathbf{w}.$$

$${\nabla }_{{\mathbf{d}}^{\prime }}\mathbf{v}=(\mathbf{I}-\mathbf{d}\otimes \mathbf{d})D\mathbf{v}\cdot {\mathbf{d}}^{\prime }=(\mathbf{v}\circ \mathbf{d}{)}^{\prime }-\left((\mathbf{v}\circ \mathbf{d}{)}^{\prime }\cdot \mathbf{d}\right)\mathbf{d},$$

$$SO(3):=\left\{\Lambda \in {\mathbb{R}}^{3\times 3},{\Lambda }^T\Lambda =\Lambda {\Lambda }^T=\mathbf{I},\det \Lambda =1\right\}.$$

$$so(3):=\left\{\hat{\mathbf{w}}\in {\mathbb{R}}^{3\times 3},\hat{\mathbf{w}}=-{\hat{\mathbf{w}}}^T\right\}.$$

$$\exp [\hat{\mathbf{\theta }}]:=\mathbf{I}+\frac{\sin \theta }{\theta }\hat{\mathbf{\theta }}+\frac{1}{2}\frac{{\sin }^2(\theta /2)}{(\theta /2{)}^2}{\hat{\mathbf{\theta }}}^2,$$

$$\frac{\mathrm{d}}{\mathrm{d}ϵ}\exp [\hat{\mathbf{\theta }}(ϵ)]=\mathrm{skew}[\mathrm{dexp}[\hat{\mathbf{\theta }}(ϵ)]{\mathbf{\theta }}^{\prime }(ϵ)]\exp [\hat{\mathbf{\theta }}(ϵ)]$$

$${\mathbb{R}}^3\cong T_{\mathbf{d}}{S}^2\oplus span(\mathbf{d}),$$

$$\mathbf{T}={\mathbf{T}}_{\perp }+{\mathbf{T}}_{\parallel },$$

$$\mathbf{\chi }[{\mathbf{d}}_0,\mathbf{d}]:=({\mathbf{d}}_0\cdot \mathbf{d})\mathbf{I}+\widehat{{\mathbf{d}}_0\times \mathbf{d}}+\frac{1}{1+{\mathbf{d}}_0\cdot \mathbf{d}}({\mathbf{d}}_0\times \mathbf{d})\otimes ({\mathbf{d}}_0\times \mathbf{d}).$$

$${\mathbf{\chi }}_{\perp }[{\mathbf{d}}_0,\mathbf{d}]=\mathbf{\chi }-\mathbf{d}\otimes {\mathbf{d}}_0,{\mathbf{\chi }}_{\parallel }[{\mathbf{d}}_0,\mathbf{d}]=\mathbf{d}\otimes {\mathbf{d}}_0.$$

$$\Lambda =\exp [\psi \hat{\mathbf{d}}]\mathbf{\chi }=\mathbf{\chi }\exp [\psi \widehat{{\mathbf{d}}_0}]$$

$${\mathbf{u}}^{\prime }\cdot \mathbf{v}=\mathbf{u}\cdot {\mathbf{v}}^{\prime }=0$$

$$\mathbf{k}={\mathbf{r}}^{\prime }\times {\mathbf{r}}^{\prime \prime },$$

$${\mathbf{u}}_{\psi }=\cos (\psi )\mathbf{u}+\sin (\psi )\mathbf{v},{\mathbf{v}}_{\psi }=-\sin (\psi )\mathbf{u}+\cos (\psi )\mathbf{v}.$$

$$\mathbf{\omega }=\mathbf{k}+{\psi }^{\prime }{\mathbf{r}}^{\prime },$$

$$\mathbf{\omega }=-({\mathbf{v}}_{\psi }\cdot {\mathbf{r}}^{\prime \prime }){\mathbf{u}}_{\psi }+({\mathbf{u}}_{\psi }\cdot {\mathbf{r}}^{\prime \prime }){\mathbf{v}}_{\psi }+{\psi }^{\prime }{\mathbf{r}}^{\prime },$$

$${\mathbf{\omega }}_{\mathbf{\chi }}={\mathbf{r}}^{\prime }\times {\mathbf{r}}^{\prime \prime }+\left(\mathbf{a}\cdot {\mathbf{r}}^{\prime \prime }\right){\mathbf{r}}^{\prime },$$

$$\mathbf{a}=-\frac{{\mathbf{d}}_0\times {\mathbf{r}}^{\prime }}{1+{\mathbf{d}}_0\cdot {\mathbf{r}}^{\prime }}.$$

$${\mathbf{\omega }}_{\Lambda }={\mathbf{r}}^{\prime }\times {\mathbf{r}}^{\prime \prime }+\left({\psi }^{\prime }+\mathbf{a}\cdot {\mathbf{r}}^{\prime \prime }\right){\mathbf{r}}^{\prime }.$$

$${\psi }^{\prime }+\mathbf{a}\cdot {\mathbf{r}}^{\prime \prime }=0,$$

$$\psi (s)=-{\int }_{s_0}^s\mathbf{a}\cdot {\mathbf{r}}^{\prime \prime }\mathrm{d}\mu ,$$

$$\Lambda (s){\mathbf{E}}_i={\mathbf{e}}_i(s).$$

$$Q:=\left\{(\mathbf{r},\Lambda ):[0,L]\to {\mathbb{R}}^3\times SO(3),\mathbf{r}(0)=\bar{\mathbf{r}},\Lambda (0)=\bar{\Lambda }\right\},$$

$$|{\mathbf{r}}_0^{\prime }(s)|=1,$$

$$ϵ:=({\mathbf{r}}^{\prime }\cdot {\mathbf{e}}_3)-1.$$

$$\mathbf{\sigma }:=({\mathbf{r}}^{\prime }\cdot {\mathbf{e}}_{\alpha }){\mathbf{e}}_{\alpha }\text{with}\alpha =1,2.$$

$$\mathbf{\kappa }:={\mathbf{e}}_3\times {\mathbf{e}}_3^{\prime }.$$

Full text

1 Introduction

Since the seminal work of Kirchhoff in 1859 on rods subjected to bending and twisting, there have been innumerable works devoted to the study of such models, starting from the classical contributions by Clebsch [1], the Cosserat brothers [2], and Love [3]. The topic, already fairly mature, is certainly rich and has motivated detailed analyses during the last fifty years [4, 5, 6, 7, 8, 9, 10]. Other recent contributions continue to research more specific aspects of the theory like its stability [11, 12], the Hamiltonian structure [13, 14], discrete models [15], relations with differential geometry [16], etc.

Despite being such a classical topic, there is a reborn interest in using and analyzing nonlinear Kirchhoff rods, since they have found applications in fields as disparate as DNA modeling [17, 18, 19], hair simulation [20, 21], cables [22, 23], Möbius bands [24], oil-well drill-strings [25], computer graphics [26], space tethers [27], vortex-jet filaments [28], knots [29, 30], climbing plants [31], catheters for medical interventions [32], etc. In many of these situations, numerical methods provide the only means for arriving to useful solutions, motivating new developments in the computational side [33, 34, 35, 36, 37].

In this work, we focus on variational formulations of Kirchhoff rods, and discuss in detail aspects related to their restriction to the transversely isotropic case. Our interest in the variational approach is mainly due to the fact that it represents the basis of Galerkin-type numerical methods, most notably finite elements. Motivated by the development of more efficient numerical methods, novel Kirchhoff formulations have been proposed in recent years [23, 35, 36, 37] that demand a careful analysis, and studying their variational structure seems a useful starting point. The latter provides a solid setting for the geometry of the configuration space and a rigorous path for obtaining the governing equations of the problem, as well as the consistent boundary and initial conditions.

Variational principles for Kirchhoff rods are, of course, not new and many references can be found where the governing equations of the rod are found from the stationarity conditions of some potential energy [8, 7], in the case of a quasistatic problem, or Hamilton’s action, in a dynamic context [13, 38]. The easiest way to formulate them is to use variational principles for geometrically exact rods with shear deformation [39, 40], and then constrain the cross section orientation to be orthogonal with respect to the tangent vector of the curve that describes the rod in the ambient space, following Kirchhoff’s hypothesis. This results in constrained optimization principles that are far from optimal from the computational point of view.

We review under which conditions the equations of a Kirchhoff rod can be obtained from a variational principle without constraints. The key for such type of formulations is to eliminate the rotation group from the configuration space of the rod [7, 23, 35, 41, 42] and replace it with a simpler set that already accounts for the Kirchhoff constraint. It will be shown that it is indeed possible to formulate such models, and provide the variational principle behind them, but only when the rod is transversely isotropic. This is a very common situation both for analysis [7, 23, 35] as well as applications [17, 30, 23, 42], so its study is of practical value. There are delicate assumptions that are sometimes implicitly made in this transformation, assumptions that are not longer valid for dynamic problems. In the latter situation, we show what other simplifications need to be done if the full rotation group is to be left out of the formulation.

In view of the previous arguments, two types of variational principles are presented first: one constrained principle for (static and dynamic) general Kirchhoff rods, and a second family of principles for rods with isotropic cross sections. Finally, we will formulate a variational principle that is valid for general rods but, when applied to transversely isotropic ones, decouples in a certain way, making it very appealing for general purpose numerical applications.

The remainder of the article is structured as follows. Section 2 presents a collection of results coming from differential geometry that will play a key role in the definition of configuration spaces. In Section 3 we describe general concepts of nonlinear rods, including deformation measures, stress resultants, and momenta. These concepts will be used throughout the rest of the article. Section 4 presents a constrained variational principle for the static and dynamic analysis of extensible and inextensible Kirchhoff rods with general cross sections. Section 5 departs from the main argument of the article to introduce and discuss parametrizations of the rotational group as composite rotations. This is the key geometrical argument behind rod models that eliminate the rotation group from their configuration space, and hence we review it carefully. With these results at hand, we discuss in Section 6 variational principles for transversely isotropic Kirchhoff rods, carefully distinguishing between the quasistatic and dynamic cases, since there are different sets of hypotheses in each case. A variational principle suitable for general and transversely isotropic rods is described in Section 7. Finally, the main results of the article are summarized in Section 8.

2 Geometrical concepts

Differential geometry plays a critical role in the description of nonlinear structural models, and rods in particular. We review in this section the geometrical concepts employed in subsequent sections.

2.1 The unit sphere

The unit sphere $S^{2}$ plays a key role in the formulation of the rod models discussed in this article. For completeness, and to clarify the notation, we summarize some of its main features, leaving to other references more detailed descriptions [43, 44, 45]. This set is a nonlinear, smooth, compact, two-dimensional manifold defined as

[TABLE]

where the dot product refers to the standard Euclidean inner product. In the context of structural mechanics, especially rods and shells, elements of $S^{2}$ are sometimes referred to as directors [46, 47] since they are used in this type of models to describe relevant geometrical directions. Special attention must be paid to the fact that this manifold possesses no special algebraic structure, specifically group-like structure.

The tangent bundle of the unit 2-sphere is the manifold

[TABLE]

Tangent vectors at a point $\bm{d}\in S^{2}$ can be alternatively be described by the relations

[TABLE]

where the symbol “$\times$” denotes the cross product between vectors in $\mathbb{R}^{3}$.

The unit 2-sphere, together with the metric inherited from $\mathbb{R}^{3}$ has the structure of a Riemannian manifold. By embedding the manifold in Euclidean three dimensional space, the covariant derivative of a smooth vector field $\bm{v}:S^{2}\to TS^{2}$ along a second vector field $\bm{w}:S^{2}\to TS^{2}$ is the vector field $\nabla_{\bm{w}}{\bm{v}}$, which evaluated at $\bm{d}\in S^{2}$ coincides with the projection of the derivative $D\bm{v}$ in the direction of $\bm{w}$ onto the tangent plane to $\bm{d}$. If $\bm{I}$ denotes the second order unit tensor and $\otimes$ the dyadic product between vectors, this projection can be expressed as

[TABLE]

In particular, if $\bm{d}:(a,b)\to S^{2}$ is a smooth one-parameter curve on the unit sphere and $\bm{d}^{\prime}$ its derivative, the covariant derivative of a smooth vector field $\bm{v}:S^{2}\to TS^{2}$ in the direction of $\bm{d}^{\prime}$ can be evaluated as

[TABLE]

which, as before, is nothing but the projection of $(\bm{v}\circ\bm{d})^{\prime}$ onto the tangent space $T_{\bm{d}}S^{2}$.

2.2 The special orthogonal group

The set of proper orthogonal tensors also plays a prominent role in rod theories and can be defined as

[TABLE]

This set possesses a group-like structure when considered with the tensor multiplication operation, and it is also a smooth manifiold, hence it is a Lie group. Its Lie algebra is the set

[TABLE]

An isomorphism exists between vectors in $\mathbb{R}^{3}$ and $so(3)$ defined as $\hat{\cdot}:\mathbb{R}^{3}\to so(3)$ such that for all $\bm{w},\bm{a}\in\mathbb{R}^{3}$, the tensor $\hat{\bm{w}}\in so(3)$ satisfies $\hat{\bm{w}}\bm{a}=\bm{w}\times\bm{a}$. The vector $\bm{w}$ is referred to as the axial vector of the skew-symmetric tensor $\hat{\bm{w}}$ and we also write $\mathrm{skew}[\bm{w}]=\hat{\bm{w}}$. The exponential map $\exp:so(3)\to SO(3)$ is a surjective application with a closed form expression given by Rodrigues’ formula

[TABLE]

with $\bm{\theta}\in\mathbb{R}^{3}$, $\theta=|\bm{\theta}|$, and $|\cdot|$ denotes the Euclidean norm. The linearization of the exponential map is simplified by introducing the map $\mathrm{dexp}:so(3)\to\mathbb{R}^{3\times 3}$ that satisfies

[TABLE]

for every $\bm{\theta}:\mathbb{R}\to\mathbb{R}^{3}$. Explicit expressions of this map, and more aspects regarding the numerical treatment of the rotation group can be found elsewhere [48, 49, 45].

2.3 Composite rotations

For any director $\bm{d}$, the three-dimensional Euclidean space can be expressed as the direct sum

[TABLE]

where $\hbox{span}(\bm{d})$ is the linear subspace spanned by $\bm{d}$. Given now two directors $\bm{d},\tilde{\bm{d}}$, we say that a second order tensor $\bm{T}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ splits from $\bm{d}$ to $\tilde{\bm{d}}$ if it can be written in the form

[TABLE]

where ${\bm{T}}_{\perp}$ is a bijection from $T_{\bm{d}}S^{2}$ to $T_{\tilde{\bm{d}}}S^{2}$ with $\ker({\bm{T}}_{\perp})=\hbox{span}(\bm{d})$, and ${\bm{T}}_{\parallel}$ is a bijection from $\hbox{span}(\bm{d})$ to $\hbox{span}(\tilde{\bm{d}})$ with $\ker({\bm{T}}_{\parallel})=T_{\bm{d}}S^{2}$. The split (11) depends on the pair $\bm{d},\tilde{\bm{d}}$ but it is not indicated explicitly in order to simplify the notation.

Let us now consider a one-parameter curve in $S^{2}$ denoted $\bm{d}_{t}$ with $t\in[0,T]$. If $\bm{d}$ is an arbitrary point on $\bm{d}_{t}$, the two dimensional space $T_{\bm{d}}S^{2}$ can be viewed, when $S^{2}$ is embedded in $\mathbb{R}^{3}$, as a tangent plane to the unit sphere at $\bm{d}$. Between $T_{\bm{d}_{0}}S^{2}$ and $T_{\bm{d}}S^{2}$ there exist, hence, infinite isomorphisms. For example, parallel transport along $\bm{d}_{t}$ provides a natural map between these two vector spaces. Another useful transformation can be obtained from the unique rotation $\bm{\chi}\in SO(3)$ that maps $\bm{d}_{0}$ to $\bm{d}$ without drill, as long as $\bm{d}_{0}\neq-\bm{d}$, and defined by

[TABLE]

The map $\bm{\chi}[\bm{d}_{0},\bm{d}]$ splits from $\bm{d}_{0}$ to $\bm{d}$ and we can define

[TABLE]

The map ${\bm{\chi}}_{\perp}[\bm{d}_{0},\bm{d}]$ is, by definition, a bijection between $T_{\bm{d}_{0}}S^{2}$ and $T_{\bm{d}}S^{2}$ which, in contrast with the one defined by parallel transport, does not depend on the whole curve $\bm{d}_{t}$ but only on the directors $\bm{d}_{0}$ and $\bm{d}$.

Given a rotation $\bm{\chi}$ that splits from $\bm{d}_{0}$ to $\bm{d}$ and any scalar $\psi\in S^{1}$ (the unit circle), the map

[TABLE]

is also a rotation that splits from $\bm{d}_{0}$ to $\bm{d}$. Geometrically, ${\bm{{\mathit{\Lambda}}}}$ maps vectors on $T_{\bm{d}_{0}}S^{2}$ to $T_{\bm{d}}S^{2}$ with a drill angle equal to $\psi$. This map will play a key role in the theory that follows.

2.4 The Bishop (natural) frame and torsion

In the current context, it is desirable to describe the cross-section orientation along the rod by means of an orthogonal frame that is somehow intrinsically related to the curve that describes the rod in the ambient space.

To introduce such a rotation field, let $\bm{r}:\left[0,L\right]\rightarrow\mathbb{R}^{3}$ be a one-parameter curve with derivative $\bm{r}^{\prime}$. Let us assume that this curve is regular, that is, $\bm{r}^{\prime}\neq\bm{0}$ everywhere on its domain. Additionally, we assume, without loss of generality, that the curve is arc-length parametrized, i.e. , $|\bm{r}^{\prime}|=1$. Now, let us consider $\bm{u},\bm{v}:[0,L]\to S^{2}$ such that $\left\{\bm{u},\bm{v},\bm{r}^{\prime}\right\}$ are mutually orthogonal vector fields along the curve. We say that $\left\{\bm{u},\bm{v},\bm{r}^{\prime}\right\}$ minimizes the rotation if the following conditions

[TABLE]

are satisfied. Let $\{\bm{u}(s_{0}),\bm{v}(s_{0}),\bm{r}^{\prime}(s_{0})\}=\{\bm{u}_{0},\bm{v}_{0},\bm{r}^{\prime}_{0}\}$ be a known value of the orthonornal triad for some $s_{0}\in(0,L)$. Using this frame as initial condition, Eqs. (15) can be integrated along $\bm{r}$, defining uniquely the vector fields $\bm{u}$ and $\bm{v}$ at all points of the curve. This curve has Darboux vector

[TABLE]

which plays an essential role in defining parallel transport.

To analyze the concept of torsion in transported frames, let us consider the rotation field $\exp[\psi\>\hat{\bm{r}}^{\prime}]$, with $\psi:[0,L]\to S^{1}$, and the rotated triad $\{\bm{u}_{\psi},\bm{v}_{\psi},\bm{r}^{\prime}\}=\exp[\psi\>\hat{\bm{r}}^{\prime}]\{\bm{u},\bm{v},\bm{r}^{\prime}\}$, with

[TABLE]

The rotated frame has Darboux vector

[TABLE]

or, equivalently,

[TABLE]

where $\psi^{\prime}$ is the torsion and $\psi$ is the torsion angle (for further details, see [50, 7]). The previous calculations show that the frame $\{\bm{u},\bm{v},\bm{r}^{\prime}\}$ – which is known as the natural or Bishop frame – is the unique one obtained by transporting $\{\bm{u}_{0},\bm{v}_{0},\bm{r}^{\prime}_{0}\}$ along the curve without torsion.

Given, as before, a known frame $\{\bm{u}_{0},\bm{v}_{0},\bm{r}^{\prime}_{0}\}$ at a point $s_{0}\in(0,L)$, one could transport it along the curve using the drill-free rotation given by Eq. (12). Ignoring, for the moment, that this map is only defined as long as $\bm{r}^{\prime}\neq\bm{r}^{\prime}_{0}$, we note that this alternatively transported frame has, in general, an induced torsion. To see this, we observe that the Darboux vector emanated from the drill-free rotation $\bm{\chi}[\bm{r}^{\prime}_{0},\bm{r}^{\prime}]$ is

[TABLE]

where

[TABLE]

The vector $\bm{a}$, responsible for the torsion of the transported frame, is non-vanishing in general, as anticipated.

More importantly, the term $\bm{a}\cdot\bm{r}^{\prime\prime}$ in the curvature $\bm{\omega}_{\bm{\chi}}$ is non-integrable, meaning that there is no scalar function such that its derivative with respect to arc-length produces torsion in the sense of Eq. (19). As a result, the drill-free map cannot be corrected via a composite map (14) yielding a frame that has no torsion. To verify this assertion, we calculate the Darboux vector of this composite map to be

[TABLE]

If the map is to have vanishing torsion, the rotation angle $\psi$ would have to satisfy

[TABLE]

or, explicitly,

[TABLE]

which clearly is a non-local quantity.

3 Geometrically exact rods

We define in this section general concepts related to the kinematics of geometrically exact rods, that is, rod models for which no approximation whatsoever is made on the size of their deformation. We describe first a very general model and then we indicate what constraints can be imposed on it to obtain theories with reduced kinematics.

3.1 General description

As customary, a rod is defined to be a three-dimensional deformable body whose length is much larger than its other two dimensions. Such a body can be described by a curve in $\mathbb{R}^{3}$ and a cross section at every point of the curve whose intersection is precisely the barycenter of the section. We will henceforth assume that these cross sections remain plane and undistorted at all possible configurations of the rod. Models that account for the distorsion of the cross section can be found elsewhere (e.g. , [8]).

Let $L$ denote the length of the rod’s centerline in the undeformed configuration, and $s\in[0,L]$ an arc-length coordinate employed to identify each point of the rod. To describe the configuration of the rod we select a fixed coordinate system with a Cartesian orthonormal basis $\{\bm{E}_{i}\}_{i=1}^{3}$. The position of the centerline point of arc-length coordinate $s$ is denoted $\bm{r}(s)$. Moreover, since the cross section of the rod is assumed to remain plane and undistorted at all configurations, two orthonormal vectors $\left\{\bm{e}_{1}(s),\bm{e}_{2}(s)\right\}$ span the cross section with coordinate $s$ at all time, and we assume that they are material vectors. Defining a third unit vector $\bm{e}_{3}(s)=\bm{e}_{1}(s)\times\bm{e}_{2}(s)$ we can uniquely identify the orientation of the cross section at $s$ by the tensor $\bm{\mathit{\Lambda}}(s)\in SO(3)$ such that

[TABLE]

To make the presentation more precise, let us consider a rod that is clamped at $s=0$. Based on the previous description, the configuration manifold of a general rod of this type is the set

[TABLE]

where the functions $\bm{r}$ and $\bm{\mathit{\Lambda}}$ are assumed to be smooth and $\bar{\bm{r}},\bar{\bm{\mathit{\Lambda}}}$ are known. Let us stress here that this configuration space is not for Kirchhoff rods, but for general models.

Among all possible configurations in $Q$, the reference configuration of the rod is defined to be $(\bm{r}_{0},\bm{\mathit{\Lambda}}_{0})$, and we choose, without any loss of generality, that the argument of these two functions coincides with the arc-length of the curve defined by $\bm{r}_{0}$, i.e. ,

[TABLE]

for all $s\in[0,L]$.

3.2 Strain measures

We study next all the possible ways in which a general rod, as defined in Section 3.1, can deform. More detailed discussions of the deformation modes of these rods can be found in [40, 8].

The first deformation measure is the axial strain, and gauges the local relative value of the rod elongation, defined as

[TABLE]

The second measure is the shear strain, and it is defined as the vector

[TABLE]

The relative change in the normal to the cross section is accounted for by the bending strain

[TABLE]

The last measure is the torsional strain defined as the scalar

[TABLE]

and accounts for the relative rotation of the cross section about the director. For convenience, these measures are often gathered in two vectors

[TABLE]

which can be related with the centerline vector and orientation rotation through the relations

[TABLE]

The vectors $\bm{\gamma}$ and $\bm{\omega}$ are spatial strain measures whose convected counterparts are, respectively,

[TABLE]

and from these, we define the convected shear and bending strains

[TABLE]

respectively. All of these measures are frame invariant under the (left) action of the special Euclidean group $SE(3)$, i.e. , the group of isometries on the configuration space $Q$. See [40, 8].

3.3 Stored energy and stress resultants

In this section we study the most general hyperelastic model for a rod, and identify the stress resultants that are work-conjugate to the strains identified in Section 3.2.

To start, let us postulate the existence of an objective stored energy function of the form $U=U(\bm{\sigma},\epsilon,\bm{\kappa},\tau;s)$. The objectivity of this function implies that

[TABLE]

for every $\bm{Q}\in SO(3)$. By selecting, in particular, $\bm{Q}=\bm{\mathit{\Lambda}}^{T}$, it follows that

[TABLE]

Hence, by expressing the stored energy in terms of convected measures — which are invariant under the action of the special Euclidean group — we guarantee its objectivity.

The differential of the stored energy function can be calculated as

[TABLE]

and we define the convected stress resultants

[TABLE]

which are work conjugate, respectively, to the shear, axial, bending, and torsional strain measures in their convected form. The spatial stress resultants are the push-forwards $(\bm{\xi},n_{\parallel},\bm{m}_{\perp},m_{\parallel})=(\bm{\mathit{\Lambda}}\bm{\mathit{\Xi}},n_{\parallel},\bm{\mathit{\Lambda}}\bm{M}_{\perp},m_{\parallel})$.

For the remainder of this article, we will focus on beam models that have no shear deformation, i.e.  $\bm{\mathit{\Xi}}=\bm{0}$. In these situations, hence, it will be unnecessary to account for the contribution to the internal energy due to shear and we will always assume that the latter is of the form $\tilde{U}=\tilde{U}(\epsilon,\bm{K},\tau;s)$.

3.4 Momenta and kinetic energy

Let us now consider a motion of the rod, that is, a time-parameterized curve in configuration space described by a pair of functions $(\bm{r}(s,t),\bm{\mathit{\Lambda}}(s,t))$ such that $(\bm{r}(\cdot,t),\bm{\mathit{\Lambda}}(\cdot,t))\in Q$ for all $t\in[0,T]$. Using a superposed dot as the notation for the derivative with respect to time, the generalized velocity of the rod is the vector field $(\dot{\bm{r}},\dot{\bm{\mathit{\Lambda}}})$ belonging, for every $t\in[0,T]$, to the tangent bundle

[TABLE]

The time derivative of the centroid position is the velocity $\bm{v}=\dot{\bm{r}}$, and the derivative of the rotation can be written as

[TABLE]

where $\bm{w}$ and $\bm{W}$ are the spatial and convected angular velocities, respectively. These fields can be split as in

[TABLE]

with

[TABLE]

The kinetic energy density of the rod is defined as

[TABLE]

The constant $A_{\rho}$ is equal to the product of the cross section area and the density of the rod’s material. The tensor $\bm{i}_{\rho}$ is the spatial inertia of the cross section. The convected inertia, constant in time, is the pullback of $\bm{i}_{\rho}$ to the reference configuration, namely,

[TABLE]

The convected inertia $\bm{I}_{\rho}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ is a second order, symmetric, positive definite tensor that splits from $\bm{E}_{3}$ to $\bm{E}_{3}$, i.e. , it is of the form

[TABLE]

where $\bm{I}_{\perp}$ maps bijectively $\mathrm{span}(\bm{E}_{1},\bm{E}_{2})$ onto itself and satisfies $\bm{I}_{\perp}\bm{E}_{3}=\bm{0}$. By pushing this inertia to the current configuration, we find that the spatial inertia admits a similar split

[TABLE]

where now $\bm{i}_{\perp}$ splits from $\bm{e}_{3}$ to $\bm{e}_{3}$. As a consequence of the structure of the inertia tensor, the rotational part of the kinetic energy density can be written in the following equivalent ways:

[TABLE]

The translational and rotational momenta of the rod are conjugate to the velocities as in

[TABLE]

and we note that we can introduce a convected version of the momentum $\bm{\pi}$ by pulling it back, as before, and defining

[TABLE]

Due to the particular structure of the inertia the momenta can also be split, as before, as in

[TABLE]

with

[TABLE]

4 The canonical variational principle for Kirchhoff rods

In this section we obtain the governing equations of Kirchhoff rods from Hamilton’s principle of stationary action in the most straightforward way.

4.1 Lagrangian

Kirchhoff rods are a subclass of the general rod models described in Section 3, with the condition of vanishing shear deformation, i.e. , $\bm{\mathit{\Sigma}=\bm{0}}$. In view of Eq. (29), this constraint is equivalent to having $\bm{r}^{\prime}$ parallel to the vector $\bm{e}_{3}$, or equivalently

[TABLE]

Instead of trying to define a configuration space for this class of rods, we will work with the general space defined in Eq. (26), and constrain the kinematics of the model using Lagrange multipliers to impose pointwise relation (53).

To make the exposition as simple as possible, we assume that the rod is clamped at the section corresponding to $s=0$. Moreover, a known conservative distributed force, per unit of reference length, $\bar{\bm{n}}:[0,L]\to\mathbb{R}^{3}$ is applied along the rod, and a conservative point force $\tilde{\bm{n}}$ is applied at the free end, namely $s=L$. Under this loading, the potential energy of the rod is

[TABLE]

Using the definitions of Section 3.4, the kinetic energy of this system can be calculated to be

[TABLE]

and the Lagrangian is just $\mathcal{L}=T-V$.

4.2 Variations of the strain measures and rates

The stationarity conditions for the action will be obtained using calculus of variations. We gather next some results that will prove necessary for the computation of the functional derivatives and, later, for the linearization of the model.

To introduce these concepts, let us consider a curve of configurations $(\bm{r}_{\iota},\bm{\mathit{\Lambda}}_{\iota})$ parametrized by the scalar $\iota$ and given by

[TABLE]

where $\delta\bm{r}:[0,L]\times[0,T]\to\mathbb{R}^{3}$ and $\delta\bm{\theta}:[0,L]\times[0,T]\to\mathbb{R}^{3}$ are arbitrary variations that satisfy

[TABLE]

The curve $(\bm{r}_{\iota},\bm{\mathit{\Lambda}}_{\iota})$ passes through the configuration $(\bm{r},\bm{\mathit{\Lambda}})$ when $\iota=0$ and has tangent at this point equal to

[TABLE]

For future reference let us calculate the variation of the derivative $\bm{\mathit{\Lambda}}^{\prime}$. To do so, let us first define the arc-length derivative of the perturbed rotation, that is,

[TABLE]

Then, the variation of $\bm{\mathit{\Lambda}}^{\prime}$ is just

[TABLE]

With the previous results we can now proceed to calculate the variations of the strain measures, as summarized in the following theorem.

Theorem 4.1**.**

The linearization of the three strain measures $(\epsilon,\bm{K},\tau)$ is

[TABLE]

Proof.

The strain measures of the one-parameter curve of configurations $(\bm{r}_{\iota},\bm{\mathit{\Lambda}}_{\iota})$ are

[TABLE]

where $\bm{e}_{i,\iota}=\bm{\mathit{\Lambda}}_{\iota}\bm{E}_{i}$. The variation of the axial strain measure is computed as follows:

[TABLE]

The variation of the bending strain is obtained from its definition employing some algebraic manipulations and expression (60) as follows:

[TABLE]

The linearization of the torsional strain follows similar steps:

[TABLE]

∎

The linearization of the rates is almost identical to the one of the strains. For convenience, we present without proof the attendant results in the following theorem:

Theorem 4.2**.**

The linearization of the three time rates $(\bm{v},\bm{W}_{\perp},w_{\parallel})$ is:

[TABLE]

4.3 Governing equations

Hamilton’s principle of stationary action states that the governing equations of the clamped-free rod are the Euler-Lagrange equations of the constrained action functional

[TABLE]

with unknown fields $(\bm{r},\bm{\mathit{\Lambda}},\dot{\bm{r}},\dot{\bm{\mathit{\Lambda}}})$ in $TQ$ and Lagrange multiplier $\bm{\eta}:[0,L]\to\mathrm{span}(\bm{E}_{1},\bm{E}_{2})$. The main result of the section is as follows:

Theorem 4.3**.**

The Euler-Lagrange equations of the action (67) are:

[TABLE]

where the stress resultants are defined as

[TABLE]

Eq.(68) can be rewritten as

[TABLE]

where $\bm{d}=\bm{e}_{3}=\frac{\bm{r}^{\prime}}{|\bm{r}^{\prime}|}$. The natural boundary conditions at $s=L$ are

[TABLE]

Proof.

The first part of the theorem follows from a systematic calculation of $\delta S$, the variation of the action, and the results of Section 3.2, and thus we omit a detailed derivation. It should be noted that $\bm{n}$ is the contact force on the cross section, which can be additively decomposed on an axial part, $\bm{n}_{\parallel}$, and the shear force $\bm{n}_{\perp}$, which appears as a result of the constraint. The contact torque $\bm{m}$ itself is the sum of the bending moment $\bm{m}_{\perp}$ and the torsional moment $\bm{m}_{\parallel}$.

Let $\bm{d}=\bm{r}^{\prime}/|\bm{r}^{\prime}|$. Assuming $\bm{r}$ is a smooth function in both $s$ and $t$, and $|\bm{r}^{\prime}|>0$, condition (68c) is equivalent to the constraint $\bm{d}=\bm{e}_{3}$. Then, $\bm{n}_{\parallel}=\frac{\partial\tilde{U}}{\partial\epsilon}\bm{d}$ and $\bm{r}^{\prime}\times\bm{n}=\bm{r}^{\prime}\times\bm{n}_{\perp}$. From Eq. (68b) if follows that

[TABLE]

Since $\bm{r}^{\prime}$ and $\bm{n}_{\perp}$ are perpendicular, the shear resultant can be found to be

[TABLE]

which can be inserted in (68a) to obtain Eq. (70a). Then, projecting Eq. (72) onto the $\bm{d}$ direction, Eq. (70b) follows. It bears emphasis that the covariant derivative ensures that the differentiation is compatible with the underlying manifold structure. ∎

Eqs. (70) are derived solely from variational arguments, including the correct boundary conditions. One possible simplification would be to ignore the contribution of the rotational inertia, leading to a model discussed by Meier [41], and obtained directly as a projection of the Simo-Reissner beam theory. Another possible simplification could be to ignore the extensibility of the beam. A variational derivation of the corresponding governing equations would be almost identical to the one presented above. However, in this case, both the axial and shear strains vanish and thus a constraint must be included that not only imposes that $\bm{r}^{\prime}$ is parallel to $\bm{e}_{3}$, as in Eq. (53), but rather that these two quantities are identical. To impose this restriction, the constraint (53) must be replaced by the stronger one

[TABLE]

Since these are three constraints that must be satisfied pointwise, the field of Lagrange multipliers that would be required to impose them would be functions from $[0,L]$ to $\mathbb{R}^{3}$. A variational principle for this class of inextensible, Kirchhoff rods has been presented by Antman [8].

4.4 Linearization

The governing equation of the Kirchhoff rod, as given by Theorem 4.3 can be linearized resulting in the Rayleigh model [51], as shown in the following result:

Theorem 4.4**.**

Consider a straight rod with constant cross section, aligned in its reference configuration with the $\bm{E}_{3}$ axis and let its centerline position be written as

[TABLE]

Assuming that the torsion angle is $\varphi$, the linearized Euler-Lagrange equations from Theorem 4.3 are:

[TABLE]

Eqs. (76a)-(76b) correspond to Rayleigh’s beam equations for bending, including rotational inertia [51]; Eqs. (77c)-(77d) model, respectively, the axial and torsional behavior. The attendant boundary conditions at $s=0$ are

[TABLE]

Lastly, the natural boundary conditions at $s=L$ are:

[TABLE]

Proof.

Let us introduce the smallness parameter $\varsigma$ and redefine the displacement fields $u$, $v$ and $w$ as $\varsigma\bar{u}$, $\varsigma\bar{v}$ and $\varsigma\bar{w}$, respectively. Similarly, the torsion angle $\varphi$ is redefined as $\varsigma\bar{\varphi}$. Such re-scaled fields allow us to reveal the order of each term in the smallness parameter ($\varsigma^{n}$ with $n=0,1,\ldots,N,\ldots,\infty$). The linearization consists in retaining only those terms up to order $n=1$ and therefore, every term for $n>1$ is then collected by the term $\mathcal{O}(\varsigma^{2})$ since no further order distinction is necessary.

Next, we employ the asymptotic expansion to linearize the governing equation and boundary conditions. For this purpose, the rotation tensor associated to the cross section of the rod can be expressed as

[TABLE]

Let us not that $\bm{\mathit{\Lambda}}$ is not an arbitrary rotation tensor, but one that guarantees that normal vector to the cross section remains parallel to $\bm{r}^{\prime}$ and rotates about that same direction an angle $\varphi$ with respect to the natural frame.

The following expansions follow directly from the definition of the nonlinear terms:

[TABLE]

Replacing these expansions in Eqs. (70)-(71), removing the terms of order $\mathcal{O}(\varsigma^{2})$, and recovering the original fields $u$, $v$, $w$ and $\varphi$, we finally obtain Eqs. (76)-(78). ∎

5 Interlude: Kirchhoff rods theories with composite rotations

The Kirchhoff rod models described in Section 4 are formulated in a configuration space (26) identical to the one employed in rods with shear. The Kirchhoff constraint is then variationally included in the action and Lagrange multipliers are added to impose it.

Several works can be found in the literature where the use of Lagrange multipliers is avoided by defining a new configuration space, one in which the constraint is already accounted for, sparing the need for Lagrange multipliers in the variational formulation (e.g. ,[34, 23, 35, 41, 38]). If $\bm{d}=\bm{r}^{\prime}/|\bm{r}^{\prime}|$, these references define the rotation of each cross section as a composite map

[TABLE]

One rotation $\bm{R}$ maps $\bm{d}_{0}$ to $\bm{d}$ and is followed by a rotation about $\bm{d}$ of angle $\psi$. Such a rotation definition requires only a smooth field $\bm{r}$ plus an additional scalar for the second rotation, avoiding a full parameterization of the rotation field and the use of Lagrange multipliers altogether.

These formulations seem very appealing, since they not only remove the rotation group from the configuration space, but avoid the use of Lagrange multipliers. They suffer however, from three drawbacks that might be important. First, this kind of composite parameterizations often has singularities. A common choice for $\bm{R}[\bm{d}_{0},\bm{d}]$ is $\bm{\chi}[\bm{d}_{0},\bm{d}]$, the unique rotation that maps $\bm{d}_{0}$ to $\bm{d}$ without drill defined in Eq. (12), that is undefined when $\bm{d}=-\bm{d}_{0}$. Such singularities might not be relevant for problems with small rotations, or incremental solutions. However, a complete, self-contained theory of Kirchhoff rods cannot be based on a singular parameterization.

The second drawback of composite parameterizations is related to the imposition of boundary conditions. Taking the derivative with respect to the arc-length of Eq. (81) the curvature $\bm{K}$ and torsional strain $\tau$ can be calculated. A careful derivation shows that, if $\bm{R}$ is obtained with expression (12), the strain $\bm{\omega}$ is of the form

[TABLE]

with $k_{\alpha}=(\bm{d}\times\bm{d}^{\prime})\cdot\bm{e}_{\alpha}$ for $\alpha=1,2$, and the torsional part of the strain is not $\psi^{\prime}$ but has and additional term $\phi$ that comes from the torsion of the drill-free frame, as discussed in Section 2.4. This new term might result in difficulties for imposing boundary conditions, since the variationally consistent Dirichlet boundary conditions for angles involve $\psi$ and a complex functions of $\bm{r}$ and $\bm{r}^{\prime}$.

One could use an intermediate rotation $\bm{R}$ whose torsion is zero and depends only on the curve $\bm{r}$, that is, Bishop’s frame $\bm{B}$. Using such a frame would remove the indeterminacy in the boundary conditions, since the strain $\bm{\omega}$ would be simply of the form

[TABLE]

with $\bm{e}_{1}=\bm{B}\bm{E}_{1},\bm{e}_{2}=\bm{B}\bm{E}_{2}$. While using this frame would solve many of the difficulties alluded to above, there is no explicit expression for Bishop’s frame depending solely on $\bm{d}_{0}$ and $\bm{d}$.

The third drawback is more subtle and has not been identified, to the authors’ knowledge, previously. When using composite frames, the angular velocity is a complex function. For a Bishop (natural) frame one might expect that the angular velocity of the cross section would be of the form

[TABLE]

but this is not the case. The component of the angular velocity in the $\bm{d}$ direction can not be calculated in closed form, and a Lagrangian can not be formulated in terms of $\bm{r},\psi$ and their derivatives.

To further illustrate this issue, let us compute the angular velocity of the composite rotation (81). It has the form

[TABLE]

where

[TABLE]

If the angular velocity is to have a simple expression of the form (84), we could try to replace the spin rotation $\exp[\psi\hat{\bm{d}}]$ in Eq. (81) by $\exp[(\psi+\xi)\hat{\bm{d}}]$, so that the new angular velocity would be of the form

[TABLE]

and the term $\bm{a}\cdot\dot{\bm{r}}^{\prime}$ could be eliminated by imposing the constraint

[TABLE]

As proven in Section 2.4, a restriction of this type is non-integrable. Moreover, choosing such a correction angle might spoil the simple form of the (spatial) curvature of the frame.

To see this more explicitly, consider the convected strain $\bm{\mathit{\Omega}}$ and angular velocity $\bm{W}$ of the composite frame whose components in the $\bm{E}_{3}$ direction are, respectively, $\psi^{\prime}$ and $\dot{\phi}$. A general result for the compatibility of strain and angular velocity [6] states that

[TABLE]

From this relation, a straightforward manipulation gives

[TABLE]

proving that, in general, $\psi\neq\phi$.

6 Variational principles for transversely isotropic Kirchhoff rods

In this section we study rods that have transversely isotropic cross sections, in the sense that the stored energy function and the in-plane inertia are of the form

[TABLE]

We note that $|\bm{K}|=|\bm{\kappa}|=|\bm{e}_{3}^{\prime}|$ is a scalar bending strain that is frame invariant. This class of Kirchhoff rods are interesting for two reasons: first, they are very common in applications and second, they admit certain simplifications in their formulation as compared with the models presented in Section 4. In particular, the quasistatic formulation admits a greatly simplified variational principle. The dynamic case, however, demands a careful consideration.

6.1 Quasistatic problems

The simple form of the stored energy function for transversely isotropic rods allows to bypass completely the use of rotations in the formulation, although a new unknown field, the twist $\psi$, is added to account for the torsional deformation. Geometrically speaking, the isotropy of the bending response avoids the need to pull-back the bending strain to the reference configuration in order to calculate the bending moment, since it is known a priori that this moment is parallel to the bending strain. Since the rotation tensor is needed for this transformation and for the computation of the torsional strain, only the latter remains necessary. In the following derivation, however, we avoid the use of the rotation altogether by introducing an additional configuration term.

The configuration space of this model, for the rod clamped at $s=0$, is the manifold

[TABLE]

and the potential energy is the functional

[TABLE]

with $\bar{\bm{n}},\bar{m}_{\parallel}$ being, respectively, known fields of forces and tangent moment per unit length, and $\tilde{\bm{n}},\tilde{m}_{\parallel}$ a known point force and a known tangent moment at the end $s=L$. In this case, the strain measures have the simple form

[TABLE]

where, as before, $\bm{d}:=\bm{r}^{\prime}/|\bm{r}^{\prime}|$. The governing equations of this model are obtained as the stationarity conditions of $V$.

6.1.1 Strain variations

The strains defined in Eq. (94) are based on the strain measures of the general rod model (cf. Section 3.2). However, since the rotation is no longer an independent field, the variations also need to be redefined. To calculate the strain variations we introduce, as in Section 4.2, a curve of perturbed configurations in $Q$

[TABLE]

where $\iota\in\mathbb{R}$, $\delta\bm{r}:[0,L]\to\mathbb{R}^{3}$ and $\delta\psi:[0,L]\to\mathbb{R}$ are arbitrary fields that satisfy

[TABLE]

The advantage of the new configuration space is apparent already, since the definition of the perturbed configurations is additive, avoiding the use of the exponential map.

The variations of the strain measures are presented in the following result.

Theorem 6.1**.**

The variations of the strains in the transversally isotropic Kirchhoff rod are:

[TABLE]

with

[TABLE]

Proof.

The variations of $\epsilon,|\bm{K}|$ and $\tau$ are obtained by a systematic application of the chain rule and the definition of the perturbed configurations (95). For convenience, we write

[TABLE]

and we define the last parenthesis to be equal to $\delta\bm{\beta}$, an arbitrary vector field on the rod, orthogonal to $\bm{d}$. The variation of $\epsilon$ then follows trivially.

To calculate the variation of the bending strain $|\bm{K}|$, it is convenient to write a one-parameter curve of directors in the form

[TABLE]

whose derivative with respect to the arc-length is

[TABLE]

Thus, the variation of $\delta\bm{d}^{\prime}$ is

[TABLE]

which, together with

[TABLE]

yields the variation of the second strain. The variation of the torsional strain is trivial. ∎

6.1.2 Quasistatic equilibrium equations

The equilibrium equations of the rod are obtained from the Euler-Lagrange equations of the potential energy.

Theorem 6.2**.**

The equilibrium equations for a transversely isotropic Kirchhoff rod are:

[TABLE]

where the stress resultants are defined as

[TABLE]

and the natural boundary conditions at $s=L$ are

[TABLE]

Proof.

The theorem follows from a systematic calculation of $\delta V$, the variation of the potential energy, which is developed next:

[TABLE]

Defining the axial force, bending, and torsion moment, respectively, as in Eq. (105) and integrating by parts, it follows that

[TABLE]

Finally, expanding $\delta\bm{\beta}$ as in Eq. (98), and integrating by parts a second time we arrive at the final expression for the variation:

[TABLE]

Using the condition $\delta V=0$ and the arbitrariness of the variations, the theorem is proved. ∎

6.2 Difficulties with dynamic problems

The variational principle of transversely isotropic Kirchhoff rods can be recovered from the general case when the rotation field is chosen to be the composition of Bishop’s frame and a rotation of angle $\psi$ around the tangent vector $\bm{d}$. As a result of the isotropy, Bishop’s frame does not play any role and its calculation can be avoided, leaving a model with only four unknown fields, namely, the three components of $\bm{r}$ and $\psi$.

The simplicity of the four-field formulation makes it very appealing for its use in numerical computations and it has been exploited in many works. Its extension to dynamical problems has also been considered, but not always in a fully correct fashion.

In a transversely isotropic rod the strain measure $\bm{\omega}$ can be expressed as

[TABLE]

and $\bm{\mathit{\Omega}}=\bm{\mathit{\Lambda}}^{T}\bm{\omega}$, when $\bm{\mathit{\Lambda}}$ is the composite rotation based on Bishop’s frame. In view of the parallelism between the strain $\bm{\omega}$ and the spatial angular velocity $\bm{w}$, it is tempting to assume that the latter must be of the form

[TABLE]

where the spin velocity satisfies $\dot{\phi}=\dot{\psi}$. This formula for the angular velocity has been used in the past [35], but it is not true, in general, as discussed in Section 5.

6.3 Variational formulation of dynamic problems

It remains to find a variational formulation for dynamic problems of transversely isotropic Kirchhoff rods, one that take advantage of the four-field formulation. As discussed before, this is not straightforward because there is no simple way to express the kinetic energy in terms of $\bm{r}$ and $\psi$. To formulate the full kinetic energy it is thus necessary to introduce the rotation field again in the variational principle.

One might keep the advantages of the four-field formulation by neglecting the kinetic energy associated to the spin around $\bm{d}=\bm{r}^{\prime}/|\bm{r}^{\prime}|$, i.e. , assuming its density to be of the form

[TABLE]

where $\bm{w}_{\perp}=\bm{d}\times\dot{\bm{d}}$. In view of the isotropy property of the beam, the latter can be written simply as

[TABLE]

This approximation on the form of the kinetic energy precludes the use of the model for problems in which a rod spins around it centerline although it might be an admissible assumption for many practical problems (cable deployment, DNA modeling, etc.) and has often been neglected tacitly (e.g. , [27]). In any case, since the field $\psi$ will appear in the governing equations of the model, but not its rate nor acceleration, the latter will admit torsional waves of infinite speed. In order to recover the structure of an hyperbolic initial boundary-value problem one might choose to regularize the kinetic energy (113) by defining it to be

[TABLE]

recognizing from the outset that this is not the exact expression for the kinetic energy but merely a convenient approximation. We present next the variational formulation and governing equations of the transversely isotropic Kirchhoff rod with the approximated kinetic energy (114). As in the general case, we use Hamilton’s principle of stationary action to derive the governing equations for the rod under consideration.

In contrast with the model presented in Section 4, the variational principle for the model under consideration does not include any constraint. The action functional is of the form

[TABLE]

where the kinetic energy $T$ is the integral over the rod of the density $k$ as given by Eq. (114) and the potential energy $V$ is defined in Eq. (93). The Euler-Lagrange equations of the action (115) are obtained using the same type of calculations as in previous models. We present without proof the main result:

Theorem 6.3**.**

The dynamic equilibrium equations of a transversely isotropic Kirchhoff rod are:

[TABLE]

where the stress resultants are defined in Eq. (105), and the natural boundary conditions at $s=L$ are

[TABLE]

We conclude this section reflecting upon whether the use of dynamic formulations for transversely isotropic Kirchhoff rods is a sensible choice or not. As discussed at length, this type of models always entails approximations on the kinetic energy that can only be justified by arguing that the torsional effects play a small role in the model. In these cases, however, one wonders if it would not be better to employ, from the outset, a rod model that completely eliminates the torsion from the equations [44].

7 A joint variational principle for general and transversely isotropic Kirchhoff rods

Sections 4 and 6 discuss, respectively, variational principles for general and transversely isotropic Kirchhoff rods. The configuration space defined for the first type of rods is too rich, and Lagrange multipliers have to be used in the variational principle to constrain the rod’s kinematics. This is in contrast with the configuration space of the second kind, which fits exactly the kinematics of the rod and makes unnecessary to employ constraints in the corresponding variational principle. On the downside, this latter formulation is only valid for transversely isotropic rods.

In this section we propose a third variational principle that is valid for rods of general cross sections, but simplifies when the model is transversely isotropic. As with the principles of Section 6, the new functional is strictly exact only for quasistatic problems. Its extension to dynamic problems requires that the part of the kinetic energy associated with the rotation of the cross section about its director be neglected.

The idea of the new developments is to use, as in previous sections, Hamilton’s principle of stationary action to obtain the equations of motion of Kirchhoff rods, but in a way that the solution of certain fields can be decoupled from the rest of unknowns when the rod is transversely isotropic. In the general case, the variational principle is completely equivalent to the one studied in Section 4 if the problem is quasistatic or the contribution from the inertia mentioned above is ignored.

The new principle is defined in terms of the unknown fields $\bm{r},\bm{\mathit{\Lambda}},\psi:[0,L]\to\mathbb{R}^{3}\times SO(3)\times\mathbb{R}$ and the Lagrange multipliers $\bm{\eta},\mu:[0,L]\to\mathbb{R}^{2}\times\mathbb{R}$. The unknown fields refer, as in previous sections, to the position of the curve of centroids, the section orientation, and the torsional angle. Based on these, the potential energy of the rod has the usual form

[TABLE]

but now the strain measures are defined as

[TABLE]

with

[TABLE]

We note that the strain measures are completely identical to the ones defined in Eqs. (28)-(31) and Eq. (94) when the Kirchhoff constraint is verified, but the arguments are different. In contrast with the model of Section 6, the full bending strain enters now the stored energy function, but the role played by the rotation $\bm{\mathit{\Lambda}}$ in the definition of $\epsilon$ and $\bm{K}$ is different than in the general model of Section 4.

Next, the kinetic energy is defined to be

[TABLE]

with $\bm{w}_{\perp}=\bm{d}\times\dot{\bm{d}}$. We recall that the rotational part of this energy can be written with convected objects using the relation:

[TABLE]

Finally, to impose Kirchhoff’s constraint and that $\psi$ be the torsional angle, we define the action of the model as the constrained functional:

[TABLE]

Barring the simplifications in the definition of the kinetic energy, the actions (123) and (67) are clearly equivalent, the only difference being the addition of a new variable $\psi$ defined by the relation $\psi^{\prime}=\bm{\mathit{\Omega}}\cdot\bm{E}_{3}$.

7.1 Strain and velocity variations

The definitions of the strains (119) and the kinetic energy density in Eq. (121) are different from the ones of Sections 4 and 6, although identical when the Kirchhoff constraint is verified. Their variations will be used in Section 7.2 and need to be calculated anew.

To this end, let us define the configuration space for this principle considering, for simplicity, once again the case of a rod clamped at the end $s=0$:

[TABLE]

Proceeding as in Sections 4 and 6, we define the one-parameter curve of configurations

[TABLE]

where $\iota\in\mathbb{R}$, and $\delta\bm{r},\delta\bm{\theta},\delta\psi,\delta\bm{\eta},\delta\mu$ being arbitrary variations with the properties

[TABLE]

vanishing at $t=0$. Then, a systematic application of the concept of linearization yields the following result, which we present without proof since it is very similar to Theorem 6.1.

Theorem 7.1**.**

The variations of the strains (119) are:

[TABLE]

with $\delta\bm{\beta}$ defined as in Eq. (98). Similarly, the variation of the convected angular velocity is given by

[TABLE]

7.2 Governing equations

The stationarity conditions of the action (123) are completely equivalent to Eqs. (70), as the following result shows:

Theorem 7.2**.**

The Euler-Lagrange equations corresponding to the action (123) are:

[TABLE]

where the stress resultants are defined in Eq. (105). The attendant natural boundary conditions at $s=L$ are:

[TABLE]

The natural boundary conditions at $s=L$ follow easily from the last integral in Eq. (132) and the arbitrariness of the variations at the free end of the beam.

Proof.

The variation of the action gives:

[TABLE]

Replacing the formulas (127)-(128) for the variations, using the property that the variations vanish at $t=0$, integrating by parts in time the first integral and in arc-length the remaning ones, we get, after some straightforward manipulations:

[TABLE]

Next, we show that the condition $\delta S=0$ gives the same differential equations of the general rod model described in Section 4. First, we note that the Lagrange multipliers impose strongly the conditions

[TABLE]

The first of these conditions provides an interpretation for the Lagrange multiplier $\mu$. The second one identifies $\psi$ as the torsion angle. The third one is equivalent to the relation $\bm{d}=\bm{e}_{3}$. We will use these relations to simplify the notation in the expressions that follow.

Since the variations $\delta\bm{r}$ are arbitrary, the first integral in Eq. (132) must vanish, and thus

[TABLE]

with $\bm{n}_{\perp}=(\bm{\mathit{\Lambda}}\bm{\eta})\times\bm{e}_{3}$. Next, we consider the second integral in Eq. (132), which must also vanish in view of the arbitrariness of the variation $\delta\bm{\theta}$. Hence

[TABLE]

Deriving the relation $0=\bm{m}_{\perp}\cdot\bm{d}$ with respect to the arc-length we get that $\bm{m}_{\perp}\cdot\bm{d}^{\prime}=-\bm{m}_{\perp}^{\prime}\cdot\bm{d}$, and it follows that

[TABLE]

Then, projecting both sides of Eq. (135) in the direction of $\bm{d}$ and using Eq. (136) we obtain:

[TABLE]

with $\bm{m}=\bm{m}_{\perp}+\mu\,\bm{d}$. Using Eq. (133), a simple manipulation of Eq. (137) yields

[TABLE]

which coincides with Eq. (129b). Projecting Eq. (136) onto the plane orthogonal to $\bm{d}$ we get

[TABLE]

from where we can obtain the closed form expression for the symbol $\bm{n}_{\perp}$:

[TABLE]

Replacing the value of $\bm{n}_{\perp}$ in Eq. (134) we arrive at the final expression:

[TABLE]

∎

Interestingly, when formulating this variational principle for a transversely isotropic rod several simplifications follow. First, by definition, the stored energy function depends on $\bm{K}$ only through its modulus, so the potential energy is not a function of the section orientation. Then, again by definition of isotropy, the inertia $\bm{I}_{\perp}$ is a multiple of the identity in $\mathrm{span}(\bm{E}_{1},\bm{E}_{2})$, and thus the kinetic energy does not depend on the section orientation either. As a result, for the transversely isotropic case, the fields $\bm{r},\psi$ can be found by applying Hamilton’s principle to the action (123) without solving for $\bm{\mathit{\Lambda}},\bm{\eta},\mu$. In a way, the problem decouples: the solution for $\bm{r},\psi$ is independent of the rest of unknowns, which can be calculated after the expressions for the former are found. In many practical applications, the latter might even be superfluous.

It bears emphasis that the action (123) is not just the action of one of the two models presented, respectively, in Sections 4 and 6, with some added terms. The key for this new variational formulation is the reparametrization of the strain measures and the approximation of the kinetic energy. In fact, it is easy to verify that the functional of Section 4 for general rods does not simplify when the model is transversely isotropic.

8 Summary

This article presents a systematic analysis of Kirchhoff rods for general and transversely isotropic cross sections. The framework is that of variational analysis that, elegantly, provides the governing and boundary, conditions of each model, and helps to identify potential pitfalls in their geometric setting.

The first result is a variational principle for general rods that can deform axially, in bending, and in torsion, but not in shear. This is a minor modification of the well-known counterpart of inextensible Kirchhoff rods and provides a clear illustration of the formulation of structural models by means of constrained variational principles.

A very common motivation for mechanicians is the development of a constraint-free theory for Kirchhoff rods. We discuss that this is just possible, in the quasistatic case, for transversely isotropic beams, and we provide the corresponding variational principle. For dynamic problems we argue that it is not possible to find an exact constraint-free theory, not even for transversely isotropic rods. For the latter case we postulate approximate models that admit a variational principle free of any constraint.

The previous results provide two alternative variational principles: one valid for any Kirchhoff rod – but constrained – and a second, constraint-free, valid for transversely isotropic rods. With a view in numerical methods, we propose a third route: a variational principle that is valid for general Kirchhoff rods but simplifies in the transversely isotropic case, so that the constraints can be solved pointwise in closed form, as long as the kinetic energy is approximated.

9 Acknowledgements

I.R. was funded by project DPI2017-92526-EXP from the Spanish Ministry of Economy, Industry and Competitiveness.