SciELO - Scientific Electronic Library Online

vol.27 número4A method to determinate the thickness control parameters in cold rolling process through predictive model via neural networksNumerical simulation of two dimensional compressible and incompressible flows índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados




Links relacionados


Journal of the Brazilian Society of Mechanical Sciences and Engineering

versão impressa ISSN 1678-5878versão On-line ISSN 1806-3691

J. Braz. Soc. Mech. Sci. & Eng. v.27 n.4 Rio de Janeiro out./dez. 2005 



Retrograde orbits perturbed by a third-body



A. F. B. A. Prado

Instituto Nacional de Pesquisas Espaciais; INPE-DMC; Av. dos Astronautas 1758; Caixa Postal 515; 12227-010 São José dos Campos, SP. Brazil




This paper develops a semi-analytical study of the perturbation caused to a spacecraft by a third body involved in the dynamics. There are several important applications for this research, such as to calculate the effect of lunar and solar perturbations on high-altitude Earth satellites. In the present research the goal is to study the evolution of retrograde orbits around the Earth. There is a special interest to see under which conditions a near-circular orbit remains near-circular. The existence of circular, equatorial and frozen orbits are also considered, The results are valid for any system of primaries by making a time transformation that depends on the masses of the bodies involved. Several plots will show the time-histories of the Keplerian elements of the orbits involved.

Keywords: Astrodynamics, third-body pertutbation, retrograde orbits, space trajectories




The majority of papers that considered the third-body perturbation problem studied the perturbation due to the Sun and the Moon in a satellite in orbit around the Earth. Kozai (1959) developed the main secular and long-period terms of the disturbing function due to the lunisolar perturbations in terms of the orbital elements of the satellite, the Sun and the Moon. This research would be later expanded by Musen, Bailie and Upton (1961) to include the parallactic term in the disturbing function. After that, Kozai (1962) studied the problem of secular perturbations in asteroids with high inclination and eccentricity, considering that they are perturbed by Jupiter, which is assumed to be in a circular orbit around the Sun. Blitzer (1959) obtained estimates of the lunisolar disturbances using methods of classical mechanics, however only for the secular terms. In the sequence, Cook (1962) used the Lagrange's planetary equations to obtain expressions for the variation of elements during a revolution of the satellite and for the rate of variation of the same elements. In that same year, Kaula (1962) derived general terms of the disturbing function for the lunisolar perturbation, using equatorial elements for the Moon, but it didn't supply a definitive algorithm for the calculations. Again Kozai (1965) approached that problem and included indirect terms of the disturbance due to the alteration of the terrestrial flattening due to those forces.

In the 1970's, that subject was studied again. Giacaglia (1973) obtained the disturbing function for the disturbance of the Moon using ecliptic elements for the Moon and equatorial elements for the satellite. Secular, long and short period terms were calculated and expressed in a closed form. Kozai (1973) developed an alternative method for the calculation of the lunisolar disturbances. The disturbing function was expressed in terms of the orbital elements of the satellite and the polar geocentric coordinates of the Sun and the Moon. The secular and long period terms are derived by numerical integration and the short period terms are obtained analytically.

In the following decade, Hough (1981) studied the effects of the lunisolar disturbance in orbits close to the inclinations 63.4° and 116.6° (critical inclinations with respect to the geopotential of the Earth) and concluded that the effects are significant in high altitudes.

All of the preceding references represent fundamental contributions in the subject and they possess an analytic focus, dedicated to the derivation of equations. In the present work a more practical approach is used with the idea of complementing the existent literature. Some recent studies have provided numeric comparisons. In particular, Prado (2002) studied the perturbation in a lunar satellite for direct orbits and the lifetimes of orbits around the natural satellites of the Solar System.

