Accepted: August 09, 2022 | Published Online: August 11, 2022

Citation: Carrera DHZ, Weber HI (2022) Stability Analysis in Rotational Systems Using the Tumbling Coefficient. Int J Astronaut Aeronautical Eng 7:061

Copyright: © 2022 Carrera DHZ, et al. This This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

The goal of this work is the investigation of the necessary conditions for the existence of tumbling in rotational motion of rigid bodies. For this purpose, we choose a methodology to simplify the study of rotational motions with great amplitude, for example free bodies in the space, allowing an extension of the analysis to non-conservative systems. The developed method can be easily applied to orbital motion of a satellite in outer space, where the angular velocities are continuously monitored, and actuators are activated to seek for rotational stability. For any rotational motion, the projection of the angular velocity along the principal axes of inertia must be known, defining completely the initial condition of motion. In this paper coordinate systems are established in order to allow a simple analytical work on the equations of motion. As a result of the methodology used in this work it is proposed the definition of a parameter, calling it tumbling coefficient, to measure the intensity of the tumbling and also the amplitude of the motion, even passing the limits of stability in the concept of Lyapunov. In the motion of bodies in space only when this tumbling coefficient has a positive value, the no tumbling condition is guaranteed. The tumbling is inherent for a rotation around the axis of intermediate moment of inertia but can also occur by action of external impulses when the conical movement characteristic of the stable attitude is de-characterized. Another relevant result of this work associates for different initial conditions of motion the areas of stability and instability in the Magnus Triangle representation. In the study of non-conservative systems for an oblate body, a sufficient condition to investigate damped motion is given by the tumbling coefficient, and this condition is satisfied by viscous torques. This paper seeks to highlight the physical understanding of phenomena and the influence of various parameters that are important in the process.

Keywords

Attitude dynamics, Spin stabilization, Satellite, Tumbling, Nonlinear dynamics, Quaternion, Stability

Introduction

All classical literature in Mechanics spends some time describing the motion of free bodies in space: there is a good number of analytical possibilities to investigate this problem. Nevertheless, there are so many aspects to handle, that it always will be possible to have a new look of some particular situation. The general motivation of this analysis is to search the limits of stable operation of a free body spinning in space with unequal principal moments of inertia spinning in space after suffering the impact of some external perturbation.

The ideal initial condition will find the body spinning about one of the axes of principal moments of inertia. if it is the intermediate axis of inertia it will be naturally unstable, but for the other two possibilities the stability will depend on the impact itself [1]. This known behavior will allow to validate the results.

Our analysis will consider the unitary vector r_s directed along the spin axis which, for no perturbations, is aligned with the angular momentum vector H calculated for the center of mass. But a misalignment between both vectors may exist due to, for example, a remaining rotation after the satellite separation from the launcher or an external perturbation [2,3]. For a moment free condition, this means that the vector r_s starts rotating around the vector H. To visualize the motion we consider a sphere (Figure 1) divided in two hemispheres [4,5], where H points to the north. The tumbling of the satellite occurs if r_s keep changing hemispheres [6].

An appropriate choice of reference frames will make all the investigation easier [1]. Coordinates for the inertial frame (S_I) consider its z axis aligned with the angular momentum vector H and the body frame (S_B) will have its z axis directed along the r_s vector, which points according to one of the principal moments of inertia (depending on the geometry being investigated). The x and y axes of (S_B) will be directed along the other principal directions of the body. The full position of (S_I) is then established with its x axis orthogonal to the plane defined by the vectors H and r_s (Figure 1).

The initial condition at t = 0 can be established just after a positive impact along the XI axis injects momentum and displaces H to a φ₀ position from Z^B aligned with the axis of principal inertia r_s the identification of its intensity will be given by angle φ₀. But it can also be at an arbitrary instant when the attitude and the angular speed are known. We will see that for the stability analysis the initial attitude must be defined, and it will also be needed to draw trajectories in space for the body’s motion [3].

Angular momentum in inertial frame, spin vector in body frame and at t = 0 in inertial frame, the normalized inertia tensor [7] for the center of mass in principal axes body frame are:

$${}^{I}H\text{ = }\left[\begin{array}{c}0\\ 0\\ h\end{array}\right];{}^{B}r{}_s\text{ = }\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right];{}^{I}r{}_{s0}\text{ = }\left[\begin{array}{c}0\\ r_{y0}\\ r_{z0}\end{array}\right]\text{ (1)}$$

$${}^{B}J\text{ = }{I}_s\left[\begin{array}{ccc}{\mu }_1& 0& 0\\ 0& {\mu }_2& 0\\ 0& 0& 1\end{array}\right]\text{where it is considered }{\mu }_1\text{ = }=\frac{I_1}{I_s}>{\mu }_2\text{ = }\frac{I_2}{I_s}\text{ (2)}$$

