Improved quaternion-based integration scheme for rigid body motion

Rotation quaternions are frequently used for describing the orientation of non-spherical rigid bodies. Their compact representation by four numbers and disappearance of numerical problems, such as gimbal lock, are reasons for using them. We describe an improvement of a predictor–corrector direct multiplication (PCDM) method for numerically integrating the rigid body equations of motion with rotation quaternions. The method only uses quaternions to describe the orientation, so no rotation matrices are needed in the implementation. A predictor–corrector approach is used to update the quaternions each time step, such that no renormalization is needed at the end of the time step. The PCDM method suggested by Zhao and Van Wachem is improved such that forces and torques are calculated at the correct time using position and orientation information at that same time. This is achieved by using a leapfrog approach in the improved version, in which the linear and angular velocities and rotation quaternions are defined at half time steps, while whole time step information of these quantities is calculated as part of the improved integration scheme. The improved PCDM scheme is compared with the original implementation for rotational kinetic energy conservation, accuracy of object orientation and angular velocity, and rate of convergence for different time steps. With the modifications that we propose, the improved method has a true second-order rate of convergence, without the need for explicit renormalization of the quaternions. Furthermore, the method is applicable to problems with position and velocity dependent torques, while still only a single force/torque evaluation is needed per time step. For objects experiencing torque, the improved PCDM method performs better than the original method, now showing a true second-order rate of convergence, and much smaller errors in the prediction of object orientation and angular velocity while still requiring only a single torque evaluation per time step.


Introduction
Rigid body motions occur on many different scales, from the movement of molecules at the atomistic scale, to colloidal particles suspended in a fluid at the mesoscopic scale, to granular particles at the macroscopic scale.At each of these scales, ever increasing computational power allows for larger systems to be considered, which calls for numerical integration schemes that are both accurate and not too costly to be used.
On the molecular dynamics scale, the rotation and orientation of non-spherical molecules is important to determine the forces acting on them [13,7].A commonly used rigid body molecular dynamics integration method is the splitting technique [6,11], which is both symplectic and time reversible.Also on the mesoscopic scale, rigid body integration schemes employing quaternions are used [10].Here, an additional complication is that the forces and torques can be a function of linear and angular velocities.This complication also arises when describing granular flows with discrete element models (DEM) [3,4], where dissipative collisional interactions are important.
When the motion of non-spherical particles is concerned, it is common to attach two frames of reference to the object: (1) a lab frame in which the particle translates and rotates and (2) a body frame which rotates with the object.Consequently, one has to define the transformation between the lab and body frame.Different descriptions are in common use to perform this transformation [14,5].Amongst the most important are Euler angles, rotation matrices and rotation quaternions.
Euler angles use a set of three angles (φ , θ , ψ) to rotate the object in a predetermined sequence.The actual implementation uses a rotation matrix, which is a function of the three Euler angles.The Euler angle approach is straightforward and simple to implement, but suffers from gimbal lock when two of the three axes of rotation are parallel to each other.
The rotation matrix approach [2,5] uses a unitary 3 × 3 matrix to describe a rotation in three dimensions.Nine components are used to describe a rotation which means that this description uses more information than is actually needed.In addition, the matrix only describes a rotation in space, when its column vectors remain normalized and orthogonal to each other.It is computationally expensive to enforce these constraints in a dynamical simulation.
Rotation quaternions use four parameters to describe the rotation.Only when a quaternion is normalized, it represents a rotation in three dimensions.Details on the properties and algebra will be omitted and can be found in Refs.[5,9].The advantages of rotation quaternions are that no gimbal lock can occur, and the fact that the normalization condition is easily computed by taking the norm of the quaternion.
Usually rotation quaternions are used in combination with rotation matrices, i.e. the actual rotation of a vector is carried out by multiplication with the rotation matrix corresponding to the current rotation quaternion.This involves more calculations than strictly necessary and increases inaccuracies due to round-off errors.It is better to only use the rotation quaternions in the update of the rigid body orientation.Zhao and Van Wachem [17] introduced a quaternion based integration scheme that uses no rotation matrices.The method is based on a predictor-corrector scheme to update the angular velocity and rotation quaternion.The quaternion update automatically conserves the norm.Therefore no renormalization is required.However, for objects which experience a non-constant torque, we will show that the accuracy of the original scheme of Ref. [17] is not second order, as claimed by the authors, but closer to first order.
It is important to have a higher order accurate integration scheme because for particle-based simulations the highest computational costs are usually associated with the evaluation of forces and torques from pair interactions.With N the number of interacting particles in the system, the costs of force/torque evaluations scale as αN β with an exponent β clearly larger than 1, and usually also with a large prefactor α.In contrast, the other operations of a typical integration scheme, such as propagating the positions based on the velocities or propagating the velocities based on the forces, can be done independently for each particle, and the computational costs therefore scale linearly with N. As a consequence, for simulations with large N it is more important to have a more accurate integration scheme which allows for larger time steps δt (thus requiring a lower number of expensive force and torque evaluations) than to have an integration scheme that uses e.g. less memory or has fewer operations linear in N per time step.
An example of an integration scheme for DEM is that of Wang et al. [16].This type of integration scheme is common, in the sense that a straightforward Taylor series expansion is used to update the quaternion to the new time level.The disadvantage of such schemes is that the norm of the quaternion is not intrinsically conserved.Recent integration schemes for rigid body dynamics [1,12] introduce a momentum and energy conservative scheme that intrinsically conserves the quaternion norm.The method in Ref. [1] relies on the Hamilton formulation of a system of differential-algebraic equations (DEAs) in which the normalization constraint is introduced via a Lagrange multiplier, which is eliminated in the final scheme.Ref. [12] extends this by also allowing for additional constraints.
In this paper we propose an improvement of the scheme of Zhao and Van Wachem [17], such that the method retains the intrinsic quaternion norm conservation, can be used in cases where the forces and torques are a function of position, orientation, linear and angular velocity.Most importantly, we will show that the improved method is second order accurate in time, allowing for much larger time steps to be used than in the original scheme.The improved method still only uses a single force and torque evaluation per time step.In Section 2 the improved integration method is introduced.Section 3 summarizes and discusses the results of validation test cases.Finally, Section 4 gives the conclusions.

