Abstract

On the Numerical Integration of Rigid Body Nonlinear Dynamics in Presence of Parameters Singularities

Marcelo A. Trindade trindade@mec.puc-rio.br

Rubens Sampaio

Department of Mechanical Engineering-PUC

22453-900 Rio de Janeiro, RJ. Brazil

rsampaio@mec.puc-rio.br

Introduction

The representation of the three-dimensional motion of rigid bodies or multibodies is, commonly, the major problem in the study of their nonlinear dynamics. This difficulty is due to complexities introduced by the representation of large rotations that the bodies may undergo during their motion. Since the modeling of multibody systems has become of great importance in areas such as aerospace and aeronautics structures design and analysis of complex robotic manipulators, the understanding of the nonlinear dynamics of general bodies under large rotations has been object of analysis and revision (Cardona, 1989, Trindade and Sampaio, 2000). The analysis of large displacements of flexible structures, such as beams (Bathe and Bolourchi, 1979, Simo and Wong, 1991) and shells (Betsch etal, 1998, Ibrahimbegovic etal, 1995), also motivated the use of finite rotations in configuration update procedures. Some computational codes were developed to treat the modeling and simulation of large, rigid or flexible, multibody systems. It was observed that one of the major problems is to choose the parametrization of the problem, since it determines the number of degrees of freedom that must be calculated, i.e., the cost of modeling and simulation.

Several parameter systems can be used to represent finite three-dimensional rotations. Applications of finite rotations in computational mechanics are reviewed by Atluri (1995), Géradin and Rixen (1995) and Trindade and Sampaio (2000). Details on the parametrization of large rotations can be found in the articles of Argyris (1982) and Stuelpnagel (1964), and in the classical text of Goldstein (1980). In the present work two rotation parameters systems are studied, namely Euler's angles and rotation vector, each of them composed of three parameters.

Theoretically three-parameter representations are not unique and present singular points where the rotation cannot be represented (Stuelpnagel, 1964). However, to the authors’ knowledge, the effect of these singularities on the numerical integration and on the representation of the rigid body motion has not yet been presented in the open literature. In such cases, the common solution is to use another parametrization composed of four parameters, like, for example, Euler's parameters. But, the use of a four parameters system increases the number of degrees of freedom (dof) by one for each body and, also, includes one algebraic equation for each body. Consequently, the global equations of motion of the system will be increased by N (number of bodies) dof and, moreover, will be a system of differential-algebraic equations, instead of a system of ordinary differential equations (ODE), which are much more difficult to solve. Géradin (1994), Rochinha and Sampaio (2000) and Simo and Wong (1991) studied the heavy top dynamics to show that the equations should be written in terms of convected operators and to analyze instabilities and dissipation of some integrators.

It is the objective of this work to present a numerical analysis of the singularities of finite rotation representation and their effect on the numerical integration of the nonlinear dynamics of rigid bodies. This is done through the analysis of the heavy top nonlinear dynamics where the rotations are represented with the Euler's angles and the rotation vector. It is shown that the system equations can be stiff depending on initial conditions and rotation parametrization, consequently, several integration algorithms, for stiff and non-stiff problems, are considered. The case where the Euler's angles approach their singular points is analyzed in detail and it is shown that, on the contrary of what may be expected, 1) in this particular problem, the numerical integration is very efficient, even more than for any other case, and 2) in spite of the uncertainty on the Euler's angles themselves, the body motion is well represented.

The numerical analysis of a rigid body nonlinear dynamics near the singular points of the finite rotation representation and the study of the effect of this type of singularity on the numerical integration are the main originalities of this work.

Rigid Body Kinematics

Let us consider the motion of a particle P of a rigid body, represented by the position of the body at two different instants of time, as showed in Figure 1.

We attach a fixed (inertial) frame E in the space and a moving frame b to the body. The operator R(t) defines the relative orientation of the frames, such that b_i = R(t) E_i (i=1,2,3), therefore the position vector of the particle P is

where x₀ is the translation of the body center of mass 0 and R(t) represents the body rotation. The velocity vector can be simply obtained by time-derivation

Solving (1) for X and replacing it in (2), we get

where the skew-symmetric matrix associated with the vector w, that we denote the spatial angular velocity, is defined as

Hence, one may write the velocity of the point P in terms of the velocity of the center of mass, the distance relative to the center of mass and the angular velocity of the body. Let us recall the following isomorphism that relates the skew-symmetric matrix and the angular velocity w,

where ´ represents the vector product and