Angular velocity, director cosine attitude matrix and angular momentum in body frame are the next equations (3), while an index 0 will define the variables at the initial condition:

$$\begin{array}{l}{}^B\omega \text{ = }\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_s\end{array}\right];{}^{B}T{}_{}^I\text{ = }\left[\begin{array}{ccc}{D}_{11}& D_{12}& D_{13}\\ D_{21}& D_{22}& D_{23}\\ D_{31}& D_{32}& D_{33}\end{array}\right];\text{ (3)}\\ {}^{B}H\text{ = }{}^{B}J{}^B\omega \text{ = }{I}_s\left[\begin{array}{c}{\mu }_1{\omega }_1\\ {\mu }_2{\omega }_2\\ {\omega }_s\end{array}\right]\text{ = }\left({}^{B}T{}^I\right){}^{I}H\text{ = }h\left[\begin{array}{c}{D}_{13}\\ D_{23}\\ D_{33}\end{array}\right]\end{array}$$

The deviation between spin axis and angular momentum axis defines the misalignment of the spin axis and is a measure for the impact [6]. It is obtained from the scalar product between these vectors, as shown by Equation (4).

$$\cos {\varphi }_0\text{ = }\frac{r_{so}^{\text{T}}H}{h}\text{ = }{D}_{330}\text{ (4)}$$

Using Tait-Bryan (Cardan) sequential rotations [2] an equivalent rotation matrix ^BT^I to the attitude matrix can be defined using the angles α, β e γ. According to the rotation sequence

$$\underset{(x,y,z)}{{\text{S}}_{\text{I}}}\stackrel{\alpha (x)}{\to }\underset{(x,y\text{'},z\text{'})}{\text{Q}}\stackrel{\beta \text{(}y\text{')}}{\to }\underset{(x\text{'}\text{'},y\text{'},z\text{'}\text{'})}{\text{R}}\stackrel{\gamma \text{(}z\text{'')}}{\to }\underset{(x\text{'}\text{'}\text{'},y\text{'}\text{'}\text{'},z\text{'}\text{'})}{{\text{S}}_{\text{B}}}$$

Results

$$\begin{array}{l}{}^{I}r{}_{s0}={}^I{T}_0^B{}^{B}r{}_s=\left[\begin{array}{c}\sin {\beta }_0\\ -\sin {\alpha }_0\cos {\beta }_0\\ \cos {\alpha }_0\cos {\beta }_0\end{array}\right]=\left[\begin{array}{c}{D}_{310}\\ D_{320}\\ D_{330}\end{array}\right]=\left[\begin{array}{c}0\\ r_{y0}\\ r_{z0}\end{array}\right];\text{ (5)}\\ \begin{array}{c}\sin {\beta }_0=0\to {\beta }_0=0\\ r_{y0}=-\sin {\alpha }_0\\ r_{z0}=\cos {\alpha }_0\end{array}\end{array}$$

A negative α₀ angle is necessary for the described impact.

There is a second solution for (5) to be considered if the body is tumbling and it is achieved for

$${}^{I}r{}_s\text{ = }{\left[\begin{array}{ccc}0& -r_{y0}& -r_{z0}\end{array}\right]}^T\text{ (6)}$$

and corresponds to angles $\alpha \text{ = }{\alpha }_0+\pi \text{and }\beta \text{ = 0}$.

Considering the initial condition (0 state)

$${}^{\text{I}}T{}^{\text{B}}\text{ = }\left[\begin{array}{ccc}\cos {\gamma }_0& -\sin {\gamma }_0& 0\\ \cos {\alpha }_0\sin {\gamma }_0& \cos {\alpha }_0\cos {\gamma }_0& -\sin {\alpha }_0\\ \sin {\alpha }_0\sin {\gamma }_0& \sin {\alpha }_0\cos {\gamma }_0& \cos {\alpha }_0\end{array}\right]\text{ (7)}$$

It is possible to relate the initial position of the gamma angle with the attitude matrix in (3). According to Equations (3) and (4) the initial attitude is defined in Tait-Bryan angles:

$${\beta }_0\text{ = }0;\cos {\alpha }_0\text{ = }\frac{I_s{\omega }_{s0}}{h};\tan \gamma \text{ = }\frac{D_{130}}{D_{230}}\text{ = }\frac{{\mu }_1{\omega }_{10}}{{\mu }_2{\omega }_{20}}\text{ (8)}$$

