Nonlinear dynamic analysis of shear- and torsion-free rods using isogeometric discretization and outlier removal

Abstract

In this paper, we present a discrete formulation of nonlinear shear- and torsion-free rods introduced by Gebhardt and Romero (Acta Mechanica 232(10):3825–3847, 2021) that uses isogeometric discretization and robust time integration. Omitting the director as an independent variable field, we reduce the number of degrees of freedom and obtain discrete solutions in multiple copies of the Euclidean space $\left( \mathbb {R}^3\right)$, which is larger than the corresponding multiple copies of the manifold $\left( \mathbb {R}^3 \varvec{\times } S^2\right)$ obtained with standard Hermite finite elements. For implicit time integration, we choose the same integration scheme as Gebhardt and Romero in (2021) that is a hybrid form of the midpoint and the trapezoidal rules. In addition, we apply a recently introduced approach for outlier removal by Hiemstra et al. (Comput Methods Appl Mech Eng 387:114115, 2021) that reduces high-frequency content in the response without affecting the accuracy, ensuring robustness of our nonlinear discrete formulation. We illustrate the efficiency of our nonlinear discrete formulation for static and transient rods under different loading conditions, demonstrating good accuracy in space, time and the frequency domain. Our numerical example coincides with a relevant application case, the simulation of mooring lines.

1 Introduction

Nonlinear rods have a plethora of applications in science and engineering, for example, in the analysis of DNA molecules [1,2,3], the dynamics of cables [4, 5], the mechanical analysis of Möbius bands [6], or the stability of elastic knots [7, 8], among others. The shear-free model of rods is based on the assumption of cross-sections that remain flat and perpendicular to the tangent vector associated with the curve that describes the rod axis [9, 10]. In the context of linear rods, the Euler-Bernoulli and Rayleigh models are well-established [11, 12]. For nonlinear rods, one of the most widely used models is the so-called Kirchhoff rod, which can be considered a generalization of the Rayleigh model [13,14,15].

In general, it is not possible to formulate the governing equation of non-shearable rods through a truly unconstrained variational statement, particularly in dynamics problems, due to the non-integrable nature of vanishing shear deformations [9, 10]. In [16], Romero and Gebhardt developed an unconstrained variational formulation for this type of rod, but rely on certain simplification hypotheses. Recently, in [17], Gebhardt and Romero introduced a new unconstrained structural model for nonlinear initially straight rods that do not exhibit shear and torsion. This model provides a variational formulation for shear- and torsion-free rods and is a special case of the static and dynamic variational principles for nonlinear Kirchhoff rods developed in [16]. Moreover, it can be considered as the non-shearable counterpart of the torsion-free beam model introduced in [18]. The main advantage of the nonlinear rod formulation [17] is that it is an unconstrained variational statement which can be employed for both static and dynamic problems. It is in its simplest representation and does not include any non-integrable constraints such as the one enforcing non-twisting conditions [17]. This relies on the two approaches employed when deriving this rod formulation in kinematic and energetic settings, particularly to model the shear- and torsion-free behavior: (i) the construction of a configuration space that completely prevents shear, and (ii) discarding the torsion contribution to potential and kinetic energy.

We are aware of the cable formulation introduced in [19] which, given the same material law and initial straight geometry, essentially models similar slender beam-like structures as the rod formulation introduced in [17], which are represented by their middle curve and are shear- and torsion-free. Nevertheless, these two formulations are two different formulations due to the following two essential differences: (i) the rod formulation [17] does not depend on the definition or choice of the Frenet-Serret triad or any triad, while this is essential for the one introduced in [19], particularly, for its strain measure; and (ii) the rod formulation [17] considers the contribution of the rod director to the kinetic energy, which is not the case when using the cable formulation [19]. Moreover, to achieve the shear- and torsion-free behavior, these formulations employ different approaches: while the formulation [19] is based on the cable kinematics, which are equivalent to the Euler-Bernoulli assumptions, the rod formulation [17] models this behavior via constructing appropriate configuration spaces. Such configuration spaces also require the straight initial geometry of the rods, while the cable formulation in [19] is valid for arbitrary initial cable geometry.

Among the several approaches to solve the governing equations of non-shearable rods, we can mention classical nodal and isogeometric finite elements [19,20,21,22,23,24]. One of the key advantages of isogeometric finite elements is the higher-order smoothness of spline basis functions, which naturally fulfills the $C^1$ continuity required by the rod formulation. As a consequence, they have broad applications in the analysis of beam and shell structures, see e.g. [25,26,27,28,29,30,31,32]. For the recently developed nonlinear rod formulation in [17], the spatial discretization scheme applied so far is the one based on nodal finite elements. It relies on cubic Hermite functions to represent the discrete rod configuration, which in turn is decomposed into nodal positions and nodal directors. We refer to this scheme hereinafter as the standard discretization scheme. The discrete solutions obtained lie in the manifold $\left( \mathbb {R}^3 \varvec{\times } S^2\right) ^n$, where n is the number of discrete nodes. This discretization scheme establishes the first attempt to numerically solve the shear- and torsion-free Kirchhoff rod.

In this paper, we investigate an alternative spatial discretization scheme in the context of isogeometric analysis (IGA) for the rod formulation of [17]. In particular, we discretize the rod configuration in terms of the position of control points, without considering the director as an independent variable field. Hence, the number of degrees of freedom can be reduced and the discrete solution lies in multiple copies of the Euclidean space $\mathbb {R}^3$ which is a larger space than the corresponding multiple copies of the manifold $\left( \mathbb {R}^3 \varvec{\times } S^2\right)$ of the standard scheme. We utilize the higher-order smoothness of spline functions that naturally fulfill the $C^1$ continuity required by the rod formulation (and beyond). We illustrate, via static benchmarks of two- and three-dimensional cantilever rods, that isogeometric discretizations and the standard scheme achieve a comparable level of accuracy. Moreover, we show for a geometrically nonlinear cantilever rod bent to a circle that the convergence behavior is comparable to an optimal convergence in the $H^2$ semi-norm as linear fourth-order problems, but shows a smaller rate in the $H^1$ semi-norm and the $L^2$ norm. We also show via this example that decreasing the continuity of spline basis functions generally does not affect the accuracy and convergence, but reduces the convergence rate in $H^2$ semi-norm when using odd polynomial degrees.

For time integration in our dynamic computations, we employ the same implicit integration scheme as [17], which is a hybrid combination of the midpoint and trapezoidal rules. This type of implicit scheme has been shown to achieve second-order accuracy, approximately preserve the energy, and exactly preserves the linear and angular momentum [33,34,35]. We show, via dynamic benchmarks of two- and three-dimensional rods, that the isogeometric discretization scheme using B-splines with $C^1$ continuity or higher is less robust than the standard one. We improve its robustness via the strong approach of outlier removal introduced in [36]. We illustrate, via an example of an unconstrained rod subjected to out-of-plane vanishing forces, that the mass term associated to the inertia is irregular. Hence, the configuration-dependent mass matrix of the studied formulation behaves irregularly and, therefore, cannot be simplified to a constant matrix. Finally, we test our rod formulation for the nonlinear behavior of swinging rods under conservative, non-conservative, and pulsating forces. Our results indicate that our discrete isogeometric scheme is an efficient tool for such nonlinear computations.

The structure of the paper is as follows: In Sect. 2, we briefly review the nonlinear rod formulation. In addition, we derive the external forces induced by a surrounding flow, considered in our numerical examples. In Sect. 3, we discuss discretization in space with isogeometric finite elements, the resulting semi-discrete formulation, and differences when compared to the standard discretization scheme based on Hermite functions. We also briefly recap the implicit time integration scheme that is applied in our transient computations. In Sect. 4, we numerically demonstrate the robustness of isogeometric discretizations for two- and three-dimensional benchmarks and improve it via a strong approach for outlier removal. In Sect. 5, we apply the isogeometric nonlinear rod formulation to a swinging rubber rod subjected to different loading conditions, which can be considered as a relevant application case for the simulation of mooring lines. In Sect. 6, we summarize our results and draw conclusions.

2 Nonlinear shear- and torsion-free rods

In this section, we briefly review the formulation of nonlinear shear- and torsion-free rods in a continuous setting introduced in [17]. We then describe and derive the external forces induced by a surrounding flow that are considered in the numerical examples of this work. We start with a brief recap of required fundamental equations and definitions in differential geometry that are later utilized for the rod formulation.

2.1 Preliminaries

Consider an arbitrary regular one-parameter curve $\varvec{\varphi }\,=\, \varvec{\varphi }(s)$ in the ambient space $\mathbb {R}^3$, where $s \,\in \, [0,L]$ is the arc-length coordinate. Since $\varvec{\varphi }$ is regular, its first derivative with respect to s, denoted as $\varvec{\varphi }^\prime$, is non-zero, i.e. $\varvec{\varphi }^\prime \,\ne \, {\textbf {0}}$. The Frenet-Serret moving frame associated with the curve $\varvec{\varphi }$ is then:

$$\begin{aligned} {\textbf {d}} \, := \, \frac{\varvec{\varphi }^\prime }{\left| \varvec{\varphi }^\prime \right| } \, , \quad {\textbf {t}} \, := \, \frac{\varvec{\varphi }^{\prime \prime }}{\left| \varvec{\varphi }^{\prime \prime }\right| } \, , \quad {\textbf {b}} \, := \, \frac{\varvec{\varphi }^\prime \, \times \varvec{\varphi }^{\prime \prime }}{\left| \varvec{\varphi }^\prime \, \times \varvec{\varphi }^{\prime \prime }\right| } \, , \end{aligned}$$

(1)

where $(\cdot )^\prime$ denotes the first derivative with respect to the arc-length s, i.e. $(\cdot )^\prime = \partial (\cdot ) / \partial \, s$, and $\left| \cdot \right| \,:\, \mathbb {R}^3 \,\rightarrow \, \mathbb {R}_{\ge \,0}$ denotes the Euclidean vector norm. We refer to ${\textbf {d}}$ as the director of the curve $\varvec{\varphi }$. We note that ${\textbf {t}}$, and thus also ${\textbf {b}}$, is ill-defined at points where $\left| \varvec{\varphi }^{\prime \prime }\right| = 0$, while the director ${\textbf {d}}$ is well-defined everywhere along the curve $\varvec{\varphi }$.

The director ${\textbf {d}}$ lives in the unit sphere $S^2 := \{ \, {\textbf {d}} \,\in \, \mathbb {R}^3 \, \vert \, {\textbf {d}} \cdot \, {\textbf {d}} \,=\, 1 \, \}$ that is a nonlinear, smooth, compact, two-dimensional manifold [37, 38]. The tangent bundle associated with $S^2$ is also a manifold, which is given by $T S^2 :=$ $\{ \, ({\textbf {d}},{\textbf {c}}) \,\in \, S^2 \,\times \,\mathbb {R}^3 \,, {\textbf {d}} \cdot \, {\textbf {c}} \,=\, 0 \}$. We recall that the covariant derivative of a smooth vector field ${\textbf {v}}\,:\,S^2 \,\rightarrow \, T S^2$ along a vector field ${\textbf {w}}\,:\,S^2 \,\rightarrow \, T S^2$ is a vector field in $T S^2$ evaluated at ${\textbf {d}}$, given by:

$$\begin{aligned} \nabla _{{\textbf {w}}} \, {\textbf {v}} \,:=\, \left( \, {\textbf {I}} \,-\, {\textbf {d}} \,\otimes \, {\textbf {d}} \,\right) \, D\,{\textbf {v}} \,\cdot \, {\textbf {w}}, \end{aligned}$$

(2)

where ${\textbf {I}}$ denotes the identity matrix, and $D\,{\textbf {v}}$ the derivative of ${\textbf {v}}$. The covariant derivative $\nabla _{{\textbf {w}}} \, {\textbf {v}}$ is the projection of $D\,{\textbf {v}}$ in the direction of ${\textbf {w}}$ onto the tangent plane at ${\textbf {d}}$ [37, 38].

For the rod formulation [17] considered in this work, we are particularly interested in the covariant derivative of $\varvec{\varphi }^\prime$ in the direction of ${\textbf {d}}^\prime$. Applying (2), this covariant derivative takes the following form [17]:

$$\begin{aligned} \nabla _{{\textbf {d}}^\prime } \, \varvec{\varphi }^\prime \,=\, \underbrace{\left( \,{\textbf {I}} - {\textbf {d}} \, \otimes \, {\textbf {d}}\, \right) }_{\mathcal {P}_{{\textbf {d}}}} \, \varvec{\varphi }^{\prime \prime } \,, \end{aligned}$$

(3)

where ${\textbf {d}}^\prime$ is computed by taking the derivative of (1), i.e., ${\textbf {d}}^\prime \,=\, \frac{1}{\left| \varvec{\varphi }^\prime \right| } \, \mathcal {P}_{{\textbf {d}}}\, \varvec{\varphi }^{\prime \prime }$. We refer to $\mathcal {P}_{{\textbf {d}}}$ as the orthogonal projection operator.

2.2 Strong and weak forms

Let the curve $\varvec{\varphi }$ now be the configuration of Kirchhoff rods, dependent on the arc-length s and time t, $\varvec{\varphi }\,=\, \varvec{\varphi }(s,\,t)$, $(s,\,t) \,\in \, [0,\,L] \, \times \, [0,\,T]$, that are initially straight, shear-, torsion-free, and transversely isotropic [17]. Next, let us Consider the following set for the rod configurations:

$$\begin{aligned} \begin{aligned}&\mathcal {D}:= \left\{ \varvec{\varphi }\, \in \, \left[ C^2 (0, \,L)\right] ^3, \; \left| \varvec{\varphi }^\prime \right| \, > \, 0, \; \right. \\&\qquad \qquad \left. \varvec{\varphi }(0,\,t) = \textbf{0}, \; \varvec{\varphi }^\prime \, (0,\,t) \, = \, {\textbf {E}}_3 \right\} \,, \end{aligned} \end{aligned}$$

(4)

where $C^2(0, \,L)$ is the space of $C^2$ continuous functions on $(0, \,L)$, ${\textbf {E}}_i$, $i=1,2,3$, are the canonical Cartesian basis of $\mathbb {R}^3$. For simplicity, we adopt here the clamped boundary condition at $s=0$.

We recall, from [17], the strong form of the equations of motion governing the space-time evolution for the Kirchhoff rod:

$$\begin{aligned} \begin{aligned}&{\textbf {n}}^\prime + \left( \, \frac{1}{\left| \varvec{\varphi }^\prime \right| } \, {\textbf {d}} \,\times \, \nabla _{{\textbf {d}}^\prime } \, {\textbf {m}} \, \right) ^\prime \, = \, \\&\qquad \qquad A_\rho \,\ddot{\varvec{\varphi }} \,+\, \left( \, \frac{1}{\left| \varvec{\varphi }^\prime \right| } \, {\textbf {d}} \,\times \, I_\rho \, \nabla _{\dot{{\textbf {d}}}} \, \dot{{\textbf {d}}} \right) ^\prime \,-\, {\textbf {f}}^{\text {ext}} , \end{aligned} \end{aligned}$$

(5)

where ${\textbf {n}}$ and ${\textbf {m}}$ are the stress measures, defined as:

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

(6)

respectively, which are conjugated with the following strain measures:

$$\begin{aligned} \varvec{\epsilon } \,:=\, \varvec{\varphi }^\prime \,-\, {\textbf {d}} \,, \qquad \varvec{\kappa } \,:=\, {\textbf {d}} \,\times \,{\textbf {d}}^\prime \,. \end{aligned}$$

(7)

Here, $A_\rho$ and $I_\rho$ are the mass per unit length and the inertia density, respectively, i.e. $A_\rho \,=\, \rho \,A$ and $I_\rho \,=\, \rho \,I$, where $\rho$ is the mass density, A the cross-section area and I the moment of inertia of the rod. ${\textbf {f}}^{\text {ext}}$ is the external generalized forces, and the dot notation in the superscript denotes the derivative with respect to time t, i.e. $\dot{(\cdot )} \,=\, \partial (\cdot ) / \partial \, t$. We note that since the director ${\textbf {d}}$ is well-defined along the rod $\varvec{\varphi }\,\in \, \mathcal {D}$ (see also (4)), as discussed in the previous subsection, the strain measures (7) are also well-defined at every point of the rod.

At time $t=0$, we require the following initial conditions:

$$\begin{aligned}&\varvec{\varphi }\,=\, \varvec{\varphi }_0 \quad&\text {on } (s,\,t) \,\in \, [0,\,L] \, \times \, [0] \,, \end{aligned}$$

(8a)

$$\begin{aligned}&\dot{\varvec{\varphi }} \,=\, {\textbf {v}}_0 \quad&\text {on } (s,\,t) \,\in \, [0,\,L] \, \times \, [0] \,. \end{aligned}$$

(8b)

Additionally, we require at all times the following boundary conditions; for instance, clamped-free ends:

$$\begin{aligned}&\text {on } (s,\,t) \,\in \, [0] \, \times \, [0,\,T]: \,, \nonumber \\&\quad \varvec{\varphi }\,=\, {\textbf {0}}\,, \qquad \varvec{\varphi }^\prime \,=\, {\textbf {E}}_3 \,, \end{aligned}$$

(9a)

