Mathematical Methods Applied to the Celestial Mechanics of Artiﬁcial Satellites
View this Special IssueResearch Article  Open Access
Jarbas Cordeiro Sampaio, Rodolpho Vilhena de Moraes, Sandro da Silva Fernandes, "The Orbital Dynamics of Synchronous Satellites: Irregular Motions in the 2 : 1 Resonance", Mathematical Problems in Engineering, vol. 2012, Article ID 405870, 22 pages, 2012. https://doi.org/10.1155/2012/405870
The Orbital Dynamics of Synchronous Satellites: Irregular Motions in the 2 : 1 Resonance
Abstract
The orbital dynamics of synchronous satellites is studied. The 2 : 1 resonance is considered; in other words, the satellite completes two revolutions while the Earth completes one. In the development of the geopotential, the zonal harmonics and and the tesseral harmonics and are considered. The order of the dynamical system is reduced through successive Mathieu transformations, and the final system is solved by numerical integration. The Lyapunov exponents are used as tool to analyze the chaotic orbits.
1. Introduction
Synchronous satellites in circular or elliptical orbits have been extensively used for navigation, communication, and military missions. This fact justifies the great attention that has been given in literature to the study of resonant orbits characterizing the dynamics of these satellites since the 60s [1–14]. For example, Molniya series satellites used by the old Soviet Union for communication form a constellation of satellites, launched since 1965, which have highly eccentric orbits with periods of 12 hours. Another example of missions that use eccentric, inclined, and synchronous orbits includes satellites to investigate the solar magnetosphere, launched in the 90s [15].
The dynamics of synchronous satellites are very complex. The tesseral harmonics of the geopotential produce multiple resonances which interact resulting significantly in nonlinear motions, when compared to nonresonant orbits. It has been found that the orbital elements show relatively large oscillation amplitudes differing from neighboring trajectories [11].
Due to the perturbations of Earth gravitational potential, the frequencies of the longitude of ascending node and of the argument of pericentre can make the presence of small divisors, arising in the integration of equation of motion, more pronounced. This phenomenon depends also on the eccentricity and inclination of the orbit plane. The importance of the node and the pericentre frequencies is smaller when compared to the mean anomaly and Greenwich sidereal time. However, they also have their contribution in the resonance effect. The coefficients , , which define the argument in the development of the geopotential can vary, producing different frequencies within the resonant cosines for the same resonance. These frequencies are slightly different, with small variations around the considered commensurability.
In this paper, the 2 : 1 resonance is considered; in other words, the satellite completes two revolutions while the Earth carries one. In the development of the geopotential, the zonal harmonics and and the tesseral harmonics and are considered. The order of the dynamical system is reduced through successive Mathieu transformations, and the final system is solved by numerical integration. In the reduced dynamical model, three critical angles, associated to the tesseral harmonics and , are studied together. Numerical results show the time behavior of the semimajor axis, argument of pericentre and of the eccentricity. The Lyapunov exponents are used as tool to analyze the chaotic orbits.
2. Resonant Hamiltonian and Equations of Motion
In this section, a Hamiltonian describing the resonant problem is derived through successive Mathieu transformations.
Consider (2.1) to the Earth gravitational potential written in classical orbital elements [16, 17] where is the Earth gravitational parameter, m^{3}/s^{2}, , , , , , are the classical keplerian elements: is the semimajor axis, is the eccentricity, is the inclination of the orbit plane with the equator, is the longitude of the ascending node, is the argument of pericentre, and is the mean anomaly, respectively; is the Earth mean equatorial radius, km, is the spherical harmonic coefficient of degree and order , and are Kaula’s inclination and eccentricity functions, respectively. The argument is defined by where is the Greenwich sidereal time, ( is the Earth’s angular velocity, and is the time), and is the corresponding reference longitude along the equator.
In order to describe the problem in Hamiltonian form, Delaunay canonical variables are introduced,, , and represent the generalized coordinates, and , , and represent the conjugate momenta.
Using the canonical variables, one gets the Hamiltonian ,
with the disturbing potential given by
The argument is defined by and the coefficient is defined by
The Hamiltonian depends explicitly on the time through the Greenwich sidereal time . A new term is introduced in order to extend the phase space. In the extended phase space, the extended Hamiltonian is given by
For resonant orbits, it is convenient to use a new set of canonical variables. Consider the canonical transformation of variables defined by the following relations:
where are the modified Delaunay variables.
The new Hamiltonian , resulting from the canonical transformation defined by (2.9), is given by
where the disturbing potential is given by
Now, consider the commensurability between the Earth rotation angular velocity and the mean motion . This commensurability can be expressed as considering and as integers. The ratio defining the commensurability will be denoted by . When the commensurability occurs, small divisors arise in the integration of the equations of motion [9]. These periodic terms in the Hamiltonian with frequencies are called resonant terms. The other periodic terms are called short and longperiod terms.
The short and longperiod terms can be eliminated from the Hamiltonian by applying an averaging procedure [12, 18]: The variables and represent the short and longperiod terms, respectively, to be eliminated of the Hamiltonian .
The longperiod terms have a combination in the argument which involves only the argument of the pericentre and the longitude of the ascending node . From (2.10) and (2.11), these terms are represented by the new variables in the following equation:
The shortperiod terms are identified by the presence of the sidereal time and mean anomaly in the argument ; in this way, from (2.10) and (2.11), the term in the new variables is given by the following equations:
The term represents the other variables in the argument , including the argument of the pericentre and the longitude of the ascending node , or, in terms of the new variables, and , respectively.
A reduced Hamiltonian is obtained from the Hamiltonian when only secular and resonant terms are considered. The reduced Hamiltonian is given by
Several authors, [11, 15, 19–22], also use this simplified Hamiltonian to study the resonance.
The dynamical system generated from the reduced Hamiltonian, (2.16), is given by
The equations of motion , , and defined by (2.17) are
From (2.18) to (2.20), one can determine the first integral of the system determined by the Hamiltonian .
Equation (2.18) can be rewritten as
and substituting (2.21) and (2.22), one obtains
Now, (2.23) is rewritten as
In this way, the canonical system of differential equations governed by has the first integral generated from (2.24):
where is an integration constant.
Using this first integral, a Mathieu transformation
can be defined.
This transformation is given by the following equations: The subscript 1 denotes the new set of canonical variables. Note that , and the is an ignorable variable. So the order of the dynamical system is reduced in one degree of freedom.
Substituting the new set of canonical variables, , , , , , , , , in the reduced Hamiltonian given by (2.16), one gets the resonant Hamiltonian. The word “resonant” is used to denote the Hamiltonian which is valid for any resonance. The periodic terms in this Hamiltonian are resonant terms. The Hamiltonian is given by
The Hamiltonian has all resonant frequencies, relative to the commensurability , where the argument is given by
with
The secular and resonant terms are given, respectively, by and .
Each one of the frequencies contained in , , is related, through the coefficients , , to a tesseral harmonic . By varying the coefficients , , and keeping / fixed, one finds all frequencies concerning a specific resonance.
From , taking, , , , , and , one gets
The Hamiltonian is defined considering a fixed resonance and three different critical angles associated to the tesseral harmonic ; the critical angles associated to the tesseral harmonic have the same frequency of the critical angles associated to the with a difference in the phase. The other terms in are considered as shortperiod terms.
Table 1 shows the resonant coefficients used in the Hamiltonian .

