On a nonlinear rod exhibiting only axial and bending deformations: mathematical modeling and numerical implementation

Abstract

In this work, we present the mathematical formulation and the numerical implementation of a new model for initially straight, transversely isotropic rods. By adopting a configuration space that intrinsically avoids shear deformations and by systemically neglecting the energetic contribution due to torsion, the proposed model admits an unconstrained variational statement. Moreover, as the natural state of the rod is the trivial one and the mechanical properties are homogeneous on the cross section, the need for pull-back and push-forward operations in the formulation is totally circumvented. These features impose, however, some smoothness requirements on the stored energy function that need to be carefully considered when adopting general constitutive models. In addition to introducing the rod model, we propose a spatial discretization with the finite element method, and a time integration with a hybrid, implicit scheme. To illustrate the favorable features of the new model, we provide results corresponding to numerical simulations for plane and three-dimensional problems that are investigated in the static and dynamic settings. Finally, and to put the presented ideas in a suitable context, we compare solutions obtained with the new model against a rod model that allows for torsion and shear.

1 Introduction

Rods undergoing stretching, shearing, bending and torsion have been subject of a very active research starting from the seminal contributions by Kirchhoff [25], Clebsch [10], Cosserat and Cosserat [11] and Love [30]. As the topic has a long tradition in applied mathematics and mechanics, there is already a plethora of excellent references dealing with this subject that certainly is very rich (see for instance the works by Antman [1,2,3], Simó [45], Singer and colleagues [24, 28, 46], Gorely and coworkers [20, 27, 31] and Audoly and collaborators [4, 6, 7], inter alia). Among all the models available in the literature, those rods that deform without shearing play a special role (see for instance O’Reilly [39] and/or Giusteri and Fried [18] for a deep exposure of the subject). As a consequence of this particular deformation behavior, the cross section, which is assumed to remain flat, is perpendicular to the tangent vector associated with the curve that describes the rod axis. Employing this hypothesis, we have, in the linear context, two classical models that are well known: the one by Euler and Bernoulli and the one by Rayleigh [23], where the latter does consider the inertial contribution due to the rotation of the cross section. Going further to the nonlinear realm, we have the generalization of the Rayleigh model, the so-called Kirchhoff rod that represents a point of departure for this work.

Non-sheareable rods have been deeply investigated from a theoretical point of view. However, there is a renewed interest with respect to their application in very dissimilar research areas such as the analysis of DNA molecules [5, 43, 44], motion prediction of human hair [7, 26], dynamics of cables [8, 12], mechanical analysis of Möbius bands [36], dynamics of drill strings [47], interactive simulations [40], bifurcation analysis of tethers [48], motion prediction of vortical filaments [13], stability of elastic knots [4, 24], clambering plants [32], morphoelastic rods [29, 37], stents and catheters [49], cooking-induced deformation of spaghetti [19] and others. Obtaining closed analytical solutions is in general not possible, mainly due to the complex nature of the governing equations. Therefore, in most of the cases the only way to obtain approximate solutions is through numerical methods. This fact has motivated the development of new computational approaches for elastic rods exhibiting no shear deformation [9, 21, 22, 34, 50,51,52].

Due to the non-integrable nature of vanishing shear deformations, it is well known that formulating the governing equation of this class of rods through a truly unconstrained variational statements is not possible, especially for dynamic problems. Moreover, we showed in a recent work [41] that the formulation of a nonlinear rod model relying on a unconstrained variational statement is in principle affordable, but some further simplifying hypotheses are still necessary. However, neither such a mathematical formulation nor its numerical implementation is available in the literature to the best of our knowledge. In connection with this, we describe as follows the main novelties of this work.

In the first place, we present a new structural model for initially straight, transversely isotropic rods, which disregards the energetic contribution of torsion and exhibits no shear deformation. Such a model relies on a configuration space which, by construction, completely prevents shearing and thus results in an unconstrained approach. As the natural state of the rod is the trivial one, the pull-back and push-forward operations are completely circumvented, and thus, no introduction of an exogenous triad, which would require a non-local constrained approach, is necessary. Discarding the torsional contribution to the potential and kinetic energies avoids the inclusion of further non-integrable constraints such as the one required to enforce the non-twisting condition [15]. Moreover, as a direct consequence of the adopted kinematic and energetic settings, this new model fully admits a variational statement, keeps the problem in its minimal irreducible representation, and is in fact an specialization of the static and dynamic variational principles for nonlinear Kirchhoff rods previously developed by the authors [41]. The proposed model can be also considered as the non-shearable counterpart of the torsion-free beam model presented in [42]. This contribution turns to be a main innovation of the current work.

In the second place, we present an strategy for the numerical implementation of the new model proposed. To this end, we consider for the spacial discretization a finite element approach that relies on a $C^1$-class approximation. It is worthwhile to mention that through this approximation, the frame-invariance and path-independence of the discrete model are identically retained. For the time discretization, we consider a hybrid combination of the midpoint and trapezoidal rules. Such an implicit hybrid combination, while being very simple, exhibits exact preservation of the linear and angular momenta as well as an excellent energy behavior and thus is good and robust enough for the purposes of the current work. The resulting discrete model can then be completely described by primal variables (there is no necessity of dual variables such as Lagrange multipliers) and thus results to be very robust with respect to conditioning issues typically exhibited by mixed primal-dual formulations. We use then the proposed numerical implementation of the new model to investigate static and dynamic, planar and three-dimensional cases with which we show its very promising features. This means that the proposed model has low computational cost, is robust and turns out to be suitable for real-time applications. In addition, and to put our ideas into a adequate context, we include results obtained with a nonlinear beam model that allows for shear and torsion. These examples will show numerically that the solution of that nonlinear beam model converges to the one of the newly proposed model for increasing shear stiffness. This contribution represents another contribution of the current work.

The remainder of this manuscript is as follows. Section 2 presents the mathematical formulation of the new model. This includes the kinematic description and the derivation of suitable deformation measures. Then, it introduces the potential and kinetic energies and, after applying a variational argument, concludes with the weak and strong forms of the governing equations. Section 3 deals with the discretization of the model in space and time. Section 4 presents numerical investigations of increasing complexity that show the promising features of the proposed model. To close the work, Sect. 5 presents concluding remarks.

2 Model description and formulation

2.1 Admissible motions

To derive the model’s governing equations of a rod that is transversely isotropic, disregards the energetic contributions of torsion, and does not exhibit shear deformation, we first need to introduce a suitable geometric setting. Instead of using a constrained theory, we select a configuration space that intrinsically disables shear deformation. In this way, and in contrast to other Kirchhoff rod theories, we will completely avoid the formulation of a primal-dual problem that requires Lagrange multipliers.

Let $\{\varvec{E}_i\}_{i=1}^3$ denote the canonical Cartesian basis of $\mathbb {R}^3$ and define $\varvec{N}\equiv \varvec{E}_1$, $\varvec{B}\equiv \varvec{E}_2$, $\varvec{D}\equiv \varvec{E}_3$. With this notation, we define the configuration manifold of torsion-free, transversely isotropic Kirchhoff rods to be the set

$$\begin{aligned} {\mathscr {Q}}:= \left\{ \varvec{\varphi }\in [C^{2}(0,S)]^3,\ |\varvec{\varphi }'(s)|>0,\ \varvec{\varphi }(0) = \varvec{0}, \ \varvec{\varphi }'(0) = \varvec{D} \right\} , \end{aligned}$$

(1)

where S is the length of the rod and we have selected clamped boundary conditions at $s = 0$ for concreteness. A motion of the rod is a curve of configurations parameterized by time and we write $\varvec{\varphi }_t(s)=\varvec{\varphi }(s,t)$ for $s\in [0,S]$ and $t\in [0,T]$. Moreover, we reduce the class of admissible motions to the ones that have as initial configuration

$$\begin{aligned} \varvec{\varphi }_0(s) = s\varvec{D}\ . \end{aligned}$$

(2)

The constant vector $\varvec{D}$ belongs to the unit sphere $S^2:=\left\{ \varvec{d}\in \mathbb {R}^3\; |\; \varvec{d}\cdot \varvec{d}= 1\right\}$, and thus, it is referred to as a director. In addition, any configuration $\varvec{\varphi }\in {\mathscr {Q}}$ defines a curve $s\mapsto \varvec{\varphi }(s)$ in ${\mathbb {R}}^3$ with Frenet–Serret triad

