Local Coordinates on Lie Groups for Half-Explicit Time Integration of Cosserat Rod Models with Constraints

Explicit Runge-Kutta methods are the gold standard of time integration methods for non-stiﬀ problems in system dynamics since they combine a small numerical eﬀort per time step with high accuracy, error control and straightforward implementation. For the analysis of beam dynamics, we couple them with a local coordinates approach in a Lie group setting to address large rotations. Stiﬀ shear forces and inextensibility conditions are enforced by internal constraints in a coarse grid discretization of a geometrically exact beam model. The resulting non-stiﬀ constrained systems are handled by a half-explicit approach that relies on the constraints at velocity level and avoids all kinds of Newton-Raphson iteration. We construct half-explicit Runge-Kutta Lie group methods of order up to ﬁve that are equipped with an adaptive step size strategy using embedded Runge-Kutta pairs for error estimation. The methods are tested successfully for a roll-up maneuver of a ﬂexible beam and for the classical ﬂying spaghetti benchmark.


Introduction
Geometric numerical methods preserve essential properties of the flow of differential equation models [11].Such strategies are important for the robust, efficient and numerically stable dynamical simulations of highly flexible structures.In space, they are used to get coarse grid discretisations that reflect the characteristic behaviour of geometrically exact beam models.In time domain, geometric integration is tailored to nonlinear configuration spaces that are typical of mechanical systems with large rotations.In the present paper, we address the latter aspect and present a novel Runge-Kutta Lie group integrator for constrained mechanical systems and its application to beam dynamics.The method is half-explicit and avoids all kinds of Newton-Raphson iterations which are a bottleneck for the efficiency of classical implicit integrators in beam analysis.Stiff shear forces are represented by internal constraints [16,17] that may be used as well to enforce a beam's inextensibility.These internal constraints are combined with a coarse grid space discretisation and result finally in a constrained system that is for typical application scenarios non-stiff in its differential part.The equations of motion of constrained mechanical systems form differentialalgebraic equations (DAEs) with Lagrange multipliers as algebraic variables that couple the constraints to the equilibrium equations for forces and moments [9,11].Following the approach of Brasey and Hairer [5], the equations of motion are transformed analytically before half-explicit time discretisation resulting in the index-2 formulation with constraint equations at the level of velocity coordinates.In that way, the velocity vector is constrained to the null space of the constraint Jacobian that is known from null space approaches like the one of Betsch [4].In the context of Lie group time integration, these velocity constraints have been used before in the RATTLie integrator [13] and in the application of generalized-α methods to the stabilized index-2 formulation of the equations of motion [2].The successful application of generalized-α and other Newmark type Lie group integrators in flexible multibody dynamics [6,7,9] relies on the specific Lie group structure of the configuration space for mechanical systems with large rotations [19,20].The Lie group setting is attractive in this context since it allows to describe large rotations and the orientation of bodies and flexible structures without any singularities [24].The spaces one came across are the space of Special Orthogonal matrices SO (3) or the space of unit quaternions S 3 , their direct or semi-direct product with R 3 and tensor products thereof.The half-explicit Runge-Kutta Lie group integrator is tested for a geometrically exact beam model that goes back to Lang and Linn [16,17].The space discretisation follows the variational integration approach of Hante [12,14] with nodal variables (q i , x i ) ∈ S 3 ⋉ R 3 .The semi-direct product of S 3 and R 3 may be considered to be isomorphic to the Special Euclidean group SE(3), taking into account that S 3 is a double covering of the Special Orthogonal group SO(3), i.e. q, −q ∈ S 3 represent the same rotation matrix R ∈ SO(3).In the work of Munthe-Kaas [21], see also Iserles et al. [15], classical Runge-Kutta methods are generalised to Lie group integrators by the help of the exponential map that defines local coordinates on the Lie group by a tangent space parametrisation.The possibility to work locally in the tangent space, which is a linear space with the well known operations among its elements, simplifies substantially the numerical solution of a system set in a Lie group.The use of the exponential map assures that the numerical solution remains in the Lie group, given that we work with local coordinates.In a general Lie group setting, the construction and implementation of higher order Runge-Kutta Munthe-Kaas methods is challenging since they involve frequent evaluations of matrix exponentials [18] and the approximation of its derivatives [21].There are no such problems in the application to SO(3) since Rodrigues' formula allows to evaluate the exponential map in closed form and a similar expres-sion may be found for its derivative [15].Park and Chung generalised this approach to the Special Euclidean group SE(3) and pointed out that any classical explicit Runge-Kutta method may be generalised in that way to a Lie group integrator on SE(3) with the same order of convergence [23].Several recent studies focus on the closed form expression of the exponential map and its derivatives, as well as their efficient evaluation [12,27,29].Additionally, one may be concerned with occurring singularities in the origin, which are solved by extended investigation of the approximation of the closed form by a Taylor polynomial [12,Appendix B].The present paper considers non-stiff integrators that are a quasi-standard in nonlinear system dynamics, such as the default integrator ode45 in MATLAB's ODE suite [25], addressing two main challenges: Non-linear configuration spaces and constraints.We aim to solve the former by introducing Lie group integrators and the latter by adapting the existing numerical method of half-explicit Runge-Kutta type [5].