The simulations done for this problem used attitude quaternions q = [ q₁ q₂ q₃ λ ]^T motivated by the frequent crossing of stability borders [8]. The initial conditions (Equation 9) are transformed to quaternions, the procedure for relate the quaternion components of the sequential rotations (Tait-Bryan) is explained in [9] and [6]. Integration of the quaternion equation result in rotation matrices along the time, allowing the visualization of the motion.

$$\begin{array}{l}{q}_0\text{ = }L\left[\begin{array}{c}\sin {\alpha }_0\cos {\gamma }_0+\sin {\alpha }_0\\ \sin {\alpha }_0\sin {\gamma }_0\\ \cos {\alpha }_0\sin {\gamma }_0+\sin {\gamma }_0\\ 1+\cos {\gamma }_0+\cos {\alpha }_0\cos {\gamma }_0+\cos {\alpha }_0\end{array}\right]\text{ (9)}\\ L\text{ = }\frac{1}{2\sqrt{1+\cos {\gamma }_0+\cos {\alpha }_0\cos {\gamma }_0+\cos {\alpha }_0}}\end{array}$$

Normalized Equations

Beside the above-mentioned parameters, there is still the kinetic energy to be considered, which will be constant together with the angular momentum if we admit the system conservative and torque free. Considering only rotation, the kinetic energy can be related to the angular momentum using Equation (3):

$$Ec\text{ = }\frac{1}{2}{\omega }^{\text{T}}H\text{ = }\frac{I_s}{2}\left({\mu }_1{\omega }_1^2+{\mu }_2{\omega }_2^2+{\omega }_s^2\right)\text{ = }\frac{h^2}{2I_s}\left(\frac{D_{13}^2}{{\mu }_1}+\frac{D_{23}^2}{{\mu }_2}+D_{33}^2\right)\text{ (10)}$$

Which suggests us to introduce a parameter K as follows:

$$\begin{array}{l}K\text{ = }\frac{2Ec{I}_s}{h^2}\text{ = }\frac{D_{13}^2}{{\mu }_1}+\frac{D_{23}^2}{{\mu }_2}+D_{33}^2\text{ (11)}\\ K\text{ = }\frac{{\mu }_1{\omega }_1^2+{\mu }_2{\omega }_2^2+{\omega }_s^2}{{\mu }_1^2{\omega }_1^2+{\mu }_2^2{\omega }_2^2+{\omega }_s^2}\end{array}$$

K represents how the perturbation distributes between kinetic energy and the square of the angular momentum [10]: when the satellite is spinning without perturbation around the chosen axis of one of the principal moments of inertia K = 1 [11]. The imposed perturbation will increase or decrease K. It also will depend only on the geometry of the satellite and of the angular velocity components in the body frame [3], as we can see in the Equation (11).

The range of values for K depend on the geometry [6,10]. We investigate three geometries:

The oblate (flattened) body: 1 > µ₁ > µ₂ also 1 < K < 1/µ₂;

The intermediate body: µ₁ > 1 > µ₂ also 1/µ₁ < K < 1/µ₂;

The prolate (stretched) body: µ₁ > µ₂ > 1 also 1/µ₁ < K < 1.

We introduce an equivalent angular velocity ν for the normalization process [11,12] based on the square of the length of the angular momentum vector (Equation 12) and introduce non-dimensional components of the angular velocity:

$$\begin{array}{l}{h}^2\text{ = }{I}_1^2{\omega }_1^2+I_2^2{\omega }_2^2+I_s^2{\omega }_s^2\text{ = }{I}_s^2{\nu }^2\text{ (12)}\\ {\overline{\omega }}_1\text{ = }\frac{{\omega }_1}{\nu };{\overline{\omega }}_2\text{ = }\frac{{\omega }_2}{\nu };{\overline{\omega }}_s\text{ = }\frac{{\omega }_s}{\nu }\end{array}$$

In this way we get two equations for the problem, representing the conservation of kinetic energy and of angular momentum:

$${\mu }_1^2{\overline{\omega }}_1^2+{\mu }_2^2{\overline{\omega }}_2^2+{\overline{\omega }}_s^2\text{ = }1\text{ (13)}$$

$${\mu }_1{\overline{\omega }}_1^2+{\mu }_2{\overline{\omega }}_2^2+{\overline{\omega }}_s^2\text{ = }\frac{2Ec}{I_s{\nu }^2}\text{ = }K\text{ (14)}$$

If we introduce $p\text{ = }{\mu }_1{\overline{\omega }}_1$, $q\text{ = }{\mu }_2{\overline{\omega }}_2$, $r\text{ = }{\overline{\omega }}_s$ then the possible solutions are obtained by the intersection of a sphere with an ellipsoid:

$$p^2+q^2+r^2\text{ = }1;\frac{p^2}{{\mu }_{1}K}+\frac{q^2}{{\mu }_{2}K}+\frac{r^2}{K}\text{ = 1 (15)}$$

