A modified HHT method for the numerical simulation of rigid body rotations with Euler parameters

In multibody dynamics, the Euler parameters are often used for the numerical simulation of rigid body rotations because they lead to a relatively simple form of the rotation matrix which avoids the evaluation of trigonometric functions and can thus save computational time. The Newmark method and the closely related Hilber–Hughes–Taylor (HHT) method are widely employed for solving the equations of motion of mechanical systems. They can also be applied to constrained systems described by differential algebraic equations. However, in the classical versions, the use of these integration schemes have a very unfavorable impact on the Euler parameter description of rotational motions. In this paper, we show analytically that the angular velocity for a rotation about a single axis under a constant moment will not increase linearly but grows slower. This effect, which does not appear for Euler angles, can be even observed if the numerical damping parameter α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\alpha $\end{document} in the HHT method is set to zero. To circumvent this problem without losing the advantage of Euler parameters, we present a modified HHT method which reduces the damping effect on the angular velocity significantly and eliminates it completely for α=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\alpha =0$\end{document}.


Introduction
In the field of multibody dynamics, the system of equations of motion is highly complex due to the mixed set of differential and algebraic equations. This kind of systems present particular characteristics complicating their numerical treatment. This paper discusses on two mathematical aspects in order to solve also highly complex multibody problems arising in practice. One aspect concentrates on the rotational parametrization in order to get a simplified structure of the equations of motion and the second aspect focuses on the numerical B W. Steiner wolfgang.steiner@fh-wels.at 1 University of Applied Sciences Upper Austria, Josef Ressel Center for Advanced Multibody Dynamics, Stelzhamerstraße 23, Wels 4600, Austria time integration procedure in order to solve the system of equations in reasonable time also for industrial applications.
Euler parameters as a special choice of quaternions are frequently used to define the orientation of the rigid bodies of a multibody system since the rotation matrix becomes quite simple in this case and the equations of motion are simplified. Compared to other parametrizations of SO(3), Euler parameters are singularity-free and avoid the introduction of trigonometric functions. These preferences contribute to reduce the computational time in a multibody simulation program, although an inner constraint has to be observed as four parameters e 0 , e 1 , e 2 , e 3 are used for describing the three rotational degrees of freedom of a rigid rotation. The equations of motion for rigid bodies in terms of the unit quaternions can be found, e.g., in pioneering work by Nikravesh [19], and as well for flexible bodies, e.g., in the work by Yoo and Haug [27] and in the references cited therein. Another advantage of Euler parameters is that the Jacobian matrices of these nonlinear equations are only moderately complicated compared to other parametrizations and can be implemented in a multibody simulation code by analytical formulas [22]. This property is essential for obtaining the adjoint equations of a multibody system which are introduced for solving optimization and inverse problems efficiently; see e.g. [15].
As one possibility for the time discretization of the governing differential algebraic system of equations in index-3 formulation, the well-known Newmark integration formulas can be chosen [18]. Since numerical damping cannot be introduced with the Newmark family integrators unless the order of accuracy is decreased, Hilber, Hughes and Taylor [13] proposed the so-called HHT α-method which can still maintain second order accuracy and stability. A variable damping depending on the value of the parameter α in the interval [− 1 3 , 0] is involved. The smaller the value of α the larger the artificial damping can be chosen in order to damp out high frequency oscillations that are of no interest. When choosing α = 0, there is no damping. However, if Euler parameters are used for the description of the generalized coordinates regarding the rotation, the Newmark formulas also influence the numerical behavior of the Euler parameters and therefore may introduce spurious damping with a high amount also in simple examples in which the system is conservative. An investigation on how the numerical energy dissipation changes with the step size is relevant.
Several authors have discussed on the problem of the artificial numerical damping of the HHT time integration procedure, see e.g. [9,12,17,20], but to the best knowledge of the authors, not in the case when Euler parameters are utilized.
One of the fundamental discussions on energy-momentum conserving algorithms has been introduced by Simo et al. [25]. Therein, various implicit time stepping algorithms are presented which preserves exactly the total angular momentum and the total energy. Various authors in the field of multibody dynamics used the Newmark formulas within the HHT α-method in order to get rid of high frequency oscillations that often occur in multibody applications due to involved constraints and flexible bodies. An approach proposed by Bauchau [4] deals with high frequencies in the multibody system with a self-stabilization algorithm which enforces the constraints in a multibody system in an index-2 formulation similar to the GGL method [11], but avoiding the introduction of additional Lagrange multipliers. Another approach concerning the stabilization of constrained mechanical systems has been presented by García Orden [10], in which the artificial energy introduced in the system can be controlled from various points of view. The group around Arnold, Cardona and Brüls [2,3] suggests a Lie group time integration of constrained systems and presents a stabilized index-2 formulation in terms of Euler parameters; see e.g. [1]. Terze et al. [26] propose a Lie group integration method for an index-1 formulation for constrained multibody systems by a reformulated Munthe-Kaas algorithm based on the work by Müller et al. [14] in which the importance of constraint satisfaction in the Lie group configuration is pointed out. Moreover, the group around Betsch and Steinmann discusses on the constrained integration of rigid body dynamics in [7], in which the main idea is to enforce the Euler parameter constraint on velocity level without resorting to a particular rotational parameterization, while in [6] especially a conserving numerical time integration in terms of quaternions is presented also following the ideas of rigid body dynamics in terms of directors, as e.g. in [5].
The present paper discusses the numerical damping effect of the HHT algorithm on the rotational motion of a rigid body if Euler parameters are used. Therefore, we consider the simple case of a rotation about a single axis under the influence of a constant moment which should lead to a linearly increasing angular velocity. We first derive an explicit formula for the numerical damping resulting from one time integration step, which shows the influence of the time step size and of the process parameter α. We will see that the dissipative effect exists even if α is set to zero.
However, the dissipation for a simple rotational motion is a severe problem of the HHT method in combination with Euler parameters which does not occur for other parametrizations of SO(3). The reason is that a constant rotation is described by oscillations of the Euler parameters which are numerically damped like "physical" oscillations of flexible structures. Whereas numerical damping is welcome in the latter case, in particular for high frequency oscillations, it is very unfavorable in the first case, e.g. for analyzing applications in rotordynamics. But since both, the HHT method and Euler parameters, have also significant advantages in multibody dynamics a solution of the outlined incompatibility would be preferable.
Hence, in the second part of this paper, we develop a modified HHT time integration scheme for the use of Euler parameters. This method is a compromise between the pure Newmark method and its canonical extension to SO(3), which was first presented in [24] and applies the Newmark formulas completely on the Lie algebra of SO(3). Since the latter approach requires the computation of the exponential map to update the rotation during each time increment, the original computational efficiency of Euler parameters is lost. Moreover, the determination of Jacobian matrices is cumbersome in this case. The key idea of the new method proposed below is to use the modified Newmark formula on velocity level only and to add the velocity constraint for Euler parameters. We show that the numerical damping effect is significantly improved and vanishes if α = 0 with this approach.