In this paper several topics related to this problem will be studied for retrograde orbits around the Earth, that is an important special case. In particular, the so called "critical angle of the third-body perturbation," which is a value for the inclination such that any near-circular orbit with inclination above this remains near-circular, is discussed in detail. The assumptions of our model are similar to the ones made in the restricted three-body problem: a) There are only three bodies involved in the system: one body with mass m0 fixed at the origin of the reference system; a massless spacecraft in a elliptic three-dimensional orbit around this body and a third body in a circular orbit around this same central body in the plane x-y; b) The motion of the spacecraft is assumed to be a three-dimensional Keplerian orbit with its orbital elements disturbed by the third body. The motion of the spacecraft is studied under the double-averaged analytical model with the disturbing function expanded in Legendre polynomials up to fourth-order. The double-average is taken over the mean motion of the satellite and over the mean motion of the disturbing body. The fact that modern computers can easily integrate numerical trajectories using complex models for the dynamics does not invalidate the use of models based in analytical approximations. The most important reason for this is that a double-averaged model, like the one shown here, can eliminate short-period periodic perturbations that appear in the trajectories. In that way, smooth curves that show the evolution of the mean orbital elements for a long time period can be constructed, which give a better understanding of the physical phenomenon studied and allow the study of long-term stability of the orbits in the presence of disturbances that cause slow changes in the orbital elements. Note that the truncated equations of motion could be numerically integrated much faster than the full equations of the restricted three-body problem. The next sections present: the mathematical model used, the study of circular, near-circular, equatorial and frozen orbits.


The Mathematical Model

This section derives the equations required by the mathematical model used during the simulations made in this research. It is assumed that the main body with mass m0 is fixed in the center of the reference system x-y. The perturbing body with mass m' is in a circular orbit with semi-major axis a' and mean motion n' (given by the expression n'2 a'3 = G[m0 + m']). The massless spacecraft is in an elliptic three-dimensional orbit with orbital elements: a (semi-major axis), e (eccentricity), i (inclination), (argument of periapsis), (longitude of the ascending node) and the mean motion is n (given by the expression n2a3 = Gm0). In this situation, the disturbing potential that the spacecraft has from the action of the disturbing body is given by:

where , G is the gravitational constant and S is the angle between the line that connects the massive central body and the perturbed body (the spacecraft) and the line that connects the massive central body and the perturbing body (the third body).

Using the traditional expansion in Legendre polynomials (assuming that r' >> r) the following expression can be found:

where Pn are the Legendre Polynomials.

For the models used in this research it is necessary to calculate the parts of the disturbing function due to P2, P3 and P4 (R2, R3 and R4, respectively). They are:

A substitution using the expression n'2 a'3 = G[m0 + m'])was made between the first and the second member of those equations. The next step is to average those quantities over the short period of the satellite as well as with respect to the distant perturbing body. The standard definition for average used in this research is , where M is the mean anomaly, which is proportional to time.

To perform the average over R2, R3 and R4 one proceeds as follows. First define the quantities a = and b = , where is the unit vector pointing from the central body to the disturbing body and and are the usual orthogonal unit vectors, functions of (i,w, W), in the plane of the satellite orbit, with pointing towards the periapsis. For the special case considered here of circular orbits for the disturbing body the following relations are available:

Using this definition and the geometry involved, it is possible to relate the angle S to the positions of the perturbing and perturbed bodies. This can be made by the equation (where f is the true anomaly of the satellite):

Then, combining this equation with equations (3) to (5), the perturbing potential becomes a function of the orbital elements of the satellite. Next, it is necessary to replace the true anomaly (f) by the eccentric anomaly (E). This is made by the well-known relations:, r/a = 1 - e cos(E). Then, the integrations required to obtain the averages are realized in terms of the eccentric anomaly, not in terms of the mean anomaly. To do that the relation dM = (1 - ecos(E))dE is also used. After this process, the following identities appear (Prado, 2002):

After using those quantities, the expressions (3) to (5) become:

Next, one takes the second average with respect to the disturbing body to eliminate the variable M'. To do this, it is necessary to held the Keplerian elements of the spacecraft constant during the process of averaging. This is possible due to the hierarchy of time scales: period of satellite << period of disturbing body << period of slow oscillations in the orbital elements. After making these assumptions the following identities, always valid for circular orbits only, are obtained:

Applying those identities in the expressions (6) to (8), the results are:


The partial derivatives required for the equations of motion are:


The next step is to obtain the equations of motion of the spacecraft. They come from the Lagrange's planetary equations in the form that depends on the derivatives of the disturbing function R with respect to the Keplerian elements.

where fi (a,e,i,w) represents the contribution of the fourth-order term and they are given by:

There are some conclusions that come directly from the equations of motion: i) the term m' is a constant that multiplies all the equations of motion, so it is equivalent to a time transformation in the system of the type µ't = t*. So, all the results based in those equations are valid for any system of primaries in a proportional time scale. The same is true for the semi-major axis (that is present in the equations in the term "n") in the second-order model; ii) the ratio gives an idea of the relative importance of the second and fourth-order terms. The importance of the fourth-order terms increase when the semi-major axis of the perturbed body increases. This importance also increases with the eccentricity of the perturbed body due to the terms that depend on the eccentricity; iii) The difference between the analytical solutions and the full numerical simulations also increases with those variables, what is confirmed by numerical integrations and by the fact that the expansions are made in terms of the eccentricity.



