Simulation of viscoelastic Cosserat rods based on the geometrically exact dynamics of special Euclidean strands

We propose a method for the description and simulation of the nonlinear dynamics of slender structures modeled as Cosserat rods. It is based on interpreting the strains and the generalized velocities of the cross sections as basic variables and elements of the special Euclidean algebra. This perspective emerges naturally from the evolution equations for strands, that are one-dimensional submanifolds, of the special Euclidean group. The discretization of the corresponding equations for the three-dimensional motion of a Cosserat rod is performed, in space, by using a staggered grid. The time evolution is then approximated with a semi-implicit method. Within this approach we can easily include dissipative effects due to both the action of external forces and the presence of internal mechanical dissipation. The comparison with results obtained with different schemes shows the effectiveness of the proposed method, which is able to provide very good predictions of nonlinear dynamical effects and shows competitive computation times also as an energy-minimizing method to treat static problems.


Introduction
The modeling and simulation of beams is of great importance in the engineering practice to analyze the configurations and stress distributions of a wide variety of mechanical structures, with sizes ranging from those of pipelines and cables to those of microactuators. When these structures are sufficiently slender or very flexible, they can undergo large displacements even in the small-strain and linear-response regime, and geometric nonlinearities must be taken into account to capture their mechanical behavior.
Since the seminal work by Simo and Vu-Quoc [1], the number of publications and numerical methods related to geometrically-exact beam models has been growing significantly. Nevertheless, given the variety of applications and the different features pertaining to each method, no universal standard is available for an efficient simulation of such models. On the other hand, it has become clear that the theory of special Cosserat rods (as presented for instance by Antman [2]) provides the optimal mathematical framework to deal with slender structures, as it comprises all of the classical beam models as special cases.
In the literature, we can find approximation schemes based on a discrete mechanical analogue for the rod, such as those by Bertails et al. [3], Bergou et al. [4], Giusteri & Fried [5], Jung et al. [6], Lang, Linn & Arnold [7], and Linn [8]. A similar structure is shared by the finite-element approaches by Borri and Bottasso [9], Cao, Liu & Wang [10], and Spillmann & Teschner [11]. Several other methods are based upon discretizing the evolution equations for the continuum rod. Many represent the rod via the position and orientation of nodal cross sections. In this way, the computation of the strains and stresses associated with twist, bending, stretching, and shearing of the rod relies on interpolation procedures that introduce some important arbitrariness in the calculations [12]. In other cases, nonlinear shape functions are used to approximate the configuration of the rod, as done by Patil & Althoff [13] and Howcroft et al. [14].
To simplify the derivation and the structure of the rod equations, a fundamental step is to view not only the strains but also the generalized velocities of the rigid cross sections as elements of the Lie algebra associated with the special Euclidean group of rigid body motions. This perspective led Simo, Marsden & Krishnaprasad [32] and Hodges [33,34] to derive the intrinsic rod equations from the variations of a Hamiltonian functional expressed solely in terms of Lie algebraic quantities. Casting the equations in a linear space such as the algebra se(3) has several computational advantages in reference to interpolation strategies and the imposition of linear constraints.
Starting from the approach of Holm & Ivanov [35], in Section 2 we derive, in a rather straightforward way, evolution equations that correspond to the mixed formulation proposed by Hodges [33] with the addition of dissipative forces (both due to internal and external viscous phenomena) and of an equation that translates an important compatibility condition on the evolution of velocities and strains, necessary to close the system of partial differential equations. Positional and rotational degrees of freedom never appear in the equations, since they are merely recovered following the evolution of the rod placement from the initial configuration as driven by the generalized velocities. We then introduce a finite-difference scheme and discretize the evolution equations on a staggered grid so as to avoid shear-locking effects. In Section 3, we show the effectiveness of our method by applying it to the solution of both static and dynamic problems that involve viscoelastic rods featuring possibly curved relaxed shapes and anisotropic cross sections.
We believe that the method presented here provides a synthesis of many of the features towards which the recent literature on geometrically exact rods is converging. In particular, the theoretical setting enjoys a significant degree of mathematical transparency, the evolution takes place in a linear space with degrees of freedom represented in the most economic way, there are no limitations on the geometry of the cross sections and on the relaxed shapes, and the local nature of the representation allows for a straightforward application of external forces and the combination, by means of boundary conditions, of several structural elements. Other intersting aspects, relevant for specific applications, are mentioned in Section 4.
The target application that we have in mind is the study of the nonlinear dynamics of viscoelastic beams, but a strongly dissipative evolution can be used also as an alternative energyminimization method to retrieve static solutions. To assess the usefulness of our method, we implemented it in the Python language and compared its performance with a selection of published results, obtained with rather different schemes, and with results produced by an established commercial software. We find that, in spite of the simplicity of our formulation and of the discretization schemes that we have adopted, the method achieves very good results in solving both static and nonlinear dynamical problems, with competitive computational times.