$$\begin{aligned}&\text {on } (s,\,t) \,\in \, [L] \, \times \, [0,\,T]: \,, \nonumber \\&\quad {\textbf {n}} \,+\, \frac{1}{\left| \varvec{\varphi }^\prime \right| } \, {\textbf {d}} \,\times \, \left( \, \nabla _{{\textbf {d}}^\prime } \, {\textbf {m}} \,-\, I_\rho \, \nabla _{\dot{{\textbf {d}}}} \, \dot{{\textbf {d}}} \,\right) \,=\, {\textbf {0}} \,, \end{aligned}$$

(9b)

$$\begin{aligned}&\frac{1}{\left| \varvec{\varphi }^\prime \right| } \, {\textbf {d}} \,\times \, {\textbf {m}} \,=\, {\textbf {0}} \,. \end{aligned}$$

(9c)

According to [17], the weak form corresponding to (5) is then:

$$\begin{aligned} \begin{aligned} \int _0^S \,&\delta \varvec{\varphi }\,\cdot \, \left( \mathcal {M}\left( \varvec{\varphi }^\prime \right) \,\hat{\nabla }_{\dot{\varvec{\varphi }}} \, \dot{\varvec{\varphi }} \,+\, \right. \\&\qquad \left. \mathcal {B}\left( \varvec{\varphi }^\prime ,\,\varvec{\varphi }^{\prime \prime }\right) ^T \, \varvec{\sigma } \,-\, {\textbf {f}}^{\text {ext}} \,\right) \, \textrm{d} \, s \,=\, 0, \end{aligned} \end{aligned}$$

(10)

where the mass operator, $\mathcal {M}$, and the linearized strain operator, $\mathcal {B}$, are given by:

$$\begin{aligned}&\mathcal {M} \,=\, \mathcal {M}\left( \varvec{\varphi }^\prime \right) \, :=\, A_\rho \, {\textbf {I}} \,+\, \nonumber \\&\qquad (\cdot )^{\prime \,T} \, I_\rho \, \frac{1}{\left| \varvec{\varphi }^\prime \right| ^2} \, \mathcal {P}_{{\textbf {d}}}\, (\cdot )^\prime \, \end{aligned}$$

(11)

$$\begin{aligned}&\mathcal {B} \,=\, \mathcal {B}\left( \varvec{\varphi }^\prime ,\,\varvec{\varphi }^{\prime \prime }\right) \,:=\, \nonumber \\&\quad \begin{bmatrix} {\textbf {I}} \,-\, \frac{1}{\left| \varvec{\varphi }^\prime \right| } \, \mathcal {P}_{{\textbf {d}}}&{} {\textbf {0}} \\ -\frac{1}{\left| \varvec{\varphi }^\prime \right| ^2} \, \left[ \varvec{\varphi }^{\prime \prime }\right] _\times \, \mathcal {H}_{{\textbf {d}}}&{} \frac{1}{\left| \varvec{\varphi }^\prime \right| } \, \left[ {\textbf {d}}\right] _\times \end{bmatrix} \; \begin{bmatrix} (\cdot )^\prime \\ (\cdot )^{\prime \prime } \end{bmatrix} \,. \end{aligned}$$

(12)

Here, $\varvec{\sigma } \,:=\, [{\textbf {n}} \quad {\textbf {m}}]^T$, $\mathcal {H}_{{\textbf {d}}}$ is the Householder operator,¹$\mathcal {H}_{{\textbf {d}}}:=\, {\textbf {I}} \,-\, 2 \, {\textbf {d}} \,\otimes \, {\textbf {d}}$, and $[{\textbf {a}}]_\times$ denotes the skew-symmetric matrix of a vector ${\textbf {a}} \,=\, \left[ a_1 \quad a_2 \quad a_3\right] ^T$, i.e.:

$$\begin{aligned}{}[{\textbf {a}}]_\times \,=\, \begin{bmatrix} 0 &{} -a_3 &{} a_2 \\ a_3 &{} 0 &{} -a_1 \\ -a_2 &{} a_1 &{} 0 \end{bmatrix} \,. \end{aligned}$$

The field covariant derivative $\hat{\nabla }_{(\cdot )} \, ( \cdot )$ is the extension of the covariant derivative (2).

3 Isogeometric discrete rod model and implicit time integration

In this section, we first discuss an alternative spatial discretization scheme of the rod formulation (10), reviewed in the previous section. We employ isogeometric discretizations, which utilize the higher-order continuity of smooth spline functions fulfilling the $C^1$-continuity required by the considered rod formulation. We show that this alternative discretization scheme yields different semi-discrete formulations, reduces the number of degrees of freedom, and leads to a larger solution space than the standard one based on nodal finite elements using cubic Hermite functions. We then numerically demonstrate, for a geometrically nonlinear rod bent to a circle, that the isogeometric discretizations show optimal convergence in $H^2$ semi-norm and a smaller convergence rate in $H^1$ semi-norm and $L^2$ norm. In addition, we briefly review the implicit time integration scheme employed in this work. We close this section with a discussion of applying the strong approach of outlier removal [36] to improve the robustness of isogeometric discretizations.

3.1 Spatial discretizations

The rod formulation (10) requires discretizations of at least $C^1$-continuity. To fulfill this, the standard discretization scheme based on nodal finite element employs cubic Hermite functions and discretizes both the nodal spatial position and nodal director as two variable fields [17]. In this work, we want to utilize the higher-order continuity of smooth spline functions that allow lower polynomial degree and one can omit the director as a variable field. Thus, we discretize the rod configuration, $\varvec{\varphi }(s,\,t) \,\in \, \mathcal {D}$, and its variation, $\delta \varvec{\varphi }(s,\,t)$, by a weighted finite sum of $m_\text {dof}$ B-splines, $N_i \,(s)$, with continuity $C^r$ and polynomial degree p [39, 40], where r is the continuity order, $1 \,\le \, r \,\le \, p-1$, as follows:

$$\begin{aligned}&\varvec{\varphi }(s,\,t) \,\approx \, \varvec{\varphi }_h \,(s,\,t) \nonumber \\&\qquad = \sum _i^{m_\text {dof}} \, N_i \, (s) \, {\textbf {x}}_i\,(t) \,=\, {\textbf {N}} \, {\textbf {q}} \,, \end{aligned}$$

(13a)

$$\begin{aligned}&\delta \varvec{\varphi }(s) \,\approx \, \delta \varvec{\varphi }_h (s) \nonumber \\&\qquad = \sum _i^{m_\text {dof}} \, N_i \, (s) \, \delta {\textbf {x}}_i \,=\, {\textbf {N}} \, \delta {\textbf {q}} \,. \end{aligned}$$

(13b)

Here, $\varvec{\varphi }_h\,=\,\varvec{\varphi }_h\,(s,\,t) \in \mathbb {R}^3$ denotes the discrete rod configuration in space, ${\textbf {x}}_i \, \in \mathbb {R}^3$ is the time-dependent position of the $i^{\text {th}}$ control point, ${\textbf {q}}\,=\,{\textbf {q}}(t) \in (\mathbb {R}^3)^{m_\text {dof}}$ is the vector of unknown time-dependent coefficients, and $\delta \varvec{\varphi }_h(s)$, $\delta {\textbf {x}}_i$, and $\delta {\textbf {q}}$ their variations, respectively. The discrete director ${\textbf {d}}$ and strain/stress measures follow directly from their definitions in (1) and (6)–(7), respectively.

Introducing (13) into the variational formulation (10), we obtain the following semi-discrete formulation:

$$\begin{aligned} \begin{aligned}&\text {Find } {\textbf {q}}(t) \,\in \, \mathbb {R}^{3\,m_\text {dof}}, t\,\in \,[0,\,T], \text { such that}: \\&\int _0^S \, \delta {\textbf {q}} \, \cdot \, \left( \, {\textbf {M}} ({\textbf {q}}) \, \nabla _{\dot{{\textbf {q}}}} \, \dot{{\textbf {q}}} \,+\, {\textbf {B}} ({\textbf {q}})^T \, \varvec{\sigma }_h \right. \\&\qquad \quad \qquad \qquad \left. -\, {\textbf {N}}^T \, {\textbf {f}}^{\text {ext}} \,\right) \, \textrm{d} \, s \,=\, {\textbf {0}} \quad \forall \, \delta {\textbf {q}} \,\in \, \mathbb {R}^{3\,m_\text {dof}} . \end{aligned} \end{aligned}$$

(14)

Here, the mass matrix ${\textbf {M}}$ and the matrix ${\textbf {B}}$, resulting from the operators (11) and (12), respectively, are:

$$\begin{aligned}&{\textbf {M}} \,=\, {\textbf {M}} ({\textbf {q}}) \,=\, \underbrace{A_\rho \, {\textbf {N}}^T \, {\textbf {I}} \, {\textbf {N}}}_{{\textbf {M}}_1} \,+\nonumber \\&\qquad \qquad \underbrace{I_\rho \, \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^2} \, \left( {\textbf {N}}^\prime \right) ^T \, {\textbf {P}}_{{\textbf {d}}}\, {\textbf {N}}^\prime }_{{\textbf {M}}_2} \\&{\textbf {B}} \,=\, {\textbf {B}}({\textbf {q}}) \,:=\, \nonumber \\&\qquad \qquad \begin{bmatrix} {\textbf {I}} \,-\, \frac{1}{\left| \varvec{\varphi }_h^\prime \right| } \, {\textbf {P}}_{{\textbf {d}}}&{} {\textbf {0}} \\ -\frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^2} \, \left[ \varvec{\varphi }_h^{\prime \prime }\right] _\times \, {\textbf {H}}_{{\textbf {d}}}&{} \frac{1}{\left| \varvec{\varphi }_h^\prime \right| } \, \left[ {\textbf {d}}_h\right] _\times \end{bmatrix} \; \begin{bmatrix} {\textbf {N}}^\prime \\ {\textbf {N}}^{\prime \prime } \end{bmatrix} \, , \end{aligned}$$

(15)

with ${\textbf {P}}_{{\textbf {d}}}\,=\, {\textbf {I}} \,-\, {\textbf {d}}_h \,\otimes \, {\textbf {d}}_h$, and ${\textbf {H}}_{{\textbf {d}}}\,=\, {\textbf {I}} \,-\, 2\,{\textbf {d}}_h \,\otimes \, {\textbf {d}}_h$. The discrete stress measures $\varvec{\sigma }_h$ are:

$$\begin{aligned} \varvec{\sigma }_h \,=\, \begin{bmatrix} {\textbf {n}}_h \\ {\textbf {m}}_h \end{bmatrix} \,&=\, \begin{bmatrix} EA \, \varvec{\epsilon }_h \\ EI \, \varvec{\kappa }_h \end{bmatrix} \nonumber \\&= \begin{bmatrix} EA \, \left( \, \varvec{\varphi }_h^\prime \,-\, {\textbf {d}}_h \,\right) \\ EI \, {\textbf {d}}_h \,\times \, {\textbf {d}}_h^\prime \end{bmatrix} \,. \end{aligned}$$

(16)

The term ${\textbf {M}} ({\textbf {q}}) \, \nabla _{\dot{{\textbf {q}}}} \, \dot{{\textbf {q}}}$ in (14), derived in [17], takes the following form:

$$\begin{aligned}&{\textbf {M}} ({\textbf {q}}) \, \nabla _{\dot{{\textbf {q}}}} \, \dot{{\textbf {q}}} \,=\, {\textbf {M}} \, \ddot{{\textbf {q}}} \,-\, \nonumber \\&\quad 2\,I_\rho \, \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^2} \, \left( {\textbf {N}}^\prime \right) ^T \, \left[ \, \frac{1}{\left| \varvec{\varphi }_h^\prime \right| } \, \left( {\textbf {d}}_h \, \cdot \, \dot{\varvec{\varphi }}_h^\prime \right) {\textbf {P}}_{{\textbf {d}}}\,+\, {\textbf {d}}_h \,\odot \, {\textbf {d}}_h \,\right] \, \dot{{\textbf {q}}} \nonumber \\&\quad \,+\, 2\,I_\rho \, \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^3} \, \left( {\textbf {N}}^\prime \right) ^T \, \left[ \, {\textbf {P}}_{{\textbf {d}}}\,\odot \, \left( \dot{\varvec{\varphi }}_h^\prime \,\otimes \, {\textbf {d}}_h \right) \,\right] \, \dot{\varvec{\varphi }}_h^\prime \,, \end{aligned}$$

(17)

where $\odot$ denotes the symmetric product between two vectors ${\textbf {a}}_1$, ${\textbf {a}}_2$, or two second-order tensors ${\textbf {A}}_1$, ${\textbf {A}}_2$, that is:

$$\begin{aligned}&{\textbf {a}}_1 \,\odot \, {\textbf {a}}_2 \,=\, \frac{1}{2} \, \left( \,{\textbf {a}}_1 \,\otimes \, {\textbf {a}}_2 \,+\, {\textbf {a}}_2 \,\otimes \, {\textbf {a}}_1 \, \right) \,, \end{aligned}$$

(18)

$$\begin{aligned}&{\textbf {A}}_1 \,\odot \, {\textbf {A}}_2 \,=\, \frac{1}{2} \, \left( \,{\textbf {A}}_1 \, {\textbf {A}}_2 \,+\, {\textbf {A}}_2^T \, {\textbf {A}}_1^T \, \right) \,. \end{aligned}$$

(19a)

Remark 1

The semi-discrete formulation (14) includes neither rotational degrees of freedom nor nodal directors. For cases involving prescribed directions such as clamped boundaries, one can employ for example the bending strip method [41] or the Nitsche’s method [42]. In this work, we strongly enforce the prescribed directions as built-in constraints in the trial and test spaces via the extraction operator (see Sect. 3.4 and Algorithm 1).

3.2 Isogeometric versus classical nodal finite elements

The isogeometric discretization scheme and the standard one based on nodal finite element using cubic Hermite functions both are based on the isoparametric concept. When the former applies cubic $C^1$ B-splines, the basis functions of these two schemes span the same function space. However, they belong to two different classes of finite element methods. While the standard scheme employs classical nodal finite elements, the isogeometric scheme is in the context of isogeometric analysis, where we deal with control points instead of element nodes [43, 44]. Moreover, the former requires two variable fields that are the nodal spatial position and the nodal director [17], while the latter only considers the positions of the control points as a variable field. Hence, the isogeometric scheme can reduce the number of degrees of freedom (dofs). In particular, a discretization with $n_e$ elements using $C^r$ B-splines of degree p leads to $3\left[ n_e \, (p-r) + r + 1\right]$ dofs². Applying B-splines with maximum continuity $C^{p-1}$ leads to $3(n_e + p)$ dofs, while the standard scheme leads to $5(n_e+1)$ dofs. For example, employing B-splines of minimum required polynomial degree of $p=2$, that are $C^1$ continuous, leads to a smaller number of dofs of $3(n_e+2)$ for the same number of elements. Using cubic $C^1$ B-splines that are in the same function space as cubic Hermite functions, however, leads to more dofs of $6(n_e+1)$ for the same number of elements. Due to the different variable fields of these two schemes, they lead to different solution spaces. Particularly, using the standard scheme results in a discrete solution in the manifold $\varvec{\varphi }_h \in \left( \mathbb {R}^3 \varvec{\times } S^2\right) ^{m_\text {dof}}$ since the director belongs to a unit sphere $S^2$. This necessary means that the standard scheme also preserves the manifold structure of the continuous rod. Isogeometric discretizations, however, lead to a solution in the Euclidean space $\varvec{\varphi }_h \in \left( \mathbb {R}^3\right) ^{n_\text {dof}}$, which is a larger space but does not ensure the underlying manifold structure of the rod.

Remark 2

We note that in general, spline basis functions have supports up to $p+1$ elements [43, 44], where p is their polynomial degree, which is larger than the support of basis functions in classical finite elements. However, this does not lead to a larger bandwidth of matrices obtained from the two discretization schemes [44, p. 92-97]. Hence, given the same number of degrees of freedom, these two schemes require similar computational cost. For a detailed comparison between the two discretization schemes, we refer to [44, p. 92-97].

Fig. 1

Convergence of the relative error between the geometrically exact beam and the nonlinear rod model [17], computed with different discretizations on a mesh of 40 elements, obtained at the last load step. (Color figure online)

We now numerically demonstrate, for two- and three-dimensional static benchmarks from [17, Sec. 5.1], that the isogeometric scheme, with a possibly smaller number of dofs, and the standard one approximately achieve the same accuracy. We consider an initially straight, transversely isotropic, clamped rod of 40 m with an axial stiffness of $EA \,=\, 100$ N, bending stiffness of $EI \,=\, 200$ Nm$^2$, and subjected to an in-plane and out-of-plane loading. We compare the response obtained with isogeometric discretizations using quadratic and cubic $C^1$ B-splines, with classical nodal finite elements using cubic $C^1$ Hermite functions, and with a geometrically exact beam model including shear and torsion deformations using linear $C^0$ Lagrange polynomials. The results of the last two approaches are provided by [17]. For the first two approaches, we discretize the rod with the same number of 40 uniform elements with the force divided into 55 uniform load steps. The number of dofs of the isogeometric scheme using quadratic, cubic B-splines, and the standard scheme is 126, 246, and 205, respectively. We have chosen a tolerance of $10^{-10}$ for the Newton–Raphson scheme, which required up to 6 iterations in all cases. We obtained visually indistinguishable snapshots when compared to the standard scheme (see [17, Fig. 1,4]) in all cases. In Fig. 1, we compare the difference between the geometrically exact beam formulation and the isogeometric beam formulation with the one between the former and the standard Hermite scheme. To this end, we plot the relative $L^2$ error between the deformed rod obtained with the geometrically exact beam and the one obtained with the isogeometric (black and green curves), and the standard scheme (blue), as a function of the shear stiffness GA employed in the geometrically exact beam model. In [17], the authors employed this convergence study to point out that the geometrically exact beam model converges to the model using their rod formulation when GA increases, since the normal directors tend to become tangent to the deformed rod axis with increasing GA. We adopt this study to illustrate that the isogeometric scheme results in the same behavior as the standard one, as observed in Fig. 1 for both the in-plane and out-of-plane loading cases, irrespective of the spline basis employed in the isogeometric scheme. We note that the constant error level between the geometrically exact beam and the nonlinear rod of [17] at large values of GA is due to shear and torsion deformations considered in the former but not in the latter.

