A Lie group variational integration approach to the full discretization of a constrained geometrically exact Cosserat beam model

In this paper, we will consider a geometrically exact Cosserat beam model taking into account the industrial challenges. The beam is represented by a framed curve, which we parametrize in the configuration space S3⋉R3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb{S}^{3}\ltimes \mathbb{R}^{3}$\end{document} with semi-direct product Lie group structure, where S3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb{S}^{3}$\end{document} is the set of unit quaternions. Velocities and angular velocities with respect to the body-fixed frame are given as the velocity vector of the configuration. We introduce internal constraints, where the rigid cross sections have to remain perpendicular to the center line to reduce the full Cosserat beam model to a Kirchhoff beam model. We derive the equations of motion by Hamilton’s principle with an augmented Lagrangian. In order to fully discretize the beam model in space and time, we only consider piecewise interpolated configurations in the variational principle. This leads, after approximating the action integral with second order, to the discrete equations of motion. Here, it is notable that we allow the Lagrange multipliers to be discontinuous in time in order to respect the derivatives of the constraint equations, also known as hidden constraints. In the last part, we will test our numerical scheme on two benchmark problems that show that there is no shear locking observable in the discretized beam model and that the errors are observed to decrease with second order with the spatial step size and the time step size.

Internal constraints may be introduced to enforce the rigid cross sections to remain perpendicular to the center line. In that way, some of the high-frequency solution components are eliminated and the full Cosserat beam model is reduced to a Kirchhoff beam model [17]. Space discretization results in equations of motion that form differential-algebraic equations (DAEs) of index 3 and may be solved efficiently by standard methods in this field [15]. The parametrization of rotations by unit quaternions helps to avoid singularities in simulation scenarios with large rotations and reduces the number of arithmetic operations per function evaluation (compared to mathematically equivalent approaches) [15]. On the contrary, the corresponding configuration space is nonlinear and normalization conditions need to be enforced during time integration at each time step if the equations of motion are embedded in higher dimensional linear spaces [30]. As an alternative to this brute force approach, the nonlinear structure of the configuration space may be addressed directly using Lie group integrators [9,19,33]. The performance of Lie group integrators may depend strongly on the specific Lie group formulation [5] with clear benefits for formulations that are based on the special Euclidean group [26,34]. For rotations being parametrized by unit quaternions, the corresponding Lie groups are multiples of the semi-direct product S 3 R 3 , see [11,19]. This semi-direct product (as well as the special Euclidean group) refers to rigid body motions in space and its application in beam theory helps to reduce the risk of locking phenomena [21,34], see also the last part of the present paper. Recently published paper [21] follows a very similar approach, while using higher order integration schemes but being restricted to the unconstrained case. In the present paper, however, the straightforward incorporation of internal and external constraints allows, e.g., enforcing the cross sections to remain normal to the centerline tangent such that transverse shearing is inhibited, see also [15,17] as well as Sect. 3.5 below.
Our research was inspired by industrial challenges that are studied in the Horizon 2020 Marie Skłodowska-Curie Innovative training network "THREAD" [10]. Twelve European research groups aim at new modeling and simulation techniques for highly flexible structures in system dynamics [3]. In fact, the applications of rod model simulation in the industries are numerous and various, ranging from cable harnesses in automotive applications and endoscopes in biomedical engineering via yarns in braiding machines to ropeway systems [10]. In industrial applications, static and quasistatic solution techniques are still dominating, but the THREAD project will also result in efficient and numerically stable methods for dynamic problems.
To this end, we consider coarse grid discretizations that preserve certain geometrical properties of the equations of motion and achieve numerical stability in that way [9]. Following the framework of variational integration in space and time [8,19], we discretize the action integral before solving the variational problem, which proves to be favorable in longterm simulation [24]. Internal or external constraints introduce extra terms in the Lagrangian that finally result in constraint equations and constraint forces in the fully discretized model that may be considered, e.g., by (nearly) energy preserving methods from molecular dynamics [1,12]. Restricting the analysis to low order methods, we get discrete equations of rather simple structure that are tailored to low or moderate accuracy requirements. Constraints are enforced in a stable way using a discontinuous approximation of the Lagrange multipliers [11] that allows enforcing the constraints and hidden constraints at the level of position and at the level of velocity coordinates simultaneously.
In this paper, we will use the well-known ingredients of Lie group structured configuration spaces, variational integrators, and numerical treatment of DAEs, to find a way of simulating thin beams efficiently and accurately at the same time. Following the approaches of Cardona, Zupan, Brüls, and others [4,6,37,38], we try to provide this investigation in a language that is close to the implementation with entities that bear a physical meaning.
The remaining part of the paper is organized as follows: In Sect. 2, we present the configuration space. We define in detail the Lie group and its operations. In particular, we explain how to use unit quaternions to describe rotation and which operations we can perform with them. The tools needed in the following study are mentioned. In Sect. 3, we derive the continuous equations of motion for the Cosserat beam model. We evaluate the expression for the kinetic and potential energy and we also consider the constraint as a different kind of energy, such that we are able to write the action integral. Using the action integral, we can apply the variational principle to obtain the continuous equations of motion. In Sect. 4, we are able to fully discretize the equations both in space and in time. We directly discretize the equations after applying the Hamilton's principle. Finally, in Sect. 5 we present our numerical experiments. We analyze how the steps in arc length and the time step sizes influence the errors.