The computational model
We propose a method for the description of the nonlinear dynamics of slender beams that is based on extensions of the SE (3)-strand equations described by Holm and Ivanov [35], with a suitable mechanical interpretation of stresses and momenta as dual to strains and velocities. Within this approach, we can easily include dissipative effects due to both the action of external forces and the presence of internal mechanical dissipation. Moreover, the relaxed shape of the rod can be arbitrarily prescribed.

Rod kinematics
We view a rod as a one-parameter family of rigid cross sections labelled by s ∈ [0, L], where L is a reference length. Each cross section is characterized by the position x of its center of mass and three orthonormal vectors d 1 , d 2 (in the cross-sectional plane) and d 3 (orthogonal to the cross section). The Euclidean transformation that translates the origin into x and rotates the Euclidean reference basis (e 1 , e 2 , e 3 ) into (d 1 , d 2 , d 3 ) is an element of the special Euclidean group SE (3) of rigid body motions. In this way, we can identify the rod at each time instant with an SE (3)-strand, that is a curve on SE (3) parametrized by s.
We take as time parameter t ∈ [0, +∞) and define the field P : [0, L] × [0, +∞) → Mat 4×3 (R) by P := (x, d 3 , d 1 , d 2 ) T , namely, each 4 × 3 matrix P(s) has in its rows the components of the vectors x, d 3 , d 1 , and d 2 at s. Then we define The twist density u 3 , curvatures u 1 , u 2 , the stretching σ 3 , and the shearing densities σ 1 and σ 2 are the six strains that describe the shape of the rod. We will identify the collection of them with the six-component strain vector U = (u 3 , u 1 , u 2 , σ 3 , σ 1 , σ 2 ). The natural (relaxed) shape of the rod is given by fixing preferred strainsŪ (s) for each value of s. The placement of the rod in the three-dimensional ambient space is given by the solution of the ODE where the prime denotes differentiation with respect to s, with the conditions P 0 (t) given at s = 0. The nonlinear relation between the strains and the placement can be seen from the definitions One of the advantages of our approach is to avoid using these relations in the simulation process.
If we now consider the motion of each cross section we find that, for any given s, the time derivative of P, denoted by a superimposed dot, is given bẏ where the spin-velocity matrix Ω has the very same structure of L, with and features three linear velocities v i and three spins w i (i = 1, 2, 3). We will identify the collection of them with the generalized velocity vector V = (w 3 , w 1 , w 2 , v 3 , v 1 , v 2 ). The definitions of velocities and spins in terms of the placement components read v i :=ẋ · d i , for i = 1, 2, 3, w 1 :=ḋ 3 · d 2 , w 2 :=ḋ 1 · d 3 , and w 3 :=ḋ 2 · d 1 .