$$\begin{aligned} \varvec{d}:= \frac{\varvec{\varphi }'}{|\varvec{\varphi }'|} \ , \qquad \varvec{n}:= \frac{\varvec{\varphi }''}{|\varvec{\varphi }''|} \ , \qquad \varvec{b}:= \frac{\varvec{\varphi }'\times \varvec{\varphi }''}{|\varvec{\varphi }'\times \varvec{\varphi }''|} , \end{aligned}$$

(3)

that satisfy $\varvec{b}= \varvec{d}\times \varvec{n}$. We note that the triad is ill-defined in those points where $|\varvec{\varphi }''|=0$.

Remark 1

In all the regular points of a configuration, we can define the map

$$\begin{aligned} \varvec{\chi }= \varvec{n}\otimes \varvec{E}_1 + \varvec{b}\otimes \varvec{E}_2 + \varvec{d}\otimes \varvec{E}_3 \end{aligned}$$

(4)

with $\varvec{\chi }:\mathbb {R}^3\rightarrow \mathbb {R}^3$ and inverse

$$\begin{aligned} \varvec{\chi }^{-1}:= \varvec{E}_1\otimes \varvec{n}+ \varvec{E}_2\otimes \varvec{b}+ \varvec{E}_3\otimes \varvec{d}=\varvec{\chi }^T \end{aligned}$$

(5)

that satisfies $\det (\varvec{\chi })=1$ and thus being an element of the Lie group SO(3). This rotation can be employed to pull-back and push-forward vectors and tensors defined at regular points of the current configuration to the reference one, and vice versa. Existing rod models in the literature employ instead a rotation without drill [8, 9, 21, 22, 33,34,35]. However, such an approach introduces a fictitious torsion and can only be applied in an incremental form; for more details see [41]. To avoid all the aforementioned problems, we argue that a rod model based on the configuration space (1) must completely circumvent the use of this kind of mappings between frames, showing in the remainder of this work that this is indeed possible.

2.2 Strain measures, their variations and the linearized strain operator

Next, we investigate all possible ways in which a rod of the type considered before can deform. The deformation modes that need to be properly captured are those related to the elongation and curvature’s change of the rod. An important fact to be noted is that due to the transversal isotropy, the strain measures are required to be independent of the orientation of the rod’s cross section. This suggests indeed that it is possible to find scalar measures for the axial and bending deformations of the rod.

Two spatial strain measures are introduced first. They are defined as