Finally, a last transformation of variables is done, with the purpose of writing the resonant angle explicitly. This transformation is defined by
So, considering (2.31) and (2.32), the Hamiltonian is found to be with constant and
Since the term is constant, it plays no role in the equations of motion, and a new Hamiltonian can be introduced,
The dynamical system described by is given by
The zonal harmonics used in (2.34) and (2.35) and the tesseral harmonics used in (2.36) to (2.41) are shown in Table 2.

The constant of integration in (2.34) to (2.41) is given, in terms of the initial values of the orbital elements, , , and , by or, in terms of the variables and ,
In Section 4, some results of the numerical integration of (2.43) are shown.
3. Lyapunov Exponents
The estimation of the chaoticity of orbits is very important in the studies of dynamical systems, and possible irregular motions can be analyzed by Lyapunov exponents [23].
In this work, “GramSchmidt’s method,” described in [23–26], will be applied to compute the Lyapunov exponents. A brief description of this method is presented in what follows.
The dynamical system described by (2.43) can be rewritten as
Introducing
Equations (3.2) can be put in the form The variational equations, associated to the system of differential equations (3.3), are given by where is the Jacobian.
The total number of differential equations used in this method is , represents the number of the motion equations describing the problem, in this case four. In this way, there are twenty differential equations, four are motion equations of the problem and sixteen are variational equations described by (3.4)
The dynamical system represented by (3.3) and (3.4) is numerically integrated and the neighboring trajectories are studied using the GramSchmidt orthonormalization to calculate the Lyapunov exponents.
The method of the GramSchmidt orthonormalization can be seen in [25, 26] with more details. A simplified denomination of the method is described as follows.
Considering the solutions to (3.4) as , the integration in the time begins from initial conditions , an orthonormal basis.
At the end of the time interval, the volumes of the dimensional produced by the vectors are calculated by where is the outer product and is a norm.
In this way, the vectors are orthonormalized by GramSchmidt method. In other words, new orthonormal vectors are calculated, in general, according to
The GramSchmidt method makes invariant the dimensional subspace produced by the vectors in constructing the new dimensional subspace spanned by the vectors .
With new vector , the integration is reinitialized and carried forward to . The whole cycle is repeated over a longtime interval. The theorems guarantee that the dimensional Lyapunov exponents are calculated by [25, 26]:
The theory states that if the Lyapunov exponent tends to a positive value, the orbit is chaotic.
In the next section are shown some results about the Lyapunov exponents.
4. Results
Figures 1, 2, 3, and 4 show the time behavior of the semimajor axis, angle, argument of perigee and of the eccentricity, according to the numerical integration of the motion equations, (2.43), considering three different resonant angles together: , , and associated to , and three angles, , , and associated to , with the same frequency of the resonant angles related to the , but with different phase. The initial conditions corresponding to variables and are defined for , , and given in Table 3. The initial conditions of the variables and are 0° and 0°, respectively. Table 3 shows the values of corresponding to the given initial conditions.