Constitutive assumptions and Euler-Poincaré equations
We can now introduce the momentum and stress fields, P and Σ, as where the symmetric positive definite matrices M(s) ∈ Mat 6×6 (R) and A(s) ∈ Mat 6×6 (R) represent, respectively, the rigid-body inertia (linear density, determined by the geometry of the cross sections) and the elastic stiffnesses at each cross section. The elastic response of the rod is here modeled with a linear function of the difference between the current strains and the preferred ones.
For the test cases considered in what follows, the rigid cross sections with surface area A are assumed to have uniform mass density. We choose d 1 , d 2 , and d 3 aligned with the principal axes of inertia of each cross section and denote by I 1 , I 2 , and I 3 the corresponding second area moments. We further assume that the elastic response does not couple different strain components. Under these assumptions, the inertia and stiffness matrices take the form where E is the Young modulus of the material and G is the shear modulus, related to E by G = E/(2 + 2ν), that involves the Poisson ratio ν. We thus see that the field Σ represents torques and forces (tensions), while P is a linear density of angular and linear momenta. We observe that both U and V represent elements of the special Euclidean algebra se (3), while P and Σ are in the dual algebra se(3) * (isomorphic to se(3) in this finite-dimensional setting).
To derive the evolution equations for the rod in terms of the evolution of the pairs (V, U ) or (P, Σ) we start from the variational approach of Holm and Ivanov [35] specialized to the quadratic Lagrangian action This involves the total kinetic energy of the rod and the total elastic energy and, taking the first variation of S, we can apply Hamilton's principle and obtain the evolution equations for the conservative dynamics of an elastic rod.
It is now important to specify what type of variations are appropriate. In fact, the fields U and V as functions of (s, t) are associated with derivatives of the rod placement P, that is in one-to-one correspondence with a strand in SE (3). The latter is a time-dependent curve g : [0, L]×[0, +∞) → SE (3) and so we need to consider variations of U and V that are constrained to be consistent with the geometric nature of the rod descriptions.

Compatibility condition
The elements of the algebra se(3) are tangent vector to SE (3) at the identity. The strain and spin-velocity matrices at (s, t) are tangent vectors at g(s, t) ∈ SE (3) generated by motions along s or t, respectively, and, in a suitable matrix representation, are given by and By taking derivatives of the forgoing expressions we obtaiṅ From the difference of these equations and, considering the equality of cross derivatives, we arrive at the compatibility conditionL − Ω = LΩ − ΩL.

Adjoint operator and its dual
In the matrix representation of se(3) we can identify commutators with adjoint operators as In the six-component vector representation in which L ∼ U and Ω ∼ V we consistently define where X k denotes the k-th component of the vector X.
We also need to compute the dual operator, namely the adjoint-transpose operator ad T , that is defined in relation to a duality pairing which, in our case, is the Euclidean scalar product in R 6 . For all X, Y ∈ se(3) and Z ∈ se(3) * ∼ = se(3), we have , from which we can infer that, in this representation, ad T X = (ad X ) T .

Constrained variations
The Euler-Poincaré evolution equations for an SE (3)-strand with Lagrangian action S can be derived by considering the constrained variations δV =Ẋ + ad V X and δU = X + ad U X for an arbitrary test field X : In fact, if we consider a variation δg of g in the group SE (3) we have X ∼ g −1 δg and An analogous computation proves the expression for δU . By assuming that the variation field X vanishes at the ends of the domain of integration, the first constrained variation of the action S gives The stationarity condition δS, X = 0 for any test field X then implieṡ that encodes the local balances of linear and angular momentum for an elastic rod in the absence of external or dissipative forces. More general boundary conditions on X can be considered for practical purposes, but the evolution equations would remain the same, and only the boundary conditions satisfied by the solutions would change. We stress that the derivation of equation (6) does not rely upon the choice of a linear constitutive equation. In fact, it preserves its form if we simply define Σ as the derivative of the elastic energy density with respect to U .

Dynamic equations for a viscoelastic rod
Now that we have derived from a variational principle the conservative evolution equations (6), we are in a position of including additional force and torque densities of a possibly dissipative nature. The simplest choices of dissipative phenomena consist in an internal dissipation that depends linearly on the time derivative of the strains and an external viscous drag that depends linearly on the generalized velocity. With these choices, the evolution equations for the dynamics of a viscoelastic rod areṖ where the second relation translates the compatibility condition (4), namely the constraint imposed on the rigid-body motion of each cross section by the fact that they should collectively move as a continuous rod. The six-component vector F represents external force and torque densities acting on each cross section, while the (symmetric positive definite) matrices D ex and D in contain the damping coefficients associated with external drag and internal mechanical dissipation, respectively. We can substitute definitions (3) in (7) and (8) and obtaiṅ so that (9) and (10) constitute a system of PDEs in the unknown (P, Σ).