A Lie group structured configuration space
In this section, we will discuss configuration spaces with Lie group structure in general, the configuration space S 3 R 3 in particular, and how we can find a Lie group structure of a configuration space. Lie groups are well researched theoretically [36]. They also have been applied in numerical integrations [14], and especially in the description of mechanical systems [4,18,26] and beams in particular [7,8,20,33]. Note that for Lie groups, we will use the notations of Brüls et al. [4][5][6]. This notation, as well as the notation for unit quaternions, is intentionally implementation oriented avoiding notational shortcuts such that it becomes very clear how formulae must be or should be implemented.
Position Let us first consider the position of a mass point in space. Usually, the set R 3 is used as a configuration space, where any x ∈ R 3 is a possible position of the mass point. Now let us instead consider x ∈ R 3 as a translation of the mass point from some reference position. Now we can apply two of these translations x, y ∈ R 3 , leading us to the concatenation x + y ∈ R 3 . It can be seen that with this operation +, the configuration space R 3 actually becomes a linear Lie group (R 3 , +) with identity element 0 ∈ R 3 and inversion x → −x.
Orientation Now we consider the orientation of a rigid body in space. Here, we will use unit quaternions in order to represent orientation, where • 2 is the Euclidean norm. As before, we consider any p ∈ S 3 to represent a rotation from some reference configuration of the rigid body. Concatenation of two unit quaternions p, q ∈ S 3 leads to the quaternion multiplication where each of the quaternions p, q ∈ S 3 are decomposed into their scalar (or real) parts p 0 , q 0 ∈ R and their vector (or imaginary) parts p, q ∈ R 3 with Note that the quaternion multiplication is also defined for non-unit quaternions in R 4 . Again, we end up with a Lie group (S 3 , * ) with identity element e = [1, 0, 0, 0] and inversion 1 This Lie group, in contrast to the linear Lie group, is not commutative. Further, we define the rotation R(p) a of a vector a ∈ R 3 with respect to a unit quaternion p ∈ S 3 as Note that antipodal unit quaternions p, −p ∈ S 3 represent the same rotation: for all p ∈ S 3 and a ∈ R 3 . 2 Rigid bodies Now we repeat the above process with a rigid body in space, that has an orientation p ∈ S 3 and a position x ∈ R 3 . Now we consider (p, x) with p ∈ S 3 , x ∈ R 3 to be a rigid body motion The concatenation of two of these rigid body motions associated with the tuples (p 1 ,

leads us to an operation
since the rotation R is linear in its vector argument and it holds R(p 1 ) R(p 2 ) a = R(p 1 * p 2 ) a. It can be easily seen that the operation • is smooth and has an identity element e = ([1, 0, 0, 0] , 0), each element having an inverse and the inversion is also smooth. Since S 3 and R 3 are clearly smooth manifolds, (S 3 R 3 , •) is a Lie group. We will call the Lie group (S 3 R 3 , •) or, by abuse of notation S 3 R 3 , because it is not a direct product of the Lie groups (S 3 , * ) and (R 3 , +) but rather a semi-direct product [11]. The dimension of S 3 R 3 is n = 6, since it has dimension six as a smooth manifold. If we would have used rotation matrices SO(3) instead of unit quaternions to characterize rotation or orientation respectively, we would have arrived at the special Euclidean Lie group (SE(3), ·).

Derivative vectors
We consider now configuration spaces G with dimension n that have a Lie group structure (G, •) with an identity element e ∈ G. If we consider a differentiable curve q(t) ∈ G, each derivative is an element of the corresponding tangent space:q(t) ∈ T q(t) G. The Lie group structure now gives us the possibility to expressq(t) aṡ where v(t) ∈ T e G is an element of the tangent space at e for all t and dL q(t) (e) is the differential of the left translation defined by L q 1 (q 2 ) = q 1 • q 2 . The tangent space T e G = g is called the Lie algebra of G. Instead of handling elements v(t) ∈ g of an abstract tangent space, we will express v(t) as the velocity vector v(t) ∈ R n , which bears a physical meaning [4,37]. The identification of g and R n is done by choosing a linear vector space isomorphism • : R n → T e G. This concept can be applied to arbitrary real variables other than time t in which case we speak of a derivative vector.
In the case of S 3 R 3 , we have and the physical meaning is that Ω ∈ R 3 is the angular velocity and U ∈ R 3 is the velocity of a rigid body described by the configuration (p, x). Note that both Ω and U are expressed with respect to the body-fixed frame {R(p) e 1 , R(p) e 2 , R(p) e 3 }.
Exponential map Now we will define the exponential map exp : Usually, the exponential map exp : T e G → G is used, but since we are using the concept of velocity vectors, exp as the concatenation of exp and the tilde operator is easier to handle. If we consider a differentiable curve v(t) ∈ R n , then the derivative of exp v(t) is given by where T : R n → R n×n is called the tangent operator. 3 In the Lie group S 3 R 3 , the exponential map is given by for [Ω , U ] = 0 and exp(0) = ([1, 0, 0, 0] , 0). Here, T S 3 is the tangent operator that is associated with the Lie group (S 3 , * ) as well as with the Lie group of rotation matrices (SO(3), ·). The tangent operator for S 3 R 3 actually coincides with the tangent operator for the Lie group (SE(3), ·), which can be found in [5]. Further, we will also need the inverse T −1 (x) = T(x) −1 of the tangent operator. The tangent operators that are not already defined here and their inverses are given in the Appendix.

Logarithm
The exponential map exp is injective in a neighborhood of the origin. 4 We call the inverse function the logarithm log which maps Lie group elements g ∈ G that are sufficiently close to the identity element e ∈ G back to R n such that exp log(g) = g. For the Lie group S 3 R 3 , the logarithm is given by for (p, x) ∈ S 3 R 3 sufficiently close to the identity element, where log S 3 is the logarithm of the Lie group (S 3 , * ) and is given in the Appendix.
Hat operator We will now introduce the hat operator •: R n → R n×n by the relation where σ a and σ b are differentiable curves in G with σ a (0) = σ b (0) = e and velocity vectors a and b at the origin, respectively. We can see that the hat operator is actually the adjoint operator formulated in the language of velocity vectors with the relation 5 In fact, the hat operator is used to write down the tangent operator as well as its inverse in a series from which the closed expressions can be derived: where B k are the Bernoulli numbers with B 1 = −1/2, see [25]. In the case of the Lie group S 3 R 3 , we have for Ω, U ∈ R 3 . Like the tangent operator, the hat operator of S 3 R 3 coincides with the hat operator of SE(3), [33].

Derivative vectors of bivariate functions
Let us now consider a function q : R 2 → G that is differentiable and has derivative vectors x(x, y), y(x, y) ∈ R n with respect to x and y, respectively. Then there is this very nice relationship between the partial derivatives of the derivative vectors: see, e.g., [5], where this was used with d/dy being the variation and x = t .
A differential operator Now we will consider functions Ψ : G → R k whose arguments are elements of a Lie group. The derivative dΨ (g) : T g G → R k for g ∈ G is then a linear function on the tangent space T g G. We can formulate dΨ (g) in the language of velocity vectors and introduce the differential operator D, where DΨ (g) is defined by for all w ∈ R k . The actual expression for DΨ (g) can often be calculated easily by considering a differentiable curve q(t) ∈ G with velocity vectors v(t) and considering We can use the differential operator D in order to calculate the derivative of the logarithm: for all q ∈ G close enough to the identity. This can be shown by differentiating exp log(q) = q and using (1). Furthermore, we can calculate the derivative of nonlinear differences by Interpolation Lastly, we will introduce an interpolation function that is a smooth mapping in ξ with constant derivative vectors log(a −1 • b): and it holds Ip(0; a, b) = a and Ip(1; a, b) = b. In the context of the linear Lie group (R n , +), where exp and log are the identity maps, Ip is just the regular linear interpolation, making Ip a generalization of linear interpolation to Lie groups. For more information on Lie group interpolation see, e.g., [34].

The continuous equations of motion
In this section, we consider the continuous equations of motions of a constrained Cosserat beam. This model was originally formulated by Simo [31] and then reformulated by Lang and Linn [15,17]. A similar approach to our was followed in [9,33], where the authors had used rotation matrices instead of unit quaternions. Celledoni et al. [7] had used quaternions but not semi-directly linked with the positions in order to describe an unconstrained Cosserat beam model.

Geometrically exact beams
We consider the Cosserat beam to be a framed curve or, in the context of our Lie group, a curve q(•, t) ∈ S 3 R 3 at each time instant t . The configuration q(s, t) = p(s, t), x(s, t) then describes the position x(s, t) ∈ R 3 and the orientation p(s, t) ∈ S 3 of the cross section at s. We assume that s represents the arc-length of the undeformed beam and that deformation of the cross section is negligible, such that it can be considered rigid.
The velocity vectors v(s, t) ∈ R 6 then contain the velocity U (s, t) ∈ R 3 and the angular velocity Ω(s, t) ∈ R 3 of the rigid cross section expressed with respect to the body-fixed frame:

v(s, t) = Ω(s, t) U (s, t) .
In the same way, we can consider the derivative vectors w(s, t) with respect to the spatial parameter s: In a similar way as before, w(s, t) contains the stress Γ (s, t) ∈ R 3 and material curvature K(s, t) ∈ R 3 expressed with respect to the body-fixed frame: Here, K * (s) ∈ R 3 and Γ * (s) ∈ R 3 are the material curvature and stress of the undeformed beam. We will summarize those two quantities to

Energies
We will now give expressions for kinetic and potential energy of the Cosserat beam model. The kinetic energy T at each time instance t is given by the integral of a kinetic energy density where M ∈ R 6×6 is a mass matrix given by where ρ > 0 is the material density, A > 0 is the area of the cross section, I 1 , I 2 , J > 0 are the principal moments of inertia of the infinitely thin cross section. Of course, this assumes that the material is homogeneous and the shape of the cross section does not change. The kinetic energy is then given by Now we give the potential energy U under the assumption of linear material behaviour. As before, the potential energy is given by the integral over the whole beam of the potential energy density Here, the matrix C ∈ R 6×6 is a diagonal matrix that contains the moduli E, G > 0 as well as the moments of inertia I 1 , I 2 , J > 0 and the area A > 0 of the cross section. The quantities A 1 > 0 and A 2 > 0 are often identified with the product of the area A with some Timoshenko shear correction factors. In practice, the products EI 1 , EI 2 , GJ, GA 1 , GA 2 and EA are determined experimentally instead of calculated from I 1 , I 2 , J, A and the experimentally obtained moduli and Timoshenko shear correction factors. The potential energy is given by

Derivation of the equations of motion
We want to derive a constrained Cosserat beam model that can be used to reduce the full Cosserat model, e.g., to a Kirchhoff beam model. In case that the quantities GA 1 and GA 2 are very large, the problem becomes very stiff. By reducing the full Cosserat model to a Kirchhoff model, we neglect the shearing and therefore remove the stiffness resulting from these large parameters. In order to introduce such constraints, we will require Ψ w(s, t) = 0 ∈ R k for all admissible s, t with some function Ψ . 6 We will later derive the equations of motion of the beam by a variational principle with an augmented Lagrangian. In order to simplify the notation, we introduce the quantity that contains the constraint-related terms. Here λ(s, t) ∈ R k are Lagrange multipliers. Now we can define the action integral with augmented Lagrangian. We introduce the generalized internal force densities The first term characterizes general internal force densities arising from bending, torsion, extension, as well as shearing and the second term characterizes generalized internal constraint force densities.
Since we are considering a beam that is free on both ends, we assume that the generalized internal force densities vanish at both ends of the beam: We call these natural boundary conditions. Now we can derive the equations of motion by Hamilton's principle [24], see also [15,18,22] with an augmented Lagrangian. As before, we use the concept of derivative vectors for the variation: The boundary conditions for the derivative vectors with respect to the variation are assumed to be Now we can interchange variation and double integration, use (2) with respect to time and space, apply partial integration and end up with, by omitting the arguments s and t , since the boundary terms of the partial integrations vanish due to the natural boundary conditions of the beam and (8). The term δw can be expressed by w and δq by using (2). The above equation has to hold for all variations δq(s, t) and δλ(s, t), and, therefore, the factors of those variations have to vanish, leading us to the continuous equations of motion for all s ∈ (0, L) and t ∈ (t 0 , t e ).

Generalizations
As in [17], external forces are added by where F (s, t, q) is a density of external generalized forces that may depend on s, t , q as well as any derivatives of q. Applying the same technique as before, we end up with the equations of motion where we have just added the term F (s, t, q) to the right hand side of the dynamic equations (9a). A viscoelastic beam could be simulated by adding damping. This can be incorporated by with a diagonal matrix D ∈ R 6×6 , containing dissipative material constants [15,17]. Here, the term 2D ·ẇ(s, t) could also be interpreted as the derivative of a dissipative power densitẏ with respect toẇ(s, t). Treating the additional integral in the same way as before, we arrive at the same equations of motion (9a), (9b), but with generalized internal forces G(s, t) = C · w(s, t) − w * (s) + 2D ·ẇ(s, t) + DΨ w(s, t) · λ(s, t).

Internal constraints
In the case of a Kirchhoff beam model, we require the cross sections to be perpendicular to the center line, so the first two components of the material stress vector Γ (s, t) have to vanish. This is equivalent to In the case of an inextensible Kirchhoff beam model, we also require that the beam remains parametrized by arc-length and therefore Γ (s, t) = 0, which is equivalent to

Equations of motion in the inertial frame
The equations of motion (9a), (9b) may not look familiar to the reader. If, however, we consider an unconstrained beam, separate equation (9a), into two equations in R 3 and formulate them with respect to the inertial frame by applying the rotation p(s, t) from q(s, t) = p(s, t), x(s, t) , we end up with the familiar equations with internal moments m(s, t) ∈ R 3 and forces f (s, t) ∈ R 3 , the velocity u(s, t) =ẋ(s, t) ∈ R 3 , the angular velocity ω(s, t) ∈ R 3 as well as the inertia tensor i(s, t) ∈ R 3×3 with respect to the inertial frame: for all a, b ∈ R 3 . We prefer, however, to look at the equations of motion in slightly more complicated form (9a), (9b) because it agrees well with our choice of configuration space and, more importantly, the mass matrix M is constant unlike the inertia tensor i(s, t).

The fully discretized equations of motion
In this section, we will concentrate on the discretization of the Cosserat beam. In order to do this, we will heavily rely on the toolbox of variational integrators [24]. A similar approach was used in [8,9,20].

Discretization of the functional space
In order to obtain fully discrete equations of motion of the constrained beam, we will reduce the number of functions that we consider in the variational principle (7). First, we consider a spatial grid 0 = s 0 < s 1 < · · · < s K = L as well as a time grid with step sizes Furthermore, we define the midpoints of the spatial grid s k+ 1 /2 = s k + s k+1 2 as well as the interior of the space-time grid cells We consider only trajectories q d that take the form , which are interpolations of the data (s k , t (n) , q (n) k ) in the Lie group [34]. The trajectory q d is of course continuous everywhere on [0, L] × [t (0) , t (e) ] and its degrees of freedom are its values at the grid points q d (s k , t (n) ) = q (n) k ∈ G for k = 0, . . . , K and n = 0, . . . , N. Notice, however, that q d is only almost everywhere differentiable: The spatial derivative does not have to exist at all spatial grid points s = s k and the time derivative does likewise not have to exist at all time grid point t = t (n) . On the contrary, inside a grid cell C (n) k , the spatial derivatives are constant along a fixed time t * and the time derivatives are constant along a fixed arc length s * . More importantly, we give names to the spatial derivative along the spatial grid points and to the time derivative along the time grid points, again using the concept of derivative and velocity vectors: with using (6), Furthermore, we define the derivative vectors with respect to space of the undeformed configuration of the beam w * k+ 1 /2 = w * (s k+ 1 /2 ).
We do not need to assume a special kind of function for the Lagrange multipliers λ, but we allow the Lagrange multipliers λ d (s, t), corresponding to the trajectories q d (s, t), to be discontinuous in time at all time grid points t = t (n) . We give the following names to their one-sided limits around t (n) in the middle of the spatial grid points: We will later see that both limits are required to respect the hidden constraints.

Derivation of the discrete equations of motion
Now we split up the action integral in order to approximate the integral of the Lagrange function over a space-time grid cell Here we have used the linearity of the integral and interchanged the order of integration in one of the terms. In order to do this we have chosen to apply the trapezoidal rule to the inner integrals and the midpoint rule to the outer integrals: Using these second-order accurate approximations, we can now define a discrete Lagrangian L d q (n) k , q (n) k+1 , q (n+1) k , q (n+1) k+1 , λ (n,+) k+ 1 /2 , λ (n+1,−) Further, we will use the abbreviation in order to keep the notation easier to read. Due to the earlier second-order approximations, we now know thatC as well as We call the sum over all discrete Lagrangians the discrete action S d (q (n) k ) k,n , (λ (n,+) k+ 1 /2 ) k,n , (λ (n,−) which depends on all discrete configurations q (n) k as well as on all one-sided limits of the Lagrange multipliers λ (n,±) k+ 1 /2 . Now we continue along the lines of variational integrators and derive the discretized equations of motion by finding a stationary point of the discrete action with the standard assumption for the variation of the endpoints in time: ). Now we will use the linearity of the sum, shift the indices and use the assumption on the variation of endpoints in time and end up with by defining the missing discrete Lagrangian derivatives as zero in order to facilitate the notation: Since the above equation has to hold for all variations, we obtain the discrete equations of motion in the form for k = 0, . . . , K and n = 1, . . . , N − 1. We can now rearrange the first equation to The left and right hand side of the equation define two canonical momenta p (n,+) k and p (n,−) k , respectively. The equation can thus be interpreted as an instance of momentum matching [24], leading to well-defined values of the momenta at the node. Here, we identify the moments by using mass matrix M multiplied with velocities v (n) k instead: This means we have This approach is called a discrete Legendre transform [24] and the transformation from momenta to velocities can be considered the application of the inverse continuous Legendre transform. Now we can consider the left hand equation and the right hand equation of (14) independently and get, by shifting n by one in the right hand equation Note, that we extend the validity of the equation (15a) to n = 0 and (15b) to n = N − 1 by arguing that we could have started already at t (−1) = t (0) − ε and ended at t (n+1) = t (n) + ε with some 0 < ε 1.
In order to obtain the derivatives of the Lagrangian, we will need the derivatives of the velocity vectors v (n+ 1 /2) k and w (n) k+ 1 /2 with respect to the q (n) k by which they are defined, see (4a), (4b) and (12a), (12b): Now, the derivatives of the discrete Lagrangian can be written as We can see that the derivatives of the discrete Lagrangian with respect to the two different Lagrange multipliers, λ (n,+) k+ 1 /2 and λ (n,−) k+ 1 /2 , both give us the same constraint equation only with shifted index n. We will drop one of the redundant equations and instead replace it by the constraint equation on velocity level This constraint (16) is often called the hidden constraint, because it can be revealed by differentiating the constraint equation on position level with respect to time. Assuming that v (n) k actually approximates the velocity vector of q (n) k with respect to time, we define Now, we can put equations (15a), (15b), (13c), (16), (12a), as well as (12b) together and obtain -by using the derivatives of the discrete Lagrangian -the full discrete equations of motion for k = 0, . . . , K and n = 0, . . . , N − 1. Here G (n,±) k+ 1 /2 are internal generalized forces given by for n = 0, . . . , N and k = 1, . . . , K − 1. All quantities in (17a)-(17e) with lower index −1, − 1 /2, K + 1 /2 or K + 1 are defined to vanish in order to facilitate notation; e.g., equation (17b) for k = 0 takes the form

Generalizations
External forces Now using the variational principle (10) instead of Hamilton's one (7), we will consider external forces. Since we already know how to treat δS , we will focus on the second term. We insert the trajectories q d , split the integration interval, approximate s, t)ds/ s k , use the trapezoidal rule for the remaining integral and finally shift indices:

s, t (n) )ds
for k = 1, . . . , K − 1 and n = 0, . . . , N and the terms vanish for lower indices less than 0 or greater than K. Now the term F (n) k will appear on the left hand side of (13a) and when this equation gets rearranged, we make sure that the products with t (n) come to the left side and the products with t (n−1) to the right side of (14). This will lead to the extra term on the right hand side of (17b), as well as on the right hand side of (17d).
Internal damping Now we will consider internal damping of the beam. In order to do this, we will use the variational principle (11). Again, we insert q d and approximate the remaining integral, since we already dealt with δS (q d , λ d ). We know that the derivative vectors of q d (s, t (n) ) are constant w (n) k+ 1 /2 for s ∈ (s k , s k+1 ) and thatẇ (n) k+ 1 /2 ≈ d/dt w d (s k+ 1 /2 , t (n) ). We split the integration range and subsequently apply the midpoint rule in space and the trapezoidal rule in time: Now we can calculate the variations of w (n) k+ 1 /2 using (4a), (4b) and end up, by shifting indices, with where all quantities with lower index less than 0 or more than K vanish. As in the case of external forces, we make sure that the extra terms that have the factor t (n) get written on the left hand side of (14) and the terms with factor t (n−1) on the right hand side of (14). We notice, however, that in the equations of motion, the transposed inverse of the tangent operator can be factored, leading to the exact same equations of motion (17a)-(17e), where the generalized internal forces now include the damping term: where all quantities with lower index less than 0 or more than K vanish.