Half-explicit Runge-Kutta Lie group integrators
In the present section, we elaborate the definition of the method.In the early 1990's, Brasey and Hairer introduced the half-explicit Runge-Kutta methods for semi-explicit index-2 DAEs in Hessenberg form [5]. Some years later, Murua [22] and Arnold [1] proposed a modification, which simplifies the order conditions of higher methods and makes the approach more efficient.In each time step, the methods start with an explicit stage and evaluate in all later stages velocity vectors in the null space of the corresponding constraint Jacobians.

Equations of motion and local coordinates
Firstly, we consider a flexible multibody system in the Lie group G := (S 3 ⋉ R 3 ) N +1 , with N denoting the number of beam edges after space discretisation [12].The elements q ∈ G have 7(N + 1) components: 4 for the orientation and 3 for the position of each frame that is attached to one of the N + 1 nodes of this discretisation.Let e be the identity element of the Lie group G.We introduce the elements of the tangent space ṽ ∈ T e G, that summarise both angular and translational velocity.The tangent space T e G at the identity element has been identified with R 6(N +1) through an isomorphism called the tilde operator • : R 6(N +1) → T e G.By the use of the previous elements, we set the equations of motion for a constrained flexible multibody system as expressed in [6] q = DL q (e) • ṽ (1a) where M(q) is the mass and inertia matrix, g(t, q, v) is the vector of external and internal forces, and B(q) is the gradient of the constraint function at the position level Φ(q) in the sense of [6, Eq. ( 9)] Equations ( 1) are the index-3 formulation of the equations of motion.The right hand side DL q (e) • ṽ of the kinematic equation (1a) evaluates at y = e the directional derivative of the left translation L q (y) = q • y in the direction of ṽ ∈ T q G q(t) = lim ϵ→0 q(t) • exp (ϵṽ(t)) − q(t) ϵ .
We can perform an index reduction by differentiating in time the constraints in Equation (1c) The index-2 formulation of the equations of motion is q = DL q (e) • ṽ (5a) It is analytically equivalent to (1) and avoids the time-consuming evaluation of curvature terms in the index-1 formulation of the equations of motion [28,30].
Given the equations of motion, the solution in local coordinates, in a neighborhood of t = t n , is where θ n (t) denotes a local parametrization of the Lie group with θ n (t n ) = 0. Differentiating (6) in time, the time derivative of q(t) is given by Comparing this equation with (3), they have to coincide for finite ϵ up to higher order terms: Given that d dϵ exp θn (t) + ϵ θn (t) with dexp − θn (t) denoting the right trivialised differential of exp , see [15, Definition 2.18], we get and ( 8) results in ṽ(t) = dexp − θn (t) θn (t).(10) For practical computations, the equivalent representation in terms of the tangent operator [6,7] proves to be favourable: The dexp operator may be evaluated in terms of iterations of the adjoint operator [15, Equation (2.44)] resulting in a series expansion of the tangent operator T, see [7], that may be summarised to closed-form expressions for all Lie groups of practical relevance in multibody dynamics, see [15,23] and the more detailed discussion in subsection 2.3 below.