Euler parameters
We start with a brief review of the Euler parameter description for rotational motions, which are e.g. derived in Nikravesh [19] or Shabana [21]. The Euler parameters e = (e 0 , e 1 , e 2 , e 3 ) T are defined by an axis of rotation v = (v x , v y , v z ) T and by the corresponding rotation angle ϕ in the following way: A rotation matrix A ∈ SO(3), which defines the orientation of a body fixed reference frame in space, can be written in the form Since the rotation vector v in Eq. (1) is a unit vector, Euler parameters must satisfy the constraint equation The respective constraint condition on velocity level reads 2e 0ė0 + 2e 1ė1 + 2e 2ė2 + 2e 3ė3 = 2e Tė = 0.
For the subsequent computations the following properties of L and e, which can easily be verified, are needed: where I 3 and I 4 denote the three and four dimensional unity matrices. The rotational equations of motion for an unconstrained rigid body parameterized in the three dimensional space, e.g. by using Euler angles, can be written as Jω +ωJω = m (8) in which J = diag(I x , I y , I z ) represents the inertia tensor in the principal axes frame for the sake of simplicity. The vector m denotes the sum of all external moments acting on the body. The matrixω describes a skew-symmetric matrix associated with the angular velocity vector ω. In terms of Euler parameters ω andω can be expressed by Eq. (3). Moreover, it can be shown thatω = 2LL T . If Eq. (8) is premultiplied with 2L T and enhanced by the constraint force 1 2 C T e λ = eλ resulting from Eq. (5), the equations of rotational motion in terms of Euler parameters are obtained: A detailed derivation of Eq. (9) is presented in [23]. Together with the constraint equations Eq. (5) it forms a set of differential algebraic equations for e(t) and the Lagrange multiplier λ(t) if the external moment m is known. The Lagrange multiplier λ(t) turns out to be zero, which can be seen after a premultiplication of Eq. (9) with e T . Since Le = 0, this results in e T eλ = 0 or λ = 0. But note that this holds only for the description of the rotational motion given by Eq. (9). If the constraint equation is used to modify Eq. (9), the Lagrange multiplier can also attain non-zero values, see [23] for details.