Necessary Conditions for Tumbling

Tumbling occurs when the spin vector r_s crosses from one hemisphere [11], defined according to H, to the other (Figure 1). The component r_z = 0 at this moment. Considering Equations (1), (2) and (3):

$$\begin{array}{l}{}^{I}r{}_s=\left[\begin{array}{c}{r}_x\\ r_y\\ r_z\end{array}\right]=\left[\begin{array}{c}{D}_{31}\\ D_{32}\\ D_{33}\end{array}\right];{}^{B}H=h\left[\begin{array}{c}{D}_{13}\\ D_{23}\\ D_{33}\end{array}\right]=I_s\left[\begin{array}{c}{\mu }_1{\omega }_1\\ {\mu }_2{\omega }_2\\ {\omega }_s\end{array}\right]\text{ (16)}\\ \Rightarrow r_z=0\Rightarrow {\omega }_s=D_{33}=0\end{array}$$

We see that the spin velocity becomes zero at this instant and the body is in the border condition to tumble. Also in (Equation (15)) r = 0 and only two parameters (p e q) must be considered to evaluate the start of the tumbling condition. An ellipse crosses a circle in 4 points (with same center) when a semi-major axis of ellipse is outside the circle and semi-minor is inside the circle. That means in Equation (15) that we look for an intersection between a circle and an ellipse [11]. There are two possible solutions for this intersection:

$$\begin{array}{l}\left\{{\mu }_{1}K>1\text{and}{\mu }_{2}K<1\right\}\text{ (17)}\\ \text{or}\\ \left\{{\mu }_{1}K<1\text{and}{\mu }_{2}K>1\right\}\end{array}$$

Both conditions will be attended introducing the tumbling coefficient in Equation (18). The acronym of this coefficient is derived from the word in Portuguese language "Capotamento" meaning tumbling.

$$Capt\text{ = }\left({\mu }_{1}K-1\right)\left({\mu }_{2}K-1\right)<0\text{ (18)}$$

This necessary condition will also be sufficient for the tumbling, occurring always when the tumbling coefficient is negative. A positive coefficient will be a guarantee that no tumbling can occur, and the motion is stable.

Some Simulations

The same body will be used for three different geometries (Figure 2), changing the spin axis form the principal axis of largest moment of inertia to the smallest. Therefore we start with an oblate body (largest inertia in z) µ₁ = 0.9 µ₂ = 0.5, then take the x axis for spin (it will be the new z axis) and get µ₁ = 1/0.9 = 1.11 µ₂ = 0.5/0.9 = 0.56 for the intermediate body an finally µ₁ = 2 µ₂ = 1.8 for the prolate body.

Figure 3 shows the tumbling coefficient for these cases, depending on the value of K. In this figure, it's observed that for the oblate body the Capt starts at a positive value, and remains positive for small impacts (K is close to 1), as the impact increases the Capt value decreases and becomes negative (K > 1.11), the rotational motion becomes unstable.

For the prolate body the Capt starts at positive value K = 1 (values close to 1 do not appear in the Figure 3 but they do exist), as the impact increases the K value decreases and the Capt value decreases also, for strong impact the rotational motion becomes unstable (K < 0.56).

In the case of intermediate body the Capt value is always negative. These results are consistent and expected as indicated in the introduction of this paper.

The intermediate case (which we know unstable [1,13]) can be investigated directly from the tumbling coefficient (Equation (18)) when we substitute the K value from Equation (14);

$$Capt\text{ = }\left({\mu }_1\left({\mu }_1{\overline{\omega }}_1^2+{\mu }_2{\overline{\omega }}_2^2+{\overline{\omega }}_s^2\right)-1\right)\left({\mu }_2\left({\mu }_1{\overline{\omega }}_1^2+{\mu }_2{\overline{\omega }}_2^2+{\overline{\omega }}_s^2\right)-1\right)<0\text{ (19)}$$

and with help of Equation (13) transform the inequality condition in:

$$Capt\text{ = }\left({\mu }_2{\overline{\omega }}_2^2\left({\mu }_1-{\mu }_2\right)+{\overline{\omega }}_s^2\left({\mu }_1-1\right)\right)\left({\mu }_1{\overline{\omega }}_1^2\left({\mu }_2-{\mu }_1\right)+{\overline{\omega }}_s^2\left({\mu }_2-1\right)\right)<0\text{ (20)}$$

For a body in the intermediate axis rotation condition, µ₁ > 1 > µ₂, Equations (18) and (20) will always be satisfied and tumbling is inherent.

