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

Authors: R. C. Domingos; R. Vilhena de Moraes; A. F. Bertachini De Almeida Prado
Journal: Mathematical Problems in Engineering
 (2008)

---

## 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.


---

## Body

## 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 (x,y) 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 m0,
and a second body with mass m′ in an elliptic orbit around this same central
body in the plane x−y.This paper is structured as follows. In Section2, 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 massm0 is fixed in the center of the reference system x−y.
The perturbing body, with mass m′,
is in an elliptic orbit with semimajor axis a′,
eccentricity e′,
and mean motion n′ (given by the expression n′2a′3=G[m0+m′]). The massless spacecraft m is in a generic two-dimensional orbit which
orbital elements are a (semimajor axis), e (eccentricity), i (inclination), ω (argument of periapsis), Ω (longitude of the ascending node), and n (mean motion) given by the expression n2a3=Gm0.
This system is showen in Figure 1. Here, G is the gravitational constant, r and r′ are the radius vectors of bodies m and m′,
and S is the angle between these radius vectors.Figure 1
Illustration of the dynamical system [16].Using the traditional expansion in Legendre
polynomials (assuming thatr′≫r), the disturbing potential is given by
[21](2.1)R=μ′G(m0+m′)r′∑n=2∞(rr′)nPn[cos(S)],where μ′=m′/(m0+m′),G is the gravitational constant, Pn are the Legendre polynomials, and S is the central elongation between the
perturbed body (the spacecraft) and the perturbing body (the second body).The termsn=1 and n=0 are not included in
(2.1). The part R2 of the disturbing function due to the second body is given by(2.2)R2=μ′n′2a22(a′r′)3(ra)2[3cos2(S)−1],where cos(S) is(2.3)cos(S)=αcos(f)+βsin(f), where α=P^⋅r^′,β=Q^⋅r^′,f is the true anomaly of the satellite, r^′ is the unit vector pointing from the central
body to the disturbing body, all of them functions of f′ and ω′ (true anomaly and argument of the periapsis of
the disturbing body, resp.). The usual orthogonal unit vectors P^ and Q^ are functions of i, ω,
and Ω in the plane of the satellite orbit, with P^ pointing toward the periapsis. For the case of
elliptic orbits the products α and β are written as(2.4)α=cos(ω)cos(Ω−f′−ω′)−cos(i)sin(ω)sin(Ω−f′−ω′),β=−sin(ω)cos(Ω−f′−ω′)−cos(i)cos(ω)sin(Ω−f′−ω′).A substitution using (2.3) was made in (2.2). Then we
have(2.5)R2=μ′n′2a22(a′r′)3(ra)2[3(α2cos2(f)+2αβcos(f)sin(f)+β2sin2(f))−1].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 is(2.6)〈F〉=12π∫02πFdM,where M 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 anomaliesf and f′ by the eccentric anomalies E and E′.
To do this task, we use some well-known relations from the Celestial Mechanics
[21] given by the
following equations:(2.7)sin(f)=(1−e2)1/21−ecos(E)sin(E),cos(f)=cos(E)−eecos(E)−1,(ra)=1−ecos(E),dM=(1−ecos(E))dE.Using those equations, we can
obtain the following relations:(2.8)〈(ra)2cos2(f)〉=(1+4e2)2,〈(ra)2sin2(f)〉=(1−e2)2,〈(ra)2cos(f)sin(f)〉=0,〈(ra)2〉=2+3e22.After using those equations, the disturbing potential
averaged over the eccentric anomaly that the spacecraft has from the action of
the disturbing body becomes(2.9)〈R2〉=3μ′n′2a24(a′r′)3[α2(1+4e2)+β2(1−e2)−(23+e2)].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(2.10)〈(a′r′)α2〉=[12+34e′2+1516e′4][cos2(ω)+cos2(i)sin2(ω)],〈(a′r′)β2〉=[12+34e′2+1516e′4][sin2(ω)+cos2(i)cos2(ω)].After performing these averages for the perturbed and
perturbing bodies, the equation〈〈R2〉〉 obtained from the double-averaged disturbing
function is(2.11)〈〈R2〉〉=K[2(3cos2(i)−1)+3e2(3cos2(i)−1)+15e2sin2(i)cos(2ω)],where(2.12)K=μ′n′2a216[1+32e′2+158e′4].The partial derivatives of〈〈R2〉〉 with respect to a,e,i,
and ω can thus be written as(2.13)∂〈〈R2〉〉∂a=an′2μ′8[1+32e′2+158e′4][2(3cos2(i)−1)+3e2(3cos2(i)−1)+15e2sin2(i)cos(2ω)],∂〈〈R2〉〉∂e=6K[e(3cos2(i)−1)+5esin2(i)cos(2ω)],∂〈〈R2〉〉∂i=3K[−2sin(2i)−3e2sin(2i)+5e2sin2(i)cos(2ω)],∂〈〈R2〉〉∂ω=−30Ke2sin2(i)cos(2ω).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(2.14)dadt=0,dedt=15μ′n′2e1−e28n[1+32e′2+158e′4]sin2(i)sin(2ω),didt=−15μ′n′2e216n1−e2[1+32e′2+158e′4]sin(2i)sin(ω),dωdt=3μ′n′28n1−e2[1+32e′2+158e′4][(5cos2(i)−1+e2)+5(1−e2−cos2(i))cos(2ω)],dΩdt=3μ′n′2cos(i)8n1−e2[1+32e′2+158e′4][5e2cos(2ω)−3e2−2],dMdt=−μ′n′28n[1+32e′2+158e′4][(3e2+7)(3cos2(i)−1)+15(1+e2)sin2(i)cos2(ω)]. Analyzing these equations of motion, it can be noticed
the following (see [13, 17]).(i)
a 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 ofn) 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 variablese,i,
and ω.(iv)
The equation for the longitude of the ascending
nodeΩ depends on the variables e,i,
and ω but it does not influence their motion.(v)
Whene=0 and/or i=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 eccentricitye of the spacecraft, we can see immediately that
it increases with e′,
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 e′ is different from zero, e increases by a scale factor of [1+(3/2)e′2+(15/8)e′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)e′2+(15/8)e′4] that is an approximation of the mean value of
the function (a′/r′)3 that is given by (1−e′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−e′2)(−3/2) (solid curve) and its mean value [1+(3/2)e′2+(15/8)e′4] (dashed line) as a function of the value of e′.
Analyzing this figure, we can see that there is a tendency to diverge when the
value e′ increases, but they have a good agreement for
small values of e′.(viii)
A pseudo-time can be introduced in the equations
by the expressiont∗=[1+(3/2)e′2+(15/8)e′4] where t is the period of the oscillation in the
circular model and t∗ is the period of oscillation when the value of e′ is considered. In this way, we have t=[1+(3/2)e′2+(15/8)e′4]−1t∗,
which means that the period of oscillations decreases when the eccentricity of
the disturbing body increases.Figure 2
Evolution of(1−e′2)(−3/2) (solid curve) and its mean value [1+(3/2)e′2+(15/8)e′4] (dashed line) as a function of the value of e′.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 bodye′,
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 ofe and i 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 inclinationi0 (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 i=0.684719203 radians or 39.23152048 degrees [13].In Prado [17] a lunar satellite with a0=0.01 (3844 km), e0=0.01, and ω0=Ω0=0° was used to show the behavior of the
evolutions of the inclination and the eccentricity with time for near-circular
orbits with i0≤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≤e′≤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 Figures7–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 i0=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.Inclination and eccentricity as a function of the time. The
perturbing body was in an elliptic orbit withe′=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)Inclination and eccentricity as a
function of the time. The perturbing body was in an elliptic orbit withe′=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)Inclination and eccentricity as a
function of the time. The perturbing body was in an elliptic orbit withe′=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)Inclination and eccentricity as a function of the time. The perturbing body was
in an elliptic orbit withe′=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)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).
(a)(b)(c)(d)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 ofe′.
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 e′ 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 partR2 of the disturbing function) is not very
accurate to study orbits with higher values of e′ and further expansions have to be made.To illustrate the behavior of the inclination and
eccentricity when the value ofi0 is around the critical value (38° ≤ i0 ≤
43°), Figure 8 shows the evolution of the
eccentricity and the inclination as a function of time for e′=0.0 and 0.3. The figures for other values of e′ 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.Inclination and
eccentricity as a function of time. The perturbing body was in elliptic orbit
withe′=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).
(a)(b)(c)(d)The results were obtained considering the full
elliptic restricted three-body problem (solid line) and double-averaged
second-order model (dashed line). The values ofi0=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 i0.
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 i0 and/or e′,
the time for reaching the critical value is reduced. For the value of i0=38° (just below the critical inclination) the
inclination and eccentricity stay close to i0 and e0,
respectively.The next simulations show results (Figure9) for
near-equatorial and circular orbits. For this case, our simulations are for a=0.01,e=0.001,i0=5°,ω0=Ω0=0°,
and 0⩽e′⩽0.6.Eccentricity as a function of time fore′=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).
(a)(b)(c)(d)From the second-order model equations, it is
easily seen the existence of stable equatorial circular orbits. When the value
ofe0 and/or i0 is small, this imply that e and i will suffer much smaller changes with time,
either increasing from or decreasing to zero. This is because their variations
with time are proportional to sin2(i) and sin(2i).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 of10−3 and, for the second-order model, it was 10−5.We presented the results fore′=0.0, 0.4, 0.5, and 0.6 because they illustrate
very well how the oscillations of e increase with high values of e′ (see the figures for e′=0.0 and 0.6 as a comparison).The next simulations (see Figure10) are made for a
lunar satellite with the following data: a=0.013 (5000 km), e=0.1, 10⩽i⩽80 degrees, e′=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.Eccentricity as a function of the
time for0.0⩽e′⩽0.6.
The results were obtained considering the double-averaged second-order model.
(a)(b)(c)(d)


## 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.



---

*Source: 763654-2008-06-11.xml*