The dynamics of Euler parameters for a rotation about a single axis
To analyze the impact of the HHT method on the dynamics of Euler parameters we consider the special case where a body rotates only about the first principal axis and is driven by a constant moment m = (M x , 0, 0) T . The rotation vector is then given by v = (1, 0, 0) T and the Euler parameters defined by Eq. (1) reduce to The L-matrix from Eq. (4) becomes The angular velocity and the angular acceleration according to Eq. (3) are given by and ω y = ω z = 0. Equation (9) for the general three dimensional case simplifies to From Eq. (6) one could readily conclude that the second terms on the left sides vanish for the exact solution. Moreover, since also λ = 0 as shown above, both Eqs. (14) and (15) are equivalent to Using Eq. (13), we see thatφ = M x /I x and, for zero initial conditions, ϕ = (M x t 2 )/(2I x ). Inserting in Eq. (10) we finally arrive at the expected result However, the properties λ = 0 and e 0ė0 + e 1ė1 = 0 which we have used for the derivation of Eq. (16), are eventually not satisfied exactly, if we apply a numerical solution procedure. Hence, we do not further simplify Eqs. (14) and (15) at this point and discuss what happens if the HHT method is employed for solving these equations numerically.

The classical HHT α-method
As one possibility for the time discretization of second order differential equations, the classical Newmark integration formulas [18] are widely employed. They relate positions, velocities and accelerations at the time step t n with the respective quantities at a later time step t n+1 = t n + h. For Euler parameters, these formulas read If e n ,ė n andë n are already known, e n+1 andė n+1 are expressed in terms ofë n+1 . In the original Newmark method, the lacking equation forë n+1 is simply the equation of motion at t n+1 . However, the Hilber, Hughes and Taylor (HHT) α-method [13] uses a weighting of the force terms in the equations of motion at the time instances t n and t n+1 . The weighting factor α, usually defined in the interval [− 1 3 , 0], introduces numerical damping into the integration scheme. It is also used to compute the parameters β in Eq. (17) and γ in Eq. (18) as When applied to Eq. (9), the HHT method claims Equations (17), (18), (20) and the constraint e T n+1 e n+1 − 1 = 0 is a set of nonlinear equations determining e n+1 ,ė n+1 ,ë n+1 and λ n+1 which, in general, must be solved numerically by using Newton's method.

