Third-Body Perturbation in the Case of Elliptic Orbits for the Disturbing Body

Abstract

This work presents a semi-analytical and numerical study of the perturbation caused in a spacecraft by a third-body using a double averaged analytical model with the disturbing function expanded in Legendre polynomials up to the second order. The important reason for this procedure is to eliminate terms due to the short periodic motion of the spacecraft and to show smooth curves for the evolution of the mean orbital elements for a long-time period. The aim of this study is to calculate the effect of lunar perturbations on the orbits of spacecrafts that are traveling around the Earth. An analysis of the stability of near-circular orbits is made, and a study to know under which conditions this orbit remains near circular completes this analysis. A study of the equatorial orbits is also performed.

1. Introduction

An extensive literature is dedicated to the question of the perturbation done by a third-body on the trajectory of a spacecraft. Sptizer [1] used the lunar theory of Hill-Brown to study the problem on the limited case of small eccentricities and inclinations between the disturbing and the disturbed bodies. Kozai [2] derived the principal secular and long-period terms of the disturbing function due to the lunisolar gravitational attractions and Musen [3] included the parallactic term in the disturbing function. Using methods of classical mechanics, only for the secular terms, Blitzer [4] got estimate for the lunisolar disturbances. In the following years, many authors studied this problem using the disturbing function, Lagrange planetary equations, and numerical approximations [5–11]. In the references presented, many results were obtained related to the perturbations of the third body on a small mass moving close to another body and showed important analytical contributions, because their results are concentrated in derivation of equations.

Recently, many works have been presented in the literature based on numerical methods. Broucke [12] calculated the general form of the disturbing function of the third body truncated after the term of second order in the expansion in Legendre polynomials. This research, published in 2003, described the problem of the third-body perturbation on a satellite in a simplified approximated model using double average over the short period of the satellite as well as with respect to the distant perturbing body [13]. The important reason for this is to eliminate the terms due to the short time periodic motion of the spacecraft and to show smooth curves for the evolution of the mean orbital elements for a long-time period. The perturbing body was in a circular orbit in the (𝑥,𝑦) plane.

Prado and Costa [14] calculated the perturbing potential for order up to four in terms of the Legendre polynomials and later they extended the calculations to consider effects of up to order eight in the expansions in Legendre polynomials [15]. In Solórzano [16], the motion of the spacecraft is studied under the single-averaged analytical model with the disturbing function expanded in Legendre polynomials up to fourth order. The single average is taken over the mean motion of the satellite to eliminate short-period perturbations that appear in the trajectories. After that, the equations of motion are obtained from the planetary equations. Prado [17] studied this problem under the double-averaged analytical model with the disturbing function expanded in Legendre polynomials up to fourth order and the full restricted problem. The double average is taken over the mean motion of the satellite and over the mean motion of the disturbing body.

We noticed in the literature that few general expressions for calculations of the perturbing body have been done for cases of elliptic orbits [18–20]. There are no formulations that explicitly include the eccentricity of the disturbing body.

In the present work, we expanded the study done by Broucke [13], and Prado [17], including the eccentricity of the perturbing body. The double-averaged analytical model with the disturbing function expanded in Legendre polynomials up to second order. Our purpose was to study the stability of orbits of a massless spacecraft in a near-circular three-dimensional orbit around a central body with mass 𝑚₀, and a second body with mass 𝑚^′ in an elliptic orbit around this same central body in the plane 𝑥−𝑦.

This paper is structured as follows. In Section 2, we present the equations of motion used for the numerical simulations. Section 3 is devoted to the analysis of the numerical results of near-circular orbits. The theory developed here is used to study the behavior of a spacecraft around the Moon, where the Earth is the disturbing body. Several plots show the time histories of the Keplerian elements of the orbits involved. Our final comments are presented in Section 4.

2. Equations of Motion