Implementation
We have formulated a high-level Algorithm 1 to solve the equations of motion (17a)-(17e). Note that the system in line 5 is linear if there are no external forces that depend on time-derivatives of q. In that case, the Newton-Raphson iteration requires exactly one step. The starting values can be chosen by quantities that are already known from the previous time step or the current time step. It is noteworthy that the Jacobians ) K−1 k=0 . 6: end for of the systems in lines 3 and 5 have band structure [2,15] if the equations and unknowns are ordered with ascending spatial index k including the constraint equations and Lagrange multipliers. Therefore, the linear systems that appear in the Newton-Raphson method can be solved very efficiently.

Numerical experiments
In this section, we will consider two different numerical benchmark problems. The first one will be a heavily damped beam that is rolled up to form a circle. The computations of the benchmark problem often give faulty results if locking is present in the numerical approximation [28]. We will see that in our case, no locking can be observed. The second benchmark will be the so-called flying spaghetti [32], where we will study the influence of the spatial step size as well as the time step size on the overall accuracy of the approximate solution obtained by the numerical scheme.
We have implemented the equations of motion derived in Sect. 4 by applying the existing numerical scheme RATTLie [12] to the spatially discretized beam model in a method of lines approach [13].

Roll-up
We consider an unconstrained Cosserat beam of length L = 10 and material properties GA 1 = GA 2 = EA = 10 4 , EI 1 = EI 2 = GI = 500, Aρ = 1, I 1 = I 2 = J = 10 as well as viscoelastic parameters D = diag(100, 100, 100, 100, 100, 100). The beam is clamped at s = 0, where we prescribe x(0, t) = 0 and p(0, t) = [1, 0, 1, 0] / √ 2 for all times t ≥ 0. To the initially straight beam x(s, 0) = sL and p(s, 0) = [1, 0, 1, 0] / √ 2, we apply a constant external torque M = [0, 2πEI 1 /L, 0] to the free end s = L. The discretization is done with an equidistant spatial grid with K = 8 and constant time steps of t = 2 −10 ≈ 9.7×10 −4 . In Fig. 1, we have shown snapshots of the center line; and, in Fig. 2, we have shown the norm of the distance between the two endpoints (p (n) 0 , x (n) 0 ) and (p (n) K , x (n) K ) as well as the norm of the velocity vectors [(Ω (n) K ) , (U (n) K ) ] . It can be seen that after a few seconds, the beam is almost at rest due to the relatively high internal damping. We observe an equilibrium where the cross sections of both ends are almost in the same position and orientation. It is known that if a beam model has shear locking, the bending stiffness is largely overestimated and the beam would not complete the ring [28]. Since this behavior cannot be observed here even for this very coarse spatial discretization, we can conclude that shear locking is not present in our beam model. The recently published paper [21] describes a similar experiment. Furthermore, the authors of [21] prove that there is no locking in their model because  of the semi-direct product structure of the unit dual quaternions used in [21], similar to the structure of S 3 R 3 used in the present paper. 7