Fig. 2

Deformed configurations of a clamped rod bent to a circle at different load steps, computed with quadratic $C^1$ B-splines ($p=2$) and a mesh of 40 elements

Fig. 3

Convergence of relative errors of the clamped rod bent to a circle computed with $C^{p-1}$ B-splines of different degrees p (left column) and with $C^r$ B-splines, $1\,\le \,r\,\le \,p-1$ (right column). The reference rate shown in figures on the right column is the convergence rate of linear fourth-order problems [45]. (Color figure online)

To study the convergence behavior of the isogeometric scheme with mesh refinement, we consider a geometrically nonlinear benchmark of planar roll-up. Figure 2a illustrates the initial rod subjected to a bending moment $M_{\text {max}}\,=\,\frac{2\,EI\,\pi }{L}$ Nm at its free end, where L is the rod initial length. We choose the same material parameters and value of L as in the static benchmarks above. We illustrate the deformed rod in a sequence of six load steps, obtained with the isogeometric scheme using quadratic $C^1$ B-splines ($p=2$) in Fig. 2. We observe that, as expected, the deformed rod closes a circle in the last load step. In the left column of Fig. 3, we plot the convergence of the relative error in $L^2$ norm, $H^1$, and $H^2$ semi-norm, obtained with $C^{p-1}$ B-splines of different polynomial degrees p, compared to the exact circle in the last load step. As a reference, we included the convergence rate of linear fourth-order problems using isogeometric discretizations based on [45]. We observe the same optimal convergence behavior obtained in $H^2$ semi-norm as the linear case for all degrees (see Fig. 3e). The error in the $H^1$ semi-norm (see Fig. 3c), however, converges with the same rate of linear cases only when using even degrees, but with one order lower when using odd degrees. Focusing on the error in the $L^2$ norm, we see that the convergence rate is smaller for all $p\,\ge \,3$ (see Fig. 3a). Furthermore, cubic $C^2$ B-splines illustrate the same convergence rate as quadratic, and quintic functions the same as quartic ones. We note that this is not the first time that the different convergence behavior between even and odd degrees has been observed. Such behavior is well-known in isogeometric collocation methods [46, 47], however, to our best knowledge, is not yet established for nonlinear problems.

Remark 3

For fourth-order problems, one can mathematically prove the optimal convergence merely in the $H^2$ norm, not in $H^1$ or $L^2$ norms (see e.g. [45, Eq. 3.13]). Convergence can be proven in $H^1$ and $L^2$ norms, generally using the Lax-Milgram lemma and the Poincaré inequality, where the latter leads to additional scaling factors associated with the mesh size, i.e. the proof is not the one of the optimal convergence according to the Lax-Milgram lemma. This necessary means that one may not expect optimal convergence in $H^1$ and $L^2$ norms for fourth-order problems.

In the second part of the convergence study, we investigate the impact of the smoothness of spline basis functions on the convergence of isogeometric discretizations. We compute the relative error obtained with $C^r$ B-splines of different polynomial degrees $p\,\ge \,3$, where $1\,\le \,r\,\le \,p-1$, and illustrate this on the right column of Fig. 3. Different markers correspond to different continuous orders r. We can see that the continuity of B-splines does not affect the convergence rate in both the $L^2$ norm and the $H^1$ semi-norm for this planar roll-up example. Focusing on the error in the $H^2$ semi-norm, reducing the continuity of B-spline basis functions decreases the convergence rate in cases of odd degrees (see the green and purple curves in Fig. 3f), however, does not affect this rate in the case of even degree ($p=4$, red curve in Fig. 3f). We note that the impact of reducing or increasing continuity of spline basis functions is not the same when studying different benchmarks and cases, see e.g. [48, 49]. Furthermore, based on empirical results, we observe that increasing or decreasing the number of quadrature points also does not affect the convergence behavior of the discretizations employed in this example. Further numerical investigations of a linear simply supported beam subjected to sinusoidal loads, using the nonlinear rod formulation [17], shows the same convergence behavior as in the linear cases for all error norms, polynomial degrees, and continuity. Hence, we suggest that the convergence behavior in the planar roll-up example results from the nonlinear behavior captured by the employed rod formulation. A mathematical error estimate for the considered rod formulation [17] is outside the scope of this work and is postponed to future work.

3.3 Time integration scheme

For time integration of our numerical examples in subsequent sections, we apply the same implicit scheme as in [17], which is a hybrid combination of the midpoint and trapezoidal rules. This implicit scheme achieves second-order accuracy, approximately preserves energy, and exactly preserves the linear and angular momenta [17, 33], which we also verify via an example of an elastic pendulum in Appendix A. We note that the employed integration scheme is based on the one introduced in [33]. In this work, to observe the occurrence and investigate the effects of all contents of the response, including the spurious high-frequency contents (see also discussions in the next section), when using isogeometric discretizations, we choose to eliminate the dissipation terms of the original scheme. We note that this choice is also to serve the comparison purpose with the results obtained with the standard nodal finite elements employed in [17]. Moreover, in [33], the authors have discussed and shown that the original scheme [33] is one possible realization of the well-established Energy-Dissipative-Momentum-Conserving method [50,51,52,53].

Remark 4

The mass matrix (15) of the considered rod formulation [17] is configuration-dependent due to its second counterpart ${\textbf {M}}_2$. Thus, an explicit time integration scheme would not be applicable.

Consider the semi-discrete formulation (14) in space evaluated at time instant $t_{n\,+\,\frac{1}{2}}\in [t_n,t_{n+1}]$:

$$\begin{aligned} \begin{aligned}&\int _0^S \, \delta {\textbf {q}} \, \cdot \, \left( \, {\textbf {M}} ({\textbf {q}}) \, \nabla _{\dot{{\textbf {q}}}} \, \dot{{\textbf {q}}} \,+\, {\textbf {B}} ({\textbf {q}})^T \, \varvec{\sigma }_h \right. \\&\qquad \qquad \qquad \left. -\, {\textbf {N}}^T \, {\textbf {f}}^{\text {ext}} \,\right) _{n+\frac{1}{2}} \, \textrm{d} \, s \,=\, {\textbf {0}} . \end{aligned} \end{aligned}$$

(19b)

We approximate the inertial term $\left( \,{\textbf {M}} ({\textbf {q}}) \, \nabla _{\dot{{\textbf {q}}}} \, \dot{{\textbf {q}}}\,\right) _{n+\frac{1}{2}}$ using an extended version of the midpoint rule as follows:

$$\begin{aligned}&\left( \,{\textbf {M}} ({\textbf {q}}) \, \nabla _{\dot{{\textbf {q}}}} \, \dot{{\textbf {q}}}\,\right) _{n+\frac{1}{2}} \,\approx \, \nonumber \\&\qquad \frac{{\textbf {M}}\left( {\textbf {q}}_{n+1}\right) \, \dot{{\textbf {q}}}_{n+1} \,-\, {\textbf {M}}\left( {\textbf {q}}_{n}\right) \, \dot{{\textbf {q}}}_{n}}{\Delta \, t} \; +\, \nonumber \\&\quad \left\{ \, 2\,I_\rho \, \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^3} \, \left( {\textbf {N}}^\prime \right) ^T \, \left[ \, {\textbf {P}}_{{\textbf {d}}}\,\odot \, \left( \dot{\varvec{\varphi }}_h^\prime \,\otimes \, {\textbf {d}}_h \right) \,\right] \, \dot{\varvec{\varphi }}_h^\prime \,\right\} _{n+\frac{1}{2}} \,, \end{aligned}$$

(20)

where $\Delta t$ is the time step. The internal term

$\left( \,{\textbf {B}} ({\textbf {q}})^T \, \varvec{\sigma }_h\,\right) _{n+\frac{1}{2}}$ is approximated by using the trapezoidal rule as follows:

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

(21)

We also approximate ${\textbf {q}}_{n+\frac{1}{2}}$ and $\dot{{\textbf {q}}}_{n+\frac{1}{2}}$ using the trapezoidal and midpoint rules, respectively, as follows:

$$\begin{aligned} {\textbf {q}}_{n+\frac{1}{2}} \,\approx \, \frac{{\textbf {q}}_{n+1} \,+\, {\textbf {q}}_n}{2} \,, \quad \dot{{\textbf {q}}}_{n+\frac{1}{2}} \,\approx \, \frac{{\textbf {q}}_{n+1} \,-\, {\textbf {q}}_n}{\Delta \, t} \,, \end{aligned}$$

(22)

and $\dot{{\textbf {q}}}_{n+1}$ as:

$$\begin{aligned} \dot{{\textbf {q}}}_{n+1} \,\approx \, \frac{2}{\Delta \, t} \, \left( \, {\textbf {q}}_{n+1} \,-\, {\textbf {q}}_n \,\right) \,-\, \dot{{\textbf {q}}}_n \,. \end{aligned}$$

(23)

Introducing the approximations (21)–(24) into (20) leads to a system of discrete nonlinear equations in space and time:

$$\begin{aligned} {\textbf {g}} \, ({\textbf {q}}_{n+1}) \,=\, {\textbf {0}} \, , \end{aligned}$$

(24)

which can be normalized and solved using, for instance, the Newton–Raphson method, for ${\textbf {q}}_{n+1}$. $\dot{{\textbf {q}}}_{n+1}$ can be then obtained using (24). The configuration-independent external forces are evaluated at time instant $t_{n\,+\,\frac{1}{2}}$, while the configuration-dependent forces induced by a surrounding flow, discussed in Appendix B, are approximated using the midpoint rule along with the approximations (23)–(24). We note that solving (25) using the Newton–Raphson method requires each term in (21)–(22) and force terms in cases of configuration-dependent forces to be linearized. We derive this and their resulting counterparts to the tangent stiffness matrix in Appendix C.

3.4 Outlier removal

Algorithm 1

Implicit time integration scheme employing the strong approach of outlier removal of [36].

Our empirical results, discussed in the next section, indicate that the excessive high-frequency contents of the response obtained with isogeometric discretizations lead to unstable computations. Based on this observation, we propose to improve the robustness of these discretizations via the strong approach of outlier removal introduced in [36], which entirely and merely removes the spurious outlier modes corresponding to the highest frequencies. The motivation to remove the outliers is based on the fact that in nonlinear dynamic analysis, the modes are coupled. Thus, when the outliers are excited at a time instance, other high-frequency modes might be excited at that time instance and later ones as well. We note that the modes mentioned here and throughout this work in the context of nonlinear analysis are the ones corresponding to a linearized problem at a time instance. Moreover, comparing our results with the one obtained with nodal finite elements employed in [17], we want to emphasize that the occurrence of outliers is purely due to the choice of the isogeometric discretization scheme, not the employed rod formulation [17] or the chosen time integration scheme.

To remove outliers, the fundamental idea of [36] is to strongly enforce sufficient higher-order continuity at the boundary in terms of additional boundary conditions. This idea originates from the observations made with a set of second- and fourth-order eigenvalue problems, where outliers occur due to the lower-order continuity at the boundary of the discrete mode shapes. In [36], the authors strongly enforce these additional constraints by constructing a new constrained subspace of the original space of B-splines basis functions. They computed the so-called extraction operator $\mathcal {C}$ that linearly combines the original basis functions, such that the resulting basis satisfies the enforced constraints. One can compute the operator $\mathcal {C}$ by finding a basis for the null space of a matrix including merely boundary constraints (see [36, Eq. (25a)]). $\mathcal {C}$ is computed for the complete space of spline basis functions but only modifies functions associated with the boundary conditions. It has the structure of a block diagonal matrix consisting of block matrices at each constrained boundary and the identity matrix (see [36, Eq. (26)]). $\mathcal {C}$ preserves the important properties of B-splines, such as non-negativity and minimum local support.

In this work, we apply the extraction operator $\mathcal {C}$ to both the right- and left-hand sides of the semi-discrete formulation (14) by multiplying with the matrix ${\textbf {C}}$ at each time step, where ${\textbf {C}}$ is the matrix expression of $\mathcal {C}$ [36]. Since $\mathcal {C}$ only depends on the chosen spline space and the boundary conditions, ${\textbf {C}}$ remains constant and can be computed once before the time integration procedure. It is associated with the property of the chosen spline space, does not depend on the choice of the time integration scheme, and thus can be applied together with other schemes than the one chosen in this work. ${\textbf {C}}$ inherits the structure of $\mathcal {C}$, i.e. ${\textbf {C}}$ is a block diagonal matrix. In general, ${\textbf {C}}$ is not a square matrix but a sparse matrix. We illustrate the implicit time integration scheme employing the outlier removal approach [36] in Algorithm 1. We note that in our computations, together with the additional boundary conditions for outlier removal, we also strongly enforce the essential boundary conditions via ${\textbf {C}}$. We apply ${\textbf {C}}$ to the global stiffness matrix ${\textbf {K}}$ before solving the matrix equation at each time step. Since ${\textbf {C}}$ is a sparse matrix, the product ${\textbf {C}}^T {\textbf {K}} \, {\textbf {C}}$ is also a sparse matrix. For the technical details of the computation of ${\textbf {C}}$, we refer to [36]. We note that there exist other approaches to remove outliers, either strongly [54] or weakly [55, 56]. For an overview of these approaches, we refer to [36, 57] and references therein.

An alternative line of approaches that also tackle spurious high-frequency modes such as outliers is to employ time integration schemes with high-frequency dissipation such as the generalized-$\alpha$ method [19, 58], the HHT-$\alpha$ method [59], the $\rho _{\infty }$-Bathe method [60]. Such schemes are robust and have found their wide applications in structural dynamics and fluid mechanics communities. In general, integration schemes with higher-frequency dissipation damp out the spurious high-frequency responses by introducing numerical damping, which can be regulated via a dissipation parameter. In this work, since we want to observe the occurrence and investigate the effect of outlier modes on the performance of isogeometric discretizations, we chose another time integration scheme that is the hybrid combination of the midpoint and trapezoidal rules and does not damp out the high-frequency modes. Furthermore, for the comparison purpose between the isogeometric discretization scheme and the standard scheme using nodal finite elements employed in [17], we employed the same time integration scheme as in [17]. We note that using schemes with high-frequency dissipation and choosing a zero dissipation parameter, i.e. zero numerical damping, also allows the occurrence of outliers. For the comparison between two discretization schemes, however, the integration scheme employed in [17] is the simplest choice for our work.

4 Robustness of isogeometric discretizations

In this section, we numerically demonstrate, for two- and three-dimensional benchmarks, that the isogeometric discretization scheme using cubic $C^1$ or $C^2$ splines is less robust than the standard one based on nodal finite elements using cubic Hermite functions. We then discuss important factors, such as the high-frequency contents of the response and round-off errors due to floating-point arithmetic, which may negatively affect the robustness of the employed discretization scheme. We show that employing the strong approach of outlier removal [36], discussed in the previous section, improves the robustness of isogeometric schemes. In addition, we discuss the influence of the configuration-dependent mass matrix (15) on the accuracy of the response.

4.1 Two- and three-dimensional benchmarks

Fig. 4

Deformed configurations of a clamped rod subjected to an in-plane loading at different time steps, computed with different discretizations and a geometrically exact beam model. (Color figure online)

Fig. 5

The energy of a clamped rod subjected to a vanishing in-plane loading, computed with different discretizations and a geometrically exact beam model. (Color figure online)

Fig. 6

Deformed configurations of an unconstrained rod subjected to a vanishing out-of-plane loading at different time steps, computed with different discretizations and a geometrically exact beam model. (Color figure online)

We consider the two dynamic benchmarks of [17, Sec. 5.2] and compare the responses obtained with isogeometric discretizations, using cubic $C^1$ and $C^2$ B-splines, against the standard scheme. In all cases, the rod is uniformly discretized into 20 elements, which leads to 126, 69, and 105 dofs using these three approaches, respectively. We also include the solution obtained with the geometrically exact beam model including shear and torsion using linear $C^0$ Lagrange polynomials as a reference, which is provided by [17].

The first example consists of an initially straight, transversely isotropic, clamped rob subjected to the following in-plane vanishing load at its free end:

$$\begin{aligned} {\textbf {F}}(t) \,=\, {\left\{ \begin{array}{ll} \frac{t}{0.5\,t_c} \, &{} {\textbf {F}}_c \,, \quad 0\,\le \,t\,\le \,0.5\,t_c \,, \\ \frac{2}{t_c}\, \left( \,t_c\,-\,t\,\right) \, &{} {\textbf {F}}_c \,, \quad 0.5\,t_c\,<\,t\,\le \,t_c \,, \\ 0 \, &{} {\textbf {F}}_c \,, \quad t \,>\, t_c \,, \end{array}\right. } \end{aligned}$$

(25)

where $t_c$ and ${\textbf {F}}_c$ are chosen to be $t_c\,=\,0.5$ s and ${\textbf {F}}_c \,=\, (0,\, 30,\, 0)$ N, respectively. Thus, the rod deforms in the xy-plane. The rod has an initial length of 10 with an initial director of $(1,\;0,\;0)$, a circular cross-section with a diameter of 0.01 m, Young’s modulus $E\,=\,2\cdot 10^{11}$ N/m$^2$, and mass density $\rho \,=\,7900$ kg/m$^3$. We choose the same time step of $\Delta \,t \,=\, 0.005$ s, a simulation time of 20 s, and a tolerance of $10^{-10}$ for the Newton–Raphson scheme, as done in [17]. Figure 4 illustrates the deformed configurations in a sequence of eleven load steps, obtained with cubic $C^1$ (black curves) and $C^2$ (green dashed curves) B-splines, with the standard scheme using cubic $C^1$ Hermite functions (blue dashed curves), and with the geometrically exact beam model using linear $C^0$ Lagrange polynomials (red dotted curves). We note that for this example, after 16.5 s, the computations using isogeometric discretizations become unstable and the Newton scheme does not converge anymore, while the one using the standard scheme remained stable during the computation time of 20 s. Nevertheless, we observe that during the first 16.5 s the isogeometric scheme using cubic $C^1$ splines and the standard one result in virtually identical responses, since their basis functions span the same space and have the same approximation power. Comparing these results to the one obtained with the geometrically exact beam, we see that the difference between these two models increases progressively in time. This may result from different employed time integration schemes for these models, as discussed in [17]. Focusing on the responses obtained with cubic $C^2$ B-splines (green dashed curves), we can see that these differ from those obtained with cubic $C^1$ B-splines at certain time steps.

Remark 5

For this benchmark, computations using the isogeometric scheme with quadratic $C^1$ B-splines become unstable already after 3.5 s. Thus, in this work, we did not apply quadratic B-splines for any dynamic benchmark.

To gain better insights, we illustrate the time history of the kinetic, potential, and total energy resulting from the aforementioned approaches in Fig. 5. We observe that, starting around 4 s, there is a phase shift between the responses obtained with cubic $C^2$ and $C^1$ B-splines (Fig. 5a,b). In particular, using cubic $C^2$ B-splines leads to a smaller phase than $C^1$ B-splines, which consequently leads to different responses (black and green curves) observed in Fig. 4. We note that for the case of cubic spline functions, our convergence study of the planar roll-up in Sect. 3.2 (see green curves in Fig. 3b,d,f), shows a slightly higher error level when using cubic $C^2$ than $C^1$ B-splines. This decrease in the approximation power of cubic splines with higher continuity for this nonlinear rod may relate to the different phases and responses observed in Fig. 4. Focusing on the responses around 16 s, before the energy, obtained with isogeometric discretizations, shoots up, indicating unstable computations, we see that high-frequency modes are excited and lead to fluctuations in the response. Due to a smaller phase, we expect that this generally occurs earlier when using cubic $C^2$ B-splines than $C^1$ splines. This necessary means that increasing the continuity of spline basis functions reduces the robustness of the corresponding discretizations. Nevertheless, before the computations become unstable, the isogeometric scheme approximately preserves the same total energy as the standard one and the geometrically exact beam model (Fig. 5c). We conclude that it is less robust than the standard one based on nodal finite elements [17]. Using discretizations with cubic splines of higher continuity reduces the phase of the responses, which may reduce their robustness and may relate to their reduced approximation power in this case of nonlinear rods [17].

Remark 6

We note that a weaker continuity than the $C^1$ continuity, the geometric $G^1$ continuity, can be applied for Kirchhoff rods (see e.g. [21, 61]) and obtained with isogeometric discretizations, which may achieve better robustness than those using $C^1$ continuity or higher. Enforcing the $G^1$ continuity is equivalent to enforcing the continuity of the unit director at the element interface. For more details of this condition and an overview of different approaches for its enforcement, we refer for instance to [21, 61] and references therein. Studying the performance of isogeometric discretizations with $G^1$ continuity for the nonlinear rod [17] and comparing them with those of $C^1$ continuity or higher is out of the scope of this work and thus is considered for future work.

Fig. 7

The energy of an unconstrained rod subjected to a vanishing out-of-plane loading, computed with different discretizations and a geometrically exact beam model. (Color figure online)

Fig. 8

Three components of the angular momentum of an unconstrained rod subjected to a vanishing out-of-plane loading, computed with different discretizations and a geometrically exact beam model. (Color figure online)

The second benchmark is an unconstrained rod of the same length and material as the clamped rod studied above. The rod has a smaller cross-section with a diameter of 0.005 m and is subjected to the vanishing load (26), with $t_c=0.5$ s and ${\textbf {F}}_c$ given by:

$$\begin{aligned}&{\textbf {F}}_c \,=\, (-30,\, -30,\, 0)\,\text {N} \quad \text {at } s\,=\,0 \,, \\&{\textbf {F}}_c \,=\, (30,\, 30,\, 0)\,\text {N} \quad \text {at } s\,=\,L \,, \\&{\textbf {F}}_c \,=\, (0,\, 0,\, -24)\,\text {N} \quad \text {at } s\,=\,L/20 \,, \\&{\textbf {F}}_c \,=\, (0,\, 0,\, 24)\,\text {N} \quad \text {at } s\,=\,19\,L/20 \,. \end{aligned}$$

Consequently, the rod deforms freely in three-dimensional space. For this benchmark, we also choose the same time step $\Delta \,t \,=\, 0.001$ s, a simulation time of 2 s, and a tolerance of $10^{-10}$ for the Newton–Raphson scheme, as done in [17]. Figure 6 illustrates the deformed configurations in a sequence of twelve load steps obtained with the aforementioned approaches. For this example, after 1.35 s the computation using cubic $C^2$ B-splines becomes unstable and the Newton scheme does not converge anymore, while that using cubic $C^1$ B-splines remains stable during the simulation time of 2 s. We have similar observations as in the case of the clamped rod above: both the isogeometric scheme using cubic $C^1$ B-splines and the standard one show the same accuracy, while the former using $C^2$ B-splines leads to distinct responses. This is also illustrated for the energy in Fig. 7 and the three components of the angular momentum in Fig. 8. We also observe a phase shift in the case of cubic $C^2$ B-splines, before the high-frequency modes affect the response and the computation becomes unstable. We see that, for this benchmark, using smoother splines of $C^2$ also leads to less robust computation than using $C^1$, as discussed in the first benchmark of a clamped rod.

4.2 Robustness improvement with outlier removal

The results discussed in the previous subsection indicate that excited high-frequency contents in the response reduce the robustness of isogeometric discretizations. To gain better insights, we perform the Fast Fourier Transformation (FFT) of the kinetic energy in the case of the clamped rod studied in the previous subsection, subjected to half of the load during a longer simulation time, $T=100$ s. To avoid noises in the transformed signal, we consider the energy as long as it is smaller than a threshold of 40 J. Figure 9 illustrates the FFT of the kinetic energy obtained with cubic $C^1$ B-splines (blue circles), where we can observe that not only one, but almost all high-frequency modes are excited.

Fig. 9

Fast Fourier transformation of the kinetic energy for the clamped rod in Fig. 4, subjected to a half of the load during a longer time period of 100 s. (Color figure online)

Table 1 Determinant of ${\tilde{{{\textbf {A}}}}}$ in Eq. (31)

Fig. 10

The energy of the clamped rod subjected to a vanishing in-plane loading in Fig. 4, computed with different basis functions, time steps, with and without outlier removal. (Color figure online)

As discussed in Sect. 3.4, we employ the strong approach of outlier removal [36] to remove the high-frequency contents in the response, such that the computations become more robust. In Fig. 9, we include the kinetic energy (green circles) with the outliers removed. We observe that almost all high-frequency contents are removed, as expected. Figure 10 shows the kinetic, potential, and total energy of the clamped rod studied in the previous subsection. Here, we present a comparison among the energy obtained with isogeometric discretizations using cubic $C^1$ and $C^2$ B-splines, employing outlier removal, (black and green curves, respectively), the energy obtained with the standard scheme [17] (blue dashed curve), and the geometrically exact beam model (red dashed curve). We observe that, as expected, the outlier removal approach improves the robustness of the isogeometric scheme. We note and recall the conclusions in [36] that the employed outlier removal approach does not affect the accuracy of the original discretizations and hence results in the same deformed configurations. We see that in the case using cubic $C^2$ B-splines (green curves), other remaining high-frequency modes are excited, as illustrated in the inset figure in Fig. 10c. Hence, we expect that in this case, the computation will become unstable at a later time, which might occur in the case using cubic $C^1$ B-splines as well. Moreover, we can see that the phase shift between cubic $C^1$ and $C^2$ B-splines is not affected by the outlier removal approach. We thus assume that this phase shift results from the different continuity of the spline basis functions, which may relate to their different approximation power, as discussed in the previous section. To investigate the effect of the time step on the robustness, we also include the results obtained with cubic $C^1$ B-splines using a smaller time step of $\Delta \,t \,=\, 0.0025$ s (purple curve) in Fig. 10. We see that decreasing the time step also improves the robustness of the discretization scheme. Focusing on the accuracy of the results, we observe that neither employing the outlier removal approach nor decreasing the time step, negatively affects the accuracy of the responses. We conclude that the robustness of the isogeometric scheme can be improved using the strong approach of outlier removal or by decreasing the time step, without compromising the obtained accuracy.

Remark 7

For the considered rods, we also studied the effect of increasing the polynomial degree p, either of $C^1$ or $C^{p-1}$ B-splines, on the responses, combining with the outlier removal approach [36]. We observe that irrespective of constant or increasing continuity, increasing p approximately leads to the same deformed configurations and energy responses before the computations become unstable. We also saw that in general, increasing p, however, results in responses containing more high-frequency contents, i.e. less robust computations at later time.

To gain a better understanding of how outlier modes and time step size affect the robustness of the discretization scheme, we consider the free vibration of a linear one-dimensional fourth-order problem such as an unconstrained Euler-Bernoulli beam. The semi-discrete equations of motion, in matrix form, are:

$$\begin{aligned} {\textbf {M}}\, \ddot{{\textbf {u}}} \,+\, {\textbf {K}}\,{\textbf {u}} \,=\, {\textbf {0}} \,, \end{aligned}$$

(26)

where ${\textbf {M}}$, ${\textbf {K}}$, and ${\textbf {u}}$ are the global mass matrix, global stiffness matrix, and the unknown displacement vector of the control points, respectively. Employing the implicit time integration scheme discussed in Sect. 3.3, the inertial and internal elastic terms are approximated at the time instant $t_{n+\frac{1}{2}}$ as follows:

$$\begin{aligned}&{\textbf {M}}\, \ddot{{\textbf {u}}}_{n+\frac{1}{2}} \,\approx \, \frac{{\textbf {M}}\, \dot{{\textbf {u}}}_{n+1}\,-\, {\textbf {M}}\, \dot{{\textbf {u}}}_{n}}{\Delta \,t}\,, \end{aligned}$$

(27)

$$\begin{aligned}&{\textbf {K}}\,{\textbf {u}}_{n+\frac{1}{2}} \,\approx \, \frac{{\textbf {K}}\,{\textbf {u}}_{n+1} \,+\, {\textbf {K}}\,{\textbf {u}}_{n}}{2}\,. \end{aligned}$$

(28)

Inserting these approximations into (27) and applying (24) to approximate $\dot{{\textbf {u}}}_{n+1}$, we obtain the following system of equations:

$$\begin{aligned} \underbrace{\begin{bmatrix} \Delta \,t \, {\textbf {K}} &{} 2\,{\textbf {M}} \\ -2\,{\textbf {I}} &{} \Delta \,t \, {\textbf {I}} \end{bmatrix}}_{{\textbf {A}}_L} \,&\begin{bmatrix} {\textbf {u}}_{n+1} \\ \dot{{\textbf {u}}}_{n+1} \end{bmatrix} \,=\, \underbrace{\begin{bmatrix} -\Delta \,t \, {\textbf {K}} &{} 2\,{\textbf {M}} \\ -2\,{\textbf {I}} &{} -\Delta \,t \, {\textbf {I}} \end{bmatrix}}_{{\textbf {A}}_R} \, \begin{bmatrix} {\textbf {u}}_{n} \\ \dot{{\textbf {u}}}_{n} \end{bmatrix} \nonumber \\&\begin{bmatrix} {\textbf {u}}_{n+1} \\ \dot{{\textbf {u}}}_{n+1} \end{bmatrix} \,=\, \underbrace{{\textbf {A}}_L^{-1} \, {\textbf {A}}_R}_{\tilde{{\textbf {A}}}} \, \begin{bmatrix} {\textbf {u}}_{n} \\ \dot{{\textbf {u}}}_{n} \end{bmatrix} \, . \end{aligned}$$

(29)

In linear cases, the stiffness and mass matrices, and thus $\tilde{{\textbf {A}}}$, do not depend on the configuration, and remain constant during the time integration. Consequently, we can compute the response at time step $t_{n+1}$ in terms of the initial conditions as follows:

$$\begin{aligned} \begin{bmatrix} {\textbf {u}}_{n+1} \\ \dot{{\textbf {u}}}_{n+1} \end{bmatrix} \,=\, \underbrace{\tilde{{\textbf {A}}} \, \ldots \, \tilde{{\textbf {A}}}}_{(n+1) \text { times}} \, \begin{bmatrix} {\textbf {u}}_{0} \\ \dot{{\textbf {u}}}_{0} \end{bmatrix} \,=\, \tilde{{\textbf {A}}}^{n+1} \, \begin{bmatrix} {\textbf {u}}_{0} \\ \dot{{\textbf {u}}}_{0} \end{bmatrix} \, . \end{aligned}$$

(30)

The requirement for (31) to have an unique solution is that the matrix $\tilde{{\textbf {A}}}$ is a convergent matrix [62]. This necessarily requires $\left| \det \left( \tilde{{\textbf {A}}} \right) \right| \, < \, 1$. Table 1 illustrates the determinant of $\tilde{{\textbf {A}}}$ computed with cubic Hermite functions, and cubic $C^1$ B-splines with and without outlier removal, using different time steps and meshes. We observe that the determinant is either smaller or larger than 1.0 with a tolerance in the range of $\left[ 10^{-13},\,10^{-15}\right]$, i.e. machine accuracy. Such a problem can be attributed, in principle, to round-off errors due to floating-point arithmetic, where the numerical deviation from 1.0 highly depends on the time step, mesh size, and basis functions, i.e. on the time integration and spatial discretization schemes. We can also see that neither using cubic Hermite functions, employing outlier removal, nor reducing the time step, guarantees that $\left| \det \left( \tilde{{\textbf {A}}} \right) \right| \, < \, 1$. Hence, there is no indication that one discretization scheme ensures the existence of the solution at an arbitrary time step $t_{n+1}$ for (31), and the other schemes do not. For nonlinear problems, since $\tilde{{\textbf {A}}}$ is configuration-dependent, it is not trivial to identify parameters and calibrate them in order to ensure $\left| \det \left( \tilde{{\textbf {A}}} \right) \right| \, < \, 1$. In addition, there might be other factors that are decisive for the existence of the solution at $t_{n+1}$, and thus for the robustness of the discretization scheme. This requires further study and investigation, which we postpone to future work.

4.3 On the influence of the configuration-dependent mass matrix

Fig. 11

The relative change of the deformed configuration of an unconstrained rod computed with different scaling factors $\alpha$. (Color figure online)

We now briefly discuss the behavior of the configuration-dependent mass matrix (15), and whether this counterpart can be omitted. We recall that the mass matrix (15), ${\textbf {M}}$, consists of two counterparts, ${\textbf {M}}_1$ and the rotational counterpart ${\textbf {M}}_2$. While ${\textbf {M}}_1$ does not depend on the current rod configuration and remains constant, ${\textbf {M}}_2$ is configuration-dependent. In most of the common nonlinear dynamic analyses, ${\textbf {M}}_2$ is considered a small perturbation and thus is neglected, which results in a constant configuration-independent mass matrix.

To investigate the effect of ${\textbf {M}}_2$ on the system response and whether ${\textbf {M}}_2$ behaves as a regular perturbation, we consider the unconstrained rod benchmark studied in Sect. 4.1, computed with the isogeometric discretization scheme using cubic $C^1$ B-splines. In particular, we scale ${\textbf {M}}_2$ by a factor $\alpha \,\in \, [0,\,1]$. The mass matrix (15) then becomes:

$$\begin{aligned} {\textbf {M}} \,=\, {\textbf {M}}_1 \,+\, \alpha \, {\textbf {M}}_2 \,. \end{aligned}$$

(31)

Figure 11 illustrates the relative change in the $L^2$ norm of the deformed rod, $\varvec{\varphi }_h$, during the simulation time of 2 s. Different colors correspond to different values of $\alpha$. We observe that this relative change increases in time and with increasing $\alpha$. This observation implies that, for the studied rod, ${\textbf {M}}_2$ does not behave as a regular perturbation, and thus should not be neglected. Mathematical proof and analysis of this counterpart of the mass matrix is out of the scope of this paper and is postponed to future works.

5 Application to a swinging rubber rod

In this section, we apply the nonlinear rod formulation [17], discretized with isogeometric discretizations, to a swinging rubber rod, which can represent for instance a sort of mooring lines or cables. Such structures are essential components, for instance, of offshore wind turbines, pumping kites, oil and gas platforms, etc. Hence, our simulation is a part and intermediate step towards the aero-hydro-elastic simulation framework accounting for such applications. We consider conservative and non-conservative forces such as gravity, forces induced by surrounding wind or water, and pulsating forces. We illustrate that our rod formulation is able to represent the rod nonlinear behavior that is a combination of elastic vibrations and rigid body oscillations around a static position and shape, which deforms differently at different force frequencies. Based on our results in the previous section, we spatially discretize the swinging rod with cubic $C^1$ B-splines ($p=3$) and improve its robustness with the strong approach of outlier removal. We start with benchmarking our approach via an example of a swinging rod under gravitational loading [63, 64].

5.1 Benchmarking

We consider an initially straight rod of length $L=1.0$ m, with a circular cross-section, a radius of $5 \,\cdot \, 10^{-3}$ m, Young’s modulus $E\,=\,5\cdot 10^{6}$ N/m$^2$, and mass density $\rho \,=\, 1100$ kg/m$^3$. The rod is subjected merely to gravitational loading [63], which we simulate with a direction of $(0,\;-1,\;0)$ while the initial director of the rod is $(1,\;0,\;0)$. Thus, the rod deforms in the xy-plane. We choose a time step of $\Delta \,t \,=\, 0.01$ s, a simulation time of 2.4 s, and a tolerance of $10^{-10}$ for the Newton–Raphson scheme.

Fig. 12

Deformed configurations of a swinging rubber rod due to gravity, computed with cubic $C^1$ B-splines ($p=3$) and outliers removed. (Color figure online)

Fig. 13

The tip displacements of the swinging rubber rod in Fig. 12. (Color figure online)

Fig. 14

The energy of the swinging rubber rod in Fig. 12. (Color figure online)

Figure 12 illustrates the deformed configuration of the studied rod at twelve time steps, and Fig. 13 shows the time evolution of the horizontal ($u_x(t)$ in blue) and the vertical displacement ($u_y(t)$ in red) at the rod free-end. We observe in Fig. 12 that the swinging rod behaves similarly to an elastic pendulum, and thus has a stable equilibrium configuration when it is aligned with the $y-$axis and its free-end is located at $(x,y,z)=(0,-1,0)$ m. Due to the highly nonlinear nature captured by the current formulation, the rod exhibits large elastic rotations and displacements in time. Comparing these results with those of [63, Fig. 8,9], we observe an apparent difference after 2.0 s, for instance, the deformed configuration at $t=2.4$ s clearly distinguishes from that in [63, Fig. 8]. We further compare the energy obtained with our approach, illustrated in Fig. 14, with that in [63, Fig. 10]. We see that the total energy and its counterparts show approximately the same time history as the reference, with an exception of around 2.4 s. While the reference and our results both consist of high-frequency contents around this time, the reference energy does not jump to a large value as ours. This is consistent with different deformed configurations after 2.0 s observed in Fig. 12. We conclude that despite the absence of outliers in our computation, our approach using cubic $C^1$ B-splines is less robust than the scheme using classical nodal finite elements with cubic Hermite functions applied in [63]. Nevertheless, it leads to the same behavior for the swinging rod, given that the computation is stable.

5.2 Dynamic response to wind forces

Fig. 15

A wind profile with constant amplitude and rotating direction around the z-axis

Fig. 16

Deformed configurations of a swinging rubber rod due to gravity and the wind profile in Fig. 15, computed with cubic $C^1$ B-splines ($p=3$) and outliers removed

Fig. 17

The tip displacements of the swinging rubber rod in Fig. 16. (Color figure online)

Fig. 18

The energy of the swinging rubber rod in Fig. 16. (Color figure online)

We now consider the same rod as in the previous benchmark, subjected to an additional wind field with a wind profile illustrated in Fig. 15. The wind velocity has constant values but changing direction along the $z-$axis, that is:

$$\begin{aligned} \begin{aligned}&{\textbf {v}}_{wind} \, (z) \,=\, v_0 \,{\textbf {d}}_{wind}\,(z) \,=\, \\&v_0 \, \left[ {\textbf {E}}_1 \, \cos \left( \beta _0 - 2\,\beta _0 \,\frac{z}{L}\right) \,+\, \right. \\&\qquad \left. {\textbf {E}}_2 \, \sin \left( \beta _0 - 2\,\beta _0 \,\frac{z}{L}\right) \right] , \end{aligned} \end{aligned}$$

(32)

Fig. 19

The mean value and the amplitude of the horizontal tip displacement $u_x$ of an aluminum swinging rod at steady state, computed with cubic $C^1$ B-splines ($p=3$) and outliers removed

where $v_0$ is the constant value of the wind velocity, ${\textbf {d}}_{wind}$ is the director of the wind velocity, $\beta _0$ is the angle of ${\textbf {d}}_{wind}$ with respect to the $y-$axis at $z=0$, and L is the initial length of the rod. For this example, we choose $v_0 = 10$ m/s, $\beta _0 = 45^{\circ }$, a simulation time of 30 s, and an initial angle of the rod with respect to the $x-$axis of $-15^{\circ }$, that is an initial rod director of ${\textbf {d}}=(\cos (\pi /12),\;0,\;-\sin (\pi /12))$. We apply the same discretization as in the previous benchmark. Figure 16 illustrates the motion sequence of the swinging rod under the considered loading. Due to the three-dimensional wind profile, the rod shows out-of-plane deformations. Focusing on the configurations at $t=8.75$ s and $t=10$ s, which are no longer visually distinguishable, we assume that the rod has approximately reached a steady state configuration. This behavior is also illustrated in the time history of the rod tip displacement, showed in Fig. 17, and of the energy of the system, showed in Fig. 18. We observe that after about 15 s, all displacement components remain approximately constant, which is consistent with approximately zero kinetic energy at this time, i.e. the rod reaches an equilibrium configuration. Comparing this with the behavior of the rod in the previous benchmark without wind force, we can see that the wind force dampens the motion of the rod to an equilibrium configuration. This damping characteristic, generally of aerodynamic forces, can be seen when it is approximated by a Taylor expansion as follows:

$$\begin{aligned} {\textbf {F}}_f = \left. \frac{\partial {\textbf {F}}_f}{\partial {\textbf {q}}}\right| _{{\textbf {q}}={\textbf {0}}}{\textbf {q}} + \left. \frac{\partial {\textbf {F}}_f}{\partial \dot{{\textbf {q}}}}\right| _{\dot{{\textbf {q}}}={\textbf {0}}}\dot{{\textbf {q}}} + \left. \frac{\partial {\textbf {F}}_f}{\partial \ddot{{\textbf {q}}}}\right| _{\ddot{{\textbf {q}}}={\textbf {0}}}\ddot{{\textbf {q}}} +..., \end{aligned}$$

(33)

Fig. 20

The amplitude of three force components induced by surrounding water at rest, integrated over the aluminium swinging rod at steady state. (Color figure online)

where the first term gives rise to the lift/thrust forces, the second term contains the so-called damping forces, and the third term is related to the added mass forces. An aerodynamic damping matrix can then be defined as ${\textbf {D}}_{\text {fluid}}=-\left. \frac{\partial {\textbf {F}}_f}{\partial \dot{{\textbf {q}}}}\right| _{\dot{{\textbf {q}}}={\textbf {0}}}$, which is a function of the free-stream velocity, among other parameters. The aerodynamic damping strongly depends on the magnitude of the free-stream velocity $V_{\infty }=\left| {\textbf {V}}_{\infty }\right|$. When $V_{\infty }$ is below a critical velocity $V_{\infty }^C$ (subcritical condition), the damping is positive and the surrounding flow absorbs the energy of the structure. When $V_{\infty }=V_{\infty }^C$ (critical condition), the damping is zero and then the surrounding flow does not absorb or supply the structure energy. When $V_{\infty }>V_{\infty }^C$ (supercritical condition), the damping becomes negative and the surrounding flow supplies energy to the structure, i.e. favoring the emergence of fluid–structure interaction instabilities such as aeroelastic flutter. For the swinging rod studied here, we are in a subcritical condition since the rod oscillation is damped out over time. Focusing on the robustness of the employed discretization scheme, we observe that the wind force also dampens out the high-frequency modes, and thus improves its robustness. We conclude that, given the subcritical condition, employing damping forces is another approach to improve the robustness, particularly of the isogeometric discretization scheme studied in this paper.

5.3 Dynamic response to pulsating forces

Fig. 21

Deformed configurations of an aluminum swinging rod at force frequencies of 0.88 Hz and 4.90 Hz, computed with cubic $C^1$ B-splines ($p=3$) and outliers removed

Fig. 22

The tip displacements of an aluminum swinging rod computed for force frequencies of 0.88 Hz and 4.90 Hz. (Color figure online)

Our last example is a long rod submerged in water, which simulates a sort of mooring lines or cables considered in aero-hydro-elastic simulations and offshore wind engineering. To this end, we consider a combination of the gravitational field, the surrounding water at rest, and a horizontal pulsating force applied at the free end of the swinging rod. We choose a sinusoidal pulsating force that is:

$$\begin{aligned} \varvec{F}\,(t) \,=\, A_F \, \sin \left( \frac{\omega _F}{2\,\pi } \, t\right) \, \varvec{E}_1 \,, \end{aligned}$$

(34)

where $A_F$ and $\omega _F$ are the force amplitude and angular frequency, respectively. For this example, we choose an amplitude $A_f \,=\, 350.0$ kN, and different frequencies $(\omega _F/2\pi )$ ranging from 0.1 Hz to 8 Hz. The surrounding water is at rest, i.e. ${\textbf {V}}_{\infty }={\textbf {0}}$ and ${\textbf {a}}_{\infty } = {\textbf {0}}$, and the mass density of the surrounding water is $\rho _f \,=\, 1000$ kg/m$^3$. In this example, we consider an aluminium rod with an initial director ${\textbf {d}}=(0,\;-1,\;0)$, a length of 250 m, a circular cross-section with a radius of 0.02 m, Young’s modulus $E\,=\,7\cdot 10^{10}$ N/m$^2$, and mass density $\rho \,=\, 2700$ kg/m$^3$. We simulate the gravitational field with a direction of $(0,\;-1,\;0)$, which allows the rod to deform only in the xy-plane. We compute this example for 1000 s and 2000 s such that a steady state is observed during the last 100 s, which is identified as the response becomes periodic in time.

Figure 19 illustrates the mean value and the amplitude of the horizontal displacement $u_x$ at the free tip of the studied swinging rod after it reaches the steady state, as a function of the force frequencies $(\omega _F/2\pi )$. We observe that the mean value of the tip displacement (Fig. 19a), which corresponds to the position, around which the rod oscillates, i.e. the static equilibrium position, jumps at force frequencies smaller than 2.0 Hz and those larger than 5.0 Hz. It means that the equilibrium position and configuration change abruptly when the force frequency changes, which indicates a nonlinear behavior of the rod in different frequency ranges. Focusing on the amplitude of the tip displacement in Fig. 19b, we can see that after a steep increase at low frequencies, it decreases with increasing force frequency and slightly jumps downwards around 4.9 Hz (Fig. 20). This relation between the amplitude and the force frequency can be explained based on the linear theory. Consider a damped harmonic oscillator of mass m, spring stiffness k, and damping coefficient b, subjected to a sinusoidal pulsating force, $f_0\sin (\Omega t)$, of frequency $\Omega$ and amplitude $f_0$. It is well-known from classical vibration theory that the response amplitude of this linear system is:

$$\begin{aligned} A = \frac{f_0}{\sqrt{m(\Omega ^2-\omega _0^2)^2+b^2\Omega ^2}} \end{aligned}$$

(35)

where $\omega _0=\sqrt{k/m}$ is the natural frequency of the system. We can see that for the cases of $\Omega < \omega _0$, increasing $\Omega$ increases the amplitude A, and for the cases of $\Omega > \omega _0$, increasing $\Omega$ decreases A. Since the studied nonlinear rod can be considered as a damped (by considering the surrounding flow) distributed-parameter elastic pendulum subjected to a pulsating force at its free end, a similar relation between the amplitude and the force frequency is expected for both linear and the nonlinear cases. This supports our observation in Fig. 19b, except for the jump at 4.9 Hz, which associates with the highest positive peak of the mean value in Fig. 19a.

To gain better insights into the rod behavior, in Fig. 21, we plot some deformed configurations at different times for 0.88 Hz and 4.9 Hz, i.e. before and close to the frequency value where the amplitude jumps. We observe that for both two frequencies, the rod behavior is a combination of elastic vibrations and rotations as a rigid body that oscillates around a static equilibrium position. While this equilibrium position is along the vertical axis in the case of 0.88 Hz (Fig. 21a), it, as well as the equilibrium configuration, changes in the case of 4.9 Hz (Fig. 21b), since the lower part of the rod is bent to a horizontal segment, and the rod then oscillates around the new deformed position. This behavior together with the amplitude jump is also observed in the time history of the tip displacement, illustrated in Fig. 22 for 0.88 Hz and 4.9 Hz. We can also see that both two displacement components oscillate over time, while their mean value jumps to a value around 100 m in the case of 4.9 Hz. This is further reflected in the relation between the pulsating force frequency and the amplitude of forces induced by the surrounding flow, illustrated in Fig. 20. In Fig. 20, we illustrate the amplitude $A_f^{\omega }$ of three components of the resulting force induced by the surrounding flow (added mass, normal drag, and tangential drag), integrated over the rod length once the rod reaches the steady state, that is:

$$\begin{aligned} A_f^{\omega }(\omega _F)=\underset{\text {steady state}}{\max }\left[ f_f^t(t)\right] -\underset{\text {steady state}}{\min }\left[ f_f^t(t)\right] \,, \end{aligned}$$

where the integrated force $f_f^t(t)$ is:

$$\begin{aligned} f_f^t(t)=\int _0^L\left| {\textbf {F}}_f(s,t)\right| \,ds \,, \end{aligned}$$

as a function of the pulsating force frequency. We observe that the amplitude of the normal and tangential drag forces also jumps approximately at 4.9 Hz. In particular, the normal drag drops, while the tangential drag jumps upwards, which is consistent with the fact that the lower part of the rod is bent to a horizontal segment, i.e. only the upper part is mainly affected by the norm drag, while the tangential drag is the dominating force in the lower part. Regarding the robustness of the applied discretization scheme, we did not obtain unstable computations and results containing high-frequency contents. We conclude that different nonlinear behaviors of the swinging rod can be represented and studied using the rod formulation [17], together with the isogeometric discretization scheme, improved using outlier removal, and with an energy-momentum preserving implicit time integration scheme. This has been shown to be a sufficiently robust approach for studying nonlinear structures such as swinging rods modeling mooring lines, which can be further investigated for complex behaviors such as parametric resonances and chaotic behavior in future works.

6 Summary and conclusions

In this paper, we explored the application of the nonlinear formulation [17] for rods that exhibit only axial and bending deformations, using isogeometric spatial discretizations. Our results illustrate different convergence rates for odd and even polynomial degrees, which is known from other isogeometric discretization methods, see e.g. [65, 66] and the references therein. They also show that the continuity generally does not affect the accuracy and convergence, except the convergence in $H^2$ semi-norm when using odd degrees, which decreases with decreasing continuity. We demonstrated computationally for studied dynamic benchmarks of two- and three-dimensional rods that isogeometric discretizations using B-splines with $C^1$ continuity or higher are less robust than the standard scheme using Hermite functions. Increasing the smoothness of cubic spline basis functions leads to a smaller phase, which may reduce the robustness and may relate to the approximation power of the corresponding discretizations with higher continuity in the case of nonlinear rods [17]. We showed that robustness can be improved via a strong approach of outlier removal [36] without compromising the accuracy. Alternatively, reducing the time step or employing forces with damping characteristics leads to more robust computations. We have shown that the robustness is closely related to round-off errors due to floating-point arithmetic. In addition, we demonstrated computationally for an unconstrained rod subjected to out-of-plane vanishing forces that the configuration-dependent mass matrix does not behave as a regular perturbation and thus cannot be simplified to a constant matrix. Lastly, we applied our nonlinear transient formulation to a swinging rubber rod subjected to gravity, forces induced by a surrounding flow such as wind and water, and a pulsating force of different frequencies. Our results also show that the isogeometric discretization scheme is robust and reliable for such an analysis.

The results presented in this work open up several directions for future works. One open question is the accuracy and convergence behavior of the discretization scheme in nonlinear problems. Another question is the irregular behavior of the configuration-dependent mass term that, to our best knowledge, is not yet proved analytically. This is particularly interesting for choosing and developing an efficient and robust time integration scheme. It is also desirable to employ our nonlinear formulation with the isogeometric discretization scheme in highly nonlinear problems with complex loads and geometries.

Appendices

Appendix A: An elastic pendulum

Fig. 23

Schematic of an elastic pendulum with two degrees of freedom

Fig. 24

Energy, the third component of angular momentum and the second precision quotient of the elastic pendulum in Fig. 23 with the initial conditions (A1). (Color figure online)

Fig. 25

Energy, the third component of angular momentum and the second precision quotient of the elastic pendulum in Fig. 23, subjected to gravity and a parabolic wind profile, with the initial conditions (A4). (Color figure online)

To benchmark our implementation regarding the time integration scheme and to verify that the chosen implicit time integration scheme achieves second-order accuracy, approximately preserves energy, and exactly preserves the linear and angular momentum [17, 33], discussed in Sect. 3.3, we consider the dynamics of a two degrees-of-freedom elastic pendulum, illustrated in Fig. 23. Its configuration space is $\mathcal {D}_{\text {pen}}=\mathbb {R}\times SO(2)$, where SO(2) is the special orthogonal group in two dimensions. For this pendulum, the two degrees of freedom are the angular rotation with respect to the vertical axis, $\theta (t)$, and the axial deformation of the spring, $\eta (t)$. Unit vectors ${\textbf {E}}_1$ and ${\textbf {E}}_2$ are the first two canonical Cartesian bases of $\mathbb {R}^3$, $L_0$, k, m, and g are the spring’s natural length, spring stiffness, mass, and gravitational acceleration, respectively. We consider a simulation time of 30 s and a time step $\Delta \,t\,=\,0.005$ s.

Table 2 Total energy and the angular momentum of an elastic pendulum illustrated in Fig. 24

We first consider the case of no gravitational acceleration ($g\,=\,0$) and no external load with the following initial conditions:

$$\begin{aligned}&\theta (t\,=\,0) = 0\,, \qquad \eta (t\,=\,0) = 0.1\,, \end{aligned}$$

(36)

$$\begin{aligned}&\dot{\theta }(t\,=\,0) = -0.5\,, \qquad \dot{\eta }(t\,=\,0) = 0.25\,. \end{aligned}$$

(A1a)

Fig. 24a–b illustrates the time evolution of the energy and the third component of the angular momentum, $j_3$. We note that for the studied pendulum, the first and second components of the angular momentum are zero, and thus are not illustrated here. We observe that both the total energy and $j_3$ are virtually constant over the simulation time. Table 2 shows their values at four snapshots. We can see that the total energy is approximately preserved, while the angular momentum is exactly preserved up to a tolerance of the machine accuracy. Figure 24c–d shows the second precision quotient of the response during the simulation time, which is computed to verify the implementation correctness of a numerical integration scheme and is defined as [67]:

$$\begin{aligned} Q_{II} \,=\, Q_{II}(t) \,=\, \frac{\left\| {\textbf {u}}_{\Delta \,t} \,-\, {\textbf {u}}_{\Delta \,t/2} \right\| }{\left\| {\textbf {u}}_{\Delta \,t/2} \,-\, {\textbf {u}}_{\Delta \,t/4} \right\| } \,, \end{aligned}$$

(A1b)

where ${\textbf {u}}\,=\,{\textbf {u}}(t)$ is a time-dependent variable, and its subscripts, $\Delta t$, $\Delta t/2$,..., denote the time step employed to compute ${\textbf {u}}$. For well-implemented time integration schemes, the precision quotient $Q_{II}$ is approximately $2^q$, where q is the order of accuracy of the integration scheme [67]. For our implicit scheme, $q=2$. We note that our choice of the second precision quotient instead of the first precision quotient is due to the lack of an analytical solution for the studied elastic pendulum [67]. We observe in Fig. 24c–d the expected value of $2^2 = 4$ for $Q_{II}$, which verifies the correctness of our implementation.

We then consider the load case consisting of the gravitational field and the following parabolic wind profile:

$$\begin{aligned}&{\textbf {v}}_w (x_1\,t) \,=\, x_1^2 \, \left( \, 1.0 \,+\, 0.1\,\sin \left( \,\frac{1}{50}\,\frac{\sqrt{k\,/\,m}}{2\,\pi }\,t\right) \,\right) \, {\textbf {E}}_2 \,. \end{aligned}$$

(A2)

We choose the following initial conditions of the pendulum:

$$\begin{aligned}&\theta (t\,=\,0) = \pi /2\,, \qquad \eta (t\,=\,0) = 0.0\,, \end{aligned}$$

(A3)

$$\begin{aligned}&\dot{\theta }(t\,=\,0) = 0.0\,, \qquad \dot{\eta }(t\,=\,0) = 0.0\,. \end{aligned}$$

(A4a)

We observe in Fig. 25a-b that the total energy and angular momentum are no longer preserved due to the presence of external forces, as expected. We can also see that after approximately 17 s, the pendulum does not contain any kinetic energy, i.e. it approximately achieves a static equilibrium state. The second precision quotient $Q_{II}$, illustrated in Fig. 25c–d, also implies a second-order accurate time integration scheme in this case. We conclude that the chosen implicit time integration scheme approximately preserves the total energy, exactly preserves the angular momentum, and is second-order accurate.

Appendix B: Simplified model of a force field induced by a surrounding flow

Fig. 26

Schematic of a wind profile. (Color figure online)

For the numerical examples included in subsequent sections, we consider forces induced by a surrounding flow such as wind and water. Figure 26 illustrates the schematic of an exemplary wind profile acting on a rod configuration. We consider a simplified model of such forces [68] that consist of three counterparts: the resulting force due to the accelerated surrounding flow by the moving rod, which we refer to as the added mass force, the tangential drag force, and the normal drag force. To this end, the resulting force ${\textbf {F}}_{f}$ per unit length at an arbitrary point of the rod is expressed as follows:

$$\begin{aligned} {\textbf {F}}_f = C_1{\textbf {a}}_N + C_2\left| {\textbf {V}}_N\right| {\textbf {V}}_N + C_3\left| {\textbf {V}}_T\right| {\textbf {V}}_T , \end{aligned}$$

(A4b)

where ${\textbf {a}}_N$ is the normal component of the relative flow acceleration with respect to the rod, ${\textbf {V}}_N$ and ${\textbf {V}}_T$ the normal and tangential component of the relative flow velocity, respectively, and the coefficients $C_1$, $C_2$ and $C_3$ are given by:

$$\begin{aligned} \begin{aligned} C_1&= \frac{1}{4}\pi \,C_M\rho _f\,{\o }^2 , \\ C_2&= \frac{1}{2}C_N\rho _f\,{\o }, \\ C_3&= \frac{1}{2}C_T\rho _f\,{\o }, \end{aligned} \end{aligned}$$

(B5)

where $\rho _f$ is the mass density of the surrounding flow, and ${\o }$ is the diameter of the cylindrical cross-section of the rod. The coefficients $C_M$, $C_N$, and $C_T$ depend on the Reynolds number and are commonly determined experimentally [68].

To determine the forces acting on a rod segment of length ds due to the added mass, normal and tangential counterparts, it is necessary to describe the relative flow acceleration and velocity in terms of the ambient motion of the fluid and the rod. Considering the kinematic description of the rod given in the previous subsection, the quantities ${\textbf {a}}_N$, ${\textbf {V}}_N$, and ${\textbf {V}}_T$ can be expressed as follows:

$$\begin{aligned} \begin{aligned} {\textbf {a}}_N&= \mathcal {P}_{{\textbf {d}}}\,{\textbf {a}}, \\ {\textbf {V}}_N&= \mathcal {P}_{{\textbf {d}}}{\textbf {V}},\quad \text {and} \\ {\textbf {V}}_T&= \left( {\textbf {d}}\otimes {\textbf {d}}\right) {\textbf {V}}, \end{aligned} \end{aligned}$$

(B6)

where ${\textbf {V}}={\textbf {V}}_{\infty }-\dot{\varvec{\varphi }}$, and ${\textbf {a}}={\textbf {a}}_{\infty }-\ddot{\varvec{\varphi }}$. The magnitude and direction of the free-stream velocity ${\textbf {V}}_{\infty }(z,t)$ is a function of the height z above the ground level (or below the sea level if ocean structures are considered) and time t. The free-stream acceleration is the time derivative of the free-stream velocity, i.e., ${\textbf {a}}_{\infty }(z,t)=\frac{\partial }{\partial t}{\textbf {V}}_{\infty }(z,t)$. Inserting (B7) into (B5) and integrating along the rod, we obtain the resultant forces as follows:

$$\begin{aligned} {\textbf {F}}_f&= C_1\,\int _{0}^{S}\mathcal {P}_{{\textbf {d}}}\,({\textbf {a}}_{\infty }-\ddot{\varvec{\varphi }})\,ds \\&\quad + C_2\int _{0}^{S}\left| \mathcal {P}_{{\textbf {d}}}\,({\textbf {V}}_{\infty }-\dot{\varvec{\varphi }})\right| \,\mathcal {P}_{{\textbf {d}}}\,({\textbf {V}}_{\infty }-\dot{\varvec{\varphi }})\,ds \\&\quad + C_3\,\int _{0}^{S}\left| \left( {\textbf {d}}\otimes {\textbf {d}}\right) ({\textbf {V}}_{\infty }-\dot{\varvec{\varphi }})\right| \left( {\textbf {d}}\otimes {\textbf {d}}\right) ({\textbf {V}}_{\infty }-\dot{\varvec{\varphi }})\,ds. \end{aligned}$$

We note that coefficients $C_i$, $i=1,2,3$, are assumed to be independent of the position along the rod.

Appendix C: Linearization

1.1 C.1 Preliminaries

We recall the semi-discrete formulation for the rod motions in Sect. 3.1, that is:

$$\begin{aligned} \begin{aligned} {\textbf {g}}_d = \int _0^S \,&\left( \, {\textbf {M}} ({\textbf {q}}) \, \nabla _{\dot{{\textbf {q}}}} \, \dot{{\textbf {q}}} \,+\, {\textbf {B}} ({\textbf {q}})^T \, \varvec{\sigma }_h \right. \\&\qquad \left. -\, {\textbf {N}}^T \, {\textbf {f}}^{\text {ext}} \,\right) \, \textrm{d}s \,=\, {\textbf {0}}. \end{aligned} \end{aligned}$$

(B7)

Analogously, the variational formulation that describes the equilibrium of the rod is:

$$\begin{aligned} {\textbf {g}}_s = \int _0^S \, \left( \, {\textbf {B}} ({\textbf {q}})^T \, \varvec{\sigma }_h \,-\, {\textbf {N}}^T \, {\textbf {f}}^{\text {ext}} \,\right) \, \textrm{d}s \,=\, {\textbf {0}}. \end{aligned}$$

(C8)

We employ a standard approach based on a Taylor expansion of (C9) and (C8) to obtain the tangent stiffness matrix associated with ${\textbf {g}}_s$ and ${\textbf {g}}_d$, respectively. To this end, Taylor’s approximation for a vector function ${\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}})$ is given by:

$$\begin{aligned}&{\textbf {g}}({\textbf {q}}+\Delta {\textbf {q}},\dot{{\textbf {q}}}+\Delta \dot{{\textbf {q}}},\ddot{{\textbf {q}}}+\Delta \ddot{{\textbf {q}}}) = \nonumber \\&\quad {\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}}) + \text {D}\,{\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}})\cdot \left( \Delta {\textbf {q}},\Delta \dot{{\textbf {q}}},\Delta \ddot{{\textbf {q}}}\right) \nonumber \\&\quad + \text {D}^2\,{\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}}):\left( \left( \Delta {\textbf {q}},\Delta \dot{{\textbf {q}}},\Delta \ddot{{\textbf {q}}}\right) \otimes \left( \Delta {\textbf {q}},\Delta \dot{{\textbf {q}}},\Delta \ddot{{\textbf {q}}}\right) \right) \nonumber \\&\quad + O\left( \left\| \Delta {\textbf {q}}\right\| ^3,\left\| \Delta \dot{{\textbf {q}}}\right\| ^3,\left\| \Delta \ddot{{\textbf {q}}}\right\| ^3,...,\left\| \Delta \dot{{\textbf {q}}}\right\| ^2\left\| \Delta \ddot{{\textbf {q}}}\right\| ,...\right) \,, \end{aligned}$$