$$Capt\text{ = }\underset{negative}{\underbrace{\underset{positive}{\underbrace{\left({\mu }_2{\overline{\omega }}_2^2\left({\mu }_1-{\mu }_2\right)+{\overline{\omega }}_s^2\left({\mu }_1-1\right)\right)}}\underset{negative}{\underbrace{\left({\mu }_1{\overline{\omega }}_1^2\left({\mu }_2-{\mu }_1\right)+{\overline{\omega }}_s^2\left({\mu }_2-1\right)\right)}}}}<0\text{ (21)}$$

Using Equation (20) it is possible, for each geometry, to draw a stability map in the plot considering the referred representation of ω₁/ω_s with ω₂/ω_s looking for the Capt = 0 line.

Figure 4 shows the stability map [12] for these cases, depending on the geometries of bodies. It is possible to verify that the stability limits of the tumbling motion, in each body, depend only on one relation between angular velocities, the stability limits depend only on ω₂/ω_s in the oblate body, and depend only on ω₁/ω_s in the prolate body; this is because, according to the geometry of the body, one component of the Equation (20) maintains a constant sign for all possible values. This fact is useful to simplify the analysis, particularly in the analysis of a non-conservative motion, as shown later.

Analyzing Figure 4 in detail, it is possible to observe the lines with the same Capt value in the stability maps [12] (Capt is constant on this lines). In the oblate body (upper frame) these lines become horizontal when passing through zero (Capt = 0), and divide two areas of blue and green colors. In the prolate body (bottom frame) these lines become vertical when passing through zero (Capt = 0), and divide two areas of orange and celestial colors. In the intermediate body (middle frame) there exists only one area, that is, the entire area has negative values of Capt, the lines show the place for one constant Capt and divide its gradient in colors.

To calculate the tumbling condition for a certain geometry it is enough to start at a certain instant of time knowing the onboard measured angular velocity which will be ω₀ = [ω₁₀, ω₂₀, ω_s0]T in the body frame (S_B). From Equation (11) we will obtain K. If the system is conservative the simulation will be simple since we use principal moments of inertia. Using Euler’s law for the body frame [14]:

$$\begin{array}{l}\frac{dH}{dt}\text{ = }\dot{H}+\omega \text{x}H\text{ = }M\text{ = }0;H\text{ = }J\omega ;\text{ (22)}\\ \Rightarrow \dot{\omega }\text{ = }{J}^{-1}H\text{x}\omega \text{ = }{J}^{-1}\left(\left(J\omega \right)\text{x}\omega \right)\end{array}$$

Substituting the J matrix in Equation (22) we obtain the differential equation for the angular velocity, which will is:

$$\left[\begin{array}{c}{\dot{\omega }}_1\\ {\dot{\omega }}_2\\ {\dot{\omega }}_s\end{array}\right]\text{ = }\left[\begin{array}{ccc}0& -{\omega }_s/{\mu }_1& {\mu }_2{\omega }_2/{\mu }_1\\ {\omega }_s/{\mu }_2& 0& -{\mu }_1{\omega }_1/{\mu }_2\\ -{\mu }_2{\omega }_2& {\mu }_1{\omega }_1& 0\end{array}\right]\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_s\end{array}\right]\text{ (23)}$$

Figure 5 present the situation of the intersection between circle and ellipse in Equation (15) with r = 0, the value of K and the motion of spin axis in the inertial frame. To obtain the attitude we use the initial quaternion given by Equation (9) with:

$$\begin{array}{l}\tan {\gamma }_0=\frac{{\mu }_1{\omega }_{10}}{{\mu }_2{\omega }_{20}};\cos {\alpha }_0=\frac{I_s{\omega }_{s0}}{h};\text{ (24)}\\ \frac{I_s^2}{h^2}=\frac{K}{{\mu }_1{\omega }_1^2+{\mu }_2{\omega }_2^2+{\omega }_s^2}\end{array}$$

And finally, the quaternion equation using the vector product formulation whit the antisymmetric matrix replacing the vector product $vx\omega$ [1], can be written as:

$$\begin{array}{l}q\text{ = }\left[\begin{array}{c}v\\ l\end{array}\right];\left[\begin{array}{c}{}^{{\text{S}}_{\text{B}}}\dot{v}\\ \dot{l}\end{array}\right]\text{ = }\frac{1}{2}\left[\begin{array}{c}\left(lE+\tilde{v}\right){\omega }_{}\\ -v^{\text{T}}{\omega }_{}\end{array}\right]\text{ (25)}\\ \Rightarrow {}^{\text{I}}T{}^{\text{B}}\text{ = }E+2l\tilde{v}+2{\tilde{v}}^2\end{array}$$

Which will give us the DCM (Direction Cosine Matrix [2]) to position the body in the inertial frame according to Figure 1. A certain impact with known value misaligns the vectors H and r, also removes the stability of the spin motion. Some examples are shown in Figure 5. One initial velocity is chosen arbitrarily for simulating the motion, to the three different aforementioned geometries.