Discretization
We discretize the evolution equations using a standard semi-implicit time integration that linearizes the evolution operator. On the other hand, some care is in order when considering the spatial discretization. We used finite differences on a staggered grid that mirrors the character of our variables: momenta and velocities are assigned to nodal cross sections, while stresses and strains are viewed as segment quantities and collocated at the midpoint of each mesh interval ( Figure 2). We introduce a partition {0 = s 0 , s 1 , . . . , s N = L} of [0, L] and consider a space-time cell for s k ≤ s ≤ s k+1 and t n ≤ t ≤ t n+1 and denote by a subscript k or k, k + 1 nodal or segment quantities, respectively. Superscripts indicate the time instant at which the quantity is computed. Boundary conditions are imposed by adding accessory cells and nodes at the two ends of the rod. In this way, it is rather straightforward to drive or fix the motion of the rod ends or to set free-end conditions. Within this scheme, equation (9) features descretized quantities that live on the nodes s k of the partition, while (10) features descretized quantities that live on the intervals of the partition.
The discretization of equation (9) reads, for inner and free nodes with 0 < k ≤ N , and the discretization of equation (10) for inner segments (0 < k < N − 1) is with the substitutions U n k,k+1 = A −1 Σ n k,k+1 +Ū k,k+1 and V n k = M −1 P n k . The choice m = n describes a fully explicit scheme, while m = n + 1 leads to a semi-implicit scheme, that becomes fully implicit only in those cases in which the nonlinear terms are exactly vanishing (purely axial, shearing or twisting deformations; bending is excluded).
Different boundary conditions can be imposed but, for the following tests, we always need the same set of conditions. At one end of the rod we prescribe the motion through a givenP 0 (t), possibly vanishing, while Σ is free; on the other end we have a given stress Σ =Σ e (t) and P free. These conditions translate into where pedices (−1, 0) and (N, N + 1) denote the accessory segments that are added outside the physical rod to impose the boundary conditions.

Numerical results
To assess the effectiveness of our method we present a series of examples and some results about computational costs. In the small-displacement regime we can make comparisons with analytical results derived from the linearized equations of motion. On the other hand, in the nonlinear large-displacement regime we will test our numerical solutions against published results on some benchmark problems. We implemented the computational model in the Python language, exploiting the scientific computing libraries numpy and scipy and the just-in-time compilation features provided by numba. We tested our implementation on a Laptop with a 1,8 GHz Intel ® Core™ i5 processor and 8 GB of 1600 MHz DDR3 memory.

Small-displacement regime: cantilever
We compared the static solution for a cantilever, clamped at one end and subject to its own weight, obtained as a long-time limit of the dissipative dynamics for a 4 m long hollow cylinder with material parameters given in Table 1. The dissipation, useful to reach the static solution, is generated by the internal dissipation matrix D in = η in diag(1, 1, 1, 1, 1, 1), with η in = 10 −4 . The spatial domain is discretized uniformly with a variable number N of intervals. With a time step of 0.01 s, we obtain convergence to a static solution (identified by a kinetic energy below 10 −12 J) within 100 steps in all cases. The results, presented in Figure 3 (s)|. We found that this error estimate scales as h 2 , where h = L/N is the size of the mesh intervals. If we consider only the linear terms of the evolution equations, this result is consistent with our discretization that employs centered finite differences. The same scaling can be observed for the absolute error on the tip displacement (Figure 3(c,d)).
We then analyzed the vibration of a cantilever with the same physical parameters given above, initialized with a small curvature (0.01 m −1 ) and in the absence of gravity and dissipation. We varied the length but kept the number of discrete segments equal to 16. We compared the fundamental frequency of the tip displacement given by the theory, ν theory , with those computed from our solution by means of a Fast Fourier Transform algorithm, ν fft . The results, presented in Table 2, show a very good match.      et al. [14]. For further comparison of the computational performance, we solved the same problem with the commercial software Abaqus ® , using quadratic beam elements B32. The problem consists in a cantilever 45-degree bend subjected to a fixed load at the tip (the effect of weight is neglected). The relaxed configuration ( Figure 5, dot-dashed curves) spans a circular arc of 45 degrees in the xy-plane. The beam has a square cross section of side 1 m, radius of curvature of 100 m and total length L = 78.54 m. The Young modulus is E = 10 7 Pa and the Poisson ratio ν = 0. The rod is initially in the horizontal plane and the load is applied in the vertical direction. The equilibrium configuration is computed for two different values of the load (300 N and 600 N). Due to the nonlinear coupling between all of the strains produced by the curved geometry of the beam, all of the components of Σ are activated ( Figure 4) and contribute to determine the equilibrium configuration ( Figure 5). The good agreement between our method and those presented in the literature can be assessed by considering the tip displacement reported in Table 3, where, for an easier comparison of the computational efficiency, we consider a discretization with 81 nodes. As shown in Table 4, while the number of degrees of freedom we use is comparatively large (as typical of local discretization schemes) the computation remains very fast. We studied the convergence of the numerical approximation by computing solutions for different values of N , using up to N = 5120 segments. Similar to what was done for the cantilever, the static solution is achieved following a dissipative dynamics. In this case, we added an external dissipation matrix D ex = η ex diag (1, 1, 1, 1, 1, 1), with η ex = 10 −1 . With a time step of 10 s the solution converges within 100 steps for each value of N . The l ∞ error on the three-dimensional rod configuration is computed as the maximum distance between images of the same abscissa (mesh node) through the midline placement, namely sup s∈  (Figure 6). The difference from the quadratic scaling observed in the small-displacement regime can be attributed to the different weight of the nonlinear terms in the solution, that entails an approximation of order h 2 for terms that are quadratic in Σ.