In this section some results are shown related to the third body perturbation problem for retrograde orbits. This section is divided in several sub-sections to show clearly several aspects of the problem.

The Circular Orbits and the Critical Inclination

Directly from the equations of motion for the averaged models it is possible to identify the existence of circular solutions for both second and fourth-order models. It means that, in the ideal case of an orbit that starts with zero eccentricity, its eccentricity remains always zero. This occurs because the right-hand side of the equation for the time derivative of the eccentricity is zero (it is a polynomial in the eccentricity with no independent term). Another property of those orbits is that the inclination is also constant for the same reason, since the time derivative of the inclination is also a polynomial in the eccentricity with no independent term.

The evolutions of these two quantities (eccentricity and inclination) are studied under the full restricted three-body problem. The results show that the circular solutions with constant inclination do not exist in this more realistic model. The eccentricity oscillates with large amplitude. The inclination remains close to constant most of the time, but from time to time it goes to the value of the critical inclination and then it returns to its initial value. This behavior is similar to the near-circular orbits shown in Figs. 1 to 3. This result is expected, because there is no physical reason to have a strong difference in the behavior of orbits with eccentricity 0.00 and 0.01. The general conclusion is that the circular solutions with constant inclination appear due to the truncation of the Legendre polynomial and are not a physical phenomenon, at least for the conditions simulated in this research.







Another important question in this problem is the existence of a critical value for the inclination between the perturbed and the perturbing bodies. This critical inclination is related to the stability of near-circular retrograde orbits. The problem is to find under what conditions a spacecraft that starts in a near-circular retrograde orbit around the main body remains in a near-circular orbit after some time. The answer for this question depends on the initial inclination i0. There is a specific critical value such that if the inclination is lower than that the eccentricity increases and the near-circular orbit becomes very elliptic. Alternatively, if the inclination is higher than this critical value the orbit stays nearly circular. So, only retrograde orbits with inclination higher than that are useful for practical purposes. The problem of near-circular orbits is very important because usually a spacecraft that is nominally in circular orbit experiences perturbations from other sources that make its eccentricity become non-zero. In the double-averaged second-order model this critical value is i = 140.7685 degree (cos2(i) = 0.60). This critical value also represents the highest inclination that allows the existence of orbits with eccentricity, inclination and argument of periapsis constant under the second-order model. The behavior of the inclination and the eccentricity with time is studied for near-circular retrograde orbits covering a large range of initial inclination (95 < i0 < 175). For those simulations the initial orbits used always have Keplerian elements a0 = 0.1 (38440 km), e0 = 0.01, w0 = W0 = 0. The initial inclination i0 vary as shown in the figures. Remember that the time is defined such that the period of the disturbing body is 2p. In that way 1000 units of time in those figures correspond to about 160 orbits of the disturbing body (in the case of the Earth-Moon system). The figures are divided in three parts only to make the observation of the time evolution easily: values of i0 above the critical value (i0 >150°), in the region of the critical value (140° < i0 < 144°) and below the critical value (i0 < 140°). The sizes of those regions are chosen to avoid unclear figures due to the different scales involved and are not subject to a detailed study. Fig. 1 shows the behavior of the inclination. Fig. 2 shows the behavior of the eccentricity. In both figures, only the results for the fourth-order model are shown. The second-order expansion gives very similar results and the restricted three-body problem gives a series of scattered points the follows around the lines shown here, similar to what is shown in Prado (2002) for direct orbits.

The figures shows that for values of the initial inclination i0 above the critical angle (i0 > 150) the eccentricity oscillates with a very small amplitude (less than 0.025) that increases rapidly when i0 decreases and the inclination remains constant. For values of i0 around the critical value (140 < i0 <144) it is possible to see that the eccentricity oscillates with a larger amplitude (about 0.30) that increases with i0. The inclination has a very characteristic behavior in this region of i0. For values of i0 slightly above the critical angle the inclination stays close to i0. For values of i0 slightly below the critical value the inclination starts at i0, increases until the critical value and then it returns to its original value i0. Those results show that the critical angle is not a sharp separation between stable and unstable near-circular retrograde orbits. For values of i0 well below the critical value (i0 < 130) the eccentricity oscillates with increasing amplitudes that go close to 1.0. The inclination keeps its characteristic behavior of starting at i0, increasing to the critical value and then returning to its original value i0. The figure also shows that this behavior repeats itself in an endless cycle. The time required to reach the critical value increases when i0 increases. This region has a gradual transition where the eccentricity oscillates with an amplitude that increases fast with i0, reaching the value 1.0 only in the case i0 = 90. Fig. 3 shows the results in the phase diagram that plots the eccentricity vs. the inclination. The region above the critical value shows straight lines. Note that the different scales of the plots show different aspects of the comparison.