Equations of motion for a rigid body
The motion of a rigid body is governed by: The superscript b refers to the body frame, while symbols without superscipt are in the lab frame.The body has a mass m, center of mass position r, center of mass velocity v, angular velocity ω ω ω, and experiences a force F and torque τ τ τ.I b is the moment of inertia tensor in the body frame, which is constant for a rigid body.
The system in Eq. ( 1) is not complete, because the orientation of the object is not specified, only its angular velocity ω ω ω.Quaternions are used for determining the orientation of the object.In principle an evolution equation could be derived for the quaternions but instead we will make use of the following property for a quaternion q (see Ref. [17] for more details): Here n + 1 is the new time step and n the old, separated by a time interval δt.Note that the angular velocity used to propagate the quaternion is evaluated at half time step n + 1 2 , which is more accurate than using the angular velocity ω ω ω n at the start of the time step.
We would like the integration method to be generally applicable to all sorts of problems in which the force and torque will not only be a function of position and orientation, but also of linear and angular velocities.Examples include DEM with soft sphere interactions [4], and dissipative particle dynamics (DPD), where one encounters forces that are not only a function of the positions but also velocities [15].Therefore we will assume that the force F and torque τ τ τ are of the following form: For N interacting particles, we implicitly assume that the above force and torque on a particular particle depend also on the positions, velocities, orientations and angular velocities of the other particles.
To achieve a highly accurate integration scheme, one requirement is to determine the forces and torques using variables r, v, q and ω ω ω, all evaluated consistently at the same time.We will show that this can lead to a true second order accurate integration scheme.

Improved predictor-corrector direct multiplication method
The improved integration scheme is based on Ref. [17].The authors call the method a predictor-corrector direct multiplication (PCDM) method.The direct multiplication part refers to the quaternion update, which is done with a multiplication of two quaternions, as given in Eq. ( 2).The quaternion multiplication ensures that the method conserves the quaternion norm intrinsically, and is the reason that no renormalization is needed.
Our improved PCDM method uses a leapfrog approach for both the linear velocities and the angular velocities.This means that v, ω ω ω and q are defined at half time steps n + 1/2, while r is defined at whole time steps.This does not mean that the whole time step information n + 1 is not available, this information is calculated as part of the improved integration scheme.Predictions for the linear and angular velocities and rotation quaternion are calculated at appropriate time levels and later used for the correction to obtain the new time step values.Hence the method is a prediction-correction scheme.
At the old time steps n and n + 1/2 the following variables are known: τ τ τ n , and ω ω ω b n .To find the values for the variables at the new time steps n + 1 and n + 3/2 the following steps are carried out: 1. Calculate the new center of mass position r n+1 using the midpoint velocity v n+ 1 2 : This is the classical leap-frog step for particles without internal rotational degrees of freedom.2. Predict the center of mass velocity, quaternion and lab frame angular velocity at time n + 1.
(a) The prediction for the center of mass velocity is: where the prime emphasizes that it concerns a prediction of the center of mass velocity.(b) For the prediction of the quaternion at time n + 1 we will first calculate a prediction for the angular velocity at three-quarters time step in the body frame: and transform this to the angular velocity at three-quarters time step in the lab frame, using the most recent available quaternion (at n + 1 2 ): Using Eq. ( 2), we then calculate the predicted quaternion at n + 1: (c) For the prediction of the lab frame angular velocity at time n + 1 we will first predict the body frame angular velocity at time n + 1: This is transformed to the lab frame angular velocity using the new quaternion prediction at n + 1: For N interacting particles, steps 1 and 2 should be performed for all particles before proceeding to step 3.This way all variables relevant for the evaluation of forces and torques are known at the same time n + 1. 3. Calculate new force and torque in the lab frame: and determine torque and angular acceleration in the in body frame: 4. Calculate the new angular velocity in the body frame, and the new rotation quaternion: Use the new angular velocity in the body frame ω ω ω b n+ 3 2 and the new quaternion q n+ 3 2 to determine the new angular velocity in the lab frame ω ω ω n+ 3 2 : Finally, calculate the new center of mass velocity:

END
Note that our algorithm retains the favourable feature that forces and torques need to be evaluated only once per time step δt.
The main improvement is that the angular acceleration in the body frame (at the end of step 3), is now calculated at n + 1 with the torque also known at n + 1.In the scheme proposed by Zhao and Van Wachem, the angular acceleration was apparently evaluated at n + 1/2.However forces and torques were only evaluated at whole time steps (n and n + 1), meaning that old time step information was used for the torques.This half-step lagging leads to a lower accuracy than the method proposed here (close to first order instead of second order in time step), as we will show in the next section.

Results and Discussions
Four test cases were carried out with the improved PCDM integration scheme.The first test case considers energy conservation after a long simulation period.Secondly, the accuracy of the solution for the angular velocities and orientation is tested.A third case considers the rate of convergence of the method with decreasing time step.Lastly, a two-particle system is considered to show that for a more realistic multi-particle coupling the method is also accurate.

Energy conservation of a symmetric spinning top
A symmetric top with principle moments of inertia I x = I y = 1100 kg m 2 and I z = 2600 kg m 2 is given an initial angular velocity of ω ω ω 0 = (1.20,0.45, −0.6) s −1 .Because no external torque is applied, the rotation is driven by the initial angular velocity.This test case considers the energy conservation after a long simulation time of 10000 s.The time step used was equal to 1×10 −2 s.
The x, y angular velocity components of the symmetric top undergo a periodic motion as is evident from Fig. 1 (a).The z component is constant, because no external torque is applied and the top is symmetric.
The rotational kinetic energy was measured with both the original PCDM method and the new method.The error in the rotational kinetic energy is very small and the same for both methods, see Fig. 1 (b).The error is the same because no external torque acts, meaning that the improved method is essentially the same as the original method, the only difference being the shifted times at which the angular velocity is determined.Both methods show a slight increase of the error over time, which is due to the relatively large time step used.A smaller time step would produce even better energy conservation over time.

Forced rotation of a sphere by a 1D torque
A sphere marked by an arrow with initial orientation (1, 0, 0) is forced to rotate by a torque given by: The sphere has a diameter of 2 m, and a density of 1100 kg m −3 .The constant A in Eq. ( 4) is equal to 1×10 5 N m.The time step for the simulation is equal to 1×10 −4 s.This test case is equal to the third test case presented by the authors of Ref. [17].The torque will result in an accelerating rotation (see Fig. 2 (a)) of the sphere in the xz-plane.Table 1 compares the orientation of the arrow on the sphere after 1 s of simulation with the theoretical result (Theory), Table 1 Comparison of the orientation of the major axis after 1 second of simulation for different integration methods.

Orientation at time
2.15% PCDM (2.91578×10 −1 , 0, −9.56547×10 −1 ) 0.63% This work (2.93444×10 −1 , 0, −9.55976×10 −1 ) 0.0065% the scalar factor method [8,17] (Scalar factor) and the original PCDM method (PCDM).The improved method has a very small error, two orders of magnitude smaller than the original PCDM method.This is related to the fact that the torque is not evaluated at half time steps, which is needed in the original scheme.Figure 2 (b) shows that the error in the only nonzero component of the angular velocity is also small, as is the absolute error of the cosine of the angle between the initial axial orientation and the current orientation, as shown in figure 2 (c).