Large-displacement regime: dynamic solution
The geometric nonlinearities that characterize the rod dynamics in the large-displacement regime have a strong influence on both the transient and steady motion behavior of slender structures. The l ∞ error on the midline placement scales linearly in h and so does the l ∞ error on the components of Σ, the actual degrees of freedom we solve for in our scheme. The difference in comparison to the small-displacement regime can be attributed to the different weight of the nonlinear terms in the solution, that entails an approximation of order h 2 for terms that are quadratic in Σ.
To provide a first check that our method is able to correctly capture such nonlinear effects we use as benchmark the data published by Howcroft et al. [14] about two tests. We consider a rather flat cantilever beam with material parameters given in Table 5. One end of the beam is clamped so that the tangent to the midline at s = 0 points always in the xy-plane. The motion of that end is driven as detailed below. The stress-free configuration features an intrinsic curvaturē  In the first test, we consider the profile of the beam in steady rotational motion around the vertical axis. This is obtained by imposing at the clamped end a non-vanishing component of the angular momentum along the vertical axis with various frequencies. The steady radial profiles obtained with N = 32 segments match very well both the benchmark data and the profile obtained by solving the static problem in a rotating frame, in which we impose the appropriate centrifugal forces (Figure 7).
In the second test, the clamped end oscillates vertically with a frequency of 4.5 Hz and a total amplitude of 0.04 m. After 12 s we superimpose to this oscillation a rotation about the vertical axis with a frequency that increases linearly for 8 s up to 8 Hz. During this ramp, the stiffening induced by the change in curvature causes an intrinsic vibrational frequency of the nonlinear system to cross 4.5 Hz, so that a resonance phenomenon can be observed. The time evolution computed with our method nicely captures oscillations and this transient phenomenon, matching very well the selected benchmark data (Figure 8).   [14] (solid lines) the radial profile of the midline for a beam (material parameters given in Table 5) with a curved relaxed configuration (dashed line) that is spinning around the vertical axis at different frequencies (1 Hz, 2 Hz, and 8 Hz). The symbols represent solutions obtained with N = 32 segments. We check the internal consistency of our method by superimposing the profile generated by a rotating beam (crosses) with the static solution in a rotating frame with centrifugal forces (squares).  [14] (thicker turquoise line). Up to time t = 12 s only the vertical oscillation at 4.5 Hz is active. At that time, the rotation is turned on with a frequency that increases linearly and reaches 8 Hz at t = 20 s. The tip oscillates always around the position obtained without the vertical oscillation (dashed red line) but the amplitude of the oscillation presents a resonance phenomenon when an intrinsic vibrational frequency of the nonlinear system, that varies due to the change in curvature, crosses 4.5 Hz. Small discrepancies between our data and the benchmark may be attributed to the residual presence in the latter of some oscillation mode different from the 4.5 Hz signal, probably due to the propagation of transient effects from the initialization of the test that, on the contrary, are completely gone in our data, as can be appreciated from the inset.

