Some Initial Conditions for Disposed Satellites of the Systems GPS and Galileo Constellations

Abstract

Through the averaged equations we revisit theoretical and numerical aspects of the strong resonance that increases the eccentricity of the disposed objects of GPS and Galileo Systems. A simple view of the phase space of the problem shows that the resonance does not depend on the semi-major axis. Therefore, usual strategies of changing altitude (raising perigee) do not work. In this problem we search for a set of initial conditions such that the deactivated satellites or upper-stages remain at least for 250 years without penetrating in the orbits of the operational satellites. In the case that Moon's perturbation is not significant, we can identify, in the phase space, the regions where eccentricity reaches maximum and minimum values so that possible risks of collision can be avoided. This is done semi-analytically through the averaged system of the problem. Guided by this idea, we numerically found the (𝜔,Ω) values of the real unaveraged problem. In particular, for the Galileo case, the theoretical results predicted in the averaged system are in good agreement with numerical results. We also show that initial inclination of the Moon plays an important role in the search of these conditions.

1. Introduction

Broadly speaking, the GPS, GLONASS, and Galileo

[1] systems are satellite constellations which were designed mainly for positioning and navigation purposes. The first members of GPS (block I) originally were designed to have inclination of 63.4 degrees with respect to the equator, distributed in three orbital planes, each one separated from 120 degrees in the longitude of the node. The altitude is 20,200 km. The GLONASS members are similar, with slightly lower altitude (19,100 km, orbital period = 11:15 h). The European GALILEO system is still in construction, and the inclination of the satellites will be 55 or 56 degrees, with altitude of 23,615 km. All the three systems have rather similar altitudes. In order to avoid risks of collision and following American Govern instructions about debris mitigation, there are some recommendations that the disposal satellites and upper-stages should be deposited at least 500 km above or below the semisynchronous orbit [2].

In a constellation of a navigation system, the members must be kept under precise requirements of functionality. However, after some time, they have to be deactivated since some level of these requirements cannot be fulfilled for long time. The destination of these deactivated objects is a problem since they must be moved into some disposal regions in order to preclude collisions with operational members of the constellation. While these vehicles can be designed a priori to transport additional propellant (at some nonnegligible cost) to be used in some planned maneuvers to insert them in the disposal regions; the same is not true for the upper-stage. In some cases (block IIF of GPS system), due to design restrictions, this upper-stage cannot be easily guided to the disposal region. It must perform several operations after the satellite is injected in the constellation. All these operations change its final parameters [3]. Since the inclination of these vehicles is near to 55 or 56 degrees, the eccentricity suffers strong variations and even an initially circular orbit can become highly eccentric so that they can cross very easily the orbit of the operational satellites. What is interesting and also problematic is the fact that the rate of growing of the eccentricity is very sensitive to the initial parameters of the disposal orbit (eccentricity, argument of the perigee, and longitude of the node). In this work, based on a theoretical framework, we present a set of initial conditions (argument of the perigee, longitude of the node) for GPS and Galileo systems such that the disposed objects can remain at the least 250 years with small eccentricity (0.01 or 0.02) without causing any risk to the operational satellites.

The above strategy of keeping small eccentricity can generate some additional problem: after some time, the disposed vehicles will accumulate and a graveyard of these objects will be created. Therefore, a risk of collisions amongst themselves is a crucial problem, since the products of these extra collisions are almost untrackable fragments that may offer more risks to the operational elements of the constellation.

According to Jenkin and Gick [4], the strategy in the opposite direction, that is, exploiting the growth of the eccentricity in order to dilute disposal orbit collision risk has some interesting points to be considered: the percentage of disposed vehicles that will reenter in the atmosphere can be increased. Another advantage observed is, although eccentricity growth strategy increases the collision risk in the constellation, that in some cases this risk can be reversed with proper choice of the initial disposal eccentricity.

In this sense, we briefly started the investigation of some initial conditions that can cause large increase of the eccentricity, for a minimum time interval, considering also different initial inclinations of the Moon’s orbit (see Appendix B). Our calculations show that the growth of the eccentricity is rather sensitive to the Moon’s position (inclination and semi-major axis).

2. Methods

2.1. Disturbing Function of the Sun

As we want to highlight some theoretical aspects, it is instructive to write the main disturbing forces in terms of the orbital elements.

In this section we obtain the averaged disturbing function of the Sun. Following the classical procedure [5], in a reference center fixed in the Earth equator, the disturbing function of the Sun is 𝑅_⊙=𝑘²𝑀_⊙||𝑟⃗−1𝑟_⊙||−𝑟||𝑟⃗⋅_⊙𝑟||^⊙3,(2.1) where 𝑀_⊙ is the mass of the Sun, 𝑘² is the gravitational constant, and 𝑟⃗, 𝑟_⊙ are position vector of the satellite and the Sun, respectively.

Expanding (2.1) in powers of (𝑟/𝑟_⊙) up to order 2, we have

𝑅_⊙=𝑘²𝑀_⊙𝑟²𝑟^3⊙−12+32cos²(𝑆),(2.2) where S is the angular distance between the satellite and the Sun. We use the classical notations a, e, I, f, 𝜔, and Ω, for semi-major axis, eccentricity, inclination, true anomaly, argument of the perigee, and longitude of the node. The same set is used for the Sun’s elements, adding the index⊙.