Flying spaghetti
Now we will apply our discretization to the flying spaghetti benchmark [32], in order to numerically observe the order of convergence in space and time. The flying spaghetti benchmark consists of an initially straight beam of length L = 10 with material parameters: GA 1 = GA 2 = EA = 10 4 , EI 1 = EI 2 = GJ = 50, Aρ = 1, and ρI 1 = ρI 2 = ρJ = 10. Both ends of the beam are free, and at the end, where s = 0, forces and moments are applied, ramping up from t = 0 to t = 2.5 and then decreasing from t = 2.5 to t = 5. A schematic description of the initial configuration and the applied moments can be found in Fig. 3. In Fig. 4, we have provided snapshots of the configuration of the beam. In this paper, we used a Kirchhoff beam model for the flying spaghetti by constraining the full Cosserat beam model as described in the above sections.
where Q (n) k ≈ q(s k , t (n) ) is the reference solution. Accordingly, the maximal absolute discrete L 2 error in v can be defined. Note that we are using here differences of Lie group elements, but only for getting a grasp of the magnitude of the errors not for computing any values of the numerical solution! We can see in Fig. 5 that the slope of the errors in q and v over the time step size t is approximately two. This means that we can observe second order convergence behaviour numerically.
In Fig. 6, we have considered the flying spaghetti with equidistant time grid with t = 2 −12 ≈ 2.5×10 −4 and varied the number of equidistant spatial grid points K. This time, we have interpolated the beam in space [34] by using piecewise Lie group interpolation (5) in q and piecewise linear interpolation in v. Then we have plotted the maximal absolute discrete L 2 error of the interpolated beam at s = 0, L/8, . . . , 8L/8. Here, we can observe that the slope of the errors over K is approximately −2, which means that we can observe second order convergence in space. The error in the velocity vectors v, however, starts to saturate for larger K. This could have two reasons: All solutions were calculated with a rather coarse time step, and, so, this error saturation could indicate that the solutions can not become a lot more precise without reducing the time step size. The second reason could be that piecewise linearly interpolating the velocity vectors v could be incompatible with the way the discretization scheme was derived, where the velocity vectors of the discretized All in all, we can observe that our discretization scheme is convergent of second order in space and time. This is what could have been expected, since we have only used approximations of second order in the derivation of the scheme, namely Lie group interpolation of first order (which is a generalization of linear interpolation) and the trapezoidal and midpoint rule for approximating integrals. A rigorous mathematical analysis is, however, needed in order to prove the claim of second order convergence.
Lastly, we have depicted the change of mechanical energy of the flying spaghetti in Fig. 7 in the conservative realm t ∈ [5,110]. As we expect from a variational discretization, there is no systematic energy drift, see, e.g., [24], where the discrete Noether theorem is shown, stating that a variational integrator has to conserve a perturbed energy.

Conclusion
We have shortly introduced the toolbox of Lie group structured configuration spaces and the concept of velocity and derivative vectors with special attention to the Lie group of rigid body motions S 3 R 3 , where rotations are parametrized using unit quaternions S 3 . Then we have described a constrained Cosserat beam model, where the configuration space S 3 R 3 was used to describe the orientation and position of each cross section of the beam. The constraints can be used to reduce the Cosserat beam model to a Kirchhoff beam model by requiring the cross sections to remain perpendicular to the center line. The continuous equations of the beam were then derived using variational principles like Hamilton's principle with an augmented Lagrangian. Furthermore, we have discretized the infinite-dimensional space of functions by only considering trajectories, which are given by a piecewise first order Lie group interpolation. Applying again the variational principles, we arrived at the fully discrete equations of motion. The derived discrete scheme was then tested on two benchmark problems. We could observe that no shear locking was present in the scheme and that it can be numerically seen to converge with second order both in space and in time.

Funding Note Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.