Forced rotation of a cylinder by a 3D torque: rate of convergence
In the third test case the rotating motion of a cylinder is considered.This case is based on test case four in Ref. [17], but not the same.The motion is forced by a non-linear torque of the following form: The cylinder with density 1100 kg m −3 has a radius of 5 × 10 −3 m and a height of 0.03 m.The error in the orientation of the cylinder is defined in the same way as in Ref. [17]: In Eq. ( 6) α error n is the error in the orientation of the cylinder at time step n.B 0 is the initial orientation and is equal to (1, 0, 0).The orientation at time step n (b n ), is obtained with the original and improved PCDM schemes.B n is the orientation found by solving the following set of kinematic equations: with τ x , τ y and τ z from equation (5).The kinematic equations are based on 3-2-1 Euler angles.The system is solved numerically in matlab with a Runge-Kutta method.With the obtained Euler angles as a function of time a rotation matrix is formed that can be used at each time step n to rotate the B 0 into B n .No analytical solution is available for the orientation when considering larger changes, i.e. all analytical approximations assume only small changes in the orientation.
To check the rate of convergence of both original and improved PCDM methods, the time step ∆t was varied between 1.0 × 10 −3 and 1.0 × 10 −8 s.The angular velocity components y and z grow rapidly as a function of time, see Fig. 3 (a).This non-linear behaviour is much more difficult to integrate correctly compared to the symmetric top or the 1D torque.The orientation error decreases with decreasing time step, as can be seen in Fig. 3 (b).Note that a decreasing time step means going to the right on the horizontal axis of the figure.Both methods converge for all time steps tested, but notice that the improved PCDM method has a second order convergence rate while the original PCDM method seems to converge slightly faster than linear.This is in agreement with the test case four in Ref. [17], where a convergence rate between linear and quadratic was reported.6) and ( 7) as a function of the log 10 of the inverse time step.The black guideline has a slope of -1 and the magenta guideline has a slope of -2.

Coupled rotation of two cylinders
The last test case considers two rotating cylinders that are in touching contact (see Fig. 4).Both cylinders are kept at their center of mass position by massless axes, such that the only motion is rotational.Cylinder 1 and 2 have an initial angular velocity of 2 s −1 and 0.4 s −1 , respectively.Their radii are 1 m for cylinder 1 and 2.4 m for cylinder 2. At the contact point between the cylinders a tangential force is acting of the following form: The tangential force in Eq. ( 9) has a friction constant η = 5.4 kg s −1 and v t,rel is the relative velocity at the contact point.At t = 0 s the cylinders are brought into touching contact, and the angular velocity will evolve over time as shown in Fig. 5 (a).For large enough times a new steady state is reached in which the first cylinder has a decreased positive angular velocity and the second cylinder has a reversed negative angular velocity.
In Fig. 5 (b) the relative error in the angular velocity of cylinder 1 is given for both original and improved PCDM integration schemes.The error for the improved method is very small for the whole simulation time, while for the original method there are larger deviations because the torque of the old time step has to be used.

Conclusions
The use of rotational quaternions to numerically integrate the equations of motion for a rigid body has several benefits.Firstly, only four numbers are required to specify the orientation of an object.More importantly, rotation quaternions allow for an accurate determination of the orientation without gimbal lock problems.The PCDM scheme has the additional benefit of not requiring renormalization of the quaternions after each time step, leading to a higher accuracy than methods that do need this renormalization.The original scheme contained an inconsistency because the angular acceleration was evaluated with forces and torques at different times.The improved PCDM integration scheme solves this by using a leapfrog approach in which linear and angular velocities and rotation quaternions are defined at half time steps.Four test cases show that the energy conservation for a torque-free body is as good as the original method, that the orientation determination is more accurate than the original PCDM scheme, and that the improved method convergence rate is second order in time, whereas the original method is in between first and second order, and actually closer to first order.For a system of coupled objects, the method also gave more accurate results compared to the original method.The higher accuracy allows for the use of larger time steps.We therefore advise to use this algorithm in simulations with large numbers of rigid molecules, rigid mesoscale assemblies or rigid granular bodies.

Fig. 1
Fig. 1 (a) Angular velocity components for the symmetric top as function of time.(b) Relative error in the rotational kinetic energy as a function of time.

Fig. 2
Fig. 2 Forced rotation of a sphere by a 1D torque.(a) Orientation as a function of time.(b) Relative error in ω y as a function of time.(c) Absolute error of the cosine of the angle between initial axial orientation and current axial orientation.

Fig. 3
Fig. 3 (a) Angular velocity components as a function of time for a forced rotation of a cylinder by a 3D torque.(b) log 10 of the error as defined in Eqs.(6) and (7) as a function of the log 10 of the inverse time step.The black guideline has a slope of -1 and the magenta guideline has a slope of -2.

Fig. 4
Fig. 4 Visualization of the two rotating cylinders.The initial angular velocities Ω 1 and Ω 2 are equal to 2.0 s −1 and 0.4 s −1 , respectively.R 1 and R 2 are the radii of the cylinders.

Fig. 5
Fig. 5 (a) Angular velocity components for cylinders 1 and 2 as a function of time.(b) Relative error in the angular velocity of the first cylinder for the original and improved PCDM integration scheme.