In order to write cos(S) in terms of orbital elements, we have (see Figure 1)

Figure 1 

Geometry of the problem.

cos(𝑆)=𝑥𝑟𝑥𝑟^⊙_⊙+𝑦𝑟𝑦𝑟^⊙_⊙+𝑧𝑟𝑧𝑟^⊙_⊙.(2.3)

Considering classical relations of the two-body problem, we write cos(S) in terms of f, 𝑓_⊙, Ω, Ω_⊙, 𝜔, 𝜔_⊙, I, 𝐼_⊙ as follows:

cos(𝑆)=14(1+cos(𝐼))1−cos𝐼^⊙cos𝑓+𝜔+𝑓_⊙+1+𝜔_⊙+Ω−Ω_⊙4(1−cos(𝐼))1+cos𝐼^⊙cos𝑓+𝜔+𝑓_⊙+1+𝜔_⊙−Ω+Ω_⊙4(1+cos(𝐼))1+cos𝐼^⊙cos𝑓+𝜔−𝑓_⊙+1−𝜔_⊙+Ω−Ω_⊙4(1−cos(𝐼))1−cos𝐼^⊙cos𝑓+𝜔−𝑓_⊙+1−𝜔_⊙−Ω+Ω_⊙2sin(𝐼)sin𝐼^⊙cos𝑓+𝜔−𝑓_⊙−𝜔_⊙−cos𝑓+𝜔+𝑓_⊙+𝜔_⊙,(2.4) or in a compact form as follows:

cos(𝑆)=𝐴𝑎+𝐵𝑏+𝐶𝑐+𝐷𝑑+𝐸𝑒,(2.5) where

𝐴=1𝑎=c4o(s1+𝑓c+os𝜔(+𝐼)𝑓)_⊙1−cos𝐼^⊙,𝐵=1+𝜔_⊙+Ω−Ω_⊙,𝑏=c4o(s1−𝑓c+os𝜔(+𝐼)𝑓)_⊙1+cos𝐼^⊙,𝐶=1+𝜔_⊙−Ω+Ω_⊙,𝑐=c4o(s1+𝑓c+os𝜔(−𝐼)𝑓)_⊙1+cos𝐼^⊙,𝐷=1−𝜔_⊙+Ω−Ω_⊙,𝑑=4co(1s−𝑓c+os𝜔(𝐼−))𝑓_⊙1−cos𝐼^⊙,𝐸=1−𝜔_⊙−Ω+Ω_⊙,𝑒=c2ossin𝑓(𝐼+)s𝜔in−𝑓𝐼^⊙_⊙,−𝜔_⊙−cos𝑓+𝜔+𝑓_⊙+𝜔_⊙.(2.6)

In order to get rid of the short period variations, we have to obtain the averaged system and the rigorous procedure is to apply the classical von-Zeipel or Hori’s method [5, 6]. In the present case, as our interest is only to examine the long-term behavior, without retrieving the contribution of the short period terms eliminated during the averaging process, the secular disturbing function can be found simply considering [7, 8]

𝑅_⊙=21𝜋^02𝜋𝑅_⊙𝑑𝑙,(2.7) where l is the mean anomaly of the satellite. Consider

𝑅^∗⊙==𝑘𝑅_⊙2𝑀_⊙𝑎²2𝑟^3⊙×32𝑃𝐴²+𝐵²+𝐶²+𝐷²+2𝐸²−2+332𝐴²𝑍cos2𝜔+2𝑓_⊙+2𝜔_⊙+2Ω−2Ω_⊙+32𝐶²+3𝑍cos2𝜔−2𝑓_⊙−2𝜔_⊙+2Ω−2Ω_⊙2𝐵²𝑍cos2𝜔+2𝑓_⊙+2𝜔_⊙−2Ω+2Ω_⊙+32𝐷²+3𝑍cos2𝜔−2𝑓_⊙−2𝜔_⊙−2Ω+2Ω_⊙2𝑍𝐸²+2𝐶𝐷cos2𝜔−2𝑓_⊙−2𝜔_⊙+32𝑍𝐸²+3𝑍−𝐸₂+2𝐴𝐵cos2𝜔+2𝑓_⊙+2𝜔_⊙+3𝑃(𝐴𝐵++𝐶𝐴𝐷𝐷)c+os𝐵𝐶2Ωc−os2Ω(2_⊙𝜔)+3𝑃−𝐸²+𝐴𝐶+𝐵𝐷cos2𝑓_⊙+2𝜔_⊙+3𝐴𝐷𝑃cos2𝑓_⊙+3𝐴𝐶𝑍cos2𝜔+2Ω−2Ω_⊙+3𝐸𝑃(−𝐴−𝐵+𝐶++2𝜔𝐷_⊙)+co2sΩΩ−−2ΩΩ_⊙⊙+3𝐸𝑃(𝐴−𝐷)cos2𝑓_⊙+2𝜔_⊙+Ω−Ω_⊙−3𝐴𝐸𝑍cos2𝜔+2𝑓_⊙+3𝐸𝑍(𝐴−𝐶)cos2𝜔+Ω−Ω_⊙+3𝐵𝐷𝑍cos2𝜔−2Ω++22Ω𝜔_⊙⊙+Ω−Ω_⊙+3𝐵𝐶𝑃cos2𝑓_⊙+2𝜔_⊙−2Ω+2Ω_⊙+3𝐸𝑍(𝐵−𝐷)cos2𝜔−Ω+Ω+_⊙3𝐸𝑃(𝐵−𝐶)cos2𝑓_⊙+2𝜔_⊙−Ω+Ω_⊙+3𝐶𝐸𝑍cos2𝜔−2𝑓_⊙−3𝐵𝐸𝑍cos2𝑓_⊙+2𝜔+2𝜔_⊙−Ω+Ω_⊙−2𝜔_⊙+Ω−Ω_⊙+3𝐷𝐸𝑍cos2𝜔−2𝑓_⊙−2𝜔_⊙−Ω+Ω_⊙,(2.8) where 𝑃=1+(3/2)𝑒², 𝑍=(5/2)𝑒².