Computational efficiency
The semi-implicit computational scheme requires the solution of a linear system with a matrix that is updated at each time step. The system matrix is banded with 47 non-vanishing diagonals. The current implementation exploits such a structure in the system solution. The construction of the matrix is indeed the most expensive part of the algorithm and we achieved an optimal memory usage and a computation time that scales linearly with the number N of mesh intervals, proportional to the degrees of freedom (Figure 9).

Conclusions
We presented a method for the simulation of the dynamics of slender structures described within the framework of special Cosserat rods. The beam is viewed as a one-parameter family of rigid cross sections. For this reason, the Lie group SE (3) of rigid-body motions and the associated Lie algebra se(3) play a fundamental role in the description of the rod kinematics and dynamics. In fact, one can identify the placement of the rod in the three-dimensional ambient space as an SE (3)-strand, namely a curve in SE (3), and deduce the intrinsic evolution equations for an elastic rod from a variational principle, computing the Euler-Poincaré equations associated with an SE (3)-invariant Lagrangian action.
As is well known from the literature, the system Lagrangian only involves elements of the algebra se(3) that describe both the generalized velocity of each cross section and the generalized strains of the rod. One usually needs to take into account kinematic relations between the two that involve translational and rotational degrees of freedom, but this geometric setting produces additional compatibility equations, involving solely generalized velocities and strains, that encode the constraint imposed on the rigid-body motion of each cross section by the fact that they should collectively move as a continuous rod.
We extend the evolution equations to incorporate dissipative effects that can be of internal origin, depending on variations of the strains over time, and external, such as the viscous drag exerted by a surrounding fluid. These terms do not require additional degrees of freedom to be expressed and are thus perfectly compatible with the proposed intrinsic framework. In our presentation, we focused on linear constitutive prescription for both the elastic and the viscous response, but the inclusion of nonlinear laws is mathematically trivial, though it may be computationally more challenging. It should be noted that the linear elastic laws allow for a formulation of the equations that is simple and minimal, in the sense that we have twelve equations for twelve unknown fields, whereas the treatment of nonlinear constitutive laws may be more manageable in a mixed formulation involving eighteen equations, six of which would be the algebraic constitutive relations.
The viscoelastic evolution equations together with the consistency relations constitute a theoretically sound and practically flexible starting point for numerical approximation schemes. We have implemented a finite-difference approximation on a staggered grid in space and a semi-implicit time-stepping that, in spite of its simplicity, shows a very good computational performance in paradigmatic tests for both static and dynamic problems. Moreover, its good scalability makes it particularly attractive for the treatment of very large structures.
Within our general setting one can impose internal constraints such as unshearbility and inextensibility with a minimal effort, since they translate into linear constraints on the set of strains. Moreover, the fast reconstruction of the rod placement, that can be performed by applying the exponential map (see A) to the generalized velocity at each time step, allows for the inclusion of position-dependent external forces that may be of relevance for the simulation of contact and other external interactions.

Acknowledgments
This research has been supported by Tenaris (www.tenaris.com).

A Exponential map in the Special euclidean setting
To reconstruct the final placement P f of a cross section that, starting from a given P i , moves for a small time interval τ with constant spin-velocity matrix Ω can be achieved by means of the matrix exponential function, that maps the element τ Ω ∈ se(3) into an element of the matrix representation of the group SE (3), namely P f = exp(τ Ω)P i .
It is well known that the accurate computation of a generic matrix exponential may be challenging, but thanks to the peculiar structure of matrices that represent elements of se(3), we can easily find that exp(τ Ω) = Id + τ Ω + 1 − cos θ θ 2 (τ Ω) 2 + θ − sin θ θ 3 (τ Ω) 3 , where θ = τ 2 (w 2 1 + w 2 2 + w 2 3 ) is the norm of the spin component of τ Ω. With this analytical relation we can keep track of the rod placement as a post-processing of the solution to the evolution equations, that may also be important in calculating position-dependent forces acting on the system. It is important to observe that the expression (16) is ill-conditioned for θ → 0 and we found it convenient to replace it with its Taylor expansion up to o(θ 6 ) for θ < 0.1.
A completely analogous formula allows to reconstruct the placement of the rod starting from P 0 at one end an iteratively computing P k+1 = exp (s k+1 − s k )L k,k+1 P k , where L k,k+1 is the matrix associated with the strains assumed constant on the segment [s k , s k+1 ].

Data Availability Statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.