It is possible to understand those behaviors by examining the second-order equations of motion, since the comparisons showed that this model adequately represents the system. The magnitude of the time derivative of the eccentricity is dependent on the term, so it increases when the eccentricity has lower values due to the presence of e, but after a value close to 0.71 is reached, it starts to decrease due to the term. Its sign is determined exclusively by sin(w) and it imposes the oscillatory behavior shown in the plots, since w has a secular variation. For the time derivative of the inclination the same analysis can be applied. The magnitude is influenced by the term, which causes the curious behavior of showing regions of almost constant inclination alternated with sharp decreases and increases. The explanation is that for lower values of the eccentricity the term sin(w) forces the time derivative to stay close to zero and the inclination remains almost constant. When the eccentricity increases and reaches values close to 1.0, the term present in the denominator forces the time derivative to increase fast, tending to infinity, which causes the fast motion of the inclination. The alternation of the sign is caused exclusively by the term sin(w), as explained before. The argument of periapsis has a secular variation, because its time derivative is always positive in the situations considered, alternating large regions of slow and short regions of fast increases. The fast increase is also explained by the term in the denominator that goes close to zero when the eccentricity approaches 1.0.

The practical application of those results is that only near-circular orbits with inclination higher than the critical value are stable in the long range, since above this value the orbit looses its characteristics of near-circularity.

The Equatorial Orbits

Another property of the second-order system that comes directly from the inspection of the equations of motion is the existence of stable equatorial orbits. It means that if an orbit starts with i0 = 180°, the inclination and eccentricity remain constant and the orbits remain in the equatorial plane. In the second-order model this property is evident from the equations of motion. If i0 = 180°, then the right-hand sides of the expressions for and are also zero, because they are proportional to sin2(i) and sin(2i), respectively. For the fourth-order model, the expression for the time derivative for the inclination has a singularity, due to the sin(i) present in the denominator. The expression for the time derivative of the eccentricity is not zero, because the coefficients C3 and C6 that appear in the equation for f1(a,e,i,w), do not vanish when the inclination is zero. The numerical integration of the full restricted three-body problem model shows the existence of equatorial solutions (zero inclination all the time) also in this more general model, as expected due to the symmetry of the system. But, opposite to the second-order model, the eccentricity has short-period oscillations with an amplitude that depends on the initial eccentricity, but it can reach values large enough to destroy completely the circularity of the orbits.

Frozen Orbits

From the equations of motion for the second-order model it is possible to detect the existence of a family of special orbits that have constant semi-major axis, eccentricity, inclination and argument of periapsis called "frozen orbits." To obtain this family it is necessary to find the solutions of the equations:, since the time derivative for the semi-major axis is always zero. To satisfy the equation for the time derivative of the eccentricity the condition is sin(2w) = 0, which implies that cos(2w) = ±1. Analyzing the solution cos(2w) = +1 in the equation for the time derivative of the argument of periapsis, it is seen that its only possible solution is e = 1.0. So, only the solution cos(2w) = -1 is taken and it implies that w = 90° or w = 270°. The vanishing of the time derivative for the inclination does not add any new condition, since sin(2w) = 0 used above implies in di/dt = 0. From the equation and the assumed solution cos(2w) = -1 one more condition is found:

This is a relation between the eccentricity and the inclination that allows frozen orbits. From this equation the condition cos2(i)<3/5 is also obtained. This condition sets a minimum value for the inclination, that is the critical value discussed before, due to the restriction e2 > 0. Those conditions are valid only for the second-order model. When submitted to the fourth-order and to the full restricted problem model a frozen orbit is destroyed and it shows periodic oscillations in all the three orbital elements. The difference between the fourth-order and the restricted models are in the amplitude of the oscillations. The amplitudes are about ten times larger when the restricted problem is used. As an example, the frozen orbit with eccentricity 0.3 and inclination 137.63934 degrees was simulated and the oscillations for the fourth-order model were: 0.01 rad in inclination, 0.03 in eccentricity and 0.08 rad in the argument of periapsis. For the restricted model the oscillations were: 0.14 rad in inclination, 0.3 in eccentricity and the argument of periapsis assumed all values in the interval [0,2p] rad. Examining the equations of motion for the fourth-order model, it is seen that f1(a,e,i,w) = f2(a,e,i,w) = 0, but f3(a,e,i,w) ¹ 0. It means that the eccentricity and the inclination keep their time derivatives null for the fourth-order model. The destruction of the frozen orbits comes from the argument of periapsis, because its time derivative is not null, so the argument of periapsis deviates from the frozen condition, which makes sin(2w) ¹ 0 and it implies that all the time derivatives deviate from zero. Another application of the frozen orbits in the second-order model is to test the accuracy of the numerical integration of this model.