Half-explicit Runge-Kutta Lie group methods
In classical Runge-Kutta methods, we use parameters a ij , b j , c i , (i, j = 1, . . ., s) to build stage vectors and numerical solution.In particular, we evaluate s stage vectors using a ij , c i and we compute the numerical solution as a linear combination of the s stage vectors weighted by the b j parameters.When using the half-explicit Runge-Kutta methods, we may introduce more stage vectors (s ≥ s), and use the notation b j = a s+1,j , (j = 1, . . ., s).In the Lie group integrator, the stage vectors are The first stage (i = 1) is explicit with For i = 2, . . ., s, we enforce the hidden constraints at velocity level (5c) for stage vectors (i = 2, . . ., s).
Eqs. (13b,13c) form a system of linear equations in terms of Λ ni , Vni , (i = 2, . . ., s), that motivate the notation half-explicit for this class of time integration methods [5]: The total number of linear systems to be solved in each time step is s − 1.Note that the half-explicit approach avoids all iterations and does not rely on some kind of Newton-Raphson method.
The updated solution at time step t = t n+1 is where d j are additional parameters for the computation of the Lagrange multipliers.Their value, as for the other parameters, is set by order conditions and by a contractivity condition to ensure zero-stability.The method has been tested up to order five.We list here the Butcher tableaux with parameters a ij , (i = 1, . . ., s + 1, j < i) for the implemented methods [1].
Table 1: Parameters of half-explicit methods up to order p = 5.
Order 5, s = 6, s = 7 In a local coordinates approach, the order analysis in [1], valid for linear spaces, appears to be applicable also for nonlinear configuration spaces.In our case, the tests are performed for the Lie groups SE(3) and SO(3) × R 3 in subsection 2.4 and for Lie group (S 3 ⋉ R 3 ) N +1 for the applications to beam dynamics.

Implementation issues
As stated in the previous subsection, the Lie groups of interest are SO(3), S 3 , their direct or semi-direct product with R 3 and tensor products thereof.Inspired by Rodrigues' formula [15] exp we consider the closed form expression [12] exp Here and later in the current section we set ω = ∥Ω∥ 2 .
Following the approach of Hante [12], we consider the semi-direct product S 3 ⋉ R 3 that introduces a coupling between the orientation and the position variables, which results in the appearance of the tangent operator in the exponential map as it is known from SE (3).
where the tangent operator T S 3 (Ω) is At the same time we need to specify the tangent operator on S 3 ⋉ R 3 where the function In the space of interest for beam analysis G = S 3 ⋉ R 3 N +1 , the exponential map and the tangent operator are diagonal block matrices, where each block is of the form ( 18) and ( 20) respectively.We would like to remark that due to the form of the expression of the operators, one has to consider alternatives for the singularities at Ω ≈ 0. In [12, Appendix B], a deep investigation is performed to obtain the expression needed to have the numerical method without loss of accuracy.In Table 2, we present a limited output of this study, focusing on the functions that are needed in the group S 3 ⋉ R 3 .

