Study of Stability of Rotational Motion of Spacecraft with Canonical Variables

Abstract

This work aims to analyze the stability of the rotational motion of artificial satellites in circular orbit with the influence of gravity gradient torque, using the Andoyer variables. The used method in this paper to analyze stability is the Kovalev-Savchenko theorem. This method requires the reduction of the Hamiltonian in its normal form up to fourth order by means of canonical transformations around equilibrium points. The coefficients of the normal Hamiltonian are indispensable in the study of nonlinear stability of its equilibrium points according to the three established conditions in the theorem. Some physical and orbital data of real satellites were used in the numerical simulations. In comparison with previous work, the results show a greater number of equilibrium points and an optimization in the algorithm to determine the normal form and stability analysis. The results of this paper can directly contribute in maintaining the attitude of artificial satellites.

1. Introduction

Stability analysis of the rotational motion of a satellite taking into account the influence of external torques is very important in maintaining the attitude to ensure the success of a space mission.

Recently, some studies on the subject have been developed, and they motivated the development of this paper. In [1, 2] is presented a study on the stability of the rotational motion of artificial satellites in an elliptic orbit, considering disturbance due to gravity gradient torque, using canonical formulations and the Andoyer canonical variables. These studies use the procedure presented in [3] to determine a normal form of the Hamiltonian up to 4th order.

In [4] is developed a numerical-analytical method for normalization of Hamiltonian systems with 2 and 3 degrees of freedom. The normal form is obtained using the Lie-Hori method [5]. The stability analysis of the system is made by the Kovalev-Savchenko theorem [6]. The most important role of this work is the results obtained analytically for the generating function of 3rd order, necessary for determining the coefficients of the normal Hamiltonian of 4th order.

Thus, the objective of this work is to optimize the stability analysis developed in [2] to determine the equilibrium points, the normal form for dynamical systems with two degrees of freedom, and the applications of the Kovalev-Savchenko theorem [6]. It will be done by applying the expressions obtained in [4] for the coefficients of the normal 4th-order Hamiltonian.

In this paper equilibrium points and/or regions of stability are established when parcels associated with gravity gradient torque acting on the satellite are included in the equations of rotational motion.

The Andoyer variables are used to describe the rotational motion of the satellite in order to facilitate the application of methods of stability of Hamiltonian systems. The Andoyer canonical variables [7] are represented by generalized moments (𝐿₁,𝐿₂,𝐿₃) and by generalized coordinates (ℓ₁,ℓ₂,ℓ₃) that are outlined in Figure 1. The angular variablesℓ₁,ℓ₂,ℓ₃ are angles related to the satellite system Oxyz (with axes parallel to the spacecraft’s principal axes of inertia) and equatorial system OXYZ (with axes parallel to the axis of the Earth’s equatorial system). Variables metrics 𝐿₁, 𝐿₂,𝐿₃ are defined as follows: 𝐿₂ is the magnitude of the angular momentum of rotation 𝐿₂⃗,𝐿₁ is the projection of 𝐿₂⃗ on the 𝑧-axis of principal axis system of inertia (𝐿₁=𝐿₂cos𝐽₂, where 𝐽₂ is the angle between the z-satellite axis and 𝐿₂⃗), and 𝐿₃ is the projection of 𝐿₂⃗ on the Z-equatorial axis (𝐿₃=𝐿₂cos𝐼₂, where 𝐼₂ is the angle between Z-equatorial axis and 𝐿₂⃗).

Figure 1 

The Andoyer canonical variables.

The nonlinear stability of equilibrium points of the rotational motion is analyzed here by the Kovalev-Savchenko theorem [6], which requires the normalized Hamiltonian up to terms of fourth order around the equilibrium points.

The equilibrium points are found from the equations of motion. With the application of the Kovalev-Savchenko theorem, it is possible to verify if they remain stable under the influence of terms of higher order of the normal Hamiltonian.

In this paper numerical simulations were made for two hypothetical groups of artificial satellites, considering them in a circular orbit and with symmetric shape in relation to their physical and geometric characteristics. The satellites are classified as medium and small sized; they have orbital data and physical characteristics similar to real satellites.

2. Equations of Motion

The Andoyer variables, defined above, are used to characterize the rotational motion of a satellite around its center of mass [7], and the Delaunay variables describe the translational motion of the center of mass of the satellite around the Earth [8].

The Delaunay variables (𝐿,𝐺,𝐻,𝑙,𝑔,ℎ) are defined as [8] 𝐿=𝑀√𝜇𝑎,𝐺=𝐿√1−𝑒²,𝐻=𝐺cos𝐼, 𝑙 is the mean anomaly, 𝑔 is the argument of perigee, ℎ is the longitude of the ascending node,𝑀 is the mass of the satellite, 𝜇 is the Earth gravitational parameter, 𝐼 is the inclination of the orbit, 𝑎 is the semimajor axis, and 𝑒 is the orbital eccentricity.

Here it is assumed that the satellites are in a circular orbit, which differs from the study presented in [2] which considers satellites in eccentric orbits. This consideration was adopted to simplify the Hamiltonian of the problem, which is extensive [1], and to facilitate the stability analysis of the equilibrium points.