Performing a second and similar average with respect to the mean anomaly of the Sun, we get

𝑅_⊙=𝑘²𝑀_⊙𝑎²4𝑟^3⊙×𝑃41−3cos²+3(𝐼)−3cos²𝐼_⊙+9cos²(𝐼)cos²𝐼_⊙2𝑍sin²+3(𝐼)−1+3cos²𝐼_⊙cos(2𝜔)2𝑃sin²+3(𝐼)sin²𝐼_⊙cos2Ω−2Ω_⊙8𝑍(1+cos(𝐼))²−3sin²𝐼_⊙cos2𝜔+2Ω−2Ω_⊙+32𝑃𝑍ssinin((𝐼𝐼))csoisn(𝐼𝐼)_⊙sin(1+𝐼_⊙cos(𝐼))cos𝐼_⊙cos2𝜔+Ω−Ω_⊙+3cos𝐼_⊙cosΩ−Ω_⊙8𝑍1+cos²+3(𝐼)²sin²𝐼_⊙cos2𝜔−2Ω+2Ω_⊙2𝑍sin(𝐼)(1−cos(𝐼))sin𝐼^⊙cos𝐼_⊙cos2𝜔−Ω+Ω_⊙.(2.9) In the above expression, the orbit of the Sun is assumed to be a Keplerian circular orbit. The elliptic case is briefly discussed in Appendix A.

2.2. Oblateness Disturbing Function

For the oblateness, the disturbing function truncated at second order of 𝑅_𝑃/𝑟 is

𝑈₂=𝑘²𝑀_𝑇𝑅^2𝑃𝑟³𝐽₂12−32sin²(𝛽),(2.10) where 𝑅_𝑃, 𝐽₂, and 𝛽 are equatorial radius of the planet, oblateness coefficient, and geocentric latitude of the satellite, respectively. If we proceed in the exact same way as we did before, we have from the geometry of Figure 1

sin(𝛽)=sin(𝐼)sin(𝑓+𝜔).(2.11) Once β is eliminated, the average of 𝑈₂ with respect to the mean anomaly of satellite is [9, 10]

⟨𝑈₂⟩=12𝜋^02𝜋𝑅𝑈₂𝑑𝑙,_𝐽2=⟨𝑈₂⟩=14𝑛²𝐽₂𝑅^2𝑃3cos²(𝐼)−11−𝑒²^−3/2,(2.12) where n is the mean motion of the satellite.

2.3. Some Special Resonances

For close satellites, usually the oblateness is the dominant part. In this case, the main frequencies of the system are given by

𝜔̇≈3𝑛𝐽²𝑅^2𝑃4𝑎²1−𝑒²²Ω̇≈−3𝑛𝐽₂5cos²(𝐼)−1,𝑅^2𝑃2𝑎²1−𝑒²²cos(𝐼).(2.13)

The ratio of these two frequencies is

Ω̇𝜔̇≈1−2c5ocsos(𝐼²)(𝐼)=𝑘.(2.14)

Note that for k = integer, we have the special resonances which do not depend on the semi-major axis. These resonances usually affect the eccentricity [8]. For 𝑘=−2, we have 2𝜔̇+Ω̇≈0 for I = 56.06° or I = 110.99°. Another classical resonance occurs when I = 63.4°, so that 𝜔̇≈0.

3. Effects of 2𝜔̇+Ω̇ and 𝜔̇ Resonances

In order to see the effects of the resonances which affect GPS and Galileo satellites, the osculating equations of a satellite will be integrated. For the moment, as disturbers, we consider only the Sun and the oblateness (the complete Cartesian equations involving the remaining disturbers will be given in Section 4). Note that the resonant conditions to be used this time are extracted from the averaged system as presented in the precedent section (I = 56.06°, I = 63.4°). Figures 2 and 3 show the effects of both resonances on the eccentricity and on the resonant angles. Note that an initial small eccentricity reaches a significant increase.

(a)

(b)

Figure 2 

Time evolution of the eccentricity (a) and the critical angle (b). Initial conditions: a = 4.805 𝑅𝑇 (30,647 km), e = 0.005, I = 56.06°, and other elements equal to zero. In the simulations, we consider only Sun and oblateness as disturbers.

(a)

(b)

Figure 3 

Same as Figure 1. Initial conditions: a = 4.7 𝑅𝑇 (≈ 29,977 km), e = 0.005, I = 63.4°, 𝜔=90∘, and other elements equal to zero.