In this section, we present the equations of motion obtained from the mathematical models used in this research. It is assumed that the main body with mass 𝑚₀ is fixed in the center of the reference system 𝑥−𝑦. The perturbing body, with mass 𝑚^′, is in an elliptic orbit with semimajor axis 𝑎^′, eccentricity 𝑒^′, and mean motion 𝑛^′ (given by the expression 𝑛^′2𝑎^′3=𝐺[𝑚₀+𝑚^′]). The massless spacecraft 𝑚 is in a generic two-dimensional orbit which orbital elements are 𝑎 (semimajor axis), 𝑒 (eccentricity), 𝑖 (inclination), 𝜔 (argument of periapsis), Ω (longitude of the ascending node), and 𝑛 (mean motion) given by the expression 𝑛²𝑎³=𝐺𝑚₀. This system is showen in Figure 1. Here, 𝐺 is the gravitational constant, 𝑟 and 𝑟^′ are the radius vectors of bodies 𝑚 and 𝑚^′, and 𝑆 is the angle between these radius vectors.

Figure 1 

Illustration of the dynamical system [16].

Using the traditional expansion in Legendre polynomials (assuming that 𝑟^′≫𝑟), the disturbing potential is given by [21]𝜇𝑅=^′𝐺(𝑚₀+𝑚^′)𝑟^′∞_𝑛=2(𝑟𝑟^′)^𝑛𝑃_𝑛[cos(𝑆)], (2.1)where 𝜇^′=𝑚^′/(𝑚₀+𝑚^′),𝐺 is the gravitational constant, 𝑃_𝑛 are the Legendre polynomials, and 𝑆 is the central elongation between the perturbed body (the spacecraft) and the perturbing body (the second body).

The terms 𝑛=1 and 𝑛=0 are not included in (2.1). The part 𝑅₂ of the disturbing function due to the second body is given by𝑅₂=𝜇^′𝑛^′2𝑎²2(𝑎^′𝑟^′)³(𝑟𝑎)²[3cos²(𝑆)−1], (2.2)where cos(𝑆) iscos(𝑆)=𝛼cos(𝑓)+𝛽sin(𝑓), (2.3) where 𝛼=𝑃⋅𝑟^′̂,𝛽=𝑄⋅𝑟^′̂,𝑓 is the true anomaly of the satellite, 𝑟^′̂ is the unit vector pointing from the central body to the disturbing body, all of them functions of 𝑓^′ and 𝜔^′ (true anomaly and argument of the periapsis of the disturbing body, resp.). The usual orthogonal unit vectors 𝑃 and 𝑄 are functions of 𝑖, 𝜔, and Ω in the plane of the satellite orbit, with 𝑃 pointing toward the periapsis. For the case of elliptic orbits the products 𝛼 and 𝛽 are written as𝛼=cos(𝜔)cos(Ω−𝑓^′−𝜔^′)−cos(𝑖)sin(𝜔)sin(Ω−𝑓^′−𝜔^′),𝛽=−sin(𝜔)cos(Ω−𝑓^′−𝜔^′)−cos(𝑖)cos(𝜔)sin(Ω−𝑓^′−𝜔^′). (2.4)

A substitution using (2.3) was made in (2.2). Then we have𝑅₂=𝜇^′𝑛^′2𝑎²2(𝑎^′𝑟^′)³(𝑟𝑎)²[3(𝛼²cos²(𝑓)+2𝛼𝛽cos(𝑓)sin(𝑓)+𝛽²sin²(𝑓))−1]. (2.5)

Now we need to average those quantities over the short period of the satellite motion as well as with respect to the distant perturbing body. The standard definition for average used is1⟨𝐹⟩=2𝜋₀^2𝜋𝐹𝑑𝑀, (2.6)where 𝑀 is the mean anomaly that is proportional to time.

The averages are realized in terms of the eccentric anomaly, and then it is necessary to replace the true anomalies 𝑓 and 𝑓^′ by the eccentric anomalies 𝐸 and 𝐸^′. To do this task, we use some well-known relations from the Celestial Mechanics [21] given by the following equations:sin(𝑓)=(1−𝑒²)^1/21−𝑒cos(𝐸)sin(𝐸),cos(𝑓)=cos(𝐸)−𝑒,(𝑟𝑒cos(𝐸)−1𝑎)=1−𝑒cos(𝐸),𝑑𝑀=(1−𝑒cos(𝐸))𝑑𝐸. (2.7)Using those equations, we can obtain the following relations:𝑟⟨(𝑎)²cos²(𝑓)⟩=(1+4𝑒²)2,𝑟⟨(𝑎)²sin²(𝑓)⟩=(1−𝑒²)2,𝑟⟨(𝑎)²𝑟cos(𝑓)sin(𝑓)⟩=0,⟨(𝑎)²⟩=2+3𝑒²2. (2.8)