Thus, assuming that the satellites have well-defined circular orbit, the main focus of the work is only aimed at stability analysis of the rotational movement of the satellites in study.

In this case the Hamiltonian of the problem is expressed in terms of the Andoyer and Delaunay variables [9] as follows: 𝐹𝐿₁,𝐿₂,𝐿₃,ℓ₂,ℓ₃,𝐿,𝐺,𝐻,𝑙,𝑔,ℎ=𝐹_𝑜𝐿,𝐿₁,𝐿₂+𝐹₁𝐿₁,𝐿₂,𝐿₃,ℓ₂,ℓ₃,𝐿,𝐺,𝐻,𝑙,𝑔,ℎ, (2.1) where 𝐹_𝑜 is the unperturbed Hamiltonian and 𝐹₁ is the term of the Hamiltonian associated with the disturbance due to the gravity gradient torque, both are described respectively by [10]𝐹_𝑜𝜇=−²𝑀³2𝐿²+121𝐶−1−12𝐴𝐿2𝐵¹²+121𝐴+1𝐵𝐿₂²+141𝐵−1𝐴𝐿₂²−𝐿₁²cos2ℓ₁,𝐹₁=𝜇⁴𝑀⁶𝐿⁶2𝐶−𝐴−𝐵2ℋ₁ℓ_𝑚,𝐿_𝑛+𝐴−𝐵4ℋ₂ℓ_𝑚,𝐿_𝑛, (2.2) where 𝑚=2,3 and 𝑛=1,2,3; 𝐴,𝐵, and 𝐶 are the principal moments of inertia of the satellite on 𝑥-axis, 𝑦-axis, and 𝑧-axis, respectively; ℋ₁ and ℋ₂ are functions of the variables (ℓ_𝑚,𝐿_𝑛), where ℓ₂ and ℓ₃ appear in the arguments of cosines. The complete analytical expression is presented in [1] for eccentric orbits. In this paper 𝐹₁ will be simplified considering the satellite in a circular orbit; its complete analytical expression is given in the Appendix.

Thus, the equations of motion associated with the Hamiltonian (2.1) are given by𝑑ℓ_𝑖=𝑑𝑡𝜕𝐹𝜕𝐿_𝑖,𝑑𝐿_𝑖𝑑𝑡=−𝜕𝐹𝜕ℓ_𝑖(𝑖=1,2,3). (2.3)

These equations are used to find the possible equilibrium points of the rotational movement.

In this study, it is also considered that the satellite has two of its principal moments of inertia equal, 𝐵=𝐴. With this relationship, the variable ℓ₁ will not be present in the Hamiltonian, reducing the dynamic system to two degrees of freedom, a necessary condition for applying the stability theorem chosen for analysis of equilibrium points.

3. Determining the Normal Form of the Hamiltonian of the Rotational Motion

3.1. Normal Form of the Hamiltonian of 2nd Order and Linear Stability Analysis

Because the variable ℓ₁ is a cyclic coordinate, the Hamiltonian (2.1) is written as 𝐹(𝐿₂,𝐿₃,ℓ₂,ℓ₃) and the equations of motion (2.3) can be written in vector form as [1]𝐰̇=𝐽𝐹_𝐰, (3.1) where 𝐰 is the state vector andF_𝐰 is the matrix of partial derivatives with respect to their respective variables given byℓ𝐰=²𝐿₂ℓ₃𝐿₃,F_𝐰=𝜕𝐹𝜕ℓ₂𝜕𝐹𝜕𝐿₂𝜕𝐹𝜕ℓ₃𝜕𝐹𝜕𝐿₃ (3.2) and J a symplectic matrix given by𝐽=0100−1000000100−10. (3.3)

Now, it is necessary to linearize the system around the equilibrium point. However, it is convenient to make a translation of the coordinates in the Hamiltonian (2.1) so that the origin coincides with the equilibrium point under study; it means that𝐿₁=𝐿_1𝑒,ℓ₂=ℓ_2𝑒+𝑞₁,𝐿₂=𝐿_2𝑒+𝑝₁,ℓ₃=ℓ_3𝑒+𝑞₂,𝐿₃=𝐿_3𝑒+𝑝₂, (3.4) where 𝐿_1𝑒,ℓ_2𝑒,𝐿_2𝑒,ℓ_3𝑒, and𝐿_3𝑒 are the coordinates of equilibrium points.

With this translation we can expand the Hamiltonian in Taylor series around the new origin that is nothing more than expand the Hamiltonian in the neighborhoods of the equilibrium point when 𝑞₁=𝑝₁=𝑞₂=𝑝₂=0. Then,𝐹ℓ_𝑚,𝐿_𝑛=^∞_𝑘=2𝐹_𝑘ℓ_𝑚,𝐿_𝑛=𝐹₂ℓ_𝑚,𝐿_𝑛+𝐹₃ℓ_𝑚,𝐿_𝑛+𝐹₄ℓ_𝑚,𝐿_𝑛+⋯, (3.5) where 𝑚=2,3, 𝑛=1,2,3, and 𝐹_𝑘 is the Hamiltonian expanded up to terms of order 𝑘, with 𝑘=2,3,4,….