Let us pay more attention to the case I = 56.06° which is the inclination of the members of the Galileo constellation. For this inclination, the dominant term in the 𝑅^∗⊙ is cos(2𝜔+Ω−Ω_⊙). Neglecting the remaining terms of 𝑅_⊙, the Hamiltonian of the problem is

𝐹=𝑅_𝐽2+𝑘²𝑀_⊙𝑎²2𝑟^3⊙𝑃81−3cos²−3(𝐼)−3cos²𝐼_⊙+9cos²(𝐼)cos²𝐼_⊙4𝑍sin(𝐼)sin𝐼^⊙(1+cos(𝐼))cos𝐼_⊙cos2𝜔+Ω−Ω_⊙.(3.1) Let us take 𝐿=𝑘²(𝑀_𝑇+𝑚)𝑎, 𝐺=𝐿√1−𝑒², 𝐻=𝐺cos(𝐼), l = mean anomaly, 𝜔=𝑔, and Ω=ℎ the set of the Delaunay variables. After a trivial Mathieu canonical transformation [11],

(𝐺,𝐻,𝜔,Ω)⟶𝑃₁,𝑃₂,𝜃₁,𝜃₂,(3.2) where

𝜃₁=2𝜔+Ω,𝑃₁=𝐺2,𝜃²=Ω,𝑃₂=𝐻−𝐺2,(3.3) then we have

𝑅_⊙=𝑘²𝑀_⊙𝑎²2𝑟^3⊙⎡⎢⎢⎢⎢⎣𝑃81−3𝑃¹+𝑃₂²4𝑃²¹−3cos²𝐼_⊙+9𝑃¹+𝑃₂²4𝑃²¹−3cos²𝐼_⊙4𝑍1−𝑃¹+𝑃₂²4𝑃²¹^1/2sin𝐼_⊙1+𝑃¹+𝑃₂𝑅_𝐽22𝑃₁cos𝐼^⊙cos𝜃₁⎤⎥⎥⎦,=14𝑛²𝐽₂𝑅^2𝑃3𝑃¹+𝑃₂²4𝑃²¹−11−𝐿²−4𝑃²¹𝐿²^−3/2.(3.4) Since Sun’s orbit is a Keplerian one, we also considered 𝜔_⊙=0, Ω_⊙=0. In (P_i, θ_i) variables, our Hamiltonian is a one degree of freedom problem, whose dynamics is very similar to the very well-known Lidov-Kozai resonance. In Figure 4, we consider an initial eccentricity 𝑒₀=0.005 and semi-major axis 𝑎=4.805 𝑅_𝑇 (30,647 km). This figure is very instructive: note that in the bottom part there is a large region where the satellite remains some finite time with very small eccentricity. These are the exactly region we are looking for. It corresponds to the region where 2𝜔+Ω≈0. On the other hand, we have the counterpart of this situation at the top of the figure: very high eccentricity, which occurs again for 2𝜔+Ω≈0. We can separate these two configurations and have a close view of these two cases. Only to confirm our reasoning, let us integrate the problem in Cartesian coordinates. We also have to decrease the effect of the Moon’s perturbation since in the averaged analysis we considered only 𝑅_⊙ and 𝑅_𝐽2. To do that, we consider convenient value for the semi-major axis. Figure 5 (initial conditions: a = 3.5 𝑅_𝑇 (≈22323 km), e = 0.005, I = 56.06°, and other elements equal to zero; Moon inclination I_L = 18.28°) shows that the minimum of eccentricity occurs when 2𝜔+Ω crosses zero from positive to negative values (decreasing direction), while maximum occurs when 2𝜔+Ω crosses zero from negative to positive values (increasing direction). It is worth noting that if the semi-major axis is high, then the effect of the Moon cannot be neglected, so that the problem is no more a one degree of freedom problem. In this case the search of the (𝜔,Ω) pair, such that eccentricity remains small, must be done integrating the complete equations of the motion.

Figure 4 

Level curves of Hamiltonian (𝑅⊙+𝑅𝐽2), showing the eccentricity variation versus resonant angle.

Figure 5 

Time evolution of the osculating eccentricity (top) and the osculating critical angle (bottom) for a disposal GPS satellite. Note that the minimum of the eccentricity occurs when 2𝜔+Ω is crossing zero in decreasing direction (from positive to negative value) while maximum occurs when 2𝜔+Ω is crossing zero, but in increasing direction (from negative to positive value). Initial conditions: a = 3.5 RT (22,324 km), e = 0.005, I = 56.06°, and other elements equal to zero; IL = 28.58°.

4. Special (𝜔,Ω) Initial Conditions for Galileo Case

In the previous section, we considered only the effects of the Sun and of the oblateness. Moreover, in the presence of the resonance, the main effects are governed by the long term variations, so that we eliminated the short period terms. From a theoretical point of view, this averaged system is quite efficient to highlight the basic dynamics that affects the eccentricity of the GPS and Galileo satellites. However, for a more complete and realistic study, we need to include more disturbers.

In this section we want to find some special initial conditions such that the satellites can remain stable for at least 250 years with very small eccentricity without causing any risk of collision to the operational elements of the constellation. The strategy to search these particular initial conditions is guided from the theoretical approach described in the previous section.

In this section we integrate the osculating elements of a disposal satellite of the Galileo system under the effect of the Sun, Moon, and the oblateness.

