Higher-Order Analytical Attitude Propagation of an Oblate Rigid Body under Gravity-Gradient Torque

A higher-order perturbation theory for the rotation of a uniaxial satellite under gravity-gradient torque demonstrates that known special configurations of the attitude dynamics at which the satellite rotates, on average, as in a torque-free state, are only the result of an early truncation of the secular frequencies of motion. In addition to providing a deeper insight into the dynamics, the higher order of the analytical solutionmakes it competitive when comparedwith the long-term numerical integration of the equations of motion.


Introduction
The torque-free rotation of artificial satellites may be perturbed by a variety of effects 1 .Due to the complexity of the force models to integrate, the problem of attitude propagation is commonly approached numerically, with a variety of available algorithms 2, 3 .Nevertheless, the problem can also be approached analytically, in which the analytical alternative is usually based on perturbation methods: the Euler-Poinsot problem is taken as the unperturbed part of the problem, and the other effects are perturbations of the torque-free rotation.
Among the perturbing torques that drive the rotational motion of an artificial satellite, a research topic of current interest is the study of the effects produced by external torques due to the satellite's interaction with the Earth's magnetic field 4, 5 .On the other hand, the gravity-gradient torque is often identified as one of the more important perturbations affecting the torque free rotation 6, 7 .Therefore, the model of a free rigid body perturbed by gravity-gradient torque appears in this literature on this matter as one of the basic, nonintegrable models used in the study of attitude propagation of artificial satellites 8, 9 .The usefulness of this model is not restricted to the case of artificial satellites and also fits the description of the rotational motion of natural satellites 10, 11 .

Mathematical Problems in Engineering
Whereas bodies whose attitude propagation may be of interest are triaxial in general, many celestial bodies have symmetry or near-symmetry rotation, which justifies the use of the uniaxial model in understanding the perturbed dynamics in this case.The uniaxial model is not limited to natural bodies, and this has also attracted the attention of aerospace engineers in the study of the attitude of artificial satellites 12 .
The uniaxial model still remains of interest in the study of the attitude propagation, either due to its direct application to actual problems 13 or because it can be taken as a truncation of a whole triaxial model as long as the departure from axisymmetry is small.In this last case, the uniaxial model is considered to be the zeroth order part in a perturbation approach in which the small triaxiality plays the role of a perturbation 14, 15 .In the case of natural celestial bodies, the work of external moments derived from the gravitational attraction of other bodies is normally negligible when compared with the kinetic energy of their torque-free rotation.This is why only first-order effects of the gravity-gradient torque are usually taken into account in the study of the rotation of most solar system bodies.However, even when the perturbation model only includes first order effects, extending the perturbation solution to higher orders provides a clearer insight into the dynamics involved.Thus, the secular terms of a first-order approach reveal special configurations of the satellite 13, 16 in which the satellite's attitude under gravity-gradient torque evolves, on average, as in the torque-free state, but with a slightly modified angular momentum.However, we will demonstrate that these special configurations do not survive when considering higher-order terms in calculating the solution.In addition, higher-order solutions would be definitely useful in increasing computational accuracy, thus extending the validity of an analytical approach for much longer intervals.For artificial satellites, the gravity-gradient effect may be much more important than in the case of natural celestial bodies, and the inclusion of second-order effects would be imperative when investigating the attitude propagation in the long term.In these case,s when the gravity-gradient second-order effects are nonnegligible, other effects such as a small triaxiality of the rigid body, or the orbit eccentricity or other external torques, may introduce observable frequencies in the rotation solution.
We will deal with the problem of the attitude propagation of a satellite under gravitygradient torque.In order to get insight into the contribution of higher-order terms in calculating the solution, we will make some simplifying assumptions and focus on the specific case of a uniaxial satellite under the disturbing gravitation of a perturber in circular orbit.The problem is solved by perturbation theory up to higher order effects in the gravity-gradient, thus extending the applicability of analytical results of 13 .The successive approximation method used by 13 now becomes too intricate for computing an analytical solution in which higher-order effects are involved.Therefore, from the beginning, we resort to the method of averaging using Lie transforms 17-19 to solve the nonintegrable perturbed problem.It deserves to be mentioned that approaching attitude dynamics problems with Lie transforms is not new and has been used from long ago in finding either analytical or semianalytical solutions of perturbed rotational motion see, for instance 6, 7, 20, 21 .
We use Deprit's method 18 , which is specifically designed for automatic machine computation which allows the perturbation theory to be easily computed to higher orders, although it is enough to compute the second-order effects introduced by the gravitygradient torque to show how the secular behavior changes with respect to that arising from lower-order truncation theories.Specifically, we will show that special configurations of the attitude dynamics at which the satellite has been established as rotating as in a torquefree state except for periodic terms are only the result of an early truncation of the secular frequencies.The new solution not only provides the required insight into the long-term dynamics, but also reveals it is highly competitive when compared with the numerical integration of the equations of motion.