Considering a initial velocity ω₀ = [1.2 2.2 1]T.

In the oblate body (µ₁ = 0.9 µ₂ = 0.5) is obtained:

K = 1.4 and Capt = -0.078, therefore, tumbling occurs, and the motion is unstable.

Results to the prolate body (µ₁ = 2 µ₂ = 1.8):

K = 0.56 and Capt = 0.001, then, non tumbling occurs, and the motion is stable.

And finally, for the intermediate body (µ₁ = 1.11 µ₂ = 0.56):

K = 1.24 and Capt = -0.117, as expected, the motion is unstable.

To identify the occurrence of tumbling with the definition of only one parameter is useful and simple. The numerical results are consistent with the reasoning developed in this work. Thus, the stability or instability of the motion can be predicted by observing only one parameter.

Stability Chards Using Magnus Triangle

Magnus triangle is a practical way to represent rigid bodies of different geometries inside of the closed area of a triangle [15]. Each side represents relations of moments of inertia and if two of these values are known the lines connecting this points to the opposite vertex will cross defining the body. The third relation will be obtained connecting the third vertex to this crossing and looking for the value at the third triangle side.

The inequality Equation (18) will be explored for different body geometries using the triangle. The initial condition is given, and we may look how bodies will behave. Then, the regions where µ₁ K > 1 and µ₂ K < 1 (or µ₁ K < 1 and µ₂ K > 1) are identified according to a positive or negative tumbling coefficient Capt as stable or unstable.

In Figure 6, only one variable was considered (µ₂ K), a curved line can be observed for the value µ₂ K = 1, dividing the triangle area in two, one white and another yellow. There are 2 dots, the dot 1 that defines an oblate body and the dot 2 that defines a prolate body.

Also shown in this Figure 6 is the way to locate the points (or bodies), defining a line that starts from the vertex and reaches the opposite side of the triangle at the correspondent variable value.

In Figure 6, Figure 7 and Figure 8, the triangle is divided in 6 areas by its barycentric lines [15], the two upper areas belonging to an oblate body (for example the dot 1), the two lower areas belonging to a prolate body (for example the dot 2), and the intermediate areas (like the butterfly) belong to an intermediate body.

With an initial condition ω₀ = [0.5 0.5 1]T Figure 7 shows these regions. When Capt < 0 there is instability, and the region is colored grey. Two dots are highlighted in the triangle: it is the same body, once with a spin axis which turns it oblate µ₁ = 0.9 and µ₂ = 0.5 (dot 1), the other changes the spin axis and it rotates as an prolate body µ₁ = 2 and µ₂ = 1,8 (dot 2). We can notice that the body in the first condition is much nearer to the stability border [16].

Figure 8 compares two initial conditions:

ω₀ = [1 1 1]T (left) and

ω₀ = [0.2 1.5 1]T (right).

The curved lines correspond to µ₁ K = 1 and µ₂ K = 1.

In the last figure, the diagram loses the symmetry, which depends on the components of the angular velocity orthogonal to the spin axis. In this case, the velocity components on the x and y axes are different.

Non-Conservative Motion of an Oblate Body

The conservative analysis done up to now will be extended considering the non-conservative motion of an oblate body, which could represent a satellite. The idea is to use the tumbling coefficient as a stability parameter, in the sequence of the motion the same analysis is repeated, and the coefficients are compared.

Oblate body

The oblate body presents 1 > µ₁ > µ₂ and, as can be seen from Figure 3, also 1 < K < 1/µ₂ (in the range of abscissa of Figure 3), i.e., µ₂ < µ₂ K < 1. Therefore, in Equation (18), the positive or negative value of Capt will be defined by the part we call C:

$$Capt\text{ = }C\left(1-{\mu }_{2}K\right)\Rightarrow C\text{ = }1-{\mu }_{1}K\text{ (26)}$$

But we also know that Capt will be positive in the range 1 > K > 1/µ₁ (Figure 3). The non-dimensional K is time dependent in the non-conservative case [6].

Parameter C has a linear dependence of K (Equation 26). The development of motion towards stability needs a C that always increases (which means that K is tending to the limit K = 1). In other words, the temporal derivative of C has to be positive to leave the instability region. We can see this well in the Figure 9.

Using the definition of K given by Equation (11), now kinematic energy and angular momentum depend on time, the equation can be derived, so:

$$C\text{ = }1-{\mu }_1\frac{2Ec{I}_s}{h^2};\dot{C}\text{ = }\frac{2I_s{\mu }_1}{h^3}\left(2Ec\dot{h}-\dot{E}ch\right)\text{ (27)}$$