The equations for the osculating elements (exact system), including Moon, are

𝑟̈⃗=−𝑘²(𝑀+𝑚)𝑟³𝑟⃗−𝑘²𝑀_⊙||𝑟𝑟⃗−⃗−𝑟𝑟^⊙_⊙||³−𝑟^⊙||𝑟_⊙||³−𝑘²𝑀_𝐿||𝑟𝑟⃗−⃗−𝑟𝑟^𝐿_𝐿||³−𝑟^𝐿||𝑟_𝐿||³+𝑃^𝐽2𝑃,(4.1)_𝐽𝑥=−𝑘²𝑀𝐽₂𝑅^2𝑃23𝑟𝑥⁵−152𝑧²𝑥𝑃𝑟⁷,(4.2)_𝐽𝑦=−𝑘²𝑀𝐽₂𝑅^2𝑃23𝑟𝑦⁵−152𝑧²𝑦𝑃𝑟⁷,(4.3)_𝐽𝑧=−𝑘²𝑀𝐽₂𝑅^2𝑃29𝑟𝑧⁵−152𝑧𝑟³⁷,(4.4) where 𝑃_𝐽2 is the acceleration due to the oblateness, whose Cartesian components are given by 𝑃_𝐽𝑥, 𝑃_𝐽𝑦, 𝑃_𝐽𝑧 [10]. The second and third terms in (3.3) are the contributions of the Sun and the Moon, respectively, and M, m, M_L are the masses of the Earth, satellite, and Moon, respectively. The position vectors of the satellite, Sun, and Moon are indicated by 𝑟⃗, 𝑟_⊙, and 𝑟_𝐿.

As we said before, we take 500 km above of the nominal altitude of the constellation. The initial elements are fixed to a = 4.805R_T (30,647 km), e = 0.005, l = 0°, and I = 56.06°. We consider two cases for the Moon’s inclination, I = 18.28° and 28.58°. We show that the initial value of the inclination is important as shown in Figures 6 and 7. In these figures, we show the pair (𝜔,Ω) such that the disposal object remains at least 250 years with eccentricity smaller than 0.01, so that there is no risk of collision with any member of the constellation. The black region corresponds to initial conditions such that the maximum eccentricity is less than 0.01. In the green region, the maximum eccentricity is less than 0.02. The two straight lines represent the exact condition 2𝜔+Ω=𝑘𝜋 (in particular we only plot the case 𝑘=0). Note that, in special, the black dots (Figure 6) are in fact located in places predicted from the previous theoretical model. For the remaining figures, the black dots are slightly shifted (upward) from the line 2𝜔+Ω=0. We believe that this is caused by the strong perturbation of the Moon. Figure 8 shows the time evolution of the eccentricity for integration whose initial conditions are obtained from Figure 6 (small square in the bottom). As expected, the eccentricity remains very low, while if we take (𝜔,Ω) outside the marked regions in Figure 6 or Figure 7, a significant increase is verified as shown in Figure 9. The initial conditions (𝜔,Ω) used in this case correspond to the red circle given in Figure 6.

Figure 6 

Black dots: represent (𝜔,Ω) values such that a satellite with a = 30,647 km remains at least 250 years with 𝑒MAX≤0.01. Green dots: the same, but 𝑒MAX≤0.02. Blue dots: curve satisfying 2𝜔+Ω=0. Moon’s inclination: IL = 18.28°. Note that most of the “stable” (black dots) (𝜔,Ω) points satisfy 2𝜔+Ω=0, 2𝜋 with Ω≈0, 𝜔=𝜋.

Figure 7 

Same of Figure 6, but now 𝐼𝐿=28.58∘. Note that the distribution of the stable (𝜔,Ω) is very sensitive under changes in 𝐼𝐿.

(a)

(b)

Figure 8 

Time evolution of the eccentricity (a) and the critical angle (b) obtained from integration of (4.1). Initial conditions: a = 4.805 𝑅𝑇, e = 0.005, I = 56.06°. Initial (𝜔,Ω): (24°, 0°): black, (23°, 2°): red, (18°, 8°): blue. These initial conditions were extracted from the red square shown in Figure 6. Initial conditions of the Moon: aL = 380, 367.2 km, eL = 0.0276, IL = 18.28°, Ω𝐿=12.11∘, 𝜔𝐿=92∘, and 𝑙𝐿=337∘. As expected, the eccentricity remains less than 0.005 for at least 250 years.

(a)

(b)

Figure 9 

Same as Figure 8 but considering initial conditions from the red small circle of Figure 6. Initial conditions: a = 4.805 𝑅𝑇 (30,647 km), e = 0.005, I = 56.06°, 𝜔= 90°, Ω= 180°. This time eccentricity grows fastly since the corresponding initial condition was not taken from the black region of Figure 6.

5. (𝜔,Ω) Conditions for GPS Case

Here let us consider the GPS system. Again we take I = 18.28 and I = 28.58 for the Moon’s inclination. As before, the importance of the Moon’s inclination is clear.

6. Tesseral and Sectorial Harmonics