The effect of the HHT method and the Euler parametrization on a simple rotation
To study the effect of the HHT method when it is used in combination with Euler parameters we consider again a single rotation about the x-axis and compute one time integration step explicitly. In this case, we can solve Eqs. (14) and (15) and set e 2 = e 3 = 0. The Newmark formulas for e 0 and e 1 read For simplicity, we characterize here the state variables at time step t n (denoted as old state) with a prime whereas the respective variables at t n+1 (denoted as new state) are not specifically marked. For the old state, we assume ϕ = 0 as the rotation angle andφ = ω as the angular velocity. Then, from Eq. (10), we obtain As Eqs. (14) and (15) should be satisfied at t n , we may conclude that After inserting Eqs. (22) and (23) in Eq. (21) we obtain the following equations for the new state: The remaining equations for the new state are obtained from specializing Eq. (20) for the rotation about the x-axis. By using also Eqs. (22) and (23) for the old state, we get Equations (24), (25), (26) and the constraint equation e 2 0 + e 2 1 = 1 determine the new state variables e 0 , e 1 ,ė 0 ,ė 1 ,ë 0 ,ë 1 and λ. The Lagrange multiplier λ can be readily eliminated. For that, we multiply Eq. (25) with e 1 , Eq. (26) with −e 0 and add both equations. By using also e 2 0 + e 2 1 = 1 and dividing by 4I x , we obtain After insertion of the Newmark formulas given in Eq. (24), this equation contains only the unknownsë 0 andë 1 . Moreover, after inserting Eq. (24) also in the constraint equation e 2 0 + e 2 1 = 1 another equation forë 0 andë 1 is obtained. The latter one is a quadratic equation from whichë 0 can be expressed byë 1 and also inserted in Eq. (27). Finally, we end up with a single fourth order polynomial equation forë 1 , which is hard to solve analytically. However, if we introduce a series expansion of the equation up to third order terms in the time step size h, the quartic equation reduces to a quadratic equation inë 1 . From its solution the remaining state variables e 0 , e 1 ,ė 0 ,ė 1 ,ë 0 can be computed just by insertion. We are especially interested in the angular velocity associated to the new state which is defined by Eq. (12). After a further series expansion in terms of h the new angular velocity reads where Eqs. (19) have also been taken into account. We can now see that the terms up to first order in h accelerate the rotation correctly according to the physical law I xω = M x . However, the succeeding terms cause numerical damping if α ∈ [− 1 3 , 0]. The damping in the leading term is proportional to ω 3 . Hence, the effect increases significantly with increasing angular velocity, but it can be eliminated by setting α = 0. However, numerical damping is still present even for α = 0, since the third order term in h is −M x /I x ω 2 h 3 /4 if α is zero. We conclude that the application of the HHT method is no good choice for simulating systems with spinning bodies if the rotations are described with Euler parameters.

A modified HHT-scheme for Euler parameters
To overcome the problem of spinning bodies outlined in the previous section, we suggest a modification of the Newmark formulas for the simulation of rotational motions which are described in Euler parameters. Modified Newmark formulas for rotational motions were already presented by Simo and Vu-Quoc in 1988 [24]. They basically introduced the rotation vector θ n transforming the rotation matrix A n at time t n into the rotation matrix A n+1 at time t n+1 . With the skew-symmetric tensorθ n associated to θ n and with the exponential map for SO (3), the relation between A n and A n+1 can be written in the following form: Herein,θ n is interpreted as an element of the Lie algebra of SO(3). Physically speaking, the vector θ n is a rotation vector in the inertial frame which is represented in the body fixed frame by the vector A T n θ n . Since the time derivatives of the body fixed rotation vector are equivalent to the body fixed angular velocity vector ω and the angular acceleration vectorω, Simo and Vu-Quoc applied the Newmark formulas to discretize the relation between θ , ω andω at the time instances t n and t n+1 in the following form: Equations (30), (31) and (29) can be interpreted as a canonical extension of the classical Newmark formulas to the structure of the rotation group SO(3). In principle, they are applicable to any parametrization of the rotation matrix. However, the practical use of the modified Newmark formulas requires the computation of the exponential map and, for solving the nonlinear equation for the new state at t n+1 , its linearization with respect to the rotation parameters.
As an alternative, which is simpler to implement, we propose to replace in the classical Newmark formulas only Eq. (18) with Eq. (31) and to still use Eq. (17) instead of Eq. (30). By inserting the relations ω = 2Lė andω = 2Lë in Eq. (31), we obtain (after dividing by 2) Since these are only three equations whereas Eq. (18) contains four equations for the four Euler parameters, we add Eq. (6) at t n+1 , i.e. e T n+1ė n+1 = 0, as the fourth equation. Hence, the modified Newmark formula replacing Eq. (18) reads The 4 × 4-matrix on the left side is orthogonal and can be inverted easily by transposing. Hence, the latter equation can be resolved forė n+1 , yieldinġ Using the general relation L T L = I 4 − ee T and adding Eq. (17), our modified Newmark formulas for the Euler parameters are summarized as follows: These formulas replace Eqs. (17) and (18) and guarantee that the angular velocity vector stays in the Lie algebra of SO (3). Equations (20), (33), (34) and the constraint e T n+1 e n+1 = 1 form a set of nonlinear equations determining e n+1 ,ė n+1 ,ë n+1 and λ n+1 . Note that, unlike Eq. (18), Eq. (34) is now nonlinear in the Euler parameters e n+1 . But like in the classical Newmark method, e n+1 can be expressed byë n+1 from Eq. (33) and after inserting in Eq. (34),ė n+1 becomes a function ofë n+1 , too. Hence, Eq. (20) and the constraint equation e T n+1 e n+1 = 1 finally contain onlyë n+1 and λ n+1 as unknowns. Like in the classical HHT method, the solution of these equations must be computed with Newton's method. In Appendix B it is shown, how the Jacobian matrix changes when the modified Newmark formulas are applied to the rotational part of the equations of motion.