The Hessian 𝑃 of the problem is calculated from the Hamiltonian expanded up to 2nd order around the equilibrium point:𝜕𝑃=²𝐹𝜕𝑞²₁𝜕²𝐹𝜕𝑞₁𝜕𝑝₁𝜕²𝐹𝜕𝑞₁𝜕𝑞₂𝜕²𝐹𝜕𝑞₁𝜕𝑝₂𝜕²𝐹𝜕𝑝₁𝜕𝑞₁𝜕²𝐹𝜕𝑝²₁𝜕²𝐹𝜕𝑝₁𝜕𝑞₂𝜕²𝐹𝜕𝑝₁𝜕𝑝₂𝜕²𝐹𝜕𝑞₂𝜕𝑞₁𝜕²𝐹𝜕𝑞₂𝜕𝑝₁𝜕²𝐹𝜕𝑞²₂𝜕²𝐹𝜕𝑞₂𝜕𝑝₂𝜕²𝐹𝜕𝑝₂𝜕𝑞₁𝜕²𝐹𝜕𝑝₂𝜕𝑝₁𝜕²𝐹𝜕𝑝₂𝜕𝑞₂𝜕²𝐹𝜕𝑝²₂. (3.6)

Thus, the system of linearized equations can be written as𝒲̇=𝐽𝑃𝒲, (3.7) where 𝒲 is the state vector in the new variables𝑞𝒲=¹𝑝₁𝑞₂𝑝₂. (3.8)

So, to find the eigenvalues of the reduced linear system, we should diagonalize the matrix 𝐽𝑃 present in (3.7) by means of the determinant det(𝜆𝐼−𝐽𝑃)=0, where 𝐼 is an identity matrix and 𝜆 are the eigenvalues of the matrix 𝐽𝑃. If these eigenvalues are pure imaginary 𝜆_𝑗=±𝑖𝜔_𝑗, the normal form of 2nd-order Hamiltonian can be expressed as [2]𝐻^∗₂𝑞_𝑗,𝑝_𝑗=²_𝑗=1𝜔_𝑗2𝑞_𝑗²+𝑝_𝑗². (3.9)

3.2. Extension of the Normal Form of the Hamiltonian to Higher Orders

Extending the process of normalization of the Hamiltonian to higher orders (𝐻₃,𝐻₄,…), the Hamiltonian normal in terms of variables 𝑞_𝑖 and 𝑝_𝑖 can be expressed by𝐻𝑞_𝑗,𝑝_𝑗=²_𝑗=1𝜔_𝑖2𝑞_𝑗²+𝑝_𝑗²+𝐻^∗₃𝑞_𝑗,𝑝_𝑗+𝐻^∗₄𝑞_𝑗,𝑝_𝑗+⋯. (3.10)

To obtain the Hamiltonian normal 𝐻^∗₃,𝐻^∗₄,…, it is necessary to use the method of Lie-Hori. According to [11], it is desirable that the normal Hamiltonian is written in complex variables before using this method.

This transformation of variables is given by𝑞_𝑗=1√2𝑥_𝑗+𝑖𝑦_𝑗,𝑝_𝑗=1√2𝑖𝑥_𝑗+𝑦_𝑗 (3.11) with inverse given by𝑥_𝑗=1√2𝑞_𝑗−𝑖𝑝_𝑗,𝑦_𝑗=1√2−𝑖𝑞_𝑗+𝑝_𝑗. (3.12)

Thus, the normal Hamiltonian of 2nd order is given by𝐻^∗₂𝑥_𝑗,𝑦_𝑗=²_𝑗=1𝑖𝜔_𝑗𝑥_𝑗𝑦_𝑗 (3.13) and normal Hamiltonian until the 4th order in complex variables is represented by𝐻^∗𝑥_𝑗,𝑦_𝑗=²_𝑗=1𝑖𝜔_𝑗𝑥_𝑗𝑦_𝑗+𝐻^∗₃𝑥_𝑗,𝑦_𝑗+𝐻^∗₄𝑥_𝑗,𝑦_𝑗+⋯, (3.14) where 𝐻^∗₃(𝑥,𝑦) and 𝐻^∗₄(𝑥,𝑦) can be expressed briefly as𝐻^∗₃𝑥,𝑦=_{|𝛼|+𝛽=3}ℎ_3,𝛼,𝛽𝑥^𝛼𝑦^𝛽,𝐻^∗₄𝑥,𝑦=_{|𝛼|+𝛽=4}ℎ_4,𝛼,𝛽𝑥^𝛼𝑦^𝛽. (3.15)(3.16)

In (3.15) and (3.16), is always assumed it 𝛼,𝛽>0. For (3.15), |𝛼|+|𝛽|=3 means that the sum of the exponents of the new variables must obey the order of the Hamiltonian 𝐻^∗₃, and, for (3.16), |𝛼|+|𝛽|=4 means that the sum of the exponents of the new variables must obey the order of the Hamiltonian𝐻^∗₄. These equations are presented entirely in [1, 4].

Thus, we have the Hamiltonian written in complex variables. It is important to note that 𝐻₂ is already normalized, and now we can apply the method of Lie-Hori to find the normal form to higher orders 𝐻^∗₃,𝐻^∗₄.