Case study
As a proof of concept we refer to the classical benchmark of the heavy top.We use the model as described in [3], where the orientation of the top is parametrized by rotation matrices in SO(3).
The heavy top has a fixed point in the origin of the inertial frame and it is free to rotate about that point.The motion under the influence of gravity is described via the position of the center of mass in the inertial frame x ∈ R 3 and the orientation of the body using a rotation matrix R ∈ SO(3), i.e. the variable q in Eq. ( 5) is q(t) = (R(t), x(t)).The data, omitting physical units, include the center of mass in the body-attached frame X = (0, 1, 0) ⊤ , the mass m = 15.0, the inertia tensor w.r.t. the center of mass J = diag(0.234375,0.46875, 0.234375) and the constant gravitational acceleration vector γ = (0, 0, −9.81) ⊤ .The constraints at position level refer to the rotation about a fixed point, which can be modelled with the constraint Φ(q) = R ⊤ x − X.As initial configuration, we set the rotation R(0) = I 3 and the angular velocity Ω(0) = (0, 150, −4.61538) ⊤ with initial position and translational velocity consistent with the holonomic constraints at position level and corresponding hidden constraints at velocity level.In Figure 1 and Figure 2, the proposed numerical methods are applied to the heavy top benchmark.In the plots, one may observe that the expected order of convergence is preserved for the test cases in the Lie group setting.Moreover the same convergence applies also to the Lagrange multipliers λ, right plot in Figure 1 and 2. The experiment solves the system in its index-2 formulation, where the constraints are at velocity level.Known limitation is the risk of linear growth of the residual in the holonomic constraints (drift-off ) and it appears to be group dependent.One may observe in Figure 3 that the choice of the Lie group determines the presence of the drift-off.When using a direct product SO(3) × R 3 , the constraints at position level are not completely satisfied and the numerical solution diverges systematically from the manifold of constraints.On the contrary, for the semidirect product SE(3) := SO(3) ⋉ R 3 the solution maintains the constraints both at position and velocity levels.That is in line with a detailed theoretical analysis  Fig. 2: Heavy top modelled in SE(3).Convergence study for the variable at position level q (left plot) and for the Lagrange multipliers λ (right plot).
for other Lie group integrators [3, Section 3.6].The drift-off effect will be further investigated in future research.It does not appear in the specific applications to beam dynamics we are going to show later in the paper.In a direct comparison with the Lie group DAE integrators that rely on the index-1 formulation of the equations of motion [28,30], the reduced risk of drift in the index-2 formulation ( 5) is a clear benefit of the proposed method ( 12)-( 14), see also [11, Theorem VII.

Error controlled variable time step size
In Section 2, we introduced a modified half-explicit Runge-Kutta method for Lie group settings.In the current section, we are going to expand the study to variable time step sizes.The advantages of using an error controlled variable step size are the reduction of the computational time obtaining at the same time a control over the error.In particular, when solving the equations with fixed time steps, we do not know how well the simulation is performing while it is running.For variable step sizes, we have to estimate the local error to establish a new time step size.
The theory behind the approach used here can be found in [10].

Methodology
The main challenge is the estimation of the local error for the position variable.Since we are working in a Lie group setting, the position variable, indicated before by q(t), cannot be used in sums and differences as an element of a linear space.The decision we make is to estimate the local error of the local parametrization θ n (t).
In this way, we are able to use the established theory of embedded Runge-Kutta formulas [10].To implement it, we need to introduce a new variable, which does not lie in the Lie group, and for which we can perform operations of sums and difference: y = (θ ⊤ , v ⊤ ).Here θ and v, as before, are the local parametrisation of the Lie group and the velocities respectively.
Local error estimate Firstly, we introduce the embedded Runge-Kutta formulas.At the same time step, we are going to evaluate two numerical solutions obtained with two Runge-Kutta methods with the same stage vectors Q ni , Θ ni , V ni but different weights b i , such that we obtain a first numerical approximation y 1 of order p and a second numerical approximation ŷ1 of order p. Usually, we have p = p − 1 or p = p + 1.The approximation y 1 is used to continue the integration, while the approximation ŷ1 is used to evaluate the estimate of the local error.Given user defined tolerances, the error indicator at time t = t n is where m = 6(N + 1) is the number of components of the solution vectors, Atol i , Rtol i , (i = 1, . . ., m) are vector valued user defined tolerances.
New time step size After the local error estimate is evaluated, the value of the error indicator err is compared with 1.If err ≤ 1, the computed step is accepted and the computation will continue for the time step t n+2 = t n+1 + h new .On the contrary, if err > 1, the computed step is rejected and a new computation for the time step t n → t n+1 will be performed using a new time step size, i.e. t n+1 = t n + h new .As one may observe, in both cases we need to compute a new value for the time step size.We are going to evaluate the new time step size in terms of the local error estimate, but we are going to use some multiplicative factors to avoid a too large change in the step size from one step to the next.In particular, given the previous time step size h, the new one will be [10] where µ = min (p, p).We observe that for err > 1, (1/err) < 1 and the time step size will surely diminish, while for err ≤ 1 the new time step size could still be smaller than the previous one due to the role of the multiplicative factor f ac, which is set to be smaller than one.Typical values of these factors are f ac = 0.8, f acmin = 0.2, and f acmax = 5.0, see [10].
4 Applications to geometrically exact beam model

Roll-up and Flying spaghetti
The numerical experiments in the present section are performed for Cosserat beams with internal constraints.In detail, the geometrically exact model of a rod describes the rod as a curve, the centreline C(s), parametrised by the curvilinear abscissa s, along which we position the cross sections with a specific orientation dictated by the rotation matrix or the unit quaternion [16,14].In Figure 4a, we may observe the model previously stated and in Figure 4b the discretisation in space is depicted.The discretisation of the beam is motivated by the staggered grid in [16], Fig. 4: Cosserat rod but only from a computational point of view.In our practical applications, both the quaternions, indicating the orientation of the cross section, and the position coordinates are set on the nodes of the discretisation.The methodology goes back to the work of Hante [12] and has been introduced in [14].We evaluate in the midpoints of the edges the spatial derivatives of the configuration variables where N is the total number of discretisation edges in which the beam is divided, s ∈ [0, L] is the curvilinear abscissa, q i = (p i , x i ) ∈ S 3 ⋉ R 3 and log denotes the inverse of the exponential map with image in R 6(N +1) .Finally, the trapezoidal and midpoint rule are used respectively to discretise the kinetic and the potential energy.The validation of the numerical method through numerical experiments is performed with two test cases.Both tests involve a constrained Cosserat rod as presented in [16], see also [14].The numerical experiments are called roll-up maneuver [9] and flying spaghetti [26], and are fully described in [12].The physical parameters involve the use of material parameters such as Young's modulus and geometric parameters as the area of the cross section or inertia values in different directions.In Table 3, we summarize these parameters, the initial conditions and the boundary conditions for both experiments.We define d and a time dependent input g(t) as The roll-up maneuver experiment considers a rod, fixed at one of the ends and subject to a moment on the other one, see Figure 5a.The flying spaghetti, already introduced by [26], on the contrary does not have any Dirichlet boundary conditions, having both ends free, but it has both a force and a momentum applied to one of its ends, see Figure 5b.For both test cases, the constraints arise from the use of the Kirchhoff model [16].In particular, the model avoids shear deformations of the cross sections, which remain perpendicular to the centreline.The constraint at position level is  where Γ is the material strain vector.Since we require (25), the first two components of Γ need to vanish.
For the experiments, we apply the half-explicit Runge-Kutta Lie group methods described in the sections above with coefficients that allow the methods to perform up to order five, see Table 1.In Figure 6, we observe the expected order of convergence for the method applied to the roll-up maneuver test case with N = 8 edges and for the flying spaghetti test case with N = 16 edges.We may notice that the convergence properties of the method in the classical setting are maintained in the Lie group setting.We now consider the solution of the two problems when using an algorithm for variable time step size.Here we firstly observe the behaviour of h, the simulation is then performed two more times, once with time step size 2h and once with time step size h/2.One may notice that the local error estimate for the solution with time step size h is about 2 4+1 times larger than the local error for the solution with time step size h/2, similarly when computing the quotient between the local error for 2h and h.In fact, the method we use to estimate the local error is of order four, and we know that where p is the order of the method.We implemented the fifth order method with variable step size and in Figure 8 we observe the results for both the roll-up maneuver and the flying spaghetti example.
In Figure 8, the simulations run in the time interval t ∈ [0, 5] and we notice that after a first time interval in which the time step size is comparable in size between the two test experiments, for t > 1s the step size increases substantially for the roll-up maneuver and remain bounded for the flying spaghetti.The difference in the size of the time step could be explained by the nature of the test case itself: In the flying spaghetti there is no equilibrium state that could be reached, while for the roll-up maneuver the closed ring position is an equilibrium point of the motion.The increase in the time step size could be then associated with the reached equilibrium status.Other observations may regard the behaviour in the change of time step size in function of the variable err.When err decreases we notice a spike in the step size sequence, since the error reaches lower values.Figure 9 shows the difference in step size history when increasing the bending stiffness for the roll-up maneuver test setup.As expected, the step size increases ).We observe that a smaller value of the stiffness will make the numerical integration smoother, in the sense that the local error will have less oscillations and the time step size will have less variations.We may notice that when increasing the stiffness, spikes in the local error will raise, which will require the decreasing of the time step by quite a large factor.Those differences are not unexpected, on the contrary we know that stiffer systems may require smaller time step size for their resolution if explicit integrators are used.