The simple rotation with the modified HHT method
We now show that the HHT method with the modified Newmark formulas preserves the angular momentum of a spinning body if α = 0. Like in Sect. 5, we consider a rotation about the x-axis, where e 2 = e 3 = 0. The old state at time t n is again defined by Eqs. (22) and (23) The new state at t n+1 is described by

Inserting in Eqs. (33) and (34) yields the relations
We first multiply Eq. (35) with e 1 , Eq. (36) with −e 0 and add both equations. This yields where we have already used Eq. (39). Equation (40) and the constraint equation e 2 0 + e 2 1 = 1 can be solved for e 0 and e 1 . The series expansion of e 0 in terms of h up to third order is given by (The second solution for e 0 has the constant part −1 and can be excluded, because e 0 must tend to +1 for h → 0 as we had e 0 = 1 for the old state.) Finally, we are again interested in the angular velocity, given by Eq. (12), which is associated to the new state. By using Eqs. (37) and (38) we obtain ω = 2(−e 1ė0 + e 0ė1 ) With e 2 0 + e 2 1 = 1 this simplifies to since some terms cancel, and by using Eq. (39) we end up with It can now be seen that our modified HHT method yields the physically exact solution ω = ω + hM x /I x if the damping parameter α is zero. In this case, the numerical integrator conserves the angular momentum. Moreover, if the damping parameter α is negative, the term where Eq. (41) has been inserted for e 0 , introduces a small numerical damping, which is of the order h 3 and proportional to ω 2 . Recall that the numerical damping of the classical HHT method which we have discussed in Sect. 5 is of the order h 2 and the leading term is proportional to ω 3 . Hence, the numerical properties of the proposed new integrator are significantly better than the properties of the classical HHT method.