3.2.1. Series Method of Lie-Hori

The normal form is obtained from the expansion of the Hamiltonian in terms of the Lie series given by [12, 13]𝐻^new=𝐻^∗+𝐻^∗+1,𝐺𝐻2!^∗+1,𝐺,𝐺𝐻3!^∗,𝐺,𝐺,𝐺+⋯, (3.17) where 𝐻^new is normal Hamiltonian obtained by canonical transformations around the equilibrium point, 𝐺=𝐺(𝑥_𝑖,𝑦_𝑖) is the generating function developed as a power series, where each degree of 𝐺_𝑛 is a sum of homogeneous polynomials of degree 𝑛, 𝐻^∗=𝐻^∗(𝑥,𝑦) is the original Hamiltonian, and {𝐻^∗,𝐺} is the Poisson brackets of order 𝑟+𝑠−2, where 𝑟 and 𝑠 represent the order of the polynomials 𝐻^∗ and 𝐺, respectively.

The Poisson bracket is defined as 𝐻^∗,𝐺=𝜕𝐻^∗𝜕𝑥𝜕𝐺−𝜕𝑦𝜕𝐻^∗𝜕𝑦𝜕𝐺𝜕𝑥. (3.18)

The new Hamiltonian 𝐻^new after the new canonical transformation in the neighborhoods of identity can be expressed as [2]𝐻^new=𝐻^new₂+𝐻^new₃+𝐻^new₄+⋯. (3.19)

So the new Hamiltonian, ordered by degree up to 4th order, can be expressed following the next form [11, 13]:2ndorder:𝐻^new₂=𝐻^∗₂,3rdorder:𝐻^new₃=𝐻^∗₃+𝐻^∗₂,𝐺₃,4thorder:𝐻^new₄=𝐻^∗₄+𝐻^∗₃,𝐺₃+1𝐻2!^∗2,𝐺₃,𝐺₃+𝐻^∗₂,𝐺₄. (3.20)(3.21)(3.22)

As 𝐻^new₂ is already given in its normal form, the next step is to calculate the generating function which leads to the normal form up to 4th-order terms. This procedure allows to obtain the minimum number of nonresonant monomials defined in [12]; thus, the generating function 𝐺_𝑛 is chosen to eliminate the complex variables in the total Hamiltonian (3.5) in which the monomials 𝑥_𝑗 and 𝑦_𝑗 have different exponents of the desired degree, leaving only the monomials that carry on the resonance of the intrinsic Hamiltonian systems. The following is a procedure to obtain the normal Hamiltonian to 4th order [1, 4].

(I) It is known that Hamiltonian systems are naturally resonant in even orders, 2nd, 4th, and so forth [10, 12], so to find the normal Hamiltonian until the 4th order, the 3rd disturbance 𝐻^new₃must first be removed from 𝐻^new, it means𝐻^new₃=0. (3.23)

Assuming 𝐻^∗₃ is briefly defined in (3.15) to eliminate terms of 3rd order, we consider the generating function 𝐺₃(𝑥,𝑦) expressed in abbreviated form as𝐺₃𝑥,𝑦=_{|𝛼|+𝛽=3}𝑔_3,𝛼,𝛽𝑥^𝛼𝑦^𝛽. (3.24)

Then it is necessary to find 𝐺₃(𝑥,𝑦) of the third-order homological equation (3.21), where𝐻^∗₂,𝐺₃=−𝐻^∗₃. (3.25)

After some algebraic manipulations, we find 𝐺₃(𝑥,𝑦) given by [4]𝐺₃𝑥,𝑦=_{|𝛼|+𝛽=3}−ℎ_3,𝛼,𝛽𝑥𝛽−𝛼,𝑊^𝛼𝑦^𝛽. (3.26)

Note that this equation is well defined if ⟨𝛽−𝛼,𝑊⟩≠0, where ⟨𝛽−𝛼,𝑊⟩ represents an inner product between 𝛽−𝛼=(𝛽₁−𝛼₁,𝛽₂−𝛼₂,…,𝛽_𝑛−𝛼_𝑛) and 𝑊=(𝜆₁,𝜆₂,…,𝜆_𝑛)=(𝑖𝜔₁,𝑖𝜔₂,…,𝑖𝜔_𝑛) which are pure imaginary eigenvalues obtained from the quadratic part of normal Hamiltonian. The complete analytical representation for ⟨𝛽−𝛼,𝑊⟩ was presented in [4].

For the resonance condition it is not accepted that ⟨𝛽−𝛼,𝑊⟩≠0, which is always true for monomials of third-order Hamiltonian when considering only the natural resonance of the Hamiltonian.

(II) This process is repeated for the polynomial of degree 4 finding 𝐺₄(𝑥,𝑦) that will eliminate some of monomials 𝐻^∗₄ in (3.16). The remaining monomials in 𝐻^new₄ are called resonant monomials [3]. Then (3.22) can be reduced as [1] 𝐺₄,𝐻^∗₂+𝐻^new₄=𝐻^∗₄+𝐻^∗₃,𝐺₃+1𝐻2!^∗2,𝐺₃,𝐺₃,𝐺₄,𝐻^∗₂+𝐻^new₄=𝐹_4. (3.27)(3.28) Equation (3.28) is called homological equation of 4th order, and it can be rewritten as𝐺₄,𝐻^∗₂=𝐹₄−𝐻^new₄,𝐺₄,𝐻^∗₂=𝑅₄. (3.29)