Hamiltonian Formulation
The Hamiltonian of the rigid-body rotation is where T is the kinetic energy of rotation of a rigid body around its center of mass 2 ω 1 , ω 2 , and ω 3 are the components of the instantaneous rotation vector in the frame of the principal axis of the body, and A ≤ B ≤ C denote the principal moment of inertia.In the case which concerns us, the gravity-gradient torque, within the accuracy of higher-order terms in the ratio of the linear dimensions of the satellite to those of the orbit, the potential energy V of the external torques is taken from MacCullagh's approximation 22 , where m is the mass of the rotating body, m 1 is the mass of the disturbing body, r is the distance between the centers of mass of both bodies, G is the gravitational constant, and D Aα 2 Bβ 2 Cγ 2 is the moment of inertia of the rigid body with respect to an axis in the direction of the line joining its center of mass with the disturbing body's, with direction cosines α, β, and γ, where the constraint α 2 β 2 γ 2 1 is applied.
In order to get insight into the contribution of higher-order terms of the analytical solution, we make some simplifying assumptions.First, we assume that the rotation does not affect the orbital motion; hence, we neglect the Keplerian term in the potential.Then, we assume circular orbital motion with r a constant.Finally, we assume that the rotating body is oblate; hence, A B. Therefore, In the Hamiltonian setting, the components of the instantaneous rotation ω 1 , ω 2 , ω 3 , and the direction cosine γ in 2.4 are assumed to be state functions of any specific set of canonical variables.Although Euler angles provide an immediate description of the attitude, the Hamiltonian of the free rigid body H T does not reveal all the symmetries of the Euler-Poinsot problem when using Euler canonical variables.In contrast, the integrable character of the free rigid body problem is evident when using Andoyer variables 23, 24 , which take advantage of the angular momentum integral so as to use the plane perpendicular to it, leaving aside the center of mass of the rigid body.The invariant plane is defined by the canonical variables that link the inertial and rotating frames, in which is the latter defined by the principal axes of the rotating body with the x and y axes which define its equatorial plane.
Using Andoyer variables λ, μ, ν, Λ, M, N , the kinetic energy takes the form where M is the modulus of the angular momentum vector, N is its projection on the z axis of the body frame, and ν is the angle encompassed by the x axis of the body frame and the axis defined by the intersection of the equatorial plane of the body and the invariant plane.The variables μ, the angle encompassed by the intersections of the invariant plane with both the inertial plane and the equatorial plane of the rigid body, λ, the angle encompassed by the x axis of the inertial plane and the intersection of the inertial and invariant planes, and Λ, the projection of the angular momentum vector on the axis perpendicular to the inertial plane, are cyclic in T showing that the torque-free rotation is a problem of one degree of freedom and, therefore, integrable.Furthermore, in the case of axisymmetric bodies A B, the kinetic energy is simply which is a trivial Hamiltonian integration problem.
In order to express the disturbing function introduced by the gravity-gradient in Andoyer variables, it is convenient to refer the inertial motion to the orbital plane.Then, after the usual rotations that match the body and inertial frames, we obtain the direction cosine γ in Andoyer variables, where t is time, n is the mean orbital motion, J arccos N/M is the inclination angle between the equatorial and invariant planes, and I arccos Λ/M is the inclination angle between the invariant and inertial planes Figure 1 .
The explicit appearance of time in the Hamiltonian is easily avoided by moving to a rotating frame at the same rotation rate of the orbital motion.Thus, we introduce the new variable λ − nt; then, shows that this change of the reference frame will also require the introduction of the Coriolis term −nΛ in the Hamiltonian.As a result, we now deal with the conservative Hamiltonian K μ, ν, , M, N, L H μ, ν, , M, N, L − nL, where L Λ is the conjugate momentum to the new variable .Therefore, the Hamiltonian of the rotation of the oblate rigid body in the Andoyer variables in the rotating frame takes the form where we use the notation κ −Gm 1 C − A / 2a 3 , a 1 1/A, a 3 1/C.Note that ν is cyclic, and, therefore, N remains constant in the perturbed problem.The equations of motion of the perturbed problem are obtained from the Hamilton equations

2.11
We must note that Andoyer variables are singular for zero inclination of the intermediate plane with respect to either the inertial or the equatorial planes of the body, or both.These singularities may be avoided when using other sets of variables 8, 25 .