After using those equations, the disturbing potential averaged over the eccentric anomaly that the spacecraft has from the action of the disturbing body becomes⟨𝑅₂⟩=3𝜇^′𝑛^′2𝑎²4(𝑎^′𝑟^′)³[𝛼²(1+4𝑒²)+𝛽²(1−𝑒²2)−(3+𝑒²)]. (2.9)

The next step is to obtain the second average with respect to the disturbing body. To do this, we considered the Keplerian elements of the spacecraft constant during the process of averaging [17]. Thus we obtain𝑎⟨(^′𝑟^′)𝛼²1⟩=[2+34𝑒^′2+15𝑒16^′4][cos²(𝜔)+cos²(𝑖)sin²𝑎(𝜔)],⟨(^′𝑟^′)𝛽²1⟩=[2+34𝑒^′2+15𝑒16^′4][sin²(𝜔)+cos²(𝑖)cos²(𝜔)]. (2.10)

After performing these averages for the perturbed and perturbing bodies, the equation ⟨⟨𝑅₂⟩⟩ obtained from the double-averaged disturbing function is⟨⟨𝑅₂⟩⟩=𝐾[2(3cos²(𝑖)−1)+3𝑒²(3cos²(𝑖)−1)+15𝑒²sin²(𝑖)cos(2𝜔)], (2.11)where𝜇𝐾=^′𝑛^′2𝑎²316[1+2𝑒^′2+158𝑒^′4]. (2.12)

The partial derivatives of ⟨⟨𝑅₂⟩⟩ with respect to 𝑎,𝑒,𝑖, and 𝜔 can thus be written as𝜕⟨⟨𝑅₂⟩⟩=𝜕𝑎𝑎𝑛^′2𝜇^′83[1+2𝑒^′2+158𝑒^′4][2(3cos²(𝑖)−1)+3𝑒²(3cos²(𝑖)−1)+15𝑒²sin²(𝑖)cos(2𝜔)],𝜕⟨⟨𝑅₂⟩⟩𝜕𝑒=6𝐾[𝑒(3cos²(𝑖)−1)+5𝑒sin²(𝑖)cos(2𝜔)],𝜕⟨⟨𝑅₂⟩⟩𝜕𝑖=3𝐾[−2sin(2𝑖)−3𝑒²sin(2𝑖)+5𝑒²sin²(𝑖)cos(2𝜔)],𝜕⟨⟨𝑅₂⟩⟩𝜕𝜔=−30𝐾𝑒²sin²(𝑖)cos(2𝜔). (2.13)

Now we need to quantify the resulting variations in these four orbital elements of the perturbed body. To obtain this, we will derive Lagrange's planetary equations. The results are given by𝑑𝑎𝑑𝑡=0,𝑑𝑒=𝑑𝑡15𝜇^′𝑛^′2𝑒√1−𝑒²38𝑛[1+2𝑒^′2+158𝑒^′4]sin²(𝑖)sin(2𝜔),𝑑𝑖=𝑑𝑡−15𝜇^′𝑛^′2𝑒²16𝑛√1−𝑒²3[1+2𝑒^′2+158𝑒^′4]sin(2𝑖)sin(𝜔),𝑑𝜔=𝑑𝑡3𝜇^′𝑛^′28𝑛√1−𝑒²3[1+2𝑒^′2+158𝑒^′4][(5cos²(𝑖)−1+𝑒²)+5(1−𝑒²−cos²(𝑖))cos(2𝜔)],𝑑Ω=𝑑𝑡3𝜇^′𝑛^′2cos(𝑖)8𝑛√1−𝑒²3[1+2𝑒^′2+158𝑒^′4][5𝑒²cos(2𝜔)−3𝑒²−2],𝑑𝑀=𝑑𝑡−𝜇^′𝑛^′238𝑛[1+2𝑒^′2+158𝑒^′4][(3𝑒²+7)(3cos²(𝑖)−1)+15(1+𝑒²)sin²(𝑖)cos²(𝜔)]. (2.14) Analyzing these equations of motion, it can be noticed the following (see [13, 17]).