Assume that 𝑅₄(𝑥,𝑦) and 𝐺₄(𝑥,𝑦) will be expressed in short form as𝑅₄𝑥,𝑦=_{|𝛼|+𝛽=4}𝑟_4,𝛼,𝛽𝑥^𝛼𝑦^𝛽,𝐺₄𝑥,𝑦=_{|𝛼|+𝛽=4}𝑔_4,𝛼,𝛽𝑥^𝛼𝑦^𝛽. (3.30)(3.31)

After some algebraic manipulations, we find 𝐺₄(𝑥,𝑦) given by [4] 𝐺₄𝑥,𝑦=_{|𝛼|+𝛽=4}𝑟_4,𝛼,𝛽𝑥𝛼−𝛽,𝑊^𝛼𝑦^𝛽, (3.32) where ⟨𝛼−𝛽,𝑊⟩≠0 represents an inner product as mentioned above.

Thus, using (3.34) in (3.29) it is possible to find 𝐻^new₄ after some algebraic manipulations.

This procedure can be repeated for any order that is required in order to find a normal form of the new Hamiltonian 𝐻^new. In this paper it is considered the 4th order of the 𝐻^new.

3.3. Normal Form of the Hamiltonian of 4th Order and Nonlinear Stability Analysis

As we get the normal form of 2nd-order Hamiltonian (3.9), we can apply the method described previously to find the Hamiltonian normal until the 4th order using the Hamiltonian in complex variables (3.14) represented by [1, 4]𝐻^∗𝑥_𝑗,𝑦_𝑗=𝑖𝜔₁𝑥₁𝑦₁+𝑖𝜔₂𝑥₂𝑦₂+𝐻^∗₃𝑥₁,𝑦₁,𝑥₂,𝑦₂+𝐻^∗₄𝑥₁,𝑦₁,𝑥₂,𝑦₂+⋯, (3.33) where 𝐻^∗₃(𝑥₁,𝑦₁,𝑥₂,𝑦₂) and 𝐻^∗₄(𝑥₁,𝑦₁,𝑥₂,𝑦₂) are obtained from expressions𝐻₃^∗𝑥₁,𝑦₁,𝑥₂,𝑦₂=_{|𝛼|+𝛽=3}ℎ_{3,𝛼1,𝛽1,𝛼2,𝛽2}𝑥₁^𝛼1𝑦₁^𝛽1𝑥₂^𝛼2𝑦₂^𝛽2,𝐻^∗₄𝑥₁,𝑦₁,𝑥₂,𝑦₂=_{|𝛼|+𝛽=4}ℎ_{4,𝛼1,𝛽1,𝛼2,𝛽2}𝑥₁^𝛼1𝑦₁^𝛽1𝑥₂^𝛼2𝑦₂^𝛽2. (3.34)

Using the method of Lie-Hori presented before, it is possible to determine the generating function 𝐺₃(𝑥₁,𝑦₁,𝑥₂,𝑦₂) in order to satisfy the homological equation (3.25), where 𝐺₃(𝑥₁,𝑦₁,𝑥₂,𝑦₂) was determined analytically in [4]. This generating function is responsible for the elimination of the 3rd-order terms of the Hamiltonian (3.33).

The generating function 𝐺₃(𝑥₁,𝑦₁,𝑥₂,𝑦₂) is used to find the terms of 4th order, 𝐹₄, and perform the separation described in (3.28), obtaining 𝐻₄ and 𝐺₄(𝑥₁,𝑦₁,𝑥₂,𝑦₂). The coefficients of the normal form 𝐻^new₄ can be calculated using the terms obtained in 𝐻₄ and the eigenvalues of 𝐻₂.

Thus, the normal Hamiltonian in complex variables separated by degree can be expressed as follows[4]: 2ndorder:𝐻^new₂=𝑖𝜔₁𝑥₁𝑦₁+𝑖𝜔₂𝑥₂𝑦₂,3rdorder:𝐻^new₃=0,4thorder:𝐻^new₄=𝛿₁₁𝑥₁𝑦₁²+𝛿₁₂𝑥₁𝑦₁𝑥₂𝑦₂+𝛿₂₂𝑥₂𝑦₂², (3.35)(3.36)(3.37) where 𝜔_𝑗(𝑗=1,2) is the imaginary part of eigenvalues associated with the matrix defined by the product of a matrix symplectic of order 4 with the Hessian of the Hamiltonian expanded in Taylor series up to 2nd order around the equilibrium point; 𝛿_𝑖𝑗 are real coefficients obtained from the combination of ℎ_{𝑘,𝛼,𝛽}𝑥^𝛼𝑦^𝛽(𝑘=3,4).