The local integration error of the modified HHT method
The objective of the present subsection is to make a general statement on the local integration error of the proposed integration method. In particular, we want to demonstrate that the error at the end of one integration step with respect to the time step size h has order 3 for the positions and order 2 for the velocities, as shown in [16] for the classical HHT method.
To obtain a series expansion of the local integration error we compute the differences of the positions e n+1 and the velocitiesė n+1 after an integration step and the Taylor series expansionẽ Therefore, the Newmark integration formula of Eq. (33) is rewritten in the equivalent form where the unknown x represents the change in the acceleration from time t n to t n+1 . In a next step, also Eq. (34) is rewritten in terms of x which is now much more complicated since x enters also e n+1 and the matrix L n+1 via e n+1 which has to be substituted from Eq. (45). However, since L is a linear function of the Euler parameter vector e due to Eq. (4), the matrix L n+1 can be computed by simply inserting Eq. (45) in Eq. (4). For our purposes, it is sufficient to consider only the linear terms in h given by After inserting in Eq. (32) we obtain the modified Newmark formula forė n+1 up to first order terms in h aṡ By using the Euler parameter identities Eqs. (6) and (7) and assuming a perfectly consistent state at time t n , it can be shown (see Appendix A) that L T n L nėn =ė n andL T n L nėn + L T n L nën =ë n .
Hence, we end up with the simpler expressioṅ If we expand the acceleration increment in the form we can see that it is sufficient to determine only the term x 0 to obtain e n+1 from Eq. (45) up to second order andė n+1 from Eq. (47) up to first order in h and, hence, to show that the position and velocity errors are of order h 3 and h 2 , respectively. The acceleration increment x is obtained from inserting Eqs. (45), (46) and (47) in Eq. (20), where the Lagrange multiplier λ n+1 = λ n + y is an additional unknown. The resulting equation for x and y has the general form F(x, y, h) = 0 if we consider the external moments m n and m n+1 as known. Inserting Eq. (48) and a similar series expansion for the Lagrange multiplier increment yields Since this equation must also hold for h → 0, we may conclude that Assuming also here that the old state at t n satisfies the constraint equation and its derivatives e T n e n − 1 = 0, e T nė n = 0,ė T nė n + e T në n = 0, exactly, the constraint equations reduces to Using our series expansion Eq. (48), we obtain from which the additional equation can be derived for x 0 . Collecting Eqs. (52) and (53) in matrix notation yields for the zero order terms x 0 and y 0 . Since the matrix on the left side is nonsingular, the solution is x 0 = 0 and y 0 = 0. Hence, the acceleration increment has the order x = O(h) and from Eqs. (45) and (47) we see that Finally, the differences between one time step of the modified HHT method in Eqs. (43) and (44), and the Taylor series given by Eqs. (55) and (56) is described by We note that it has been shown for general Lie group integrators that the order of convergence can even be three for positions and two for velocities [2]. However, in our context, we aim for an easy to implement integrator which has at least equal convergence properties as the classical HHT method but avoids the artificial damping effect which occurs when Euler parameters are used. As the analysis of a simple rotation in Sect. 7 demonstrates, the order of the integration error may even be higher if the damping parameter α is zero. Nevertheless, in practice, the introduction of numerical damping is unavoidable in complex industrial examples because otherwise the integrator usually fails to converge.

The use of the modified HHT method in complex multibody systems
In multibody systems, Eq. (9) is embedded in the differential algebraic equations of motion with the general form where M(q) is the mass matrix, C(q, t) the vector of constraint equations and C T q (q) its Jacobian with λ as the vector of Lagrange multipliers. Q(q, q, t) denotes the vector of gyroscopic and generalized forces. In the globally redundant description, the vector q ∈ R N contains translational, rotational and flexible components for every body of the system. The Euler parameters e i of body i are components of q representing its rotational degrees of freedom.
In the classical HHT method, the equations of motions are discretized by the formulas When using the modified HHT method, the components of q with the Euler parameters have to be replaced by Eqs. (33) and (34) whereas the remaining components of q can still be discretized with Eqs. (59) and (60). Moreover, in the HHT algorithm, the acceleration terms in the equations of motion are evaluated at t n+1 and all other terms are weighted in the interval [t n , t n+1 ] by setting Here, all terms with subscript (. . . ) n+1 are evaluated at time t n+1 and terms with subscript (. . . ) n are evaluated at time t n . After expressing q n+1 andq n+1 by Eqs. (59) and (60), and by Eqs. (33) and (34) for Euler parameters, respectively, the following equations have to be solved forq n+1 and λ: C(q n+1 , t n+1 ) =: C n+1 .
For that purpose, a Newton-Raphson scheme is commonly applied. The increments q (k) and λ (k) for an update stepq (k+1) n+1 =q (k) n+1 + q (k) , λ (k+1) = λ (k) + λ (k) are obtained from the linearization of the residua R n+1 and C n+1 with respect toq n+1 and λ given by The Jacobian matrix J (k) n+1 computes as where ( ) q = ∂/∂q and ( )q = ∂/∂q. The two matrices J 0,n+1 and J 1,n+1 contain the derivatives of q n+1 andq n+1 with respect toq n+1 according to Eqs. (59) and (60), and to Eqs. (33) and (34) for Euler parameters. Since Eqs. (59) and (33) are equivalent, we obtain the generally valid formula where I N is the N × N unity matrix. However, the computation of the Jacobian matrix is more sophisticated now. Let the vector of generalized coordinates be defined by where r i ∈ R 3 denote the translational degrees of freedom, e.g. the reference point coordinates of body i, e i its Euler parameters describing the axis orientation of the floating reference frame and f i ∈ R N f,i its flexible degrees of freedom. For integrating r i and f i Eq. (60) can be applied, whereas for e i Eq. (34) is used. Hence, J 1,n+1 is a block diagonal matrix of the form where I k always denotes the k-dimensional unity matrix. For computing the derivatives ∂ė n+1 /∂ë n+1 , Eq. (33) has to be first inserted in Eq. (34). The resulting formula reads where the skew-symmetric matrix is generated from the components of the vector w = L nėn + h(1 − γ )L nën which is given by the old state at t n and fixed during the Newton-Raphson iteration. The derivation of Eq. (69) is included in Appendix B. Note that if the classical Newmark formula Eq. (60) is applied also to the Euler parameters, the Jacobian matrix from Eq. (67) is simply given by J 1,n+1 = hγ I N . However, the additional computational effort for the assembly of J 1,n+1 does not lead to an increase of the overall computational time. In all our numerical examples we even observed rather lower computational times as with the classical Newmark formula which can be explained by the faster convergence of the Newton-Raphson process.