(i)𝑎 is constant during the integration.(ii)All the results based on these equations are valid for any semimajor axis (i.e., present in the equations in terms of 𝑛) and system of primaries in a proportional time scale.(iii)Here we have a system of three simultaneous ordinary differential equations. They contain essentially the variables 𝑒,𝑖, and 𝜔.(iv)The equation for the longitude of the ascending node Ω depends on the variables 𝑒,𝑖, and 𝜔 but it does not influence their motion.(v)When 𝑒=0 and/or 𝑖=0 there are no variations on the inclination and eccentricity, and the orbit remains circular and/or planar. These circular solutions with constant inclination appear due to the truncation of the expansion of the disturbing function and are not a physical phenomenon. Although no variations for the eccentricity and inclination are obtained, in a real case (full restricted three-body problem), the circular solutions with constant inclination do not exist. The eccentricity can oscillate with large amplitude that depends on the value of the initial eccentricity.(vi)For the eccentricity 𝑒 of the spacecraft, we can see immediately that it increases with 𝑒^′, which can be explained by the decrease of the minimum distance between the primaries. Thus, the perturbations for an orbiting spacecraft are maximum when the secondary body is near the pericenter of its orbit. Analyzing the equations, we have that, when 𝑒^′ is different from zero, 𝑒 increases by a scale factor of [1+(3/2)𝑒^′2+(15/8)𝑒^′4]. It is also interesting to notice that, when the eccentricity increases, there is a decrease in the inclination.(vii)In every equation, appears the term [1+(3/2)𝑒^′2+(15/8)𝑒^′4] that is an approximation of the mean value of the function (𝑎^′/𝑟^′)³ that is given by (1−𝑒^′2)^(−3/2). This approximation can reduce the accuracy of the numerical results. In Figure 2 there are two curves that show the evolutions of (1−𝑒^′2)^(−3/2) (solid curve) and its mean value [1+(3/2)𝑒^′2+(15/8)𝑒^′4] (dashed line) as a function of the value of 𝑒^′. Analyzing this figure, we can see that there is a tendency to diverge when the value 𝑒^′ increases, but they have a good agreement for small values of 𝑒^′.(viii)A pseudo-time can be introduced in the equations by the expression 𝑡^∗=[1+(3/2)𝑒^′2+(15/8)𝑒^′4] where 𝑡 is the period of the oscillation in the circular model and 𝑡^∗ is the period of oscillation when the value of 𝑒^′ is considered. In this way, we have 𝑡=[1+(3/2)𝑒^′2+(15/8)𝑒^′4]⁻¹𝑡^∗, which means that the period of oscillations decreases when the eccentricity of the disturbing body increases.

Figure 2 

Evolution of (1−𝑒′2)(−3/2) (solid curve) and its mean value [1+(3/2)𝑒′2+(15/8)𝑒′4] (dashed line) as a function of the value of 𝑒′.

According to these studies, in the elliptic restricted three body problems, the evolutions in time of the orbital elements of the satellite depend on its initial conditions, on the eccentricity of the disturbing body 𝑒^′, and on the mass ratio 𝜇^′. Considering that the magnitude of the relative forces on the satellite is subjected to change with the separation between the central and the disturbing bodies, so the orbit of the satellite can have different characteristics depending on its distance from the central body. If the distance between the central body and the satellite increases, the gravitational force of the central body decreases with the square of this distance, so the perturbations of the third body become more important. The increase of the distance between the central body and the satellite may cause regions of stable orbits quasiperiodic or chaotic. So, escape or collision of the satellite may occur.

3. Numerical Simulations

In this section, we show the effects of nonzero initial eccentricities for the perturbing body. We numerically investigate the variations of 𝑒 and 𝑖 for a spacecraft within the elliptic restricted three-body problem Earth-Moon spacecraft. The spacecraft is in an elliptic three-dimensional orbit around the Moon and its motion is perturbed by the Earth.