However, the Kovalev-Savchenko theorem requires that the normal Hamiltonian is in real variables, and with the applications of the transformation of variables (3.12) the normal form Hamiltonian can been expressed as follows:2ndorder:𝐻^new₂=𝜔₁2𝑞₁²+𝑝₁²+𝜔₂2𝑞₂²+𝑝₂²,3rdorder:𝐻^new₃=0,4thorder:𝐻^new₄=𝛿₁₁𝑞₁⁴+𝑝₁⁴+2𝑞₁²𝑝₁²+𝛿₁₂𝑞₁²𝑞₂²+𝑞₁²𝑝₂²+𝑝₁²𝑞₂²+𝑝₁²𝑝₂²+𝛿₂₂𝑞₂⁴+𝑝₂⁴+2𝑞₂²𝑝₂², (3.38)(3.39)(3.40) where 𝛿_𝑖𝑗 are the coefficients of the Hamiltonian normal 4th order which are expressed analytically in [4].

4. Kovalev-Savchenko Theorem

To use the stability theorem, the normal Hamiltonian of the problem is necessary as it was discussed above.

Considering the normal Hamiltonian ℋ^o, an analytics function of coordinates (𝑞_𝜈) and generalized moments (𝑝_𝜈) to a fixed point 𝑃 is expressed byℋ^o=²_𝜈=1𝜔^𝑜_𝜈2ℛ_𝜈+²_𝜈,𝜐=1𝛿^𝑜_𝜈,𝜐4ℛ_𝜈ℛ_𝜐+𝑂₅, (4.1) where 𝑂₅ represents higher-order terms; 𝜔^o_𝜈 is the imaginary part of eigenvalues associated with the matrix defined by the product of a 4th-order matrix symplectic with the Hessian of the Hamiltonian expanded in Taylor series up to 2nd order around the equilibrium point; 𝛿^o_𝜈,𝜐 depend on the eigenvalues 𝜔^o_𝜈 and the coefficients of the Hamiltonian expanded in Taylor series of 3rd and 4th order around the equilibrium point and they are presented analytically in Formiga [4], and ℛ_𝑚=𝑞_𝑚²+𝑝_𝑚²with𝑚=1,2. (4.2)

The stability analysis is performed here by the Kovalev-Savchenko theorem [6] which ensures that the motion is Lyapunov stable if the following conditions are satisfied.(i)The eigenvalues of the reduced linear system are pure imaginary ±𝑖𝜔^𝑜₁ and ±𝑖𝜔^𝑜₂.(ii)The condition 𝑘₁𝜔^𝑜₁+𝑘₂𝜔^𝑜₂≠0. (4.3) is valid for all 𝑘₁ and 𝑘₂ integers satisfying the inequality 𝑘₁+𝑘₂≤4. (4.4)(iii)The determinant 𝐷^𝑜 must satisfy the inequality 𝐷^𝑜=−𝛿^𝑜₁₁𝜔^𝑜₂²−2𝛿^𝑜₁₂𝜔^𝑜₁𝜔^𝑜₂+𝛿^𝑜₂₂𝜔^𝑜₁²≠0, (4.5) where 𝛿^o_𝜈,𝜐 are the coefficients of the normal 4th-order Hamiltonian.

This theorem says that a Hamiltonian reduced in its normal form up to 4th order, in the absence of the resonance condition of the eigenvalues associated and if the condition (4.5) is satisfied, it is guaranteed the existence of tori invariant in a neighborhood small enough of equilibrium position [14–16].

5. Computational Algorithm for the Normal Form of the Hamiltonian of Rotational Motion and Analysis of Stability

In order to synthetize the process to analyze the stability of equilibrium points, the logical sequence of the algorithm is now presented.

5.1. Stability Analysis of Equilibrium Points

See Figure 2.

Figure 2 

Flowchart representative of the stability analysis of equilibrium points.

6. Numerical Simulations

As already mentioned, we consider two types of satellites: medium sized (MP), which has similar orbital characteristics with the American satellite PEGASUS [17], and small sized (PP), which has similar orbital characteristics of the Brazilian data collection satellite SCD-1 and SCD-2 [18]. All the numerical simulations were developed using the software MATHEMATICA.

Table 1 shows a quantitative summary of the equilibrium points found and their stability, according to the criteria of Kovalev-Savchenko theorem [6]. It can be observed that, when the first condition is not satisfied, the eigenvalues associated with the matrix 𝐽𝑃 are real or are not pure imaginary. This equilibrium point is not linearly stable.

Table 1 

Quantitative summary of the classification of equilibrium points.

There were found totals of 60 equilibrium points for the MP satellite, 50 equilibrium points to the PP-1 satellite and 50 equilibrium points for the PP-2 satellite. Tables 2, 3 and 4, show two equilibrium points found in the simulations, one being Lyapunov stable and the other unstable, to satellites MP, PP1 and PP2, respectively.

Table 2 

Satellite MP: 𝐴=𝐵=3,9499∗10⁻¹kg⋅(km)², 𝐶=1,0307∗10⁻¹kg⋅(km)².

Table 3 

Satellite PP-1: 𝐴=𝐵=9.20∗10⁻⁶kg⋅(km)², 𝐶=13∗10⁻⁶kg⋅(km)².

Table 4 

PP-2: 𝐴=𝐵=28,7∗10⁻⁴kg⋅(km)², 𝐶=36,9∗10⁻⁴kg⋅(km)².