Conclusion
The time integrator we introduced in the present paper shows favourable properties: Having the structure of a half-explicit integrator, it does not need the use of a Newton-Raphson method to solve possible nonlinear systems of equations, and being designed to maintain the solution on the Lie group, it is advantageous for the evaluation of the solution of mechanical systems subject to large rotations.We observed that the application on the Lie group does not influence the stability behaviour and the classical orders of convergence are preserved in the corresponding Lie group integrator.
On the other hand, we extended the study to the control of the local error.The implementation of a variable time step size is not only helpful to reduce the computing time, but also to control the local error.The classical fifth order explicit Runge-Kutta method of Dormand and Prince [8] has been made available to constrained system on Lie groups.

Fig. 1 :
Fig.1: Heavy top modelled in SO(3) × R 3 .Convergence study for the variable at position level q (left plot) and for the Lagrange multipliers λ (right plot).

Fig. 8 :
Fig.8: The blue lines represent the test data for the roll-up test with N = 8 edges, the red lines for the flying spaghetti with N = 16 edges.The upper plot represents the trend in the variable err that is the estimate of the local error against the user defined tolerance Atol = 10 −8 and Rtol = 10 −6 , the lower plot is the new time step size after each accepted step.

Fig. 9 :
Fig. 9: The blue line represents the solution for the roll-up test with bending stiffness C K = 5 • 10 2 , the orange line the solution with bending stiffness C K = 5 • 10 4 .