(C9)

where $\text {D}^i(\cdot )$, $i=1,2,...$, is a $(i+1)$-order tensor of type (0, i), and ( : ) denotes the double-contraction tensor operation. Assuming that the higher-order terms are negligible, we obtain:

$$\begin{aligned}&{\textbf {g}}({\textbf {q}}+\Delta {\textbf {q}},\dot{{\textbf {q}}}+\Delta \dot{{\textbf {q}}},\ddot{{\textbf {q}}}+\Delta \ddot{{\textbf {q}}}) \approx \nonumber \\&\qquad {\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}}) + \partial _{{\textbf {q}}}\,{\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}})\cdot \Delta {\textbf {q}} \nonumber \\&\qquad +\partial _{\dot{{\textbf {q}}}}\,{\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}})\cdot \Delta \dot{{\textbf {q}}} + \partial _{\ddot{{\textbf {q}}}}\,{\textbf {g}}({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}})\cdot \Delta \ddot{{\textbf {q}}}\,, \end{aligned}$$

(C10)

where $\partial _{{\textbf {q}}}(\cdot )$, $\partial _{\dot{{\textbf {q}}}}(\cdot )$, and $\partial _{\ddot{{\textbf {q}}}}(\cdot )$ denote partial derivatives with respect to ${\textbf {q}}$, $\dot{{\textbf {q}}}$, and $\ddot{{\textbf {q}}}$, respectively. Introducing (C9) into (C11) yields:

$$\begin{aligned}&{\textbf {g}}_s({\textbf {q}}+\Delta {\textbf {q}}) \approx {\textbf {g}}_s({\textbf {q}}) + \nonumber \\&\quad \partial _{{\textbf {q}}}\,\left[ \int _0^S \left( \, {\textbf {B}} ({\textbf {q}})^T \, \varvec{\sigma }_h \,-\, {\textbf {N}}^T \, {\textbf {f}}^{\text {ext}} \,\right) \, \textrm{d}s\right] \cdot \Delta {\textbf {q}} \approx \nonumber \\&\quad {\textbf {g}}_s({\textbf {q}})+\left[ \int _0^S \partial _{{\textbf {q}}}\left( {\textbf {B}} ({\textbf {q}})^T \, \varvec{\sigma }_h - {\textbf {N}}^T \, {\textbf {f}}^{\text {ext}}\right) \, \textrm{d}s\right] \cdot \Delta {\textbf {q}}\,, \end{aligned}$$

(C11)

and (C8) into (C11):

$$\begin{aligned}&{\textbf {g}}_d({\textbf {q}}+\Delta {\textbf {q}},\dot{{\textbf {q}}}+\Delta \dot{{\textbf {q}}},\ddot{{\textbf {q}}}+\Delta \ddot{{\textbf {q}}}) \approx {\textbf {g}}_d({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}}) + \nonumber \\&\qquad \partial _{{\textbf {q}}}\,\left[ ...\right] \cdot \Delta {\textbf {q}} + \partial _{\dot{{\textbf {q}}}}\,\left[ ...\right] \cdot \Delta \dot{{\textbf {q}}} + \partial _{\ddot{{\textbf {q}}}}\,\left[ ...\right] \cdot \Delta \ddot{{\textbf {q}}} \nonumber \\&\quad \approx {\textbf {g}}_d({\textbf {q}},\dot{{\textbf {q}}},\ddot{{\textbf {q}}}) + \left[ \int _0^S\partial _{{\textbf {q}}}\left( ...\right) \textrm{d}s\right] \cdot \Delta {\textbf {q}} \nonumber \\&\quad +\left[ \int _0^S\partial _{\dot{{\textbf {q}}}}\left( ...\right) \textrm{d}s\right] \cdot \Delta \dot{{\textbf {q}}} + \left[ \int _0^S\partial _{\ddot{{\textbf {q}}}}\left( ...\right) \textrm{d}s\right] \cdot \Delta \ddot{{\textbf {q}}}\,. \end{aligned}$$

(C12)

In tensor notation, derivatives with respect to a vector can be expressed as³$\partial _{{\textbf {q}}}(\cdot ) = \partial _{q^i}(\cdot )\otimes {\textbf {G}}_i$, where ${\textbf {q}}=(q^1,q^2,...,q^{3m_\text {dof}})^T$, and $\{{\textbf {G}}_1,...,{\textbf {G}}_{3m_\text {dof}}\}$ is an orthonormal basis for $\mathbb {R}^{3m_\text {dof}}$. We consider the three following tangent matrices resulting from (C12) and (C13):