To guarantee a positive result for increasing values of C we need

$$2Ec\dot{h}-\dot{E}ch>0\Rightarrow \frac{\dot{E}c}{2Ec}<\frac{\dot{h}}{h}\text{ (28)}$$

The derivatives can be calculated using the definitions (matrix J is diagonal):

$$\begin{array}{l}Ec\text{ = }\frac{1}{2}{\omega }^{\text{T}}H\text{ (29)}\\ \dot{E}c\text{ = }\frac{1}{2}\left({\dot{\omega }}^{\text{T}}H+{\omega }^{\text{T}}J\dot{\omega }\right)\text{ = }\frac{1}{2}\left({\dot{\omega }}^{\text{T}}H+{\dot{\omega }}^{\text{T}}J\omega \right)\text{ = }{\omega }^{\text{T}}M\end{array}$$

$$h^2\text{ = }{}^{I}H{}^{\text{T}}{}^{I}H\Rightarrow \dot{h}h\text{ = }{}^{I}H{}^{\text{T}}\frac{d\left({}^{I}H\right)}{dt}\text{ = }{}^{I}H{}^{\text{T}}{}^{I}M\text{ (30)}$$

Replacing the derivatives in (28):

$$\begin{array}{l}2Ec{H}^{\text{T}}M-h^2{\omega }^{\text{T}}M>0\text{ (31)}\\ \Rightarrow {\left(2EcH-h^2\omega \right)}^{\text{T}}M>0\end{array}$$

Any applied torque on this oblate body (spinning on the axis of largest inertia) should obey this inequality to improve stability. A control law is possible trying to maximize de rate of change of C.

One could look for the resulting motion due to viscous damping. Using the torque as

$$M\text{ = }\left[\begin{array}{c}{a}_1{\omega }_1\\ a_2{\omega }_2\\ a_s{\omega }_s\end{array}\right]\text{ (32)}$$

To simplify the analysis, the coefficients on all degrees of freedom will be considered equal [14]: (a = a1 = a2 = as) and has to be negative for a damping condition, then:

$$\begin{array}{l}{\left(2EcH-h^2\omega \right)}^{\text{T}}M\text{ = }a{\left(2EcH-h^2\omega \right)}^{\text{T}}\omega \text{ = (33)}\\ a\left(2Ec{H}^{\text{T}}\omega -h^2{\omega }^{\text{T}}\omega \right)\text{ = }a\left(4E{c}^2-h^2{\omega }^2\right)\text{ = }\\ a\left(2Ec-\omega h\right)\left(2Ec+\omega h\right)>0\end{array}$$

Then,

$$a\left(2Ec-\omega h\right)>0\to 2Ec-\omega h<0\text{ (34)}$$

$$\omega h>2Ec\text{ = }{\omega }^{\text{T}}H\text{ = }\omega h\cos \delta \to 1>\cos \delta \text{ (35)}$$

This inequality is always satisfied when there is a misalignment of vectors ω e H. The loss of energy of the oblate body will lead to a stabilization around the spin axis of largest moment of inertia.

Prolate body

The prolate body presents µ₁ > µ₂ > 1 and, as can be seen from Figure 3, also 1/µ₁ < K < 1 (in the range of abscissa of Figure 3), i.e., 1 < µ₁ K < µ₁. Therefore, in Equation (18), the positive or negative value of Capt will be defined by the part we call B:

$$Capt\text{ = }\left({\mu }_{1}K-1\right)B\Rightarrow B\text{ = }{\mu }_{2}K-1\text{ (36)}$$

But we also know that Capt will be positive in the range 1/µ₂ > K > 1 (Figure 3). The non-dimensional K is time dependent in the non-conservative case.

Parameter B has a linear dependence of K (Equation 36). The development of motion towards stability needs a B that always increases (which means that K is tending to the limit K = 1). In other words, the temporal derivative of B has to be positive to leave the instability region.

Using the definition of K given by Equation (11), now kinematic energy and angular momentum depend on time, the equation can be derived, so:

$$B\text{ = }{\mu }_2\frac{2Ec{I}_s}{h^2}-1;\dot{B}\text{ = }\frac{2I_s{\mu }_2}{h^3}\left(\dot{E}ch-2Ec\dot{h}\right)\text{ (37)}$$

To guarantee a positive result for increasing values of B we need

$$\dot{E}ch-2Ec\dot{h}>0\Rightarrow \frac{\dot{E}c}{2Ec}>\frac{\dot{h}}{h}\text{ (38)}$$

The derivatives can be calculated using the definitions (matrix J is diagonal) in Equation 29 and 30. Replacing the derivatives in (38):