Figures 5, 6, 7, and 8 show the time behavior of the semimajor axis, angle, argument of perigee and of the eccentricity for two different cases. The first case considers the critical angles , , and , associated to the tesseral harmonic , and the second case considers the critical angles associated to the tesseral harmonics and . The angles associated to the , , , and , have the same frequency of the critical angles associated to the with a different phase. The initial conditions corresponding to variables and are defined for , , and given in Table 4. The initial conditions of the variables and are 0° and , respectively. Table 4 shows the values of corresponding to the given initial conditions.

Analyzing Figures 5–8, one can observe a correction in the orbits when the terms related to the tesseral harmonic are added to the model. Observing, by the percentage, the contribution of the amplitudes of the terms , , and , in each critical angle studied, is about 1,66% up to 4,94%. In fact, in the studies of the perturbations in the artificial satellites motion, the accuracy is important, since adding different tesseral and zonal harmonics to the model, one can have a better description about the orbital motion.
Figures 9, 10, 11, and 12 show the time behavior of the semimajor axis, angle, argument of perigee and of the eccentricity, according to the numerical integration of the motion equations, (2.43), considering three different resonant angles together; , , and associated to and three angles , , and associated to . The initial conditions corresponding to variables and are defined for , , and given in Table 5. The initial conditions of the variables and are 0° and , respectively. Table 5 shows the values of corresponding to the given initial conditions.