One may also write these relations in the inertial frame E, through the definition of a convected angular velocity W

Representation of Finite Rotations

To write the equation of motion of the body, one must represent its kinematics in terms of some chosen degrees of freedom. Therefore, the intrinsic operator R is now parametrized using two set of parameters, namely the rotation vector and the Euler's angles. The first one, defined by Y = fn represents the rotation as a unique rotation-angle of f in the direction of a unit vector n, that is, the axis of rotation. Hence, the rotation-angle is the modulus of the rotation vector f=||Y||. Then, the rotation operator can be represented by (Cardona, 1989, Trindade and Sampaio, 2000)

where is the skew-symmetric matrix corresponding to the rotation vector Y. Using (4), the angular velocity w can be written in terms of parameter Y as (for details see (Trindade and Sampaio, 2000))

defining the tangent rotation operator T(Y), in terms of the parameter Y, as

Notice that R(Y) and T(Y) are not singular, for ||Y||® 0, since it is easy to show (Trindade, 1996) that lim_{||Y||® 0}R(Y) = I and lim_{||Y||® 0}T(Y) = I, where I is the identity operator. Hence, Y=0 is not a singularity.

The Euler's angles represent the rotation as a sequence of three elementary plane rotations, one of a precession angle y in E₃ direction, denoted as R_y(y), one of a nutation angle q in E₁ direction (Figure 2), denoted as R_q(q), and one of a spin angle j in b₃ direction, denoted as R_j(j). Therefore, one may assemble these three rotations such that the global rotation operator is written as

or, in explicit form, in the basis b_i Ä b_j,

where cz= cosz and sz = sinz, (z = y,q,j). It is worthwhile to analyze the extraction of the Euler's angles from the rotation matrix R, which may be achieved using the following relations

r_ij being the ij-th element of the matrix R. From the last relations, one may observe that this inversion is not unique when q = 0,p, since in these singular points

and only the nutation angle q and the sum or subtraction of the precession and spin angles (y±j) can be determined, but not the precession angle nor the spin separately. This is the well-known inversion singularity of the Euler's angles, although, in this case, it is clear that spin and precession angles alone have no meaning and also no use.

Thereafter, using the definition of (4), one may write the expression of the angular velocity in terms of the Euler's angles derivatives vector as

where the tangent rotation operator T (y,q,j), in terms of the Euler's angles, is written, in the basis b_iÄb_j, as

Notice that det(T)= -sin q , hence for the singular case when q = 0,p, det(T)=0 so that the inverse of T is not defined.

System Modeling

The system to be analyzed, shown in Figure 2, is the heavy top with a fixed point (Sampaio and Trindade, 1998, Trindade and Sampaio, 1999). The top is characterized by a mass m and principal moments of inertia I₁, I₂ and I₃; the distance between the center of mass of the top and the fixed point (origin of the inertial frame) is L. The top is subject to the action of gravity only.

The components of the rotation vector (y₁, y₂, y₃), with respect to the frame b, are chosen as generalized coordinates and, instead of their derivatives, the components (w₁,w₂,w₃) of the angular velocity w (9) in the basis (b₁, b₂, b₃) are assumed to be the generalized velocities. Then, velocity and acceleration of the center of mass must be written in terms of the generalized velocities, leading to

The angular acceleration, being defined as the time-derivative of the angular velocity in b_iÄb_j, is simply expressed in terms of w₁, w₂, w₃ as

From (17) and (18), the inertial forces and moments are

where J is the body's moment of inertia. The gravity force F_g=-mg E₃, which is constant in the inertial frame E, can be written in the body frame b, using the relation b_i = R(t) E_i, as

Consequently, this force is time-variant in the body frame. Then, using Kane's procedure (Kane and Levinson, 1985, Trindade, 1996), the partial V_p (p=1,2,3) and angular partial velocities W_p are defined as the derivative of the linear and angular velocities, respectively. These may be written, in terms of the generalized velocities w₁, w₂, w₃ as

The generalized inertial forces are defined as the scalar product between the partial and angular partial velocities and the inertial forces and moments so that

Similarly, the generalized active forces are defined as the scalar product between the partial and angular partial velocities and the external forces and moments leading to

The sum of these inertial and external generalized forces (23) and F_p (24) leads to three nonlinear scalar equations of motion representing the dynamics of the heavy top. They may be written, in the body frame b, as

where T (Y) is defined by equation (10). One may notice that this method allows to write all inertial forces in the inertial frame, that is, in the case of several interconnected bodies, each body can be treated separately. Only forces that are fixed in the inertial frame E need to be transported to the body frame b.