$$\begin{array}{l}{h}^2{\omega }^{\text{T}}M-2Ec{H}^{\text{T}}M>0\text{ (39)}\\ \Rightarrow {\left(h^2\omega -2EcH\right)}^{\text{T}}M>0\end{array}$$

Any applied torque on this oblate body (spinning on the axis of largest inertia) should obey this inequality to improve stability. A control law is possible trying to maximize de rate of change of B.

One could look for the resulting motion due to viscous damping. Using the torque as

$$M\text{ = }\left[\begin{array}{c}{a}_1{\omega }_1\\ a_2{\omega }_2\\ a_s{\omega }_s\end{array}\right]\text{ (40)}$$

To simplify the analysis, the coefficients on all degrees of freedom will be considered equal [14]: (a = a₁ = a₂ = a_s) and has to be negative for a damping condition, then:

$$\begin{array}{l}{\left(h^2\omega -2EcH\right)}^{\text{T}}M\text{ = }a{\left(h^2\omega -2EcH\right)}^{\text{T}}\omega \text{ = (41)}\\ a\left(h^2{\omega }^{\text{T}}\omega -2Ec{H}^{\text{T}}\omega \right)\text{ = }a\left(h^2{\omega }^2-4E{c}^2\right)\text{ = }\\ a\left(\omega h-2Ec\right)\left(2Ec+\omega h\right)>0\end{array}$$

Then,

$$a\left(\omega h-2Ec\right)>0\to \omega h-2Ec<0\text{ (42)}$$

$$\omega h<2Ec\text{ = }{\omega }^{\text{T}}H\text{ = }\omega h\cos \delta \to 1<\cos \delta \text{ (43)}$$

This inequality is never satisfied when there is a misalignment of vectors ω e H. The loss of energy of the prolate body will lead to a destabilization around the spin axis of largest moment of inertia. The motion should stabilize on another axis of rotation.

Intermediate body

The intermediate body presents µ₁ > 1 > µ₂ and, as can be seen from Figure 3, also 1/µ₁ < K < 1/µ₂ (in the range of abscissa of Figure 3), i.e., 1 < µ₁ K and µ₂ K < 1.

$$Capt\text{ = }\underset{+}{\underbrace{\left({\mu }_{1}K-1\right)}}\underset{-}{\underbrace{\left({\mu }_{2}K-1\right)}}\text{ = }\underset{+}{\underbrace{-C}}\underset{-}{\underbrace{B}}\text{ = }-\underset{-}{\underbrace{C}}\underset{-}{\underbrace{B}}<0\text{ (44)}$$

$$\frac{\partial Capt}{\partial t}\text{ = }-(\dot{C}B+C\dot{B})=\underset{+}{\underbrace{-B}}\dot{C}\underset{+}{\underbrace{-C}}\dot{B}\text{ (45)}$$

$$\frac{\partial Capt}{\partial t}\text{ = }\underset{+}{\underbrace{-B}}\underset{+}{\underbrace{\dot{C}}}\underset{+}{\underbrace{-C}}\underset{-}{\underbrace{\dot{B}}}\text{ (46)}$$

The result indicates that the body will tend to remain unstable [13]. Similar to the prolate body, the motion can stabilize on another axis of rotation.

Conclusion

The main purpose of this work, investigating the border of the tumbling motion of a body rotating in space, was achieved by an adequate definition of the coordinate systems, simplifying the mathematical analysis and diminishing the variables for the initial condition of the equations. The definition of a tumbling coefficient reduced the stability analysis to one parameter: if its value is positive the satellite does not tumble, otherwise it tumbles.

The tumbling condition does not depend on a specific attitude and will be defined through angular velocities. The inertial frame will be conveniently defined at any moment considered for initial condition and the angular velocity starts the tumbling analysis.

The normalized equations of motion describing any rotational behavior of the body will allow establishing the necessary conditions for an unstable behavior. These are necessary and sufficient conditions for the tumbling occurrence. Analytical expressions defining stability are established and stable and unstable regions for the consequent rotations are identified. This work provides useful tools to predict whether the rotational motion will remain stable or become unstable.

The analysis was verified by numerical results for different inertial configurations showing that the use of a single parameter is sufficient to predict tumbling.

Magnus stability triangle is a valuable tool to synthesize the stability conditions for bodies of any geometry. Stability areas are defined, borders are analytically established and conditions that favor stability are identified.

In the analysis of the non-conservative motion of an oblate body, the tumbling coefficient can be reduced to one relevant variable. Its derivative in time will result in a sufficient condition for tumbling and may help establishing the required maneuver to stabilize the motion. The amplitude of the resulting motion will be inversely proportional to this variable. Using this criterion for a damped motion, it is possible to show that there is always convergence to a rotational stability around the axis of larger inertia.

---