$$\begin{aligned} {\textbf {K}}_S&= \partial _{q^i}({\textbf {B}}({\textbf {q}})^T{\textbf {a}})\vert _{{\textbf {a}}=\varvec{\sigma }_h}\otimes {\textbf {G}}_i \nonumber \\&\quad + {\textbf {B}}({\textbf {q}})^T\partial _{q^i}(\varvec{\sigma }_h)\otimes {\textbf {G}}_i\,, \nonumber \\ {\textbf {K}}_F&= \partial _{q^i}({\textbf {N}}^T{\textbf {f}}^{\text {ext}})\otimes {\textbf {G}}_i + \partial _{\dot{q}^i}({\textbf {N}}^T{\textbf {f}}^{\text {ext}})\otimes {\textbf {G}}_i \nonumber \\&\quad + \partial _{\ddot{q}^i}({\textbf {N}}^T{\textbf {f}}^{\text {ext}})\otimes {\textbf {G}}_i\,, \nonumber \\ {\textbf {K}}_{M}&=\partial _{q^i}({\textbf {M}} ({\textbf {q}})\,\nabla _{\dot{{\textbf {q}}}}\dot{{\textbf {q}}})\otimes {\textbf {G}}_i \nonumber \\&\quad + \partial _{\dot{q}^i}({\textbf {M}} ({\textbf {q}})\,\nabla _{\dot{{\textbf {q}}}}\dot{{\textbf {q}}})\otimes {\textbf {G}}_i \nonumber \\&\quad + \partial _{\ddot{q}^i}({\textbf {M}} ({\textbf {q}})\,\nabla _{\dot{{\textbf {q}}}}\dot{{\textbf {q}}})\otimes {\textbf {G}}_i\,, \end{aligned}$$

(C13)

where ${\textbf {a}}\in \mathbb {R}^{3m_\text {dof}}$ is a constant vector, ${\textbf {K}}_S$, ${\textbf {K}}_F$ and ${\textbf {K}}_M$ are the so-called static tangent stiffness matrix, force tangent stiffness matrix and mass tangent stiffness matrix, respectively.

1.2 C.2 Static tangent stiffness matrix

The static tangent stiffness matrix ${\textbf {K}}_S$ consists of two terms: the first term resulting from the derivative of the operator ${\textbf {B}}({\textbf {q}})$ (geometric stiffness), and the second arising from the derivative of vector $\varvec{\sigma }_h$ (elastic stiffness), i.e., ${\textbf {K}}_S={\textbf {K}}_S^g+{\textbf {K}}_S^e$. We first recall the discrete stress and strain measures from Sect. 3.1, which are:

$$\begin{aligned} \begin{aligned} \varvec{\sigma }_h = \left\{ \begin{array}{ll} {\textbf {n}}_h \\ {\textbf {m}}_h \end{array}\right\}&= \left[ \begin{array}{ll} EA\,\varvec{\epsilon }_h \\[0.3em] EI\,\varvec{\kappa }_h \end{array}\right] \\&= \underbrace{\left[ \begin{array}{ll} EA\,{\textbf {I}}_3 &{} {\textbf {0}}_3 \\ {\textbf {0}}_3 &{} EI\,{\textbf {I}}_3 \end{array}\right] }_{{\textbf {C}}(E,A,I)}\,\left\{ \begin{array}{ll} \varvec{\epsilon }_h \\ \varvec{\kappa }_h \end{array}\right\} . \end{aligned} \end{aligned}$$

(C14)

Considering the discrete linearized strain operator ${\textbf {B}}({\textbf {q}}) = \begin{bmatrix} {\textbf {B}}_{11} &{} {\textbf {0}} \\ {\textbf {B}}_{21} &{} {\textbf {B}}_{22} \end{bmatrix}$ in (16), we compute the derivative of the sub-operators ${\textbf {B}}_{11}$, ${\textbf {B}}_{21}$, and ${\textbf {B}}_{22}$, as follows:

$$\begin{aligned} \partial _{q^i}{} & {} \left[ {\textbf {B}}_{11}^T({\textbf {q}}){\textbf {a}}\right] \vert _{{\textbf {a}}={\textbf {n}}_h} \nonumber \\{} & {} = \partial _{q^i}\left. \left[ {\textbf {N}}^{\prime T}\left( {\textbf {I}}-\left| \varvec{\varphi }_h^\prime \right| ^{-1} {\textbf {P}}_{{\textbf {d}}}\right) {\textbf {a}}\right] \right| _{{\textbf {a}}={\textbf {n}}_h} \nonumber \\{} & {} = {\textbf {N}}^{\prime T}{\textbf {A}}_1\,\partial _{q^i}(\varvec{\varphi }_h^\prime ), \nonumber \\ \partial _{q^i}{} & {} \left[ {\textbf {B}}_{21}^T({\textbf {q}}){\textbf {a}}\right] \vert _{{\textbf {a}}={\textbf {m}}_h} \nonumber \\{} & {} = \partial _{q^i}\left. \left[ {\textbf {N}}^{\prime T}\left( \left| \varvec{\varphi }_h^\prime \right| ^{-2}{\textbf {H}}_{{\textbf {d}}}[\varvec{\varphi }_h^{\prime \prime }]_{\times } \right) {\textbf {a}}\right] \right| _{{\textbf {a}}={\textbf {m}}_h} \nonumber \\{} & {} = {\textbf {N}}^{\prime T}{\textbf {A}}_2\,\partial _{q^i}(\varvec{\varphi }_h^\prime ) + {\textbf {N}}^{\prime T}{\textbf {A}}_3\,\partial _{q^i}(\varvec{\varphi }_h^{\prime \prime }), \nonumber \\ \partial _{q^i}{} & {} \left[ {\textbf {B}}_{22}^T({\textbf {q}}){\textbf {a}}\right] \vert _{{\textbf {a}}={\textbf {m}}_h} \nonumber \\{} & {} = \partial _{q^i}\left. \left[ {\textbf {N}}^{\prime \prime T}\left( \left| \varvec{\varphi }_h^\prime \right| ^{-1}[{\textbf {d}}_h]_{\times } \right) {\textbf {a}}\right] \right| _{{\textbf {a}}={\textbf {m}}_h} \nonumber \\{} & {} = {\textbf {N}}^{\prime \prime T}{\textbf {A}}_4\,\partial _{q^i}(\varvec{\varphi }_h^{\prime }), \end{aligned}$$

(C15)

where:

$$\begin{aligned} {\textbf {A}}_1{} & {} = \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^{2}}\left[ 2{\textbf {d}}_h\odot {\textbf {n}}_h - 3({\textbf {d}}_h\cdot {\textbf {n}}_h){\textbf {d}}_h\otimes {\textbf {d}}_h \right. \nonumber \\{} & {} \quad \left. + \left( {\textbf {d}}_h \cdot {\textbf {n}}_h \right) \, {\textbf {I}} \right] , \nonumber \\ {\textbf {A}}_2{} & {} = -\frac{2}{\left| \varvec{\varphi }_h^\prime \right| ^{3}}\left[ 2\left( [\varvec{\varphi }_h^{\prime \prime }]_{\times }{\textbf {m}}_h\odot {\textbf {d}}_h\right) \right. \nonumber \\{} & {} \quad \left. + ([\varvec{\varphi }_h^{\prime \prime }]_{\times }{\textbf {m}}_h\cdot {\textbf {d}}_h)(2{\textbf {H}}_{{\textbf {d}}}-{\textbf {I}})\right] , \nonumber \\ {\textbf {A}}_3{} & {} = -\frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^{2}}{\textbf {H}}_{{\textbf {d}}}\,[{\textbf {m}}_h]_{\times }, \nonumber \\ {\textbf {A}}_4{} & {} = \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^{2}}[{\textbf {m}}_h]_{\times }\,{\textbf {H}}_{{\textbf {d}}}. \end{aligned}$$

(C16)

Applying the spatial approximation of $\varvec{\varphi }(s,t)$ introduced in (13), we obtain:

$$\begin{aligned} \begin{aligned} \partial _{q^i}(\varvec{\varphi }_h^{\prime })&= \partial _{q^i}\left( \frac{d}{ds}({\textbf {N}}{\textbf {q}}) \right) \\&= {\textbf {N}}^{\prime } \partial _{q^i}\left( {\textbf {q}}\right) = {\textbf {N}}^{\prime }{\textbf {G}}_i, \\ \partial _{q^i}(\varvec{\varphi }_h^{\prime \prime })&= {\textbf {N}}^{\prime \prime }{\textbf {G}}_i. \end{aligned} \end{aligned}$$

(C17)

Introducing (C18) into (C16), we obtain the first counterpart of ${\textbf {K}}_S$, ${\textbf {K}}_S^g$, that is:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_S^g&= \partial _{q^i}({\textbf {B}}({\textbf {q}})^T{\textbf {a}})\vert _{{\textbf {a}}=\varvec{\sigma }_h}\otimes {\textbf {G}}_i \\&= \left[ {\textbf {N}}^{\prime T}({\textbf {A}}_1+{\textbf {A}}_2){\textbf {N}}^{\prime } + {\textbf {N}}^{\prime T}{\textbf {A}}_3\,{\textbf {N}}^{\prime \prime } \right. \\&\qquad \left. + {\textbf {N}}^{\prime \prime T}{\textbf {A}}_4{\textbf {N}}^{\prime }\right] \, \overbrace{{\textbf {G}}_i \otimes {\textbf {G}}_i}^{{\textbf {I}}_{3m_\text {dof}}} \\&= \begin{bmatrix} {\textbf {N}}^{\prime T}&{\textbf {N}}^{\prime \prime T} \end{bmatrix}\,\begin{bmatrix} {\textbf {A}}_1+{\textbf {A}}_2 &{} {\textbf {A}}_3 \\ {\textbf {A}}_4 &{} {\textbf {0}}_3\end{bmatrix}\,\begin{bmatrix} {\textbf {N}}^{\prime } \\ {\textbf {N}}^{\prime \prime } \end{bmatrix}, \end{aligned} \end{aligned}$$

(C18)

where ${\textbf {I}}_{3m_\text {dof}}$ is the identity matrix of dimension $3m_\text {dof}\times 3m_\text {dof}$.

The second counterpart, ${\textbf {K}}_S^e$, employing (C15), can be computed as follows:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_S^e&= {\textbf {B}}({\textbf {q}})^T\partial _{{\textbf {q}}}(\varvec{\sigma }_h) \\&= {\textbf {B}}({\textbf {q}})^T\partial _{{\textbf {q}}}\left( {\textbf {C}}\left\{ \begin{array}{ll} \varvec{\epsilon }_h \\ \varvec{\kappa }_h \end{array}\right\} \right) . \end{aligned} \end{aligned}$$

(C19)

We recall $(\varvec{\delta \epsilon },\delta \varvec{\kappa })^T=\mathcal {B}(\varvec{\varphi }^{\prime },\varvec{\varphi }^{\prime \prime })\, \delta \varvec{\varphi }$. Thus, the derivative $\partial _{{\textbf {q}}}(\varvec{\epsilon }_h,\varvec{\kappa }_h)^T$, directly related to the discretized strain operator ${\textbf {B}}({\textbf {q}})$, is:

$$\begin{aligned} \begin{aligned} \partial _{{\textbf {q}}}\left( \left\{ \begin{array}{ll} \varvec{\epsilon }_h \\ \varvec{\kappa }_h \end{array}\right\} \right) ={\textbf {B}}({\textbf {q}})\,\underbrace{\partial _{{\textbf {q}}}({\textbf {q}})}_{{\textbf {I}}} . \end{aligned} \end{aligned}$$

(C20)

${\textbf {K}}_S^e$ then takes the following form:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_S^e&={\textbf {B}}({\textbf {q}})^T{\textbf {C}}\,{\textbf {B}}({\textbf {q}}), \end{aligned} \end{aligned}$$

(C21)

where the constitutive matrix ${\textbf {C}}$ is assumed to be constant. The static tangent stiffness matrix ${\textbf {K}}_S$ is then:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_S&= {\textbf {K}}_S^g + {\textbf {K}}_S^e = \\&\begin{bmatrix} {\textbf {N}}^{\prime T}&{\textbf {N}}^{\prime \prime T} \end{bmatrix}\,\begin{bmatrix} {\textbf {A}}_1+{\textbf {A}}_2 &{} {\textbf {A}}_3 \\ {\textbf {A}}_4 &{} {\textbf {0}}_3\end{bmatrix}\,\begin{bmatrix} {\textbf {N}}^{\prime } \\ {\textbf {N}}^{\prime \prime } \end{bmatrix} + {\textbf {B}}({\textbf {q}})^T{\textbf {C}}\,{\textbf {B}}({\textbf {q}}). \end{aligned} \end{aligned}$$

(C22)

1.3 C.3 Force tangent stiffness matrix

We now derive the force tangent stiffness matrix ${\textbf {K}}_F$ that is required for our computations in this work. We consider two different types of forces: (i) follower forces for two-dimensional study cases, and (ii) forces induced by a surrounding flow described in Appendix B.

1.3.1 C.3.1 Follower forces in two-dimensional studies

In two-dimensional cases, the motion of the rod is confined in a single plane at all times (see Fig. 27). Thus, we can define a lumped follower force as follows:

$$\begin{aligned} {\textbf {f}}_f = f_0\,({\textbf {E}}_2\times {\textbf {d}}_h), \end{aligned}$$

(C23)

where $f_0$ is the magnitude of the force and is assumed to be constant, and ${\textbf {E}}_2$ is the unit normal to the plane of motion. According to (C14), the tangent stiffness matrix ${\textbf {K}}_F$ resulting from ${\textbf {f}}_f$ is:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_F&= \partial _{q^i}\left( {\textbf {N}}^T f_0\,({\textbf {E}}_2\times {\textbf {d}}_h)\right) \otimes {\textbf {G}}_i \\&=\frac{f_0}{\left| \varvec{\varphi }_h^\prime \right| }{\textbf {N}}^T \, [{\textbf {E}}_2]_{\times }\,{\textbf {P}}_{{\textbf {d}}}{\textbf {N}}^\prime \,\underbrace{{\textbf {G}}_i\otimes {\textbf {G}}_i}_{{\textbf {I}}_{3m_\text {dof}}} \\&=\frac{f_0}{\left| \varvec{\varphi }_h^\prime \right| }{\textbf {N}}^T \, [{\textbf {E}}_2]_{\times }\,{\textbf {P}}_{{\textbf {d}}}{\textbf {N}}^\prime . \end{aligned} \end{aligned}$$

(C24)

Fig. 27

Schematic for a two-dimensional follower force

We note that since the configuration space of Kirchhoff rods employed in this work requires only the specification of $\varvec{\varphi }(s,t)$ and the director ${\textbf {d}}$, a modeling of three-dimensional follower forces is generally not possible.

1.3.2 C.3.2 Force induced by a surrounding flow

As mentioned in Appendix B, we consider a simplified model for a force induced by a surrounding flow, that consists of three counterparts: an added mass force, a tangential drag force, and a normal drag force. In this subsection, we derive the tangent stiffness matrix corresponding to each counterpart. To simplify the involved algebra for this derivation, we consider the resulting expression of the semi-discrete formulation after employing the implicit time integration scheme described in Sect. 3.3. We note that the derivation of the tangent stiffness matrices can be alternatively performed with the equations of motion and/or the semi-discrete equations, and thus is independent of the numerical scheme applied.

The chosen implicit scheme in this work (see also Sect. 3.3) evaluates ${\textbf {g}}_d$ at the time instant $t_{n+\frac{1}{2}}$, and thus also the force term, i.e. $({\textbf {N}}^T{\textbf {f}}^{\text {ext}})_{n+\frac{1}{2}}$. The term associated with the force is then:

$$\begin{aligned} \int _0^{S}&\left[ ({\textbf {N}}^T{\textbf {F}}_f)_{n+\frac{1}{2}}\right] \textrm{d}s = C_1\int _0^{S}\left[ ({\textbf {N}}^T{\textbf {P}}_{{\textbf {d}}}{\textbf {a}}_{h})_{n+\frac{1}{2}}\right] \textrm{d}s \\&+ C_2\int _0^{S}\left[ ({\textbf {N}}^T\left| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_{h}\right| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_{h})_{n+\frac{1}{2}}\right] \textrm{d}s \\&+ C_3\int _0^{S}\left[ ({\textbf {N}}^T\left| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_{h}\right| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_{h})_{n+\frac{1}{2}}\right] \textrm{d}s\,, \end{aligned}$$

where ${\textbf {a}}_{h}={\textbf {a}}_{\infty }-\ddot{\varvec{\varphi }}_h$ is the discrete relative acceleration of the surrounding flow, and ${\textbf {V}}_{h}={\textbf {V}}_{\infty }-\dot{\varvec{\varphi }}_h$ the discrete relative velocity.

The added mass force, denoted as ${\textbf {F}}^{am}$, is a function of the free-stream, rod acceleration, and rod configuration. Employing the midpoint rule and a standard Taylor expansion up to first order, we obtain the following approximation for ${\textbf {F}}^{am}$ at the time instance $t_{n+\frac{1}{2}}$:

$$\begin{aligned} \begin{aligned}&{\textbf {F}}^{am}_{n+\frac{1}{2}}({\textbf {q}}_{n+1}+\Delta {\textbf {q}},\dot{{\textbf {q}}}_{n+1}+\Delta \dot{{\textbf {q}}}) \\&\qquad \approx {\textbf {F}}^{am}_{n+\frac{1}{2}}({\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n+1}) \\&\qquad + \partial _{{\textbf {q}}_{n+1}}{\textbf {F}}^{am}_{n+\frac{1}{2}}({\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n+1})\cdot \Delta {\textbf {q}} \\&\qquad + \partial _{\dot{{\textbf {q}}}_{n+1}}{\textbf {F}}^{am}_{n+\frac{1}{2}}({\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n+1})\cdot \Delta \dot{{\textbf {q}}}, \end{aligned} \end{aligned}$$