Using the same procedure and the rotation and tangent rotation operators, in terms of the Euler's angles, one may write the equations of motion, using these parameters, as

where T (y,q,j) is defined by the equation (16). It is worthwhile to notice that the use of the angular velocities as variables leads to simpler equations compared to those using the derivative of the parameters.

Simulation of Heavy Top Behaviors

Although the equations of motion showed in the previous section are written in terms of the components of the spatial angular velocity w, to take advantage of the simplicity of writing inertial forces and moments in the body frame, one must write the equations in terms of the convected angular velocity W to integrate the equations (Géradin, 1994, Simo and Wong, 1991), since the body frame varies with time and one cannot sum the components of the spatial angular velocity for different times. This can be obtained using (17).

For the integration of the equations of motion, we used the MATLAB ODESuite (Shampine and Reichelt, 1997), which contains several ODE numerical integration algorithms, for stiff and non-stiff problems. Stiff problems are those where some of their variables present fast dynamics and others slow dynamics. Whereas, for non-stiff problems, all variables present the same type of dynamics. The methods for non-stiff problems are ODE45 (explicit Runge-Kutta formula of 4^th /5^th order), ODE23 (explicit Runge-Kutta formula of 2^nd /3^rd order), ODE113 (variable order formula of Adams, Bashforth and Moulton) and ADAMS (Adams predictor-corrector formula of SIMULINK); whereas the methods for stiff problems are ODE15s (implicit Numerical Differentiation Formula developed by Klopfenstein and Shampine), ODE23s (implicit method based on a modified Rosenbrock formula) and GEAR (Gear's predictor-corrector formula of SIMULINK).

The properties of the heavy top are: m = 5kg, L = 1m, I₁= I₂= 0.8kg m² and I₃= 1kg m². In the first simulation, shown in Figure 3, we consider that the top has an initial nutation angle q₀ = p/9rad and an initial spin velocity = 50rad/s, all the others initial conditions being zero; henceforth we give only the non-zero initial conditions. In this case, the two representations present the same results, but the rotation vector leads to a slower integration as we see in Table 1. This can be justified by the fact that: (i) the equations of motion using the rotation vector are more complex than the ones using Euler's angles and (ii) the components of the rotation vector present very high frequencies of oscillation, which forces a small time step in the integrator. The behavior of the rotation vector components can be observed by the behavior of its norm f in Figure 3.

In the integration of the equations written with Euler's angles, the frequencies of oscillation of the components of the convected angular velocity are much higher (fast dynamics) than the one of the Euler's angles (slow dynamics) which characterizes a stiff problem, as it can be seen in Figure 3 (nutation, angular velocity). The ratio between x and y components of angular velocity and nutation angle frequencies was evaluated through fast-fourier transform (FFT). For the present case, the dynamics of these two components of the angular velocity is 6.1250 times faster than that of the nutation angle. In addition, the magnitude of W₁ and W₂ is very small compared to that of the z component W₃ (1.75%), since the top presents a large spin velocity. That is why the x and y components of the angular velocity are less observable in Figure 3. On the contrary, the frequencies of oscillation of the components of the rotation vector (fast dynamics) are of the same order than those of the angular velocity (fast dynamics), which characterizes a non-stiff problem. Therefore, we can observe that the method ODE15s is faster than ODE45 for Euler's angles but slower for the rotation vector.

This first simulation represents the most usual behavior of the heavy top. The top presents a periodic nutation, as shown in Figure 3 (´ q), with a frequency higher than that of the precession, as we can see in Figure 3 (Position of CM), which is in agreement with the approximation given by Goldstein (1980) for the fast top.

In the second simulation, shown in Figure 4, we consider that the top has an initial nutation angle q₀ = p/9 rad, an initial spin velocity of = 50rad/s and an initial precession velocity of = -10rad/s. This simulation represents the second type of motion of the heavy top. In this case, the amplitudes of motion are greater, leading to a more complex overall motion. The top almost reaches the vertical position (q = 0). However, the integration is faster for almost all algorithms and for both parameters systems. One may notice that the initial precession velocity is opposite to the natural precession, due to the relation spin/gravity, leading to the loops observed in the trajectory of CM (Figure 4). This also leads to a higher nutation frequency, so that the ratio between x and y angular velocity components and nutation frequencies is smaller compared to that of the first case (5.0833, for the present case, versus 6.1250, for the previous one). Also, the relative amplitudes of W₁ and W₂ compared to the third component W₃ increase for this case, reaching 12.79%. Thus, they become more observable in the angular velocity plot (Figure 4).

In the third simulation, shown in Figure 5, we consider that the top has an initial nutation angle q₀ = p/9 rad, an initial spin velocity of = 50rad/s and an initial precession velocity of = 1rad/s. In this case, the simulation is similar to the first one. However, as we see in Figure 5, those initial conditions lead to a different and unusual behavior of the top. That is, the precession is almost linear in time and the nutation angle and the spin velocity are almost constant. In this case, the top is said to precess about the vertical axis (Goldstein, 1980). This may be explained by the fact that the initial precession velocity is on the sense of the natural precession and, thus, the gravity effect on the top is less noticeable. This also leads to a smaller amplitude of x and y angular velocity components compared to the z component, in particular 0.81%, which justifies the term ‘precession about the vertical axis’. The ratio between the first two angular velocity components and nutation frequencies is 6.2500, that is comparable to that of the first case.

In the fourth simulation, shown in Figure 6, we consider that the top has an initial nutation angle q₀=p/9 rad and a small initial spin velocity of = 0.01rad/s. As the initial spin velocity is very small, the top is not stabilized in the top position and it falls like a pendulum, hence approaching the singular point q = p . In this case, the integration of the equations of motion is much faster when compared to the other cases (more than four times for the Euler's angles and twenty times for the rotation vector), as one may see in Table 1. It is, however, slower when using the Euler's angles, because (i) as the determinant of T(y,q,j) approaches zero the system of ODEs becomes singular and (ii) the frequencies of the rotation vector become very small (Figure 6). This singular case will be examined in detail in the next section, although, Table 1 proves already that, for this example, the singular case does not disturb the integration and even improves it. In this case, the algorithm ODE15s was slower than ODE45, as shown in Table 1, for both Euler's angles and rotation vector because, using both parameters, the resulting problems are non-stiff. In fact, the FFT analysis of the responses for this case show that the dynamics of x and y angular velocity components is as fast as that of the nutation angle. However, for Euler's angles, the method ODE15s was much slower than ODE45 (58 seconds for ODE15s against 38 seconds for ODE45), due to the fact that for each instant the parameters approach a singularity, the algorithm must strongly reduce the time step, increasing the number of failed steps, that is, the iterations that do not satisfy the error condition and hence must be repeated. Consequently, the methods for which each step is expensive, such as ODE15s and ODE113, present slower integration. The methods ADAMS and GEAR had the maximum time step limited in order to present satisfactory graphic output.

Numerical Integration with Euler's Angles: The Singular Case

The objective of this section is to analyze the integration of the equations of motion for the special singular case (Figure 6). From the analysis of the representation of finite rotations in terms of the Euler's angles presented before in this paper, we know that this representation is theoretically singular for q = 0,p. In fact, this theoretical singularity may be explained by two different approaches. The first approach to the singularity of the representation is to consider the inverse problem, that is, from a given rotation matrix R, extract the Euler's angles. As shown in (13) and (14), this inversion is not defined when q = 0,p. In fact, we can extract the angle q and either the sum y + j or subtraction y - j, of y and j but not the two of them. If we could, then y and j would be determined. Therefore, we can expect that in a numerical extraction of these angles, some discontinuity in the angles y and j may occur. However, as the sum or the subtraction y ± j can be correctly extracted, the possible discontinuity will be the same for the two angles. The second approach is to think about the tangent operator T in terms of Euler's angles. The tangent operator T(y,q,j) relates the angular velocity w and the time derivatives of the Euler's angles ,,. To integrate the equations of motion (26), we must, in some sense, invert the matrix T and, as it can be seen in (16), the matrix T is singular for q = 0,p. Therefore, we can expect that when q approaches one of these values, the integration of the equations will be disturbed in some sense. Although these two approaches are equivalent in a theoretical sense, the numerical problems involved are different. In the case of simulation of the three-dimensional movement of a rigid body using Euler's angles parametrization, it is important to know if the possible discontinuity, that can occurs in the singular points, can affect the representation of the motion.

The case under consideration is obtained by setting a small initial spin velocity, an initial nutation angle near p and an initial nutation velocity that forces the nutation angle to pass through the singular value p. The only non vanishing initial conditions are q₀ = 3.14rad; = 0.01rad/s; = 0.5rad/s.

Although the problem was simulated with all the algorithms presented in a previous section, the results showed here corresponds to that with better performance which was, for the case presented, the Runge-Kutta formula of 4^th / 5^th order (ODE45).

As the initial spin velocity is very small, we can expect that the top movement will be similar to that of a pendulum with equivalent initial conditions. Moreover, as the initial angle and velocity are also small, we can expect the system to be almost linear. This means that the nutation angle q is expect to have an almost sinusoidal response. Also, the spin j and precession y angles are expect to have linear responses.

In the first simulation, presented in Figure 7, we considered an error tolerance of 10^-6. We can observe that the spin angle j presents a discontinuous behavior. These jumps are of magnitude p, reflecting the error in the extraction of the angles. The same jumps are observed for the precession angle. Each jump of the spin j and precession y angles is accompanied by a change in the sign of the nutation velocity q, as we see in the results. The jumps in j and y can be explained by the singularity of the tangent operator T. Indeed, from (16), the derivatives of the Euler's angles may be written in terms of the material angular velocity components (W₁, W ₂,W₃) as

so that lim_q®p = ± ¥ and lim_q®p = ± ¥, which explains the fast increase/decrease in these two angles. We observe that, when there are no discontinuities for the spin angle j, the nutation angle q presents the expected behavior, that is, an almost sinusoidal behavior. Moreover, from Figure 7, the jumps of magnitude p in the spin and precession angles are associated to the peaks in the nutation angle. This may be explained by the following analysis of the second line of last equation for j*=j ± p, where

Hence, one may notice that an increase/decrease of p in the spin angle j leads to a change in the sign of the nutation velocity q.

Next, we execute the same integration changing the error tolerance. Since the step size is variable in this method, changing the error tolerance will force the method to adjust its step size to achieve the given tolerance. As we may observe in Figures 8, 9 and 10, this change in the parameters of the method modifies also the location of the discontinuities. Moreover, there is no visible pattern in the changes. This means that, in fact, no improvement is achieved by decreasing the error tolerance. However, we observe that for all cases (Figures 7, 8, 9, 10) the energy is preserved and, moreover, the angular velocity is the same. Notice that, although the y component of the angular velocity (W₂ ) seems to be unstable in Figures 7-10, it is not. The increase of its amplitude is due to the top precession that changes continuously the plane of its ‘3D pendulum’ behavior. Hence, the top starts as a pendulum with almost only W₁ component and, then, due to precession, the ‘pendulum’ plane rotates about the z-axis leading to a reduction of W₁ component and corresponding increase of W₂ component. Consequently, when the pendulum plane has rotated of an angle p/2, the x angular velocity component vanishes and the y one reaches its maximum. The periodicity of this phenomena is evidently related to the frequency of precession, which in the present case is very small compared to the nutation frequency (0.73%).

It is interesting to see if the rotation matrix is well represented in all cases, hence its elements were evaluated for each instant of time in terms of the Euler's angles, using (12). However, the dynamics of each element of the rotation matrix is not easy to understand. Therefore, the position of the center of mass is evaluated from the rotation matrix since P(t)=R(t)[0 0 L]^T . Figures 7, 8, 9, 10 show that the position of the center of mass is the same for all cases, meaning that the discontinuities in the Euler's angles do not affect the evaluation of the rotation matrix and, consequently, of the position of the center of mass. Therefore, despite the singularities in the representation of finite rotations in terms of the Euler's angles, one gets the correct representation of the motion.

Conclusions

A comparison between the Euler's angles and the rotation vector to represent finite rotations in nonlinear dynamics of rigid bodies has been presented, through the analysis of the heavy top simulation. All classes of possible motions of the heavy top were shown (Goldstein, 1980). Numerical results show that the parameters system used to represent finite rotations plays a major role in the numerical integration of the equations of motion.

It was shown that Euler's angles lead, generally, to faster integration than the rotation vector. The only exception was the Euler's angles singular case, although this case yields the fastest integration for both parameter systems. It was also shown that the rotation vector leads to complex equations of motion and, consequently, to slower integration in any case. Also, the results show that, although the singularities of the Euler's angles disallow the correct representation of the angles themselves, the rigid body motion can always be represented. The methodology presented here has been used to simulate three-dimensional multibody systems (Trindade, 1996), where the importance of the representation of finite rotations increases, since the greater the size of the system the greater is the cost of simulating it. Moreover, in a complex multibody system, it is very difficult to know if and when any of the bodies will approach a singularity. One may conclude that if we are interested in the representation of the rigid body motion and the knowledge of the parameters themselves is less important, the Euler's angles is always a better option, even near singularities.

Manuscript received: March 2000. Technical Editor: Átila P. S. Freire.

Publication Dates

History