For the Lyapunov stable equilibrium points of Tables 2, 3, and 4 respectively, Table 5 shows the values of the angles 𝐼₂, 𝐽₂, the rotation speed 𝜔, and the rotation period 𝑇. These values characterize the nonexistence of singularities in these points; it means that the angles𝐼₂ and 𝐽₂ are not null or close to zero. By Figure 1, when these angles 𝐼₂ and 𝐽₂ are null or close to zero, the Andoyer variables 𝑙₁, 𝑙₂, and 𝑙₃ are indeterminate, because it is difficult to determine the intersection between the involved planes in the definitions of these variables. This analysis was performed for all equilibrium points found in the simulations. Table 6 shows the same conditions for the unstable points of Tables 2, 3, and 4.

Table 5 

Analysis of possible singularities of the Anloyer variables for satellites MP, PP1 of and PP2.

Table 6 

Analysis of possible singularities of the Anloyer variables for satellites MP, PP1 and PP2.

7. Conclusion

In this paper we presented a semianalytical stability of the rotational motion of artificial satellites, considering the influence of gravity gradient torque for symmetric satellite in a circular orbit. Applications in numerical simulations, performed with the software MATHEMATICA, were made to two types of satellites: medium (MP) and small (PP).

Initially the points of equilibrium were determined using the physical, orbital, and attitude characteristics of each satellite. Then the algorithm for stability analysis was applied and it was obtained 2 stable equilibrium points for the satellite MP, 7 stable points for the satellite PP-1, and 10 stable points for the satellite PP-2.

For the satellite MP it was gotten only two equilibrium points because this satellite has similar characteristics to the satellite PEGASUS, which is tumbling [17]. For satellites PP-1 and PP-2 were obtained many other equilibrium points, but most of them were discarded, because they lead to the Andoyer variables, a condition of uniqueness (it means that the angles 𝐼₂ and 𝐽₂ are null or close to zero, and the Andoyer variables 𝑙₁, 𝑙₂, and 𝑙₃ are indeterminate).

It can be observed that a larger number of stable equilibrium points were determined in comparison with the results of [1, 2], which show the stability analysis of the rotational motion with the gravity gradient torque, but with the satellite in an elliptic orbit.

An optimization was done in the algorithm of determining the normal form of the Hamiltonian and the stability analysis, using expressions obtained by [4] for the coefficients of the normal 4th-order Hamiltonian. The introduction of these expressions enabled a more effective stability analysis for equilibrium points in comparison with the results of [1, 2]. It is possible to say that the numerical simulations have become less laborious allowing analysis of data in more numbers.

This paper presents results that can directly contribute in maintaining the attitude of artificial satellites. Once the regions of stability are known for the rotational motion, a smaller number of maneuvers to maintain the desired attitude can be accomplished. In this case, a fuel economy can be generated to the satellite with propulsion system, increasing the satellite’s lifetime.

Appendix