Perturbation Approach
In general, one can resort to the classical double averaging method to find the secular terms of the disturbing function see 26 , for instance .In our case, the averaging would remove the fast angles μ and from 2.9 .However useful the classical double averaging may be in finding the relevant long-term evolution of a dynamical system, in the case of the Hamiltonian equation 2.9 , it would be limited to providing the known first order terms in the pertinent literature 13, equation 22 .
In order to reach the higher orders required in this work, the perturbation solution is computed analytically by Lie transforms using Deprit's algorithm 18 .This method is based on solving the homological equation L 0 W n K n − H n , which, in general, is a partial differential equation where the term H n comes from previous computations, the term K n of the new Hamiltonian is chosen at our convenience, and L 0 is the Lie derivative, such that now the term W n of the generating function can be solved.
We constrain the solution to the case when the spin rate of the satellite is much faster than its orbital rate, thus precluding the problem of resonances, and then we place the Coriolis term at the first order.Besides, we assume that the gravity-gradient torque is a second-order effect.Then, the Hamiltonian equation 2.9 is ordered as K m≥0 K m,0 /m!, where 3.1 This order allows us to proceed stepwise to eliminate first the rotation angle μ by means of a Lie transform and then the angle by means of a second Lie transform.Moreover, this splitting of averages has the added advantage that the homological equation of the method can be solved by quadrature in both canonical transformations.
Note that if the spin and orbital rates were of the same order, the Lie derivative of the homological equation would be the operator which involves the solution of a partial differential equation, contrary to a quadrature, for computing the generating function of the canonical transformation.As shown with the generating function W just below 22 of 13 , this way of proceeding explicitly shows the denominators that would be small when close to resonances, thus preventing the convergence of the series solution.
To avoid the explicit appearance of square roots in the automatic evaluation of Poisson brackets required by the method, it is easier for us to handle circular functions of the inclination angles I and J as state functions of the Andoyer variables.Their nonvanishing partial derivatives with respect to the Andoyer variables are

Averaging of μ
The first canonical transformation μ, ν, , M, N, L → μ , ν , , M , N, L has the effect of averaging the Hamiltonian over the rotation angle μ .The Lie derivative is now the operator and the homological equation can be solved by a simple quadrature to compute the generating function of the canonical transformation.
The Hamiltonian in the new variables is K m≥0 K 0,m /m!, which, up to the fourth order approximation, that only depends on one angle, .As a consequence, up to the truncation order, the averaging makes M constant, and, therefore J .

Elimination of
A second canonical transformation μ , ν , , M , N, L → μ , ν , , M , N, L is now computed so as to eliminate .Now, the Lie derivative is and again the homological equation is solved by quadrature.From this second transformation we obtain the double averaged Hamiltonian K m≥0 K m /m! given by 32a 1 M 3 cos I 3 − 178 cos 2 J 207 cos 4 J − 3 5 − 126 cos 2 J 153 cos 4 J cos 2 I 27κ 3 32M 3 n 2 1 − 3 cos 2 J 3 3 − 5 cos 2 I cos I .

3.9
We note that the special configurations mentioned in 13 -critical inclinations such as cos 2 I cos 2 J 1/3 or cos I cos J 0 in which the rigid body evolves, on average, as in torque-free rotation, although at a slightly different rotation rate from the unperturbed caseare exactly the result of the order of approximation used, whose perturbation solution only considered first-order effects on the gravity-gradient perturbation which are equivalent to the second-order of our present approach by Lie transforms .These special configurations no longer exist when truncating the perturbation approach at higher orders.Thus, while special inclinations such as cos 2 I cos 2 J 1/3 are preserved up to the third-order truncation of 3.9 , where and the frequencies of the averaged problem correspond to a torque-free rotation state, this is no longer true for the case cos I cos J 0, where show that λ is no longer fixed, on average, and suffers from a small precessional motion.Furthermore, the unperturbed-type state at special configurations cos 2 I cos 2 J 1/3 is also destroyed at the fourth order of the theory, where 4M 3 a 1 .

3.12
These higher-order terms may explain the observed behavior in Figure 6 of 13 .