Up to now, we have not considered tesseral and sectorial harmonics. Since GPS satellites have orbital period near to 12:00 h, the inclusion of such harmonics must be examined when drawing the figures of Section 5. While the numerical values of the coefficients of these harmonics are very small compared to the zonal harmonics, due to the 1n  :  2𝛾 resonance (𝛾 is the rotation velocity of the Earth), this contribution could cause some nonnegligible effects. We show very briefly these perturbations. There are several models of the geopotential, including most sophisticated recursive formulae to generate very high-order (JGM-3 [12], EGM96S, EGM96 [13], etc.). In this work, we do not need high order model, so that we consider only few terms. The disturbing function for the general geopotential can be written in the form [14, 15]

𝑉=𝑘²𝑀_𝑇𝑟+𝑘²𝑀_𝑇𝑟−^∞_𝑛=2𝑅𝑟^𝑃^𝑛𝐽_𝑛𝑃_𝑛0(sin(𝜙))+^∞_𝑛𝑛=2_𝑚=1𝑅𝑟^𝑃^𝑛𝐽_𝑛𝑚cos𝑚𝜆−𝜆_𝑚𝑛𝑃_𝑛𝑚(sin(𝜙))(6,.1) where 𝐽_𝑛𝑚,𝜆_𝑛𝑚 are numerical coefficients and P_nm are the associated functions of Legendre.

Figure 12 describes the fundamental axes of the following reference: the potential (6.1) is referred to the (x, y, z) which is an equatorial system fixed on the Earth; therefore, it rotates with respect to (X, Y, Z) which is an inertial system, as follows:

𝑟⃗∶the position vector of the satellite.I:inclination.ϕ:geocentric latitude.λ:longitude.Ω:longitude of node.

According to Figure 12, some simple geometrical relations can be written as

cos(2𝜆)=𝑥𝑦𝑥=₂=𝑧𝑟𝑟c=cooss𝑟s((i𝜙𝜙n))(cs𝜙oin)s,((𝜆𝜆)),,−𝑦²𝑥²+𝑦²,sin(2𝜆)=2𝑥𝑦𝑥²+𝑦²cos₂,(𝜆)=𝑥²𝑥²+𝑦²,sin²(𝜆)=𝑦²𝑥²+𝑦²𝑥₂,+𝑦²=𝑟²cos²(𝜙).(6.2)

Let us define 𝜑_𝑛𝑚=(𝑘²𝑀_𝑃/𝑟)(𝑅_𝑃/𝑟)^𝑛𝐽_𝑛𝑚cos𝑚(𝜆−𝜆_𝑛𝑚)𝑃_𝑛𝑚(sin(𝜙)).

Therefore,

𝜑₂₂=3𝑘²𝑀_𝑃𝑅^2𝑃𝑟⁵𝐽₂₂𝑥²−𝑦²cos2𝜆₂₂+2𝑥𝑦sin2𝜆₂₂.(6.3)

Proceeding in the similar way, we obtain

𝜑₃₂=15𝑘²𝑀_𝑃𝑅^3𝑃𝑟⁷𝜑𝐽₃₂𝑥²−𝑦²𝑧cos2𝜆₃₂+2𝑥𝑦𝑧sin2𝜆₃₂₃₃,=15𝑘²𝑀_𝑃𝑅^3𝑃𝑟⁷𝐽₃₃𝑥𝑥²−3𝑦²cos3𝜆₃₃+𝑦3𝑥²−𝑦²sin3𝜆₃₃.(6.4)

Note that the zonal terms (m = 0) were already considered before. In Figure 12 the system (x, y, z) is fixed on the Earth while (X, Y, Z) is an inertial system, so that (x, y, z) rotates with respect to (X, Y, Z). Therefore, to have 𝜑_𝑛𝑚 referred to (X, Y, Z), we consider the trivial rotation:

⎛⎜⎜⎝𝑥𝑧𝑦⎞⎟⎟⎠=⎛⎜⎜⎝−cosisn0(𝜃(𝜃))csoins0((𝜃𝜃))010⎞⎟⎟⎠⎛⎜⎜⎝𝑋𝑍𝑌⎞⎟⎟⎠,(6.5) where 𝜃=𝛾𝑡+𝜃₀, 𝛾 = 2𝜋/day.

Therefore considering only 𝜑₂₂, 𝜑₃₂, and 𝜑₃₃, the force to be added in (4.1) is

𝑃_𝜑=grad_𝑋𝑌𝑍𝜑₂₂+𝜑₃₂+𝜑₃₃.(6.6)

Once we have introduced XYZ system, it is clear that several angular combinations like 2𝜃−𝜆 appear in the above 𝜑_𝑛𝑚, when (6.3) and (6.4) are expressed in the classical orbital elements. Due to the nominal semi-major axis of the GPS satellites, 2𝜃−𝜆 will generate long-term variations due to their proximity of the 2  :  1 exact resonance. In fact, the orbital period of GPS satellite is about 12:00 h.

To express 𝜑_𝑛𝑚 in terms of the orbital elements is trivial; however, the best way to see the effects is to keep these terms in Cartesian coordinates without any expansion in eccentricity or inclination, as the former increases to high values, while the second is essentially high from the beginning.

The effects of these additional terms are shown in Figure 13. The interaction of several arguments coming from 𝜑₂₂,𝜑₃₂,𝜑₃₃, and also those related to the lunar disturbing function, give rise to new resonant combinations. It is interesting that the presence of different 𝜑_𝑛𝑚 cause significative differences in the behavior of the eccentricity, however, in the beginning (up to t = 370 years), all the curves are almost coincident. In particular, our numerical experiments show that if eccentricity remains small, the effects of 𝜑_𝑛𝑚 are not significant.