The disturbance Hamiltonian 𝐹₁, due to the gravity gradient torque, for a satellite in a circular orbit with two of its principal moments of inertia equal, 𝐵=𝐴, is given by [19] 𝐹₁=𝜇⁴𝑀⁷𝐿⁶𝐶−𝐴𝑀32𝐿₁𝐿₂²−12×−12+38𝐻1+𝐺²+𝐿₃𝐿₂²𝐻−3𝐺²𝐿₃𝐿₂²−38sin[2𝐼]sin2𝐼₂cosℎ−ℓ₃−38sin²[𝐼]sin²𝐼₂cos2ℎ−2ℓ₃−3𝐻161−3𝐺²sin2𝐼₂sin2𝐽₂cosℓ₂+3𝐻161−3𝐺²sin²𝐼₂sin²𝐽₂cos2ℓ₂+3𝐿161−³𝐿₂𝐿1+2³𝐿₂sin[2𝐼]sin2𝐽₂cosℎ−ℓ₃−ℓ₂+3𝐿161+³𝐿₂𝐿1−2³𝐿₂sin[2𝐼]sin2𝐽₂cosℎ−ℓ₃+ℓ₂+3𝐿161−³𝐿₂sin²[𝐼]sin𝐼₂sin2𝐽₂cos2ℎ−2ℓ₃+ℓ₂−3𝐿161−³𝐿₂sin²[𝐼]sin𝐼₂sin2𝐽₂cos2ℎ−ℓ₃+2ℓ₂−3𝐿161−³𝐿₂sin[2𝐼]sin𝐼₂sin²𝐽₂cosℎ−ℓ₃+2ℓ₂+3𝐿161+³𝐿₂sin[2𝐼]sin𝐼₂sin²𝐽₂cosℎ−ℓ₃−2ℓ₂−3𝐿321−³𝐿₂²sin²[𝐼]sin²𝐽₂cos2ℎ−2ℓ₃+2ℓ₂−3𝐿321+³𝐿₂²sin²[𝐼]sin²𝐽₂cos2ℎ−2ℓ₃−2ℓ₂+32𝐿₁𝐿₂²−1238sin²[𝐼]𝐿1−3³𝐿₂²−3×cos2𝑙+2𝑔8𝐻1−𝐺sin[𝐼]sin𝐼₂×cos2𝑙+2𝑔−ℎ−ℓ₃−3𝐻161+𝐺²sin²𝐼₂×cos2𝑙+2𝑔+2ℎ−2ℓ₃−3𝐻161−𝐺²sin²𝐼₂cos2𝑙+2𝑔−2ℎ−2ℓ₃+932sin²[𝐼]sin2𝐼₂sin2𝐽₂cos2𝑙+2𝑔+ℓ₂+932sin²[𝐼]sin2𝐼₂sin2𝐽₂cos2𝑙+2𝑔−ℓ₂−3𝐻161+𝐺𝐿1−³𝐿₂𝐿1+2³𝐿₂sin[𝐼]sin2𝐽₂×cos2𝑙+2𝑔+ℎ−ℓ₃+ℓ₂−3𝐻161+𝐺𝐿1+³𝐿₂𝐿1−2³𝐿₂sin[𝐼]sin2𝐽₂×cos2𝑙+2𝑔+ℎ−ℓ₃−ℓ₂+3𝐻161−𝐺𝐿1−³𝐿₂𝐿1+2³𝐿₂sin[𝐼]sin2𝐽₂×cos2𝑙+2𝑔−ℎ−ℓ₃−ℓ₂+3𝐻161−𝐺𝐿1+³𝐿₂𝐿1−2³𝐿₂sin[𝐼]sin2𝐽₂×cos2𝑙+2𝑔−ℎ−ℓ₃+ℓ₂+3𝐻321+𝐺²𝐿1−³𝐿₂sin𝐼₂sin2𝐽₂×cos2𝑙+2𝑔+2ℎ−2ℓ₃−ℓ₂−3𝐻321+𝐺²𝐿1+³𝐿₂sin𝐼₂sin2𝐽₂×cos2𝑙+2𝑔+2ℎ−2ℓ₃+ℓ₂+3𝐻321−𝐺²𝐿1−³𝐿₂sin𝐼₂sin2𝐽₂×cos2𝑙+2𝑔−2ℎ−2ℓ₃+ℓ₂−3𝐻321−𝐺²𝐿1+³𝐿₂sin𝐼₂sin2𝐽₂×cos2𝑙+2𝑔−2ℎ−2ℓ₃−ℓ₂−932sin²[𝐼]sin²𝐼₂×sin²𝐽₂cos2𝑙+2𝑔+2ℓ₂−932sin²[𝐼]sin²𝐼₂sin²𝐽₂cos2𝑙+2𝑔−2ℓ₂+3𝐻161+𝐺𝐿1−³𝐿₂sin[𝐼]sin𝐼₂sin²𝐽₂×cos2𝑙+2𝑔+ℎ−ℓ₃+2ℓ₂−3𝐻161+𝐺𝐿1+³𝐿₂sin[𝐼]sin𝐼₂sin²𝐽₂×cos2𝑙+2𝑔+ℎ−ℓ₃−2ℓ₂−3𝐻161−𝐺𝐿1−³𝐿₂sin[𝐼]sin𝐼₂sin²𝐽₂×cos2𝑙+2𝑔−ℎ−ℓ₃−2ℓ₂+3𝐻161−𝐺𝐿1+³𝐿₂sin[𝐼]sin𝐼₂sin²𝐽₂×cos2𝑙+2𝑔−ℎ−ℓ₃+2ℓ₂−3𝐻641−𝐺²𝐿1−³𝐿₂²sin²𝐽₂×cos2𝑙+2𝑔+2ℎ−2ℓ₃+2ℓ₂+3𝐻641−𝐺²𝐿1+³𝐿₂²sin²𝐽₂×cos2𝑙+2𝑔+2ℎ−2ℓ₃−2ℓ₂+3𝐻641+𝐺²𝐿1−³𝐿₂²sin²𝐽₂×cos2𝑙+2𝑔−2ℎ−2ℓ₃−2ℓ₂+3𝐻641+𝐺²𝐿1+³𝐿₂²sin²𝐽₂×cos2𝑙+2𝑔−2ℎ−2ℓ₃+2ℓ₂, (A.1) where: 𝐴=𝐵 and 𝐶 are the principal moments of inertia on 𝑥-axis, 𝑦-axis and 𝑧-axis respectively, 𝑀 is the mass of the satellite, 𝜇 is the gravitational constant,𝐼 is the inclination of the orbit, 𝐼₂ is the inclination of the plane of angular momentum with the plane of the equator, and 𝐽₂ is the inclination of the principal plane with the plane of the angular momentum.

In (A.1) the generalized moments (𝐿₁,𝐿₂,𝐿₃) are implicit in some terms, using the definitions of generalized moments present in the introduction and trigonometric properties, they can be explicit replacing.sin𝐼₂=𝐿1−³²𝐿₂²,sin²𝐼₂𝐿=1−³²𝐿₂²,sin2𝐼₂=2𝐿₃𝐿₂𝐿1−³²𝐿₂²,sin𝐽₂=𝐿1−¹²𝐿₂²,sin²𝐽₂𝐿=1−¹²𝐿₂²,sin2𝐽₂=2𝐿₁𝐿₂𝐿1−¹²𝐿₂². (A.2)

Acknowledgment

This present work was supported by CAPES.