An interesting question that appears in this problem is what happens to the stability of near-circular orbits. The answer for this question depends on the initial inclination 𝑖₀ (inclination between the perturbed and perturbing bodies). If this inclination is above the critical value, the orbit becomes very elliptic and the spacecraft may escape. In opposition, if the inclination is below the critical value, the orbit remains near circular [17]. In the double-averaged second-order model the critical inclination is 𝑖=0.684719203 radians or 39.23152048 degrees [13].

In Prado [17] a lunar satellite with 𝑎₀=0.01 (3844 km), 𝑒₀=0.01, and 𝜔₀=Ω₀=0^∘ was used to show the behavior of the evolutions of the inclination and the eccentricity with time for near-circular orbits with 𝑖₀≤80^∘. The perturbing body was in circular orbit.

For our numerical integrations, the initial conditions of Prado [17] were used. The eccentricity of the perturbing body in the range 0≤𝑒^′≤0.6 was examined. It means that we generalized the Earth-Moon system in order to measure the effects on the eccentricity of the primaries. The time is defined such that the period of the disturbing body is 2𝜋. Therefore, the integrations were performed for 2000 canonical units, what correspond to 320 orbits of the disturbing body (in this case of the Earth-Moon system).

The results are summarized in Figures 7–6. The curves show the evolution of the inclination and the eccentricity as a function of the time, respectively. The results were obtained considering the full elliptic restricted three-body problem and a double averaged second-order model. In general, the behavior of the eccentricity and inclination are according to the expected. The eccentricity oscillates with large amplitude, reaching the value 0.97 only in the case 𝑖₀=80^∘. The inclination remains close to constant most of the time, but from time to time it decreases to the value of the critical inclination and then returns to its initial value. The minimum in inclination occurs at the same time as the maximum in eccentricity. As expected, the averaged model presented smooth curves.

(a)

(b)

(c)

(d)

Figure 4 

Inclination and eccentricity as a function of the time. The perturbing body was in an elliptic orbit with 𝑒′=0.1. The results were obtained considering the full elliptic restricted three-body problem (above) and double-averaged second-order model (below).

(a)

(b)

(c)

(d)

Figure 5 

Inclination and eccentricity as a function of the time. The perturbing body was in an elliptic orbit with 𝑒′=0.3. The results were obtained considering the full elliptic restricted three-body problem (above) and double-averaged second-order model (below).

(a)

(b)

(c)

(d)

Figure 6 

Inclination and eccentricity as a function of the time. The perturbing body was in an elliptic orbit with 𝑒′=0.5. The results were obtained considering the full elliptic restricted three-body problem (above) and double-averaged second-order model (below).

(a)

(b)

(c)

(d)

Figure 7 

Inclination and eccentricity as a function of the time. The perturbing body was in an elliptic orbit with 𝑒′=0.6. The results were obtained considering the full elliptic restricted three-body problem (above) and double-averaged second-order model (below).

(a)

(b)

(c)

(d)

Figure 3 

Inclination and eccentricity as a function of the time. The perturbing body was in a circular orbit. The results were obtained considering the full elliptic restricted three-body problem (above) and double-averaged second-order model (below).

By focusing on the value of the eccentricity of the perturbing body, our results show differences between the results of the numerical simulations considering the full problem and the double averaged second-order model. In the two models, the amplitudes of the inclination and the eccentricity do not suffer significant changes with the increase of 𝑒^′. In the second-order model, the time required for the inclination to reach its critical value decreases when e' increases. In this case, the critical value is reached a greater number of times for a given time, when compared to the results of the full model. The numerical results of the full problem showed that there is no strong effects with the increase of 𝑒^′ on the inclination and eccentricity of the spacecraft. It can be seen that the behavior of the eccentricity and the inclination are similar for all the cases. Thus, what predominates in this process is the distance of the spacecraft to the secondary body. The spacecraft suffers a small gravitational perturbation of the secondary body and the effects of such perturbation are not significant.

The comparison between the results above shows that the second-order model (only the part 𝑅₂ of the disturbing function) is not very accurate to study orbits with higher values of 𝑒^′ and further expansions have to be made.