Figures 14 and 15 show that the (Ω, 𝜔) initial conditions are similar to those figures of Section 4. Note that differences when 𝜑_𝑛𝑚 are included are almost negligible as expected according to what we learned from Figure 13. However, in the strategy of exploiting the growth of the eccentricity, the initial (Ω, 𝜔) can be sensitive, depending on the time integration and on the number of the harmonics considered in the geopotential. In Appendix B, we show similar figures, where (Ω, 𝜔) values are the initial conditions that cause fast increase of the eccentricity in less than 250 years.

7. Conclusion

With the averaged equations, we showed the dynamics of the 2𝜔+Ω resonance. The reason of the increase of the eccentricity is essentially due to this resonance which does not depend on the value of the semi-major axis. Therefore, any change of the semi-major axes (raising the perigee) of the decommissioned object will not remove from the resonance. After showing the existence of some initial conditions in the (𝜔,Ω) domain where the eccentricity can remain very small based on the averaged simplified model, we used the complete set of equations to search this pair in (𝜔,Ω) plane. The importance of the Moon’s inclination becomes clear as shown in Figures 6, 7, 10, and 11. We obtained these initial values for GALILEO and GPS systems. For completeness, we also derived a first-order averaged system in the eccentricity of the third body (Appendix A). Then several additional resonances appear although their effects are not so relevant for the navigation system. The search of the (𝜔,Ω) pair for the maximum increase of the eccentricity can be done straightforward following the same procedure we used for minimum eccentricity. For completeness, in the disturbing function of the geopotential we also included terms coming from 𝐽₂₂, 𝐽₃₂, and 𝐽₃₃. However, their contributions for the first 370 years are not significant. In a separate paper, we intend to show the corresponding Figures 10 and 11 considering the second strategy of exploiting the increase of the eccentricity. Our experiments show that the effect of 𝐽₂₂, 𝐽₃₂, and 𝐽₃₃ terms becomes more visible only after some hundred years. Therefore, their contributions in Figures 14 and 15 are almost negligible. However, considering long-time integration, their effects and interaction with solar and lunar perturbations become important.

Figure 10 

Same of Figure 6, but now a = 26,060 km.

Figure 11 

Same of Figure 10, but now IL = 28.58°.

Figure 12 

Geometry of the problem.

(a)

(b)

Figure 13 

Time evolution of eccentricity of GPS satellite. Initial conditions: a = 27,059.74 km, e = 0.005, 𝜔= 58°, Ω= 154°. The initial conditions of the Moon correspond to the epoch August 1, 2001. The contribution of individual and all 𝜑𝑛𝑚 is shown through different colors. The differences appear only after 370 years.

Figure 14 

Same of Figure 10, including 𝜑22, 𝜑32, and 𝜑33.

Figure 15 

Same of Figure 14, but now IL = 28.58°.

Figure 16 

Black dots represent (𝜔,Ω) values such that a satellite with a = 26,559.74 km reaches 𝑒≥0.6 in 250 years. Green dots: the same, but eccentricity reaches 𝑒≥0.5 in 250 years. Blue dots: curve satisfying 2𝜔+Ω=270. Moon’s inclination: IL = 18.28°.

Figure 17 

Same of Figure 16, but now IL = 28.58°.

(a)

(b)

Figure 18 

Time evolution of eccentricity (a) and critical angle (b) of a GPS satellite with a = 26,559.74 km whose eccentricity reaches the value of 0.6 after 200 years. Initial conditions: e = 0.005, I = 56.06°, 𝜔= 0°, Ω= 270°, IL = 18.28°.

Appendices

A.  Solar Disturbing Function (Up to First Order in 𝑒)

Here we give the complete expression of the averaged disturbing function up to first order in eccentricity of the third body:

𝑅^1⊙=𝑘²𝑀^3⊙𝑎³2𝑎^3⊙𝑒_⊙×92𝑃𝐴²+𝐵²+𝐶²+𝐷²+2𝐸²−23cos𝑙^⊙+94𝐴²+9𝑍cos2𝜔+𝑙_⊙+2𝜔_⊙+2Ω−2Ω_⊙4𝐴²𝑍cos2𝜔+3𝑙_⊙+2𝜔_⊙+2Ω−2Ω_⊙+94𝐵²+9𝑍cos2𝜔+𝑙_⊙+2𝜔_⊙−2Ω+2Ω_⊙4𝐵²𝑍cos2𝜔+3𝑙_⊙+2𝜔_⊙−2Ω+2Ω_⊙+94𝐶²+9𝑍cos2𝜔−3𝑙_⊙−2𝜔_⊙+2Ω−2Ω_⊙4𝐶²𝑍cos2𝜔−𝑙_⊙−2𝜔_⊙+2Ω−2Ω_⊙+94𝐷²+9𝑍cos2𝜔−3𝑙_⊙−2𝜔_⊙−2Ω+2Ω_⊙4𝐷²𝑍cos2𝜔−𝑙_⊙−2𝜔_⊙−2Ω+2Ω_⊙+94𝑍𝐸²+9+2𝐶𝐷cos2𝜔−3𝑙_⊙−2𝜔_⊙4𝑍𝐸²+2𝐶𝐷cos2𝜔−𝑙_⊙−2𝜔_⊙+94𝑍𝐸²+9+2𝐴𝐵cos2𝜔+𝑙_⊙+2𝜔_⊙4𝑍𝐸²+2𝐴𝐵cos2𝜔+3𝑙_⊙+2𝜔_⊙+92𝑍−𝐸²+9+𝐴𝐷+𝐵𝐶cos𝑙_⊙−2𝜔2𝑍−𝐸²+𝐴𝐷+𝐵𝐶cos𝑙_⊙+2𝜔+92𝑃−𝐸²+9+𝐴𝐶+𝐵𝐷cos𝑙_⊙+2𝜔_⊙2𝑃−𝐸²+𝐴𝐶+𝐵𝐷cos3𝑙_⊙+2𝜔_⊙+92𝑃(𝐴𝐵+𝐶𝐷)cos𝑙^⊙+9−2Ω+2Ω_⊙2𝑃(𝐴𝐵+𝐶𝐷)cos𝑙^⊙+2Ω−2Ω_⊙+92𝐴𝐶𝑍cos𝑙^⊙+9−2𝜔−2Ω+2Ω_⊙2𝐴𝐶𝑍cos𝑙^⊙+2𝜔+2Ω−2Ω_⊙+92𝐴𝐷𝑃cos𝑙^⊙+9+2𝜔_⊙+2Ω−2Ω_⊙2𝐴𝐷𝑃cos3𝑙^⊙+2𝜔_⊙+2Ω−2Ω_⊙+92𝐸𝑃(𝐴−𝐷)cos𝑙^⊙+9+2𝜔_⊙+Ω−Ω_⊙2𝐸𝑃(𝐴−𝐷)cos3𝑙^⊙+2𝜔_⊙+Ω−Ω_⊙+92𝐸𝑃(−𝐴−𝐵+𝐶+𝐷)cos𝑙^⊙+9−Ω+Ω_⊙2𝐸𝑃(−𝐴−𝐵+𝐶+𝐷)cos𝑙^⊙+Ω−Ω_⊙+92𝐸𝑍(𝐴−𝐶)cos𝑙^⊙+9−2𝜔−Ω+Ω_⊙2𝐸𝑍(𝐴−𝐶)cos𝑙^⊙+2𝜔+Ω−Ω_⊙−92𝐴𝐸𝑍cos2𝜔+𝑙^⊙−9+2𝜔_⊙+Ω−Ω_⊙2𝐴𝐸𝑍cos2𝜔+3𝑙^⊙+2𝜔_⊙+Ω−Ω_⊙+92𝐵𝐶𝑃cos𝑙^⊙+9+2𝜔_⊙−2Ω+2Ω_⊙2𝐵𝐶𝑃cos3𝑙^⊙+2𝜔_⊙−2Ω+2Ω_⊙+92𝐵𝐷𝑍cos𝑙^⊙+9−2𝜔+2Ω−2Ω_⊙2𝐵𝐷𝑍cos𝑙^⊙+2𝜔−2Ω+2Ω_⊙+92𝐸𝑃(𝐵−𝐶)cos𝑙^⊙+9+2𝜔_⊙−Ω+Ω_⊙2𝐸𝑃(𝐵−𝐶)cos3𝑙^⊙+2𝜔_⊙−Ω+Ω_⊙+92𝐸𝑍(𝐵−𝐷)cos𝑙^⊙+9−2𝜔+Ω−Ω_⊙2𝐸𝑍(𝐵−𝐷)cos𝑙^⊙+2𝜔−Ω+Ω_⊙−92𝐵𝐸𝑍cos𝑙^⊙−9+2𝜔+2𝜔_⊙−Ω+Ω_⊙2𝐵𝐸𝑍cos3𝑙^⊙+2𝜔+2𝜔_⊙−Ω+Ω_⊙+92𝐶𝐸𝑍cos2𝜔−3𝑙^⊙+9−2𝜔_⊙+Ω−Ω_⊙2𝐶𝐸𝑍cos2𝜔−𝑙^⊙−2𝜔_⊙+Ω−Ω_⊙+92𝐷𝐸𝑍cos2𝜔−3𝑙^⊙+9−2𝜔_⊙−Ω+Ω_⊙2𝐷𝐸𝑍cos2𝜔−𝑙^⊙−2𝜔_⊙−Ω+Ω_⊙.(A.1)

Note that all the cosines in the above relations have 𝑙_⊙ in the argument. Since the inclination of satellite is fixed about 55° or 56°, each possible resonance can occur for one particular value of the semi-major axis. However, considering the nominal altitude of the GPS and Galileo satellites, and since 𝑒_⊙ is small, none of the above combinations of angles is significant.

B. The Eccentricity Increasing Strategy

Following the same strategy to obtain Figures 10, 11, 14, and 15, we obtain similar figures, but now the pair (𝜔, Ω) represents the initial condition of an orbit whose eccentricity attains 𝑒≥0.5 in less than 250 years. As before, all orbits start with initial a = 26,559.74 km, e = 0.005, and 𝐼=56.06^∘. The maximum eccentricity usually is reached after t = 200 years.

Acknowledgments

The authors thank CNPQ, FAPESP, and FUNDUNESP. An anonymous referee is gratefully thanked for very useful comments.