Effects of the Semi-Major Axis

The next step is to study the effects of the initial semi-major axis in the behaviour of the trajectories. Figure 4 shows the variation in the inclination for trajectories starting with initial inclination equals to 120 degrees and semi-major axis assuming the values 0.02, 0.06, 0.10, 0.12 and Fig. 5 shows the evolution of the eccentricity for the same trajectories.






This paper develops a mathematical model to study the third-body perturbation in a retrograde orbit: the double-averaged expanded to the fourth-order in Legendre polynomials. The results show in detail the behavior of retrograde orbits with respect to its initial inclination and the role of the critical inclination in the stability of near-circular orbits. They show that this critical value is a transition region where the eccentricity has an oscillation that increases in amplitude. It has also shown the existence of equatorial solutions (but with the eccentricity oscillating) and the non-existence of circular solutions under the dynamics given by the restricted three-body problem. The "frozen orbits" found in the double-averaged second-order model are studied in the fourth-order and the full restricted models. It is shown that they have their Keplerian elements disturbed by a small periodic oscillation in the fourth-order model and a large periodic oscillation in the restricted problem.



The author is grateful to CNPq (National Council for Scientific and Technological Development) - Brazil for the contract 300221/95-9 and to FAPESP (Foundation to Support Research in São Paulo State) for the contract 2000/14769-4. The author also acknowledges Roger Broucke, whose class notes provided the basis for the double-averaging method.



Blitzer, L., 1959, "Lunar-Solar Perturbations of an Earth Satellite," American Journal of Physics, Vol. 27, No. 9, pp. 634-645.        [ Links ]

Cook, G. E., 1962, "Luni-Solar Perturbations of the Orbit of an Earth Satellite," The Geophysical Journal, Vol. 6, No. 3, pp. 271-291.        [ Links ]

Giacaglia, G. E. O., 1973, "Lunar Perturbations of Artificial Satellites of the Earth," Smithsonian Astrophysical Observatory, Special Report 352, Cambridge, MA..        [ Links ]

Hough, M.E., 1981, "Orbits Near Critical Inclination, Including Lunisolar Perturbations," Celestial Mechanics, Vol. 25, No. 2, pp. 111-136.        [ Links ]

Kaula, M. W., 1962, "Development of the Lunar and Solar Disturbing Functions for a Close Satellite," Astronautical Journal, Vol. 67, pg. 300.        [ Links ]

Kozai, Y., 1959, "On the Effects of the Sun and the Moon upon the Motion of a Close Earth Satellite," Smithsonian Astrophysical Observatory, Special Report 22, Cambridge, MA..        [ Links ]

Kozai, Y., 1962, "Secular Perturbations of Asteroids with High Inclination and Eccentricity," The Astronomical Journal, Vol. 67, No. 9, pp. 591-598.        [ Links ]

Kozai, Y., 1965, "Effects of the Tidal Deformation on the Motion of Satellites," Publications of the Astronomical Society of Japan, Vol. 17, pp. 395-402.        [ Links ]

Kozai, Y., 1973, "A New Method to Compute Lunisolar Perturbations in Satellite Motions," Smithsonian Astrophysical Observatory, Special Report 349, Cambridge, MA..        [ Links ]

Musen, P., Bailie, A., and Upton, E, 1961, "Development of the Lunar and Solar Perturbations in the Motion of an Artificial Satellite," NASA-TN, D494.        [ Links ]

Prado, A.F.B.A., 2002, "Third-Body Perturbation in Orbits Around Natural Satellites". To be published in the Journal of Guidance, Control and Dynamics.        [ Links ]



Paper accepted July, 2005.



Technical Editor: Atila P. Silva Freire.

Creative Commons License Todo o conteúdo deste periódico, exceto onde está identificado, está licenciado sob uma Licença Creative Commons