To illustrate the behavior of the inclination and eccentricity when the value of 𝑖₀ is around the critical value (38^° ≤ i₀ ≤ 43^°), Figure 8 shows the evolution of the eccentricity and the inclination as a function of time for 𝑒^′=0.0 and 0.3. The figures for other values of 𝑒^′ are not presented here, but show similar behaviors. The results showed that the behavior of the eccentricity and inclination are very similar on the two models.

(a)

(b)

(c)

(d)

Figure 8 

Inclination and eccentricity as a function of time. The perturbing body was in elliptic orbit with 𝑒′=0.0 and 0.3. The results were obtained considering the full elliptic restricted three-body problem (solid line) and double-averaged second-order model (dashed line).

The results were obtained considering the full elliptic restricted three-body problem (solid line) and double-averaged second-order model (dashed line). The values of 𝑖₀=41^∘ and 43^∘ (above the critical inclination) give results that showed that the inclination oscillates with very small amplitude (≈0.01 rad), decreases until the critical value, and then returns to its value 𝑖₀. The eccentricity oscillates with significant amplitude (about 0.30). Again, when the eccentricity reaches its maximum, the inclination reaches its minimum for the case of initial condition near the critical value. For the effect of increasing 𝑖₀ and/or 𝑒^′, the time for reaching the critical value is reduced. For the value of 𝑖₀=38^∘ (just below the critical inclination) the inclination and eccentricity stay close to 𝑖₀ and 𝑒₀, respectively.

The next simulations show results (Figure 9) for near-equatorial and circular orbits. For this case, our simulations are for 𝑎=0.01,𝑒=0.001,𝑖₀=5^∘,𝜔₀=Ω₀=0^∘, and 0⩽𝑒^′⩽0.6.

(a)

(b)

(c)

(d)

Figure 9 

Eccentricity as a function of time for 𝑒′=0.0, 0.4, 0.5, and 0.6. The results were obtained considering the full elliptic restricted three-body problem (red line) and double-averaged second-order model (black line).

From the second-order model equations, it is easily seen the existence of stable equatorial circular orbits. When the value of 𝑒₀ and/or 𝑖₀ is small, this imply that 𝑒 and 𝑖 will suffer much smaller changes with time, either increasing from or decreasing to zero. This is because their variations with time are proportional to sin²(𝑖) and sin(2𝑖).

The numerical integration of the full problem and the second-order model showed constant inclination all the time. The variation of the eccentricity on the two models had very small amplitudes. The oscillations for the full model (red line) have amplitude in the order of 10⁻³ and, for the second-order model, it was 10⁻⁵.

We presented the results for 𝑒^′=0.0, 0.4, 0.5, and 0.6 because they illustrate very well how the oscillations of e increase with high values of 𝑒^′ (see the figures for 𝑒^′=0.0 and 0.6 as a comparison).

The next simulations (see Figure 10) are made for a lunar satellite with the following data: 𝑎=0.013 (5000 km), 𝑒=0.1, 10⩽𝑖⩽80 degrees, 𝑒^′=0.0 to 0.6. The color scale indicates the eccentricity of the perturbing body. The main effect of the eccentricity of the perturbing body is to reduce the period of the oscillation of the eccentricity of the perturbed body.

(a)

(b)

(c)

(d)

Figure 10 

Eccentricity as a function of the time for 0.0⩽𝑒′⩽0.6. The results were obtained considering the double-averaged second-order model.

4. Conclusions

This paper showed an analytical expansion to study the third-body perturbation for the case where the perturbing body is in an elliptical orbit. It followed the same steps used in the literature for the circular case that is based in the expansion of the perturbing function in polynomials of Legendre. The behavior of the trajectories is similar to the circular case, showing an oscillatory behavior for the eccentricity. In general, the second-order approximation shows that the main effect of the eccentricity of the perturbing body is to reduce the period of the oscillation of the eccentricity of the perturbed body. This reduction is proportional to the eccentricity of the perturbing body. The amplitude of the oscillation is not changed. Better accuracy can be attained with larger expansions.

Acknowledgment

The authors wish to express their appreciation for the support provided by CNPq under the contracts nos. 155029/2006-1, 305147/2005-6, 308111/2006-0 and FAPESP under the contract no. 06/04997-6.