Analyzing Figures 1–12, one can observe possible irregular motions in Figures 1–4, specifically considering values for m^{2}/s and m^{2}/s, and, in Figures 9–12, for m^{2}/s and m^{2}/s. These curves will be analyzed by the Lyapunov exponents in a specified time verifying the possible regular or chaotic motions.
Figures 13 and 14 show the time behavior of the Lyapunov exponents for two different cases, according to the initial values of Figures 1–4 and 9–12. The dynamical system involves the zonal harmonics and and the tesseral harmonics and . The method used in this work for the study of the Lyapunov exponents is described in Section 3. In Figure 13, the initial values for , , and are m^{2}/s and m^{2}/s, and , respectively. In Figure 14, the initial values for , , and are m^{2}/s and m^{2}/s, and , respectively. In each case are used two different values for semimajor axis corresponding to neighboring orbits shown previously in Figures 1–4 and 9–12.
Figures 13 and 14 show Lyapunov exponents for neighboring orbits. The time used in the calculations of the Lyapunov exponents is about 150.000 days. For this time, it can be observed in Figure 13 that , corresponding to the initial value km, tends to a positive value, evidencing a chaotic region. On the other hand, analyzing the same Figure 13, , corresponding to the initial value km, does not show a stabilization around the some positive value, in this specified time. Probably, the time is not sufficient for a stabilization in some positive value, or , initial value km, tends to a negative value, evidencing a regular orbit. Analyzing now Figure 14, it can be verified that , corresponding to the initial value km, tends to a positive value, it contrasts with , initial value km. Comparing Figure 13 with Figure 14, it is observed that the Lyapunov exponents in Figure 14 has an amplitude of oscillation greater than the Lyapunov exponents in Figure 13. Analyzing this fact, it is probable that the necessary time for the Lyapunov exponent , in Figure 14, to stabilize in some positive value is greater than the necessary time for the in Figure 13.
Rescheduling the axes of Figures 13 and 14, as described in Figures 15 and 16, respectively, the Lyapunov exponents tending to a positive value can be better visualized.
5. Conclusions
In this work, the dynamical behavior of three critical angles associated to the 2 : 1 resonance problem in the artificial satellites motion has been investigated.
The results show the time behavior of the semimajor axis, argument of perigee and eccentricity. In the numerical integration, different cases are studied, using three critical angles together: , , and associated to and , , and associated to the .
In the simulations considered in the work, four cases show possible irregular motions for m^{2}/s, m^{2}/s, m^{2}/s, and m^{2}/s. Studying the Lyapunov exponents, two cases show chaotic motions for m^{2}/s and m^{2}/s.
Analyzing the contribution of the terms related to the , it is observed that, for the value of m^{2}/s, the amplitudes of the terms , , and are greater than the other values of . In other words, for bigger values of semimajor axis, it is observed a smaller contribution of the terms related to the tesseral harmonic .
The theory used in this paper for the 2 : 1 resonance can be applied for any resonance involving some artificial Earth satellite.
Acknowledgments
This work was accomplished with support of the FAPESP under the Contract no. 2009/007355 and 2006/049976, SP Brazil, and CNPQ (Contracts 300952/20082 and 302949/20097).
References
 M. B. Morando, “Orbites de resonance des satellites de 24h,” Bulletin of the American Astronomical Society, vol. 24, pp. 47–67, 1963. View at: Google Scholar
 L. Blitzer, “Synchronous and resonant satellite orbits associated with equatorial ellipticity,” Journal of Advanced Robotic Systems, vol. 32, pp. 1016–1019, 1963. View at: Google Scholar
 B. Garfinkel, “The disturbing function for an artificial satellite,” The Astronomical Journal, vol. 70, no. 9, pp. 688–704, 1965. View at: Google Scholar
 B. Garfinkel, “Tesseral harmonic perturbations of an artificial satellite,” The Astronomical Journal, vol. 70, pp. 784–786, 1965. View at: Publisher Site  Google Scholar
 B. Garfinkel, “Formal solution in the problem of small divisors,” The Astronomical Journal, vol. 71, pp. 657–669, 1966. View at: Publisher Site  Google Scholar
 G. S. Gedeon and O. L. Dial, “Alongtrack oscillations of a satellite due to tesseral harmonics,” AIAA Journal, vol. 5, pp. 593–595, 1967. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 G. S. Gedeon, B. C. Douglas, and M. T. Palmiter, “Resonance effects on eccentric satellite orbits,” Journal of the Astronautical Sciences, vol. 14, pp. 147–157, 1967. View at: Google Scholar
 G. S. Gedeon, “Tesseral resonance effects on satellite orbits,” Celestial Mechanics, vol. 1, pp. 167–189, 1969. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. T. Lane, “An analytical treatment of resonance effects on satellite orbits,” Celestial Mechanics, vol. 42, pp. 3–38, 1988. View at: Google Scholar
 A. Jupp, “A solution of the ideal resonance problem for the case of libration,” The Astronomical Journal, vol. 74, pp. 35–43, 1969. View at: Publisher Site  Google Scholar
 T. A. Ely and K. C. Howell, “Longterm evolution of artificial satellite orbits due to resonant tesseral harmonics,” Journal of the Astronautical Sciences, vol. 44, pp. 167–190, 1996. View at: Google Scholar
 D. M. Sanckez, T. Yokoyama, P. I. O. Brasil, and R. R. Cordeiro, “Some initial conditions for disposed satellites of the systems GPS and galileo constellations,” Mathematical Problems in Engineering, vol. 2009, Article ID 510759, 22 pages, 2009. View at: Publisher Site  Google Scholar
 L. D. D. Ferreira and R. Vilhena de Moraes, “GPS satellites orbits: resonance,” Mathematical Problems in Engineering, vol. 2009, Article ID 347835, 12 pages, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. C. Sampaio, R. Vilhena de Moraes, and S. S. Fernandes, “Artificial satellites dynamics: resonant effects,” in Proceedings of the 22nd International Symposium on Space Flight Dynamics, São José dos Campos, Brazil, 2011. View at: Google Scholar
 A. G. S. Neto, Estudo de Órbitas Ressonantes no Movimento de Satélites Artificiais, Tese de Mestrado, ITA, 2006.
 J. P. Osorio, Perturbações de Órbitas de Satélites no Estudo do Campo Gravitacional Terrestre, Imprensa Portuguesa, Porto, Portugal, 1973.
 W. M. Kaula, Theory of Satellite Geodesy: Applications of Satellites to Geodesy, Blaisdel, Waltham, Mass, USA, 1966.
 A. E. Roy, Orbital Motion, Institute of Physics Publishing Bristol and Philadelphia, 3rd edition, 1988.
 P. H. C. N. Lima Jr., Sistemas Ressonantes a Altas Excentricidades no Movimento de Satélites Artificiais, Tese de Doutorado, Instituto Tecnológico de Aeronáutica, 1998.
 P. R. Grosso, Movimento Orbital de um Satélite Artificial em Ressonância 2:1, Tese de Mestrado, Instituto Tecnológico de Aeronáutica, 1989.
 J. K. S. Formiga and R. Vilhena de Moraes, “Dynamical systems: an integrable kernel for resonance effects,” Journal of Computational Interdisciplinary Sciences, vol. 1, no. 2, pp. 89–94, 2009. View at: Google Scholar
 R. Vilhena de Moraes, K. T. Fitzgibbon, and M. Konemba, “Influence of the 2:1 resonance in the orbits of GPS satellites,” Advances in Space Research, vol. 16, no. 12, pp. 37–40, 1995. View at: Publisher Site  Google Scholar
 F. Christiansen and H. H. Rugh, “Computing lyapunov spectra with continuous gramschmidt orthonormalization,” Nonlinearity, vol. 10, no. 5, pp. 1063–1072, 1997. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 L. QunHong and T. JieYan, “Lyapunov exponent calculation of a twodegreeoffreedom vibroimpact system with symmetrical rigid stops,” Chinese Physics B, vol. 20, no. 4, Article ID 040505, 2011. View at: Publisher Site  Google Scholar
 I. Shimada and T. Nagashima, “A numerical approach to ergodic problem of dissipative dynamical systems,” Progress of Theoretical Physics, vol. 61, no. 6, pp. 1605–1616, 1979. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 W. E. Wiesel, “Continuous time algorithm for Lyapunov exponents,” Physical Review E, vol. 47, no. 5, pp. 3686–3697, 1993. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2012 Jarbas Cordeiro Sampaio et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.