(C25)

where:

$$\begin{aligned}&{\textbf {F}}^{am}_{n+\frac{1}{2}}({\textbf {q}}_{n},{\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n},\dot{{\textbf {q}}}_{n+1}) \nonumber \\&\quad = C_1\int _0^S{\textbf {N}}^T\left\{ \frac{1}{2}\left[ ({\textbf {P}}_{{\textbf {d}}}{\textbf {a}}_{\infty })_n+({\textbf {P}}_{{\textbf {d}}}{\textbf {a}}_{\infty })_{n+1}\right] \right\} \textrm{d}s \nonumber \\&\quad - C_1\int _0^S{\textbf {N}}^T\left\{ \frac{1}{\Delta t}\left[ ({\textbf {P}}_{{\textbf {d}}}\dot{\varvec{\varphi }}_h)_{n+1}-({\textbf {P}}_{{\textbf {d}}}\dot{\varvec{\varphi }}_h)_{n}\right] \right\} \textrm{d}s \,. \end{aligned}$$

(C26)

Considering (C26) and (C27), the corresponding tangent stiffness matrix per unit of length is:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_F^{am} = C_1{\textbf {N}}^T&\left[ \frac{1}{2}\underbrace{\partial _{q_{n+1}^i}({\textbf {P}}_{{\textbf {d}}}{\textbf {a}}_{\infty })_{n+1}}_{{\textbf {f}}_1} \right. \\&\quad - \frac{1}{\Delta t}\underbrace{\partial _{q_{n+1}^i}({\textbf {P}}_{{\textbf {d}}}\dot{\varvec{\varphi }}_h)_{n+1}}_{{\textbf {f}}_2} \\&\quad \left. - \frac{1}{\Delta t}\underbrace{\partial _{\dot{q}_{n+1}^i}({\textbf {P}}_{{\textbf {d}}}\dot{\varvec{\varphi }}_h)_{n+1}}_{{\textbf {f}}_3}\right] \otimes {\textbf {G}}_i . \end{aligned} \end{aligned}$$

(C27)

We assume that both the magnitude and direction of the force vary with the altitude (or depth for marine applications). For the sake of simplicity, we consider ${\textbf {a}}_{\infty }$ as a function of the vertical coordinate z and time t, i.e., ${\textbf {a}}_{\infty }(z,t)$, where z depends on the current configuration of the rod, that is $z=\varvec{\varphi }_h\cdot {\textbf {E}}_3$. We then obtain:

$$\begin{aligned} {\textbf {f}}_1&= \left\{ \underbrace{-\frac{1}{\left| \varvec{\varphi }_h^\prime \right| }\left[ ({\textbf {a}}_{\infty }\cdot {\textbf {d}}_h){\textbf {H}}_{{\textbf {d}}}+ {\textbf {d}}_h\otimes {\textbf {a}}_{\infty }\right] _{n+1}}_{{\textbf {K}}_1}{\textbf {N}}^{\prime } \right. \nonumber \\&\quad \left. + \underbrace{\left[ {\textbf {P}}_{{\textbf {d}}}(\partial _z{\textbf {a}}_{\infty }\otimes {\textbf {E}}_3) \right] _{n+1}}_{{\textbf {K}}_2}{\textbf {N}}\right\} \,{\textbf {G}}_i \,, \nonumber \\ {\textbf {f}}_2&= \underbrace{-\frac{1}{\left| \varvec{\varphi }_h^\prime \right| }\left[ (\dot{\varvec{\varphi }}_h\cdot {\textbf {d}}_h){\textbf {H}}_{{\textbf {d}}}+ {\textbf {d}}_h\otimes \dot{\varvec{\varphi }}_h \right] _{n+1}}_{{\textbf {K}}_3}{\textbf {N}}^{\prime }\,{\textbf {G}}_i \,, \nonumber \\ {\textbf {f}}_3&= \underbrace{\left[ {\textbf {P}}_{{\textbf {d}}}{\textbf {N}} \right] _{n+1}}_{{\textbf {K}}_4}\,{\textbf {G}}_i \,. \end{aligned}$$

(C28)

Introducing (C29) into (C28), and recalling that ${\textbf {G}}_i\otimes {\textbf {G}}_i = {\textbf {I}}_{3m_\text {dof}}$ and $\Delta \dot{{\textbf {q}}}=\frac{2}{h}\Delta {\textbf {q}}$, we obtain:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_F^{am}&= C_1{\textbf {N}}^T \left[ \frac{1}{2}{\textbf {K}}_1{\textbf {N}}^{\prime } + \frac{1}{2}{\textbf {K}}_2{\textbf {N}} \right. \\&\left. - \frac{1}{\Delta t}{\textbf {K}}_3{\textbf {N}}^{\prime } - \frac{2}{\Delta t^2}{\textbf {K}}_4{\textbf {N}}\right] . \end{aligned} \end{aligned}$$

(C29)

The normal drag force, denoted as ${\textbf {F}}^{cn}$, is a function of the free-stream, rod velocity, and rod configuration. The normal drag force of the surrounding flow, evaluated at $t_{n+\frac{1}{2}}$, is:

$$\begin{aligned} \begin{aligned}&{\textbf {F}}^{cn}_{n+\frac{1}{2}}({\textbf {q}}_{n},{\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n},\dot{{\textbf {q}}}_{n+1}) \\&\quad = C_2\int _0^S{\textbf {N}}^T\left\{ \frac{1}{2}\left[ (\left| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h\right| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h)_n \right. \right. \\&\qquad \qquad \qquad \quad \left. \left. +(\left| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h\right| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h)_{n+1}\right] \right\} \textrm{d}s . \nonumber \end{aligned} \end{aligned}$$

(C30)

Analogously to the tangent stiffness matrix corresponding to the added mass force derived above, the tangent stiffness matrix per unit of length associated with the normal drag force is:

$$\begin{aligned} \begin{aligned} {\textbf {K}}_F^{cn}&= C_2{\textbf {N}}^T\left[ \frac{1}{2}{\textbf {K}}_1{\textbf {N}} - \frac{1}{2}{\textbf {K}}_2{\textbf {N}}^{\prime } - \frac{1}{\Delta t}{\textbf {K}}_3{\textbf {N}}\right] , \end{aligned} \end{aligned}$$

(C31)

where:

$$\begin{aligned} {\textbf {K}}_1&= \left| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h\right| (2{\textbf {I}}-{\textbf {P}}_{\text {u}}){\textbf {P}}_{{\textbf {d}}}(\partial _z{\textbf {V}}_{\infty }\otimes {\textbf {E}}_3)\,, \\ {\textbf {K}}_2&= \frac{\left| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h\right| }{\left| \varvec{\varphi }_h^\prime \right| }(2{\textbf {I}}-{\textbf {P}}_{\text {u}})\left[ ({\textbf {V}}_h\cdot {\textbf {d}}_h){\textbf {H}}_{{\textbf {d}}}+ {\textbf {d}}_h\otimes {\textbf {V}}_h\right] \,, \\ {\textbf {K}}_3&= \left| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h\right| (2{\textbf {I}}-{\textbf {P}}_{\text {u}}){\textbf {P}}_{{\textbf {d}}}, \\ {\textbf {P}}_{\text {u}}&= {\textbf {I}} - \frac{1}{\left| {\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h\right| ^2}({\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h)\otimes ({\textbf {P}}_{{\textbf {d}}}{\textbf {V}}_h)\,. \end{aligned}$$

The tangential drag force, denoted as ${\textbf {F}}^{ct}$, is also a function of the free-stream, rod velocity, and rod configuration. The tangential drag force of the surrounding flow, evaluated at $t_{n+\frac{1}{2}}$, is:

$$\begin{aligned}&{\textbf {F}}^{ct}_{n+\frac{1}{2}}({\textbf {q}}_{n},{\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n},\dot{{\textbf {q}}}_{n+1}) \\&\quad = C_3\int _0^S{\textbf {N}}^T\left\{ \frac{1}{2}\left[ \left( \left| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right) _n \right. \right. \\&\quad + \left. \left. (\left| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h)_{n+1}\right] \right\} \textrm{d}s\,. \end{aligned}$$

Analogously, the tangent stiffness matrix per unit of length associated with the tangential drag force is:

$$\begin{aligned} {\textbf {K}}_F^{ct}&= C_3{\textbf {N}}^T\left[ \frac{1}{2}{\textbf {K}}_1{\textbf {N}} + \frac{1}{2}{\textbf {K}}_2{\textbf {N}}^{\prime } - \frac{1}{\Delta t}{\textbf {K}}_3{\textbf {N}}\right] \,, \end{aligned}$$

where:

$$\begin{aligned} {\textbf {K}}_1&= \left| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right| \left[ 2{\textbf {I}}-{\textbf {P}}_{\text {v}}\right] ({\textbf {d}}_h\otimes {\textbf {d}}_h)(\partial _z{\textbf {V}}_{\infty }\otimes {\textbf {E}}_3)\,, \\ {\textbf {K}}_2&= \frac{\left| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right| }{\left| \varvec{\varphi }_h^\prime \right| }\left( 2{\textbf {I}}-{\textbf {P}}_{\text {v}}\right) \left[ ({\textbf {V}}_h\cdot {\textbf {d}}_h){\textbf {H}}_{{\textbf {d}}}+ {\textbf {d}}_h\otimes {\textbf {V}}_h\right] \,, \\ {\textbf {K}}_3&= \left| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right| (2{\textbf {I}}-{\textbf {P}}_{\text {v}})\left[ {\textbf {d}}_h\otimes {\textbf {d}}_h\right] \,, \\ {\textbf {P}}_{\text {v}}&= {\textbf {I}} - \frac{1}{\left| ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right| ^2}\left[ ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right] \otimes \left[ ({\textbf {d}}_h\otimes {\textbf {d}}_h){\textbf {V}}_h\right] \,. \end{aligned}$$

1.4 C.4 Mass tangent stiffness matrix

Lastly, we derive the mass tangent stiffness matrix, ${\textbf {K}}_M$, corresponding to the inertial contribution ${\textbf {M}}({\textbf {q}})\nabla _{\dot{{\textbf {q}}}}\dot{{\textbf {q}}}$ of the mass operator. Similarly to the tangent stiffness matrix corresponding to the force induced by a surrounding flow in the previous subsection, we consider the resulting expression of the semi-discrete formulation after employing the implicit time integration scheme described in Sect. 3.3. Considering (21) together with the midpoint rule formulas, we obtain the following approximation for the inertia term of the mass operator, evaluated at $t=n+\frac{1}{2}$, which is:

$$\begin{aligned}&{\textbf {f}}^{I}_{n+\frac{1}{2}}({\textbf {q}}_{n},{\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n},\dot{{\textbf {q}}}_{n+1}) \nonumber \\&\quad = \int _0^S\left\{ \frac{1}{\Delta t}\left[ \left( {\textbf {M}}\dot{{\textbf {q}}}\right) _{n+1} - \left( {\textbf {M}}\dot{{\textbf {q}}}\right) _n\right] \right. \nonumber \\&\quad + \frac{1}{2}{\textbf {N}}^{\prime T}\left[ 2\,I_{\rho } \left( \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^3}[{\textbf {P}}_{{\textbf {d}}}\odot (\dot{\varvec{\varphi }}_h^\prime \otimes {\textbf {d}}_h)]\dot{\varvec{\varphi }}_h^\prime \right) _n \right. \nonumber \\&\qquad \left. \left. + 2\,I_{\rho }\left( \frac{1}{\left| \varvec{\varphi }_h^\prime \right| ^3}[{\textbf {P}}_{{\textbf {d}}}\odot (\dot{\varvec{\varphi }}_h^\prime \otimes {\textbf {d}}_h)]\dot{\varvec{\varphi }}_h^\prime \right) _{n+1}\right] \right\} \textrm{d}s \,. \end{aligned}$$

(C32)

Employing the standard Taylor expansion up to first order for this term leads to:

$$\begin{aligned} \begin{aligned} {\textbf {f}}^{I}_{n+\frac{1}{2}}&({\textbf {q}}_{n+1}+\Delta {\textbf {q}},\dot{{\textbf {q}}}_{n+1}+\Delta \dot{{\textbf {q}}}) \\&\quad \approx {\textbf {f}}^{I}_{n+\frac{1}{2}}({\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n+1}) \\&\quad + \partial _{{\textbf {q}}_{n+1}}{\textbf {F}}^{I}_{n+\frac{1}{2}}({\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n+1})\cdot \Delta {\textbf {q}} \\&\quad + \partial _{\dot{{\textbf {q}}}_{n+1}}{\textbf {F}}^{I}_{n+\frac{1}{2}}({\textbf {q}}_{n+1},\dot{{\textbf {q}}}_{n+1})\cdot \Delta \dot{{\textbf {q}}} . \end{aligned} \end{aligned}$$

(C33)

Introducing (C32) into (C33), we obtain the following counterparts associated with the inertia contribution of the mass operator:

$$\begin{aligned}&{\textbf {K}}_1 = \partial _{q^i_{n+1}}({\textbf {M}}({\textbf {q}})\dot{{\textbf {q}}})_{n+1}\otimes {\textbf {G}}_i \\&\quad = -\frac{I_{\rho }}{\left| \varvec{\varphi }_h^\prime \right| ^3}{\textbf {N}}^{\prime T}\left\{ ({\textbf {d}}_h\cdot \dot{\varvec{\varphi }}_h^{\prime })(2{\textbf {H}}_{{\textbf {d}}}-{\textbf {I}}) \right. \\&\quad \left. \, + {\textbf {d}}_h\otimes \dot{\varvec{\varphi }}_h^{\prime } + 2\,\dot{\varvec{\varphi }}_h^{\prime }\otimes {\textbf {d}}_h \right\} {\textbf {N}}^{\prime } \,, \\&{\textbf {K}}_2 = \partial _{\dot{q}^i_{n+1}}({\textbf {M}}({\textbf {q}})\dot{{\textbf {q}}})_{n+1}\otimes {\textbf {G}}_i = {\textbf {N}}^T{\textbf {M}}({\textbf {q}}_{n+1}){\textbf {N}} \,, \\&{\textbf {K}}_3 = 2\,I_{\rho }{\textbf {N}}^{\prime T} \partial _{q^i_{n+1}}\left( \frac{[{\textbf {P}}_{{\textbf {d}}}\odot (\dot{\varvec{\varphi }}_h^\prime \otimes {\textbf {d}}_h)]\dot{\varvec{\varphi }}_h^\prime }{\left| \varvec{\varphi }_h^\prime \right| ^3}\right) _{n+1}\otimes {\textbf {G}}_i \\&\quad = \frac{-I_{\rho }\,{\textbf {N}}^{\prime \, T}}{\left| \varvec{\varphi }_h^\prime \right| ^4} \left\{ 2(\dot{\varvec{\varphi }}_h^\prime \cdot {\textbf {d}}_h) \left[ 4 \dot{\varvec{\varphi }}_h^\prime \odot {\textbf {d}}_h + (\dot{\varvec{\varphi }}_h^\prime \cdot {\textbf {d}}_h)[3{\textbf {H}}_{{\textbf {d}}}-2{\textbf {I}}]\right] \right. \\&\quad - \left. \dot{\varvec{\varphi }}_h^\prime \otimes \dot{\varvec{\varphi }}_h^\prime + (\dot{\varvec{\varphi }}_h^\prime \cdot \dot{\varvec{\varphi }}_h^\prime )[{\textbf {I}} - 2{\textbf {H}}_{{\textbf {d}}}] \right\} \,{\textbf {N}}^{\prime } \,, \\&{\textbf {K}}_4 = 2\,I_{\rho }{\textbf {N}}^{\prime T} \partial _{\dot{q}^i_{n+1}}\left( \frac{[{\textbf {P}}_{{\textbf {d}}}\odot (\dot{\varvec{\varphi }}_h^\prime \otimes {\textbf {d}}_h)]\dot{\varvec{\varphi }}_h^\prime }{\left| \varvec{\varphi }_h^\prime \right| ^3}\right) _{n+1}\otimes {\textbf {G}}_i \\&\quad = 2\frac{I_{\rho }}{\left| \varvec{\varphi }_h^\prime \right| ^3}{\textbf {N}}^{\prime T} \left\{ (\dot{\varvec{\varphi }}_h^\prime \cdot {\textbf {d}}_h){\textbf {H}}_{{\textbf {d}}}+ {\textbf {d}}_h\otimes \dot{\varvec{\varphi }}_h^\prime \right. \\&\qquad \qquad \qquad \left. + 2{\textbf {P}}_{{\textbf {d}}}\odot (\dot{\varvec{\varphi }}_h^\prime \otimes {\textbf {d}}_h) \right\} \,{\textbf {N}}^{\prime } \,. \end{aligned}$$

Recalling that $\Delta \dot{{\textbf {q}}}=\frac{2}{h}\Delta {\textbf {q}}$, we then obtain the tangent stiffness matrix per unit of length associated with the inertia term of the mass operator as follows:

$$\begin{aligned} {\textbf {K}}_M&=\left[ \frac{1}{\Delta t}{\textbf {K}}_1 + \frac{2}{\Delta t^2}{\textbf {K}}_2 + \frac{1}{2}{\textbf {K}}_3 + \frac{1}{\Delta t}{\textbf {K}}_4\right] \,. \end{aligned}$$