$$\begin{aligned} \varvec{\epsilon } := ( |\varvec{\varphi }'|-1)\, \varvec{d} \ , \qquad \varvec{\kappa } := \varvec{d}\times \varvec{d}' \end{aligned}$$

(6)

and gauge the axial and bending deformations, respectively. For future reference, let us note that the derivative of the director can also be related to the covariant derivative on the sphere

$$\begin{aligned} \varvec{d}'= \frac{1}{|\varvec{\varphi }'|} \nabla _{\varvec{d}'}\varvec{\varphi }'= \frac{1}{|\varvec{\varphi }'|} (\varvec{I}-\varvec{d}\otimes \varvec{d})\varvec{\varphi }''. \end{aligned}$$

(7)

The bending deformation can be given the following geometric interpretation: given the curve $s\mapsto \varvec{d}(s)$ on an oriented surface, here $S^2$, the curvature can be split into a normal component and a tangent one. Since $\varvec{d}(s)$ lives in $S^2$, the normal component has always unit modulus, i.e., the intrinsic curvature of the sphere. The tangent component is precisely $\varvec{\kappa }$, as defined in Eq. (6), and measures how $\varvec{d}$ fails to be geodesic [38].

Next, we strive to obtain the equations of motion of the rod through a variational argument. Prior to this calculation, we need the expressions for the strain measure variations. These are summarized in the following lemma.

Lemma 1

Let $\delta \varvec{\varphi }$ be an admissible variation of the rod configuration. Then, the variations of the two strain measures are:

$$\begin{aligned} \delta \varvec{\epsilon }&= \left( \varvec{I} - \frac{1}{|\varvec{\varphi }'|}(\varvec{I}-\varvec{d}\otimes \varvec{d})\right) \delta \varvec{\varphi }' , \end{aligned}$$

(8a)

$$\begin{aligned} \delta \varvec{\kappa }&= \frac{-1}{|\varvec{\varphi }'|^2} \left( 2 \nabla _{\varvec{d}'}\varvec{\varphi }'\odot \varvec{d}+ (\varvec{d}\cdot \varvec{\varphi }'')(\varvec{I}-\varvec{d}\otimes \varvec{d}) \right) \delta \varvec{\varphi }' + \frac{1}{|\varvec{\varphi }'|} (\varvec{I}-\varvec{d}\otimes \varvec{d}) \delta \varvec{\varphi }'', \end{aligned}$$

(8b)

with

$$\begin{aligned} \varvec{a}\odot \varvec{b} := \frac{1}{2} ( \varvec{a}\otimes \varvec{b} + \varvec{b}\otimes \varvec{a}) \end{aligned}$$

(9)

being the symmetrized outer product of two vectors $\varvec{a},\varvec{b}\in {\mathbb {R}}^3$.

Proof

The proof is based on a systematic application of the rules of differentiation. Details are hence omitted.

$\square$

It is convenient at this stage to introduce the linearized strain operator $\mathscr {B}(\varvec{\varphi }',\varvec{\varphi }'')$ such that

$$\begin{aligned} \begin{Bmatrix} \delta \varvec{\epsilon } \\ \delta \varvec{\kappa } \end{Bmatrix} = \mathscr {B}(\varvec{\varphi }',\varvec{\varphi }'')\; \delta \varvec{\varphi }\ , \end{aligned}$$

(10)

and defined as

$$\begin{aligned} \begin{aligned} \mathscr {B}(\varvec{\varphi }',\varvec{\varphi }''):=&\begin{bmatrix} \varvec{I} - \frac{1}{|\varvec{\varphi }'|}(\varvec{I}-\varvec{d}\otimes \varvec{d})&{}\quad \varvec{0}\\ \frac{-1}{|\varvec{\varphi }'|^2} \left( 2 \nabla _{\varvec{d}'}\varvec{\varphi }'\odot \varvec{d}+ (\varvec{d}\cdot \varvec{\varphi }'')(\varvec{I}-\varvec{d}\otimes \varvec{d}) \right) &{}\quad \frac{1}{|\varvec{\varphi }'|} (\varvec{I}-\varvec{d}\otimes \varvec{d}) \\ \end{bmatrix} \begin{Bmatrix} {\mathscr {D}}_1 \\ {\mathscr {D}}_2 \end{Bmatrix},\\ \end{aligned} \end{aligned}$$

(11)

where we introduced the differential operators

$$\begin{aligned} {\mathscr {D}}_\alpha = \frac{\partial ^\alpha }{\partial s^\alpha } \qquad \alpha =1,2\ . \end{aligned}$$

(12)

2.3 Stored energy and stress measures

Next, we investigate the hyperelastic models for initially straight, torsion-free, non-shearable, transversely isotropic rods and identify the stress measures that are work-conjugate to the strain measures previously introduced.

To this end, we assume the existence of an objective stored energy density function, ${\tilde{\mathscr {U}}}:{\mathbb {R}}^3\times {\mathbb {R}}^3\rightarrow {\mathbb {R}}$ that is taken to be identical for all cross sections, for simplicity. The objectivity conditions requires that for every pair of strains $\varvec{\epsilon },\varvec{\kappa }$ and every $\varvec{Q}\in SO(3)$

$$\begin{aligned} {\tilde{\mathscr {U}}}( \varvec{\epsilon }, \varvec{\kappa }) = {\tilde{\mathscr {U}}}( \varvec{Q} \varvec{\epsilon }, \varvec{Q} \varvec{\kappa })\ . \end{aligned}$$

(13)

A standard argument employed to guarantee the objectivity of this energy function is to introduce two other stored energy functions $\bar{{\mathscr {U}}}_a,\bar{{\mathscr {U}}}_b:{\mathbb {R}}^+\cup \{0\}\rightarrow {\mathbb {R}}$ and postulate

$$\begin{aligned} {\tilde{\mathscr {U}}}( \varvec{\epsilon } , \varvec{\kappa }) = \bar{{\mathscr {U}}}_a (\epsilon ) + \bar{{\mathscr {U}}}_b(\kappa ) \qquad \mathrm {with} \qquad \epsilon = | \varvec{\epsilon } | , \ \kappa = | \varvec{\kappa }| \ . \end{aligned}$$

(14)

These considerations facilitate the formulation of stored energy functions because all that is required is to postulate two scalar functions defined in the non-negative real numbers. These potentials cannot be arbitrary either. In order to guarantee the existence of solutions, elastic problems impose certain coerciveness conditions on their stored energy functions. Without claiming that it is enough to ensure existence in the finite strain range, we adopt the following quadratic stored energy functions borrowed from the small-strain rod theories

$$\begin{aligned} \bar{{\mathscr {U}}}_a(\epsilon ) = \frac{1}{2}\, EA\, \epsilon ^2\ , \qquad \bar{{\mathscr {U}}}_b(\kappa ) = \frac{1}{2}\, EI\, \kappa ^2\ , \end{aligned}$$

(15)

where E is the Young modulus of the rod material, A is cross section area and I the inertia with respect to an axis contained in the cross section.

With these definitions at hand, we can define the potential energy functional of the rod as

$$\begin{aligned} \mathscr {U}(\varvec{\varphi }) = \int _0^S{\tilde{\mathscr {U}}}( \varvec{\epsilon }, \varvec{\kappa }) \,\mathrm {d} s. \end{aligned}$$

(16)

Hence, its variation follows systematically as

$$\begin{aligned} \begin{aligned} D\mathscr {U}(\varvec{\varphi }) \cdot \delta \varvec{\varphi }&= \int _0^S\delta {\tilde{\mathscr {U}}}(\varvec{\epsilon },\varvec{\kappa }) \,\mathrm {d}s\\&= \int _0^S\left( \frac{\bar{{\mathscr {U}}}_a'(\epsilon )}{\epsilon } \varvec{\epsilon } \cdot \delta \varvec{\epsilon } + \frac{\bar{{\mathscr {U}}}_b'(\kappa )}{\kappa } \varvec{\kappa } \cdot \delta \varvec{\kappa } \right) \,\mathrm {d}s\\&= \int _0^S\delta \varvec{\varphi }\cdot (\mathscr {B}(\varvec{\varphi }',\varvec{\varphi }'')^T\varvec{\sigma })\,\mathrm {d}s, \end{aligned} \end{aligned}$$

(17)

where the stress measures, energy-conjugate to $\varvec{\epsilon }$ and $\varvec{\kappa }$,

$$\begin{aligned} \varvec{\sigma }:= \begin{Bmatrix} \varvec{n} \\ \varvec{m} \end{Bmatrix} = \begin{Bmatrix} \displaystyle \frac{\bar{{\mathscr {U}}}_a'(\epsilon )}{\epsilon } \varvec{\epsilon }\\ \displaystyle \frac{\bar{{\mathscr {U}}}_b'(\kappa )}{\kappa } \varvec{\kappa }\\ \end{Bmatrix} \end{aligned}$$

(18)

naturally arise. The stress measure $\varvec{n}$ can be identified with the axial force, parallel to $\varvec{d}$ by definition, whereas the bending moment $\varvec{m}$ is contained in the cross section of the rod and has the direction of the bending strain. In fact, for the stored energy function (14) it follows that the stress measures are

$$\begin{aligned} \varvec{n} = EA\, \varvec{\epsilon } \ , \qquad \varvec{m} = EI\, \varvec{\kappa }\ , \end{aligned}$$

(19)

extending the classical results for the axial force and bending moment in a rod. In general, we emphasize that the stored energy functions $\bar{{\mathscr {U}}}_a$ and $\bar{{\mathscr {U}}}_b$ must be such that

$$\begin{aligned} \lim _{\epsilon \rightarrow 0} \frac{\bar{{\mathscr {U}}}_a'(\epsilon )}{\epsilon } = c_\epsilon>0 \qquad \mathrm {and} \qquad \lim _{\kappa \rightarrow 0} \frac{\bar{{\mathscr {U}}}_b'(\kappa )}{\kappa } = c_\kappa >0\ .\end{aligned}$$

(20)

In fact, this condition can be used to define the axial and bending stress measures to be always identically zero when the strains $\varvec{\epsilon }$ and $\varvec{\kappa }$ vanish, respectively, as automatically done by the definitions in Eq. (19).

2.4 Kinetic energy and its variation

Next, we define a kinetic energy density of the form

$$\begin{aligned} {\tilde{\mathscr {T}}}({\dot{\varvec{\varphi }}}, {\dot{\varvec{d}}}; s):= \frac{1}{2}A_{\varrho }|{\dot{\varvec{\varphi }}}|^2+\frac{1}{2}I_{\varrho }|{\dot{\varvec{d}}}|^2 = \frac{1}{2}{\dot{\varvec{\varphi }}}\cdot \mathscr {M}(\varvec{\varphi }'){\dot{\varvec{\varphi }}}, \end{aligned}$$

(21)

where $A_{\varrho }$ is the mass per unit length and $I_{\varrho }$ is the inertia density around any axis contained in the rod section (by the assumption of isotropy, all rotational inertias are identical). Based on the previous relation, the mass operator $\mathscr {M}(\varvec{\varphi }')$ must be

$$\begin{aligned} \mathscr {M}(\varvec{\varphi }'):= A_{\varrho }\varvec{I}+\mathscr {D}_1^TI_{\varrho }\frac{1}{|\varvec{\varphi }'|^{2}}(\varvec{I}-\varvec{d}\otimes \varvec{d})\mathscr {D}_1. \end{aligned}$$

(22)

On this basis, we can define the kinetic energy functional as

$$\begin{aligned} \mathscr {T}(\varvec{\varphi }) = \int _0^S{\tilde{\mathscr {T}}}({\dot{\varvec{\varphi }}}, {\dot{\varvec{d}}}; s)\,\mathrm {d}s\end{aligned}$$

(23)

and its systematic variation is as follows:

$$\begin{aligned} \begin{aligned} \delta \mathscr {T}(\varvec{\varphi })&= \int _0^S\delta {\tilde{\mathscr {T}}}({\dot{\varvec{\varphi }}}, {\dot{\varvec{d}}}; s)\,\mathrm {d}s\\&=-\int _0^S\delta \varvec{\varphi }^{a} \mathscr {M}(\varvec{\varphi }')_{ab}\ddot{\varvec{\varphi }}^{b} + {\hat{\Gamma }}(\mathscr {M}(\varvec{\varphi }'))_{abc} {\dot{\varvec{\varphi }}}^{b} {\dot{\varvec{\varphi }}}^{c})\,\mathrm {d}s\\&=-\int _0^S\delta \varvec{\varphi }\cdot \mathscr {M}(\varvec{\varphi }'){\hat{\nabla }}_{{\dot{\varvec{\varphi }}}}{\dot{\varvec{\varphi }}}\,\mathrm {d}s, \end{aligned} \end{aligned}$$

(24)

where ${\hat{\nabla }}_{{\dot{\varvec{\varphi }}}}{\dot{\varvec{\varphi }}}$ denotes the field covariant derivative and a, b, and c are indices taking the values 1,2, and 3, for which the convention of repeated indices applies. The term $\mathscr {M}(\varvec{\varphi }'){\hat{\nabla }}_{{\dot{\varvec{\varphi }}}}{\dot{\varvec{\varphi }}}$ contains not only the contribution $\mathscr {M}(\varvec{\varphi }')\ddot{\varvec{\varphi }}$, but also the contribution ${\hat{\Gamma }}(\mathscr {M}(\varvec{\varphi }'))({\dot{\varvec{\varphi }}},{\dot{\varvec{\varphi }}})$. The corresponding connection coefficients are uniquely defined by

$$\begin{aligned} 2\,{\hat{\Gamma }}(\mathscr {M}(\varvec{\varphi }'))_{abc} := \partial _{a}\mathscr {M}(\varvec{\varphi }')_{bc}+\partial _{b}\mathscr {M}(\varvec{\varphi }')_{ac}-\partial _{c}\mathscr {M}(\varvec{\varphi }')_{ab}, \end{aligned}$$

(25)

such that ${\hat{\nabla }}_{{\dot{\varvec{\varphi }}}} \mathscr {M}(\varvec{\varphi }')= \varvec{0 }$.

Equation (24) follows from the fact that

$$\begin{aligned} \int _0^T\!\!\delta {\tilde{\mathscr {T}}}({\dot{\varvec{\varphi }}}, {\dot{\varvec{d}}}; s)\,\mathrm {d}t&= \int _0^T\!\!\delta \left( \frac{1}{2}{\dot{\varvec{\varphi }}}\cdot \mathscr {M}(\varvec{\varphi }'){\dot{\varvec{\varphi }}}\right) \!\!\,\mathrm {d}t\end{aligned}$$

(26a)

$$\begin{aligned}&= \int _0^T\!\!\!\delta \varvec{\varphi }\cdot \left( \frac{1}{2}\mathscr {D}_1^T(\partial _{\varvec{\varphi }'}({\dot{\varvec{\varphi }}}\cdot \mathscr {M}(\varvec{\varphi }'){\dot{\varvec{\varphi }}})-\mathscr {M}(\varvec{\varphi }')\ddot{\varvec{\varphi }}-\dot{\mathscr {M}}(\varvec{\varphi }'){\dot{\varvec{\varphi }}})\right) \!\!\,\mathrm {d}t\end{aligned}$$

(26b)

$$\begin{aligned}&= -\int _0^T\!\! \delta \varvec{\varphi }^{a}(\mathscr {M}(\varvec{\varphi }')_{ab}\ddot{\varvec{\varphi }}^{b}+{\hat{\Gamma }}(\mathscr {M}(\varvec{\varphi }'))_{abc}{\dot{\varvec{\varphi }}}^{b}{\dot{\varvec{\varphi }}}^{c})\,\mathrm {d}t\end{aligned}$$

(26c)

for vanishing $\delta \varvec{\varphi }$ at $t=0$ and $t=T$. In addition, we have that

$$\begin{aligned} \dot{\mathscr {M}}(\varvec{\varphi }')= -\mathscr {D}_1^TI_{\varrho }(2 \frac{1}{|\varvec{\varphi }'|^{3}} (\varvec{d}\cdot {\dot{\varvec{\varphi }}}')(\varvec{I}-\varvec{d}\otimes \varvec{d})+ \frac{1}{|\varvec{\varphi }'|^{2}}({{\dot{\varvec{d}}}}\otimes {\varvec{d}}+{\varvec{d}}\otimes {{\dot{\varvec{d}}}}))\mathscr {D}_1\end{aligned}$$

(27)

and

$$\begin{aligned} \frac{1}{2}\mathscr {D}_1^T(\partial _{\varvec{\varphi }'}({\dot{\varvec{\varphi }}}\cdot \mathscr {M}(\varvec{\varphi }'){\dot{\varvec{\varphi }}})) = -2\mathscr {D}_1^T\left( I_{\varrho }\frac{1}{|\varvec{\varphi }'|^{3}} (\varvec{I}-\varvec{d}\otimes \varvec{d})\odot ({\dot{\varvec{\varphi }}}'\otimes \varvec{d}){\dot{\varvec{\varphi }}}'\right) , \end{aligned}$$

(28)

which are required to compute ${\hat{\Gamma }}(\mathscr {M}(\varvec{\varphi }'))({\dot{\varvec{\varphi }}},{\dot{\varvec{\varphi }}})$.

3 Variational derivation of the equations of motion

Hamilton’s principle of least action states that the governing equations of a mechanical system are the Euler–Lagrange equations of the action functional

$$\begin{aligned} \mathscr {S}= \int _0^T\mathscr {L}(\varvec{\varphi }',{\dot{\varvec{\varphi }}},\varvec{\varphi }'',{\dot{\varvec{\varphi }}}';t)\,\mathrm {d}t= \int _0^T\left( \mathscr {T}(\varvec{\varphi }',{\dot{\varvec{\varphi }}},{\dot{\varvec{\varphi }}}';t)-\mathscr {U}(\varvec{\varphi }',\varvec{\varphi }'';t)\right) \,\mathrm {d}t \end{aligned}$$

(29)

whose solutions $\varvec{\varphi }\in {\mathscr {Q}}$ render the action stationarity. The main result is summarized next:

Theorem 1

The (weak) Euler–Lagrange equations corresponding to the functional (29) are:

$$\begin{aligned} \int _0^S \delta \varvec{\varphi }\cdot \left( \mathscr {M}(\varvec{\varphi }'){\hat{\nabla }}_{{\dot{\varvec{\varphi }}}}{\dot{\varvec{\varphi }}}+\mathscr {B}(\varvec{\varphi }',\varvec{\varphi }'')^T\varvec{\sigma }\right) \,\mathrm {d}s= 0 . \end{aligned}$$

(30)

Proof

The proof of the theorem follows from the systematic calculation of $\delta \mathscr {S}$, namely,

$$\begin{aligned} \delta \mathscr {S}= \int _0^T\delta \mathscr {L}(\varvec{\varphi }',{\dot{\varvec{\varphi }}},\varvec{\varphi }'',{\dot{\varvec{\varphi }}}';t)\,\mathrm {d}t\end{aligned}$$

(31)

with

$$\begin{aligned} \delta \mathscr {L}(\varvec{\varphi }',{\dot{\varvec{\varphi }}},\varvec{\varphi }'',{\dot{\varvec{\varphi }}}';t) = \int _0^S\delta \varvec{\varphi }\cdot \left( \mathscr {M}(\varvec{\varphi }'){\hat{\nabla }}_{{\dot{\varvec{\varphi }}}}{\dot{\varvec{\varphi }}}+\mathscr {B}(\varvec{\varphi }',\varvec{\varphi }'')^T\varvec{\sigma }\right) \,\mathrm {d}s, \end{aligned}$$

(32)

for $\delta \varvec{\varphi }$ satisfying $\delta \varvec{\varphi }(s,0)=\delta \varvec{\varphi }(s,T)=\varvec{0 }$. $\square$

From the variational equation (30), we can recover the strong form of the equations of motion of the Kirchhoff rod as presented next.

Theorem 2

The strong form of the equations of motion is

$$\begin{aligned} \varvec{n}'+ \left( \frac{1}{|\varvec{\varphi }'|}\varvec{d}\times \nabla _{\varvec{d}'}\varvec{m}\right) ' = \dot{\varvec{p}}+\left( \frac{1}{|\varvec{\varphi }'|}\varvec{d}\times \nabla _{{\dot{\varvec{d}}}}\varvec{\pi }\right) ' , \end{aligned}$$

(33)

with the stress measures $\varvec{n}$ and $\varvec{m}$ introduced above and the momenta $\varvec{p}$ and $\varvec{\pi }$ defined as $A_{\varrho }{\dot{\varvec{\varphi }}}$ and $I_{\varrho }{\dot{\varvec{d}}}$, respectively. Moreover, for the clamped-free case, we have the concomitant boundary conditions

$$\begin{aligned}&\varvec{\varphi }= \varvec{0},\ \varvec{\varphi }'=\varvec{D}\qquad&\mathrm {at}&\quad s=0\ , \end{aligned}$$

(34a)

$$\begin{aligned}&\varvec{n}+\frac{1}{|\varvec{\varphi }'|}\varvec{d}\times \left( \nabla _{\varvec{d}'}\varvec{m} - \nabla _{{\dot{\varvec{d}}}}\varvec{\pi }\right) = \varvec{0 }, \ \frac{1}{|\varvec{\varphi }'|}\varvec{d}\times \varvec{m}=\varvec{0 }\qquad&\mathrm {at}&\quad s=S, \end{aligned}$$

(34b)

and the initial conditions

$$\begin{aligned} \varvec{\varphi }(s,0) = \varvec{\varphi }_0(s)\ , \quad {\dot{\varvec{\varphi }}}(s,0) = \varvec{v}_0(s) \qquad \mathrm {for} \quad s\in [0,S]\ . \end{aligned}$$

(35)

Proof

The proof of the theorem can be found in Romero and Gebhardt [41], and thus, we omit further details. $\square$

4 Discretization

The initial boundary value problem defined by Eqs. (33)–(35) is difficult to solve in general and often one has to resort to numerical methods to obtain approximations to the solution. In order to obtain the latter, the governing equations of the rod have to be discretized both in space and time. In this section, we provide the details of both discretizations, starting from a finite element spatial projection and completing the process with a mixed midpoint/trapezoidal time integration.

4.1 Spatial discretization

We start by defining a finite element spatial semi-discretization. For that, consider an approximation of the rod configuration to be of the form

$$\begin{aligned} \varvec{\varphi }(s,t) \approx \varvec{\varphi }_h(s,t) = \sum _{a=1}^{\nu } \left( N_a^x(s)\, \varvec{x}_a(t) + N_a^d(s)\, \varvec{d}_a(t) \right) = \sum _{a=1}^{2\nu } {N}_a(s)\, \varvec{q}_a(t) \ , \end{aligned}$$

(36)

where $\nu$ is the number of nodes, $N_a^x,N_a^d,{N}_a:[0,S]\rightarrow {\mathbb {R}}$ are $C^1$ shape functions, $\varvec{x}_a,\varvec{d}_a,\varvec{q}_a:[0,T]\rightarrow {\mathbb {R}}^3,S^2, \mathbb {R}^3\times S^2$ are the nodal positions, directors, and generalized coordinates, functions that must be twice differentiable in time. For convenience, we introduce the matrix $\varvec{N}$ of shape functions and the global nodal vector $\varvec{q}(t)$ such that

$$\begin{aligned} \varvec{\varphi }_h(s,t) = \sum _{a=1}^{2\nu } {N}_a(s)\; \varvec{q}_a(t) = \varvec{N}(s)\;\varvec{q}(t) . \end{aligned}$$

(37)

The nodal values of the director determine the discrete director field at these points. At interior points of the elements, this field can be calculated with the expression

$$\begin{aligned} \varvec{d}_h(s,t)= \frac{\varvec{\varphi }_h'(s,t)}{|\varvec{\varphi }_h'(s,t)|}\ . \end{aligned}$$

(38)

Using the interpolations (37) and (38), the weak form (30) can be approximated as indicated next:

Lemma 2

The governing equations of motion, semi-discretized in space, are:

$$\begin{aligned} \int _0^S \left( \varvec{M}(\varvec{q})\nabla _{{\dot{\varvec{q}}}} {\dot{\varvec{q}}}+ \varvec{B}(\varvec{q})^T\varvec{\sigma }\right) \,\mathrm {d}s= \varvec{0 } \end{aligned}$$

(39)

with the concomitant initial conditions $\varvec{q}(0)=\varvec{q}_0$ and ${\dot{\varvec{q}}}(0)={\dot{\varvec{q}}}_0$.

Proof

Using the Galerkin approximation (36) of the configuration, the inertial contribution of Eq. (30) can be rewritten as

$$\begin{aligned} \delta \varvec{\varphi }_h \cdot \mathscr {M}(\varvec{\varphi }_h') {\hat{\nabla }}_{{\dot{\varvec{\varphi }}}_h} {\dot{\varvec{\varphi }}}_h = \delta \varvec{q}\cdot \varvec{M}(\varvec{q})\nabla _{{\dot{\varvec{q}}}} {\dot{\varvec{q}}}, \end{aligned}$$

(40)

with

$$\begin{aligned} \varvec{M}(\varvec{q})\nabla _{{\dot{\varvec{q}}}} {\dot{\varvec{q}}}= \varvec{N}^T\!\! \mathscr {M}(\varvec{\varphi }_h') \varvec{N}\ddot{\varvec{q}}+ \varvec{N}^T\!\! \dot{\mathscr {M}}(\varvec{\varphi }_h') \varvec{N}{\dot{\varvec{q}}}- \frac{1}{2}\varvec{N}'^T \! ( \partial _{\varvec{\varphi }_h'} ({\dot{\varvec{\varphi }}}_h \cdot \mathscr {M}(\varvec{\varphi }_h') {\dot{\varvec{\varphi }}}_h)) . \end{aligned}$$

(41)

The right-hand side of (40) includes not only the contribution $\varvec{M}(\varvec{q})\ddot{\varvec{q}}$ that corresponds to the first term of (41), but also the contribution $\varvec{\Gamma }(\varvec{q})({\dot{\varvec{q}}},{\dot{\varvec{q}}})$ that corresponds to the second and third term. The corresponding connection coefficients are uniquely defined by

$$\begin{aligned} 2\Gamma _{abc} := \partial _{a}M_{bc}+\partial _{b}M_{ac}-\partial _{c}M_{ab} \end{aligned}$$

(42)

such that $\nabla _{a}(M_{bc}) = \varvec{0 }$. For the term coming from the potential energy, we have, using again Galerkin’s approximation (36), that

$$\begin{aligned} \delta \varvec{\varphi }_h\cdot \mathscr {B}^T(\varvec{\varphi }'_h,\varvec{\varphi }''_h)\, \varvec{\sigma }= \delta \varvec{q}\cdot \varvec{B}^T(\varvec{q})\,\varvec{\sigma }, \end{aligned}$$

(43)

with

$$\begin{aligned} \varvec{B}(\varvec{q}):= & {} \begin{bmatrix} \varvec{I} - \frac{1}{|\varvec{\varphi }_h'|}(\varvec{I} - \varvec{d}_h\otimes \varvec{d}_h) &{}\quad \varvec{0 }\\ \frac{-1}{|\varvec{\varphi }_h'|^2} \left( 2 \nabla _{\varvec{d}_h'}\varvec{\varphi }_h'\odot \varvec{d}_h + (\varvec{d}_h\cdot \varvec{\varphi }_h'')(\varvec{I}-\varvec{d}_h\otimes \varvec{d}_h) \right) &{}\quad \frac{1}{|\varvec{\varphi }_h'|}(\varvec{I} - \varvec{d}_h\otimes \varvec{d}_h) \end{bmatrix}\begin{bmatrix} \varvec{N}'\\ \varvec{N}''\end{bmatrix} . \end{aligned}$$

(44)

Collecting the last two results, we obtain

$$\begin{aligned} \int _0^S \delta \varvec{\varphi }_h\cdot \left( \mathscr {M}(\varvec{\varphi }_h') {\hat{\nabla }}_{{\dot{\varvec{\varphi }}}_h} {\dot{\varvec{\varphi }}}_h + \mathscr {B}(\varvec{\varphi }',\varvec{\varphi }'')^T\varvec{\sigma }\right) \,\mathrm {d}s= \delta \varvec{q}\cdot \int _0^S\left( \varvec{M}(\varvec{q})\nabla _{{\dot{\varvec{q}}}}{\dot{\varvec{q}}}+\varvec{B}(\varvec{q})^T\varvec{\sigma }\right) \,\mathrm {d}s. \end{aligned}$$

(45)

Since, according to Eq. (30), this expression must vanish for arbitrary configurations, the integral on the right-hand side must itself be zero and the proof follows. $\square$

To complete the spatial semidiscretization, it remains to select an appropriate shape function so as to define the finite element space. A standard choice for $C^1$ approximation functions is based on Hermite polynomials defined on a reference domain $[-1,+1]$ that are later pushed to each point of the rod. To do that, we consider the $\nu -1$ line elements defined by the $\nu$ nodes. Element e has two nodes with labels $I=e$ and $J=e+1$. A matrix of element shape functions can be defined as follows:

$$\begin{aligned} \varvec{n}_e(\xi ) = \begin{bmatrix} H_{x}^{-}(\xi )\;\varvec{I}&\frac{S_e}{2} H_{d}^{-}(\xi )\;\varvec{I}&H_{x}^{+}(\xi )\;\varvec{I}&\frac{S_e}{2} H_{d}^{+}(\xi )\;\varvec{I}\end{bmatrix}, \end{aligned}$$

(46)

where $\xi \in [-1,+1]$, $S_e$ is the length of element e and we have used the standard Hermite polynomials

$$\begin{aligned} \begin{aligned} H_{x}^{-}(\xi )&= +\frac{1}{4}(2+\xi )(1-\xi )^2 ,\\ H_{d}^{-}(\xi )&= +\frac{1}{4}(1+\xi )(1-\xi )^2 ,\\ H_{x}^{+}(\xi )&= +\frac{1}{4}(2-\xi )(1+\xi )^2 ,\\ H_{d}^{+}(\xi )&= -\frac{1}{4}(1-\xi )(1+\xi )^2 . \end{aligned} \end{aligned}$$

(47)

With these interpolation functions, a map $\varvec{\Phi }_e:[-1,+1]\rightarrow {\mathbb {R}}^3$ can be defined with the expression

$$\begin{aligned} \varvec{\Phi }_e(\xi ) = H_{x}^-(\xi )\, \varvec{x}_{0,I} + H_{x}^+(\xi )\, \varvec{x}_{0,J} + \frac{S_e}{2} H_{d}^-(\xi )\, \varvec{d}_{0,I} + \frac{S_e}{2} H_{d}^+(\xi )\, \varvec{d}_{0,J}\ , \end{aligned}$$

(48)

with $\varvec{x}_{0,I},\varvec{x}_{0,J}, \varvec{d}_{0,I}, \varvec{d}_{0,J}$ being the position and directors of the element nodes in the reference configuration, respectively. The map $\varvec{\Phi }_e$ can be used to evaluate the global shape function matrix $\varvec{N}$ or any other function defined on the reference arc-length interval [0, S]. To see that, note that the mapping $\xi \mapsto \varvec{\Phi }_e(\xi )$ defines the reference placement of the rod and a bijective application $\Sigma _e:[-1,+1]\rightarrow [s_{e},s_{e+1}]$, with $s_e$ corresponding to the s coordinate of node e, as

$$\begin{aligned} \Sigma _e(\xi ) = s_e + \int _{-1}^\xi |\partial _\xi \varvec{\Phi }_e(\xi )| \,\mathrm {d} \xi \ . \end{aligned}$$

(49)

Then, the restriction of $\varvec{N}$ to element e can be evaluated as

$$\begin{aligned} \left. \varvec{N} \right| _e (\Sigma _e(\xi )) = \varvec{n}_e(\xi )\ . \end{aligned}$$

(50)

Moreover, the integral of any function $f\circ \varvec{\varphi }$ defined on [0, S] can be approximated with the integral of $f\circ \varvec{\varphi }_h$ and then evaluated as

$$\begin{aligned} \int _0^S (f\circ \varvec{\varphi })(s) \,\mathrm {d} s \approx \sum _{e=1}^{\nu -1} \int _{-1}^{+1} (f\circ \varvec{\varphi }_h \circ \Sigma _e) (\xi ) \left| \partial _{\xi }{\varvec{\Phi }_e}\right| \,\mathrm {d} \xi \ . \end{aligned}$$

(51)

This kind of integrals, or rather their numerical quadrature, will be employed in Sect. 5 to evaluate the balance equation (30) once it is discretized in time.

4.2 Time integration

Next, we describe the integration algorithm adopted to compute the dynamic simulations presented in this paper. To such an end, let us consider the governing equations, which are semi-discretized in space, evaluated at time instant $t_{n+\frac{1}{2}}\in [t_n, t_{n+1}]$, i.e., 

$$\begin{aligned} \int _0^S\left( \left( \varvec{M}(\varvec{q})\nabla _{{\dot{\varvec{q}}}}{\dot{\varvec{q}}}\right) _{n+\frac{1}{2}} + \left( \varvec{B}(\varvec{q})^T\varvec{\sigma }\right) _{n+\frac{1}{2}} - \varvec{Q}^{\text {ext}}_{n+\frac{1}{2}} \right) \,\mathrm {d}s= \varvec{0 }\ , \end{aligned}$$

(52)

where for sake of completeness, we include external generalized forces, namely $\varvec{Q}^{\text {ext}}_{n+\frac{1}{2}}$.

In a recent publication [16], we showed that the hybrid combination of the midpoint and trapezoidal rules can render robust second-order, implicit integration methods that at least exactly preserve the linear and angular momenta. Therefore, this is the approach that we are going to pursue in the following. To approximate the inertial term $\left( \varvec{M}(\varvec{q})\nabla _{{\dot{\varvec{q}}}}{\dot{\varvec{q}}}\right) _{n+\frac{1}{2}}$, we employ an extended version of the midpoint rule, i.e., 

$$\begin{aligned} \left( \varvec{M}(\varvec{q})\nabla _{{\dot{\varvec{q}}}}{\dot{\varvec{q}}}\right) _{n+\frac{1}{2}} \approx \frac{\varvec{M}\left( \varvec{q}_{n+1}\right) {\dot{\varvec{q}}}_{n+1}-\varvec{M}\left( \varvec{q}_{n}\right) {\dot{\varvec{q}}}_{n}}{h}- \frac{1}{2}\left. \partial _{\varvec{\xi }}\left( {\dot{\varvec{q}}}_{n+\frac{1}{2}}\cdot \varvec{M}(\varvec{\xi })\cdot {\dot{\varvec{q}}}_{n+\frac{1}{2}}\right) \right| _{\varvec{\xi }=\varvec{q}_{n+\frac{1}{2}}}, \end{aligned}$$

(53)

where the first term approximates $\varvec{M}(\varvec{q})\ddot{\varvec{q}}+{\dot{\varvec{M}}}(\varvec{q}){\dot{\varvec{q}}}$, i.e., the time derivative of the momentum $\varvec{M}(\varvec{q}){\dot{\varvec{q}}}$, and the second term approximates the term arising from the derivative of the kinetic energy with respect to the vector of generalized coordinates. In addition, we consider the classical approximation formulas

$$\begin{aligned} \varvec{q}_{n+\frac{1}{2}}\approx \frac{\varvec{q}_n+\varvec{q}_{n+1}}{2}~~~\text {and}~~~{\dot{\varvec{q}}}_{n+\frac{1}{2}}\approx \frac{\varvec{q}_{n+1}-\varvec{q}_{n}}{h} . \end{aligned}$$

(54)

Meanwhile, to approximate the internal term $\left( \varvec{B}(\varvec{q})^T\varvec{\sigma }\right) _{n+\frac{1}{2}}$, we employ the trapezoidal rule, i.e., 

$$\begin{aligned} \left( \varvec{B}(\varvec{q})^T\varvec{\sigma }\right) _{n+\frac{1}{2}} \approx \frac{1}{2} \left( \varvec{B}\left( \varvec{q}_{n}\right) ^T\varvec{\sigma }_n+\varvec{B}\left( \varvec{q}_{n+1}\right) ^T\varvec{\sigma }_{n+1} \right) \ . \end{aligned}$$

(55)

Finally, collecting all pieces together, we can solve the fully discretized governing equations at time instant $t_{n+\frac{1}{2}}$ to find the discrete state of the system at time instant $t_{n+1}$, i.e., $\varvec{q}_{n+1}$ and ${\dot{\varvec{q}}}_{n+1}$. The discrete state of the system at time instant $t_n$ is known, namely $\varvec{q}_n$ and ${\dot{\varvec{q}}}_n$ are data of the problem. Moreover, ${\dot{\varvec{q}}}_{n+1}$ is computed as $\frac{2}{h}(\varvec{q}_{n+1}-\varvec{q}_{n})-{\dot{\varvec{q}}}_{n}$ once $\varvec{q}_{n+1}$ has been determined through an iterative method, for instance, the exact Newton–Raphson method.

Remark 2

The extension of the integration method to include holonomic constraints, e.g., for the length preservation of nodal directors or for setting a bundle of rods in hairlock modeling among others, requires to include the constraint equation at time instant $t_{n+1}$, $\varvec{g}(\varvec{q}_{n+1})=\varvec{0 }$, and to add the discrete restrictions forces at time instant $t_{n+\frac{1}{2}}$, $\partial _{\varvec{\xi }}(\varvec{g}(\varvec{\xi }))^T\mid _{\varvec{\xi }=\varvec{q}_{n+\frac{1}{2}}}\varvec{\lambda }_{n+\frac{1}{2}}$ to the fully discretized balance equations, which completes the discrete problem statement. Specifically, in the current numerical implementation, we consider an additive update of the nodal directors as elements of $\mathbb {R}^3$ that are constrained to preserve their unit length. Nevertheless, a multiplicative actualization in $S^2$, without constraints, is also possible for the current work as done in [42]. The latter are two alternative implementations of the time integration, with and without constraints. We stress, however, that the model itself is formulated without any constraint.

Fig. 1

Cantilever rod subject to an in-plane loading: sequence of deformed configurations

Fig. 2

Cantilever rod subject to an in-plane loading: deformed configuration of several TS rods with increasing values of GA versus the deformed configuration computed with the new formulation

5 Results

In this section, we present four numerical examples chosen to show the performance of the discretized rod models previously proposed. The first two examples belong to the realm of statics, whereas that the two remaining ones are dynamic simulations. In addition, and as a further reference, we also present the results obtained with a model for geometrically exact beams that includes both torsion and shear deformations, which we refer to as the TS rod. Details on this particular formulation are to be found in [14, 17]. Such a model asymptotically converges, for the case of slender beams, i.e., $L/\sqrt{I/A} \gg 1$, with vanishing shear deformations to the Kirchhoff rod model. We note that the aim of this section is not to provide a detailed comparison between models; however, this helps to put the proposed ideas into a suitable context and to provide further insights on the properties of the new formulation developed. Such a detailed comparative analysis falls outside the present scope and thus, is not carried out.

5.1 Statics

Here, we consider an initially straight, transversely isotropic rod with a length of $40~\mathrm {m}$ that initially lies along the direction $\varvec{D}= (1, 0, 0)$. The cross-sectional properties are $EA = 100~\mathrm {N}$ and $EI = 200~\mathrm {Nm^2}$. As boundary conditions, we consider that one end is fully fixed and the other one is free to move. For estimating the shear stiffness of the TS rod, we assume $\nu = 0.3$ and thus, $GA = 38.5~\mathrm {N}$. The rod is uniformly discretized into 40 two-node finite elements. The node corresponding to the fully fixed end is labeled 1 and the one corresponding to the free end has label 41. Finally, for all numerical examples the loading was applied, for plotting purposes, in 55 uniformly distributed load levels and a relative error-based tolerance of $10^{-10}$ has been set for the Newton iterations, which also applies to the examples dealing with dynamic problems.

5.1.1 Cantilever beam subject to an in-plane loading

This first example investigates the rod’s response to a single load given by $(-1.5, 0.5, 0)~\mathrm {N}$, which is applied to the node corresponding to the free end. Figure 1 shows a sequence of twelve successive deformed configurations. The results obtained with the model proposed are drawn in blue, whereas the results obtained with the TS model are plotted in red. The fully fixed end representing the clamped boundary conditions is indicated in green. The results reveal that the two-dimensional mechanical behavior of both models is, as expected, similar, but not identical. To gain further insight regarding how the solution of the TS rod model converges to the solution obtained with the proposed model as GA increases, we present Fig. 2. For higher values of GA, the deformed configuration approaches the deformed configuration computed with the new formulation. As we are already dealing with a very slender rod, the values of GA do not need to take high values to clearly demonstrate the convergent behavior. We can observe, moreover, that for small values of GA, the normal directors are not tangent to the deformed rod axis. These normal directors tend to become tangent to the deformed rod axis as GA increases. In our formulation, the directors are always unitary and tangent to the rod axis, even in the fully discrete setting. In addition, Fig. 3 shows the convergence behavior of the solution as a function of GA, which is consistent with the linear theory; see Sect. 5.1.3. To further reduce the relative error, it is necessary to employ finer meshes, but the trend is always the same and thus, this is omitted. For practical reasons nevertheless, it is also important to highlight that for a fixed mesh as GA increases, the condition number of the tangent stiffness matrix in the TS model increases, complicating the Newton–Raphson convergence.

Fig. 3

Cantilever rod subject to an in-plane loading: convergence behavior of the solution obtained with the TS rod model

Fig. 4

Cantilever rod subject to an out-of-plane loading: sequence of deformed configurations

5.1.2 Cantilever rod subject to a out-of-plane loading

This second example considers two loads given by $(-1.5, 0.5, 0)~\mathrm {N}$ and $(1.5, -0.5, 0.5)~\mathrm {N}$, applied to the node corresponding to the free end and to the central node (located at the mid span of the rod), respectively. Figure 4 shows a sequence of twelve successive deformed configurations. For sake of consistency, we follow the same color convention as in the previous example, and it will be retained in the remainder of this work. Once again, the results reveal that the convergence behavior in the three-dimensional case is as well consistent with the linear theory; see Sect. 5.1.3. As done in the previous example, we present Fig. 5 that shows how the deformed configuration obtained with the TS rod model approaches the deformed configuration computed with the proposed model for increasing values of GA. Figure 6 shows, moreover, the convergence behavior of the solution as a function of GA. These results confirm the observations made in the previous example.

Fig. 5

Cantilever rod subject to an out-of-plane loading: deformed configuration of several TS rods with increasing values of GA vs. the deformed configuration computed with the new formulation

Fig. 6

Cantilever rod subject to an out-of-plane loading: convergence behavior of the solution obtained with the TS rod model

5.1.3 Discussion on some key properties of the model proposed

The static setting represents a suitable framework to investigate some key properties of the model proposed and its spatial discretization with the finite element method. These include the convergence of the solution with the mesh refinement as well as the convergence of the solution of the TS rod to the solution obtained with new model. Figure 7 shows the relative energy norm for a cantilever rod subject to an out-of-plane tip moment, i.e., a pure bending case. The tip moment applied is equal to $EI \pi / 2 L$ and in the exact solution the rod bends to a quarter of circle. To investigate the convergence behavior, we consider meshes consisting of 5, 10, 20, 40, 80 and 160 finite elements. As expected, the convergence behavior exhibited is optimal and consistent with those results reported in [34] for such a standard benchmark case.

Figures 3 and 6 show a saturation in the relative difference between the solution obtained with the current model and the one resulting from the TS rod. To elucidate this behavior, we perform a comparison in the linear setting of the Euler–Bernoulli and the Timoshenko beam models referred to as E–B and T models, correspondingly. To this end, let us consider a straight beam with geometrical and material properties $L = 40~\mathrm {m}$ and $EI = 200~\mathrm {Nm^2}$ to which a tip load $p = 0.05~\mathrm {N}$ is applied. Again, the meshes considered consist of 5, 10, 20, 40, 80 and 160 finite elements. Moreover, for the T model, GA is varied in the range from $10^{-2}$ to $10^6~\mathrm {N}$. Figure 8 presents the results, and from it we can draw interesting observations. The first one, that is precisely what we wanted to investigate, is the saturation for moderate to large GA that is observed in the nonlinear setting. The second and expected one is that the saturation value reduces as the number of elements increases. The third one is that the curves corresponding to all discretizations are indistinguishably in the range of relatively small GA. In the log–log plot, these are proportional to $(\text {number of elements})^{-1}$ and lie practically on the same curve. The four and last observation is that the relative error shows an abruptly downwards jump for a critical value of GA, when approached from both sides. Such a critical value of GA is mesh dependent and is higher for finer discretizations. This behavior has not been observed for the nonlinear case. We have to stress, however, that the nonlinear case considers axial and bending deformations simultaneously. Thus, it could happen that the presence of the axial deformation mitigates such a phenomenon due to stress stiffening of the rod, or simply that high-order effects mask this behavior.

The first, third and fourth observations can be explained by considered the same structure modeled, this time, with a single finite element. For the E–B model, we have the following stiffness matrix, load vector and displacement solution.

$$\begin{aligned} K_{\mathrm {E-B}} = \frac{EI}{L^3} \begin{bmatrix} 12 &{} &{} -6L \\ &{} &{} \\ -6L &{} &{} 4 L^2 \end{bmatrix}, \quad f_{\mathrm {E-B}} = \begin{Bmatrix} p\\ \\ 0 \end{Bmatrix}, \quad u_{\mathrm {E-B}} = \frac{p\,L}{EI} \begin{Bmatrix} \frac{L^2}{3}\\ \\ \frac{L}{2} \end{Bmatrix}. \end{aligned}$$

(56)

Similarly, for the T model, we have the following stiffness, load vector and displacement solution:

$$\begin{aligned} K_{\mathrm {T}} = \begin{bmatrix} \frac{GA}{L} &{} &{} -\frac{GA}{2} \\ &{} &{} \\ -\frac{GA}{2} &{} &{} \left( \frac{GA\,L}{4}+\frac{EI}{L}\right) \end{bmatrix}, \quad f_{\mathrm {T}} = \begin{Bmatrix} p\\ \\ 0 \end{Bmatrix}, \quad u_{\mathrm {T}} = \frac{p\,L}{EI} \begin{Bmatrix} \frac{GA\,L^2+4 EI}{4 GA}\\ \\ \frac{L}{2} \end{Bmatrix}. \end{aligned}$$

(57)

Now, we can define the relative difference for the bending deflection in terms of the first components of both solutions as

$$\begin{aligned} \mathrm {reldiff}(GA;EI,L) = \frac{|u_{\mathrm {T}}-u_\mathrm {E-B}|}{|u_\mathrm {E-B}|}= \frac{|GA\, L ^2-12 EI|}{4\,GA \, L^2}. \end{aligned}$$

(58)

The relative difference saturates up from moderate values of GA and whose asymptotic value is $\frac{1}{4}$. For GA equal to $\frac{12EI}{L^2}$, the relative difference reaches the minimum, i.e., the zero value. A log–log plot of this function for $GA\in {\mathbb {R}}_{>0}$ shows qualitatively the same behavior observed in Fig. 8 for the multi-element cases. Therefore, we can conclude that the saturation behavior observed in the nonlinear case is fully consistent with the linear theory.

Fig. 7

Cantilever rod subject to an out-of-plane tip moment (pure bending): convergence behavior of the solution with respect to the number of elements

Fig. 8

Cantilever rod subject to in-plane loading (linear setting): convergence behavior of the solution obtained with the T model with respect to the E–B model

Fig. 9

Clamped-free rod subject to a vanishing in-plane loading: motion sequence

5.2 Dynamics

Here, we consider two initially straight, transversely isotropic rods of length $10~\mathrm {m}$, initially aligned with $\varvec{D}= (1, 0, 0)$. The chosen material properties are $E = 2.0\!\times \!10^{11}~\mathrm {Pa}$ and $\rho = 7900~\mathrm {kg/m}^3$. Again, we assume $\nu = 0.3$ for the TS rod. The rod is uniformly discretized into 20 two-node finite elements. The loads that will be specified later are scaled with a load intensity function that starting at $t = 0~\mathrm {s}$ increases linearly from zero to one, reaching the maximum at $t = 0.5 t_\mathrm {c}~\mathrm {s}$. After that instant, it decreases linearly from one to zero, reaching zero at $t = t_\mathrm {c}~\mathrm {s}$. For $t > t_\mathrm {c}~\mathrm {s}$, the load intensity function is identically zero. For both examples considered, we assume $t_\mathrm {c} = 0.5~\mathrm {s}$. To integrate the equation of motion, we employ a hybrid combination of the midpoint and trapezoidal rules. To integrate the equations of motion corresponding to the TS rod model, we employ the energy-momentum integration described in [14, 17].

5.2.1 Clamped-free rod subject to a vanishing in-plane loading

This third example investigates a rod with diameter $0.01~\mathrm {m}$. The node corresponding to the clamped end is the first one, and the one corresponding to the free end is the last one. A single load given by $(0, 30, 0)~\mathrm {N}$ is multiplied by the load intensity function and applied to node 21. Finally, the time step considered for the simulation is $\Delta t = 0.005~\mathrm {s}$. Figure 9 shows a motion sequence. We observe that the rod moves upward reaching a configuration with the seemly maximum deformation on the upper half plane and then returns from that first apparent extreme configuration. Then the rod moves downward reaching a configuration with the seemly maximum deformation on the lower half plane and then returns from that second apparent extreme configuration. Finally, the rod moves upward snapping from the lower to the upper half plane. At the beginning, the solutions obtained with both models are practically indistinguishable, but after some time some subtle differences become noticeable. Although the responses are not identical, they remain very close and the agreement is indeed very good. For a long-time simulation, the observed differences will become larger, and then, the responses will not match any longer.

Such a behavior is to be attributable not only to the different features of the formulations, but to the different features of the time integration schemes adopted as well. Figure 10 presents the time evolution of the potential, kinetic and total energies. In addition, we include the energies computed with the TS rod model. The latter ones are indicated in black with thinner lines, following the same styling to ensure an immediate identification. The integration method adopted for the current model is not able to identically preserve energy for the load-free system. Nevertheless, the total energy shows an acceptable behavior even provided that for $t>10~\mathrm {s}$, some subtle differences are observed for both the potential and kinetic energies. Such differences are consistent with the observations made on the corresponding motion sequence. Finally, the kinetic and potential energies show some moderate-frequency content that is not present in the reference solution. This can be presumably attributed to the simple integration method, which is well known for lacking robustness when compared to more reliable structure-preserving integration methods. However, a detailed explanation of such a behavior requires further investigations that are outside the scope of this work.

Fig. 10

Clamped-free rod subject to a vanishing in-plane loading—energy

5.2.2 Flying rod subject to vanishing out-of-plane loading

The last example considers a rod with diameter $0.005~\mathrm {m}$. Due to the absence of any kind of support, the rod is then free to move in space. A first pair of loads given by $(-30, -30, 0)~\mathrm {N}$ and $(30, 30, 0)~\mathrm {N}$ is applied to the nodes 1 and 21, whereas a second pair given by $(0, 0, -24)~\mathrm {N}$ and $(0, 0, 24)~\mathrm {N}$ is applied to the nodes 2 and 20, respectively. As before, the loads are multiplied by the load intensity function and the time step chosen for the simulation is $\Delta t = 0.001~\mathrm {s}$. Figure 11 shows a motion sequence. We observe that the rod deforms and rotates in space, but the central node remains always located at the position $(0, 0, 0)~\mathrm {m}$. Since the applied forces have a vanishing resultant the linear momentum is zero at all times and it is not plotted. The results of both models are in a good agreement. As in the previous case, for a long-term simulation, the observed differences will become larger and then, the responses will not match any longer. Nevertheless, this is not an unexpected behavior. Figures 12 and 13 present the time evolution of the three components of the angular momentum and the time evolution of the potential, kinetic and total energies, respectively. As before, the same quantities computed with the TS rod are presented in black with thinner lines, following the same styling. We observe, as expected, that the components of the angular momentum are exactly preserved after the loads vanish. The energy behavior is, as in the previous case, reasonably good and its value is almost constant. Although some small differences can be observed, the quality of the results is good and are, to a good extend, similar to the reference. This suggests that the proposed model is working correctly.

Fig. 11

Flying rod subject to a vanishing out-of-plane loading: sequence of motion

Fig. 12

Flying rod subject to a vanishing out-of-plane loading: angular momentum

Fig. 13

Flying rod subject to a vanishing out-of-plane loading: energy

6 Conclusion

This article presents the mathematical formulation and the numerical implementation of a new model for initially straight, transversely isotropic rods, that disregards the energetic contribution of torsion and exhibits no shear deformation.

Since the model disregards torsion and relies on a configuration space that precludes shearing, it is an unconstrained formulation. Let us recall that, in standard rod models, the elimination of the torsional and shear deformations can only be imposed by means of non-integrable restrictions that cannot be encompassed by a standard variational framework. In contrast, the model proposed admits a variational statement, as seen in its derivation from Hamilton’s principle, and is kept in its minimal irreducible representation. This is the first novelty of the present work.

We also present an strategy for the numerical discretization of the proposed model. It relies on finite elements for the spatial discretization, and a hybrid combination of the midpoint and trapezoidal rules for the time discretization. Through several examples, including static and dynamic cases, we show that the model has promising features due to its simplicity, low computational cost, and robustness. This is the second novelty of this work.

In a few words, the article combines theoretical modeling and numerical discretization for a model that, to the authors’ knowledge, has never been discussed in the literature. The robustness and low computational cost of the numerical solution presented makes it ideally suited for real-time applications of rods, for example, in multibody dynamics, polymer simulations, hair dynamics, etc. Specific applications such as these will be tackled in future works.