Numerical Comparisons
So as to evaluate the performance of the analytical solution, we take a fictitious Earth satellite with moments of inertia A B 400 kg km 2 , C 600 kg km 2 , which we assume to be rotating at a rate of 1 rotation per minute.The satellite is assumed to be in a circular orbit of the MEO region with a semimajor axis a 13 000 km.
In the first test case, we take the satellite in a high-inclination orbit with I 70 deg and assume that it is rotating close to the axis of maximum inertia with J 10 −3 .Besides, internal units such that the initial value of the modulus of the momentum is M 1 are chosen for the integration; the other initial conditions are μ ν 0 and 1 radian.For these initial conditions, the nonaveraged equations of motion, 2.11 , are integrated for 30 orbital periods, a time interval in which the satellite completes more than 11 000 rotation cycles.Results of the numerical integration are then compared with the analytical attitude propagation.In the latter case, the initial conditions must be transformed to the double averaged phase space, resulting in μ 0 0.005824457466, ν 0 −0.005932985721, 0 1.0003166358, M 0 1.0000000025, and L 0 0.34164643181.Results of the comparison between the analytical solution and the numerical integration of the nonaveraged equations are presented in Figure 2, where the left-hand column presents the differences between the numerical integration and the analytical propagation of the secular terms of the fourth order truncation defined by the frequencies in 3.9 .The right-hand column of Figure 2 shows the differences obtained when using the full fourth-order theory, which, in addition to the secular terms propagation, includes the third-order transformation equations for recovering the medium period terms related to the averaging over , and the fourth-order transformation equations which allow us to recover the short-period terms related to the averaging of μ.
As shown in the plotting on the left of Figure 2, the periodic errors have a noticeable amplitude when compared with the secular trend.The amplitude of the errors due to periodic terms is of about 2 arc minutes for λ, although this increases up to about one degree for μ and ν.When the full fourth order analytical theory is used right-hand column of Figure 2 , the periodic errors are confined to very small values, revealing a secular trend in the order of 10 −10 radians per cycle in angular-variables errors, even though the amplitude of the periodic errors of μ and ν mask this linear trend to some extent.The computation of higher orders in the perturbation approach should improve the behavior of the analytical theory for both the secular and periodic terms.
The second test case is for a satellite of the same physical characteristics and initial configuration, except that now we take cos I cos J 1/3, one of the special configurations of the lower-order theories.The propagation of these initial conditions in the nonaveraged model shows that this configuration is very close to an unperturbed state.As shown in Figure 3, the time history of the difference between μ and the unperturbed state μ μ 0 a 1 M t seems to consist of only periodic terms, as well as what happens to the difference between ν and the unperturbed analog ν ν 0 − a 1 − a 3 Nt.In fact, there is a small linear trend of the order of 3.1 • 10 −7 t for μ and 7.7 • 10 −7 t for ν, that is masked by the amplitude of periodic oscillations.In contrast, the linear trend of the differences, in spite of being of the order of 2.9 • 10 −8 t, is better appreciated in the evolution of λ because of the notably smaller amplitude of the periodic oscillations.To highlight this difference, we superimposed a linear fit to the differences λ − λ 0 to their time history presented in the inferior plotting in Figure 3, thus revealing a clear departure from this fit represented by the straight dashed white line from zero value.
The initial state in the double-averaged phase-space is now μ 0 0.00030577890713, ν 0 −0.0005329981843, 0 0.9999991669, M 0 1.0000016387, and L 0 0.577352484, which are the initial conditions that feed the analytical solution.The comparison between the numerical integration of the nonaveraged equations, and the perturbation solution is presented in Figure 4.As in the previous example, we note that the propagation of the secular terms alone introduces important periodic errors left-hand column of Figure 4 .In contrast, when recovering the periodic terms removed from in averaging process, the errors are reduced to quite acceptable values.As presented in the right-hand column of Figure 4, the momenta obtained from the fourth-order analytical theory are mainly affected by very small periodic errors, whilst the secular ones seem to be negligible.On the other hand, the angular variables are affected by both periodic and secular errors that are apparent in their time histories.Nevertheless, the secular errors grow at the very low rates of just a few microarc seconds during an orbital period for λ and ν, and tens of microarc seconds during an orbital period in the case of μ.

Conclusions
The rotation of an oblate rigid body under gravity-gradient torque is a nonintegrable problem that may be approached by perturbations.Lower-order approaches to the solution that only consider first-order effects in the gravity-gradient perturbation are normally considered enough in some applications.However, analytical solutions that consider the effects of the gravity-gradient torque up to the second-order may be required in the engineering problem of attitude propagation.Carrying the perturbation approach up to this higher-order provides a complete insight into the long term dynamics, and the perturbation solution continues being competitive when compared with the numerical integration of the equations of motion.The use of Andoyer variables resulted crucial for the perturbation approach used, because their use notably facilitates the solution of the partial differential equations that provide the generating functions of the Lie transforms.Finally, since modern perturbation methods are designed for automatic computation, the analytical solution may be easily extended to even higher orders if required.

Figure 1 :
Figure 1: Relations between the Andoyer angles and the inclination angles I and J.

Figure 2 :
Figure 2: Errors of the fourth-order analytical solution versus the numerical integration for initial conditions μ ν 0, λ 1, I 70 deg, J 10 −3 rad.Left : only the secular terms of the analytical solution are considered.Right : the analytical solution includes secular and periodic terms.The notation δ is used for relative errors, whereas Δ means absolute differences.

2 jFigure 4 :
Figure4: Errors of the fourth-order analytical solution versus the numerical integration for initial conditions μ ν 0, λ 1, cos I cos J 1/3 rad.Left : only the secular terms of the analytical solution are considered.Right : the analytical solution includes secular and periodic terms.The notation δ is used for relative errors, whereas Δ means absolute differences.