Numerical example
The software uses Euler parameters for the description of the rotational motions and can utilize the HHT method both in the classical and modified form presented in this paper.

Heavy top
As a first numerical example we consider the heavy top (cf. Brüls and Cardona [8]), which has found much interest in mechanics and serves as a benchmark problem for our modified HHT method. Figure 1 shows the schematics of the heavy top in the three dimensional space. This numerical test indicates second order convergence for the angular velocity. Figure 3 shows the velocity coordinates of the center of mass defined in the global coordinate system over the time, while in Fig. 4 the angular velocity about the x-axis in the body fixed coordinate system is visualized. In order to proof the accuracy of the modified HHT method, the energy of the heavy top example is plotted over the time; see Fig. 5. It shows the influence of the damping parameter α on the total energy dissipation.

V8-crank drive
As a second numerical example we consider a V8 crank drive which was modeled in our multibody simulation software FreeDyn. The system is composed of 17 rigid bodies, the crank shaft, 8 conrods and 8 pistons. Each rigid body has three translational degrees of freedom for the description of the position in space and for the description of the orientation  Euler parameter have been used. The bodies are connected by algebraic constraint equations such that we end up with a set of differential algebraic equations with 119 differential equations and 118 constraint equations. The mechanical system is shown in Fig. 6. The crank drive is accelerated about the x-axis by a torque on the crank shaft depicted in Fig. 7. Figure 8 shows the angular velocity of the crank shaft computed with both methods. The damping parameter has been set to α = −0.3 and h = 1.0 × 10 −4 has been chosen for the time step size. It can be seen clearly that the classical integrator does not yield a reasonable result for this example since the angular velocity does not increase after some time due to the numerical damping, which is not in accordance with the physical laws of motion. However, the modified HHT method accelerates the crank shaft nearly correct.
For the present example, the convergence behavior of the HHT method with the modified Newmark formulas is studied numerically. For this study the damping parameter α = −0.3. In Fig. 9 the asymptotic behavior of the global error in the angular velocity of the crankshaft for h → 0 is visualized in terms of the integrated error e ω := Moreover, the influence of the damping parameter α on the solution has been analyzed. Therefore, the present example has been solved with the HHT method with the modified Newmark formulas and a time step size h = 1.0 × 10 −4 . In Fig. 10 the asymptotic behavior of the global error in the angular velocity of the crankshaft for α → 0 is visualized. It can be seen that the error is a function of the damping parameter α to the power of 2.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A
For the Euler parameters, several relations like (7)