# Effects of the Eccentricity of a Perturbing Third Body on the Orbital Correction Maneuvers of a Spacecraft

Authors: R. C. Domingos; A. F. B. A. Prado; V. M. Gomes
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The fuel consumption required by the orbital maneuvers when correcting perturbations on the orbit of a spacecraft due to a perturbing body was estimated. The main goals are the measurement of the influence of the eccentricity of the perturbing body on the fuel consumption required by the station keeping maneuvers and the validation of the averaged methods when applied to the problem of predicting orbital maneuvers. To study the evolution of the orbits, the restricted elliptic three-body problem and the single- and double-averaged models are used. Maneuvers are made by using impulsive and low thrust maneuvers. The results indicated that the averaged models are good to make predictions for the orbital maneuvers when the spacecraft is in a high inclined orbit. The eccentricity of the perturbing body plays an important role in increasing the effects of the perturbation and the fuel consumption required for the station keeping maneuvers. It is shown that the use of more frequent maneuvers decreases the annual cost of the station keeping to correct the orbit of a spacecraft. An example of an eccentric planetary system of importance to apply the present study is the dwarf planet Haumea and its moons, one of them in an eccentric orbit.


---

## Body

## 1. Introduction

Advances in astronomy have been showing a large number of planets outside the solar system, some of them in orbits with large eccentricities. In that sense, it would be interesting to know how often and how expensive would it be to keep the most important orbital elements (inclination and eccentricity) inside defined limits for a satellite orbiting one of those planets and perturbed by the central star, or even for a satellite orbiting a moon of those planets. The Solar System also has some examples of bodies with large eccentricity, like the dwarf planet Haumea and one of its moons [1]. This body belongs to the Kuiper belt and has a moon (Namaka) with eccentricity of 0.249, which is near the limit where the method shown here is valid. So, there are many possible applications of the ideas developed here. The study is made using the eccentricity of the perturbing body as a parameter, because our goal is to explore how much this eccentricity affect the orbit of the satellite. It is possible to see that, for some values of the eccentricity, no orbital maneuvers are required for longer time, in the order of decades.An orbital maneuver is a procedure performed with the goal of modifying the orbit of a satellite by using a propulsive system fixed to the spacecraft, which is able to modify the velocity of the spacecraft and, consequently, its orbit. This procedure can be used when inserting a satellite in its final orbit from a previous parking orbit where the launcher was able to deliver it, or to correct the deviations from the nominal orbit due to the action of perturbation forces. One of the first results in this field is the so-called Hohmann transfer [2]. It considers the transfer with minimum variation of velocity between two coplanar and circular orbits. Later, the three-impulsive concept appeared in the literature, in independent researches developed by Hoelker and Silber [3] and Shternfeld [4]. In this maneuver, a bielliptical three-impulsive transfer between two circular orbits has a lower variation of velocity than the one required by the Hohmann transfer if the ratio between the initial and final radius of the orbits is greater than 15.48. In the range from 11.94 to 15.58, the more economical maneuver can use two or three impulses, depending on the value of the distance where the second impulsive is applied (see Marec [5]).The next step is to consider whether a low thrust can be used to maneuver the satellite. In this situation, the propulsion system applies a finite force for a nonnegligible time. Many researches using this model are based directly on the “Primer-Vector” theory, which was developed by Lawden [6, 7]. Other strategies to solve this problem are based on Pontryagin’s maximum principle. Biggs [8, 9] and Gomes and Prado [10] show some results using this technique in researches that are the basis of the procedure used in the present paper when the low thrust maneuvers are considered.The present paper makes use of both types of maneuvers cited above to correct the perturbations from a third-body in an elliptical orbit. It is important to note that the amplitudes of the maneuvers performed here are small, because it is only required to make the satellite to return to its nominal orbit, correcting the effects of the third-body perturbation. Thus, an initial and a final orbit around the Earth are given and the goal is to search for solutions on how to transfer the satellite between those two orbits, which presents a low cost in terms of fuel consumption. In that sense, this is an orbital maneuver between two close near circular orbits. Previous researches made in this problem can be found by Marec [11], Edelbaum [12], Fernandes and Golfetto [13], Fernandes and Carvalho [14], and Fernandes [15].In this way, the present research concentrates on studying the time evolution of a satellite perturbed by a third-body, for times of the order of a few tens of years, under the single- and double-averaged models, to compare their results with the full elliptic restricted three-body problem, that will be our “full model” to represent the reality. These first results give an idea of the differences between those models in this time span, which is chosen to cover times larger than the usual spacecraft missions around the Earth. Those larger times can include missions like a permanent basis orbiting an exoplanet in an elliptic orbit around a central star or a “temporary” spacecraft orbiting asteroids in eccentric orbits. A spacecraft placed in orbit around Haumea would also benefit from those previous studies. Missions like those ones have a good potential to require duration of the order of 20–30 years. Another result of this research is to find the times when it is necessary to perform orbital maneuvers to send the satellite back to its nominal orbit, thus correcting the variations introduced by the presence of the third-body. These calculations are based on the instants that the satellite reaches one of the accepted limits of the mission for the deviations in eccentricity and/or inclination. Since the corrections in the inclination are made by out-of-plane maneuvers and the corrections in eccentricity are made by planar maneuvers, both corrections are assumed to be made independently of each other. Two different approaches are used for both situations, in order to compare the cost involved: an impulsive maneuver that is easier to calculate and is usually considered as a first model for almost all mission designers and a low thrust transfer that has the advantage of using less fuel to achieve the same final orbit, due to the higher value of the specific impulse, but requires a more complex implementation of the hardware used to control the system.


## 2. Mathematical Models

In order to study the maneuvers to control eccentricity and inclination due to the perturbations caused by a third-body in an elliptical orbit, the assumptions used to develop the single- and double-averaged analytical models are the same ones of the restricted elliptic three-body problem. It is assumed that a planet with massm
0 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
′. The satellite with mass m that can be considered negligible when compared to the other masses involved is in a generic orbit whose orbital elements are a (semi-major axis), e (eccentricity), i (inclination), ω (argument of periapsis), and Ω (right ascension of the ascending node). Its mean motion, obtained from its semimajor axis, is called n. The reason to develop averaged models is that they can be used to show the behavior of the orbit in longer time periods, by removing short term oscillations that may not be important in some studies. Using this technique, it is possible to obtain analytical equations of motion that can be analyzed to show some characteristics of the motion and that can also be integrated to give fast results in terms of the long period evolution of the orbital elements.The equations of motion of a satellite perturbed by a third-body using single- and double-averaged techniques have been studied in detail by Domingos et al. [16, 17]. In these works, the authors also considered an elliptic orbit for the perturbing body to study its effects on the orbit of the satellite. Here, we briefly summarize the equations of motion and then perform numerical integrations of those equations with the goal of comparing the models from the point of view of the effects in the orbital maneuvers required by the mission. Different initial eccentricities (0.0, 0.1, and 0.2) for the perturbing body are considered. Those values are based on Domingos et al. [17], where it is shown that the averaged models presented here have good accuracy up to those values for the eccentricity of the primaries. The system considered here has the same physical data (masses and distances) of the Earth-Moon system, with the additional flexibility of allowing the eccentricity of the orbit of the secondary body around the primary to grow up to 0.2, in order to study the effects of this parameter on the behavior of the system. These hypotheses are made to allow us to obtain results that are valid for other real and hypothetical systems where the primaries of the system have an eccentricity of this level and similar mass parameter. Some early researches studying the third-body perturbation can be found in Cook [18] and Smith [19], both of them considering a satellite around the Earth perturbed by the Sun and the Moon, but not using average techniques and taking into account short term oscillations of the orbital parameters. Smith [19] considers elliptical orbits for the perturbing body, while Cook [18] makes the hypothesis of circular motions for the perturbing bodies.


### 2.1. Single-Averaged Model

In Domingos et al. [16], Section 2, a complete review and more details about the development of the disturbing potential and the average process are shown.The equations of motion are used as in Domingos et al. [17] to compute the effects of the perturbing body on the satellite. Using an expansion in Legendre’s polynomials (considering r
′
≫
r), the second order part of the disturbing potential, considering an average over the eccentric anomaly of the satellite, is given by
(1)
〈
R
2
〉
=
μ
′
a
2
n
′
2
2
(
a
′
r
′
)
3
×
{
(
1
+
3
2
e
2
)
[
3
2
(
α
2
+
β
2
)
-
1
]
iiiiii
+
15
4
(
α
2
-
β
2
)
e
2
}
,
where r
′
=
a
′
(
1
-
e
′
2
)
/
(
1
+
e
′
cos
⁡
f
′
), cos
⁡
f
′
=
cos
⁡
M
′ + e
′
(
cos
⁡
2
M
′
-
1
) + (
9
/
8
)
e
′
2
(
cos
⁡
3
M
′
-
cos
⁡
M
′
) + ⋯, and sin
f
′
=
sin
M
′ + e
′
sin
2
M
′ + e
′
2
(
(
9
/
8
)
sin
3
M
′
-
(
7
/
8
)
sin
M
′
) + ⋯.When considering elliptic orbits for the perturbing body, the values ofα and β are given by α
=
cos
⁡
ω
cos
⁡
D
-
cos
⁡
i
sin
ω
sin
D, β
=
-
sin
ω
cos
⁡
D
-
cos
⁡
i
cos
⁡
ω
sin
D, and D
=
Ω
-
f
′
-
ω
′. The mean anomaly of the perturbing body M
′ is given by M
′
=
M
o
′
+
n
′
t. Thus, the variations in the orbital elements of the perturbed body are obtained. To do this, we derived Lagrange’s planetary equations that describe the variations of the mean orbital elements of the spacecraft. They are given by
(2)
d
e
d
t
=
-
(
1
-
e
2
)
1
/
2
n
a
2
e
∂
〈
R
〉
∂
ω
,
d
i
d
t
=
cos
⁡
i
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
ω
-
1
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
Ω
,
d
Ω
d
t
=
1
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
i
,
d
ω
d
t
=
-
cos
⁡
i
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
i
+
(
1
-
e
2
)
1
/
2
n
a
2
e
∂
〈
R
〉
∂
e
.The semimajor axis is constant since the mean anomalyM was eliminated in the perturbing function. Following this idea, we numerically integrate the equation of motion (3), obtained from Lagrange’s planetary equations. The results are shown below:
(3)
d
e
d
t
=
K
15
4
μ
′
n
′
2
e
1
-
e
2
n
×
[
sin
2
ω
(
cos
⁡
2
D
-
cos
⁡
2
i
sin
2
D
)
iiiii
(
cos
⁡
2
D
)
-
cos
⁡
i
cos
⁡
2
ω
sin
2
D
]
,
d
i
d
t
=
K
1
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
3
4
μ
′
n
′
2
a
2
×
{
e
2
[
-
5
cos
⁡
2
ω
sin
2
D
cos
⁡
2
i
3
2
iiiiiiiii
+
3
2
(
-
sin
2
D
+
sin
2
D
cos
⁡
2
i
)
iiiiiiiii
-
5
2
sin
D
cos
⁡
2
ω
(
1
+
cos
⁡
2
i
)
]
iiiiii
+
(
sin
D
-
sin
2
D
cos
⁡
2
i
)
iiiiii
3
2
-
5
cos
⁡
2
D
sin
2
ω
cos
⁡
i
}
,
d
Ω
d
t
=
K
3
4
μ
′
n
′
2
a
2
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
×
{
(
1
+
3
2
e
2
)
(
-
sin
2
D
sin
2
i
)
iiiii
+
[
cos
⁡
2
i
cos
⁡
2
ω
sin
2
D
+
sin
i
sin
2
D
sin
2
ω
]
iiiii
×
(
5
2
e
2
)
}
,
d
ω
d
t
=
K
3
4
μ
′
n
′
2
×
{
-
a
2
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
iiiii
×
[
(
1
+
3
2
e
2
)
(
-
sin
2
D
sin
2
i
)
iiiiiiiiii
+
(
5
2
e
2
)
(
cos
⁡
2
i
sin
2
D
cos
⁡
2
ω
)
iiiiiiiiii
+
1
2
sin
2
i
sin
2
D
sin
2
ω
]
iiiii
+
1
-
e
2
n
[
5
(
cos
⁡
2
D
-
cos
⁡
2
i
sin
2
D
)
cos
⁡
2
ω
6
2
iiiiiiiiiiiiiiiiiiiiiii
-
5
cos
⁡
i
sin
2
D
sin
2
ω
iiiiiiiiiiiiiiiiiiiiiii
+
a
2
[
(
e
2
)
]
1
/
2
6
2
(
cos
⁡
2
D
+
cos
⁡
2
i
sin
2
D
)
]
}
,
where K
=
(
1
+
3
e
′
cos
⁡
f
′
+
3
e
′
2
cos
⁡
2
f
′
)
/
(
1
-
6
e
′
+
15
e
′
2
).This set of equations has singularities for equatorial orbits, because the inclination appears in the denominator of some of them. There is no singularity for small eccentricities, because it does not appear in the denominator. So, it is necessary to check the results very carefully against the other models, which have no singularities, when near equatorial orbits are considered. In any case, exactly equatorial orbits were not used in the present research.


### 2.2. Double-Averaged Model

Using (1) and considering an average over the eccentric anomaly of the perturbing body (Domingos et al. [16]), the equation obtained for the double-averaged disturbing function is
(4)
〈
〈
R
2
〉
〉
=
K
′
[
2
(
3
cos
⁡
2
i
-
1
)
+
3
e
2
(
3
cos
⁡
2
i
-
1
)
iiiii
+
15
e
2
sin
2
i
cos
⁡
2
ω
]
.
Now, it is possible to quantify the resulting variations in the orbital elements of the satellite. To obtain this result, again Lagrange’s planetary equations are derived using the perturbing potential shown by (4). The results are given by
(5)
d
e
d
t
=
K
′
15
μ
′
n
′
2
e
1
-
e
2
8
n
sin
2
(
i
)
sin
(
2
ω
)
,
d
i
d
t
=
K
′
-
15
μ
′
n
′
2
e
2
16
n
1
-
e
2
sin
(
2
i
)
sin
(
ω
)
,
d
ω
d
t
=
K
′
3
μ
′
n
′
2
8
n
1
-
e
2
[
(
5
cos
⁡
2
(
i
)
-
1
+
e
2
)
iiiiiiiiiiiiiiiiiiiiii
+
5
(
1
-
e
2
-
cos
⁡
2
(
i
)
)
cos
⁡
(
2
ω
)
]
,
d
Ω
d
t
=
K
′
3
μ
′
n
′
2
cos
⁡
(
i
)
8
n
1
-
e
2
[
5
e
2
cos
⁡
(
2
ω
)
-
3
e
2
-
2
]
,
d
M
d
t
=
-
K
′
μ
′
n
′
2
8
n
[
(
3
e
2
+
7
)
(
3
cos
⁡
2
(
i
)
-
1
)
iiiiiiiiiiiiiiiii
+
15
(
1
+
e
2
)
sin
2
(
i
)
cos
⁡
2
(
ω
)
]
,
where K
′
=
[
1
+
(
3
/
2
)
e
′
2
+
(
15
/
8
)
e
′
4
].This set of equations has no singularities, so it is important to use it to verify the results obtained by the previous method, in particular when near equatorial orbits are considered.


### 2.3. The Impulsive Transfers to Correct the Eccentricity

In present study, we have considered the time evolution and maneuvers to correct only inclination and eccentricity, because they were considered the most important elements that impact the orbit of the satellite (see [13, 14]). As shown by the mathematical models used here, the semimajor axis is not affected by third-body perturbations, except for short periodic oscillations of small amplitudes. The argument of periapsis and the right ascension of the ascending node do not affect very much the mission under analysis, since they do not change the regions of the celestial body that is covered by the satellite. Those variables have a secular perturbation that causes a circular motion for both angles (Solórzano and Prado [20]).To return the eccentricity of an orbit to its original value and keeping the same semimajor axis, one possibility is to make a bi-impulsive transfer to send the satellite from an orbit O2, with a given semi-major axis “a” and eccentricity “e
2”, to an orbit O1, with the same semimajor axis “a” but with eccentricity “e
1”. Another possibility would be to perform a single impulsive maneuver at the point where the orbits intersect each other. The two-impulsive maneuver will be used, because it uses the apsis of the orbits, which are better in terms of fuel consumption, and also because splitting the maneuver in two parts makes the task of correcting errors in the impulses applied easier. Since the semimajor axis is not changed by the third-body perturbations, in terms of average, it is also desired to perform orbital maneuvers that do not change this parameter. So, when correcting the eccentricity, this is the only orbital parameter that needs to be changed by the maneuver. Figure 1 shows this situation. A good solution is a Hohmann type transfer that uses the most distant apoapsis. So, it is necessary to apply the first impulse, with magnitude Δ
V
1, to the satellite when it is at the apoapsis of the orbit O2 that is a point where the distance from the center of the Earth is r
a
2. This impulse puts the satellite in a transfer orbit O
T that is an ellipse tangent to the apoapsis of O2 and to the perigee of O1, where the distance from the center is r
p
1. When the satellite reaches the periapsis of the orbit O1, the second impulse, with magnitude Δ
V
2, is applied to the satellite and this impulse puts the satellite in the final desired orbit O2. This second impulse is applied opposite to the direction of the motion of the satellite, in order to remove energy from the satellite, thus returning its semimajor axis to the original value.Figure 1
Geometry of the bi-impulsive transfer.To obtain the magnitude of the impulses, it is possible to start with the vis viva equation (6), where v is the velocity of the satellite, r is the distance between the satellite and the center of the Earth, and μ is the gravitational parameter of the Earth, whose numerical value used here is 398600 km3/s2:
(6)
v
=
2
μ
r
-
μ
a
.Based on that equation, the magnitude of the first impulse is the difference between the velocity of the satellite when it is at the pointr
a
2 in the orbit O2 and the velocity when the satellite is at the same point but in the orbit O
T. The equation is
(7)
Δ
V
1
=
2
μ
a
(
1
-
e
2
)
-
μ
a
T
-
2
μ
a
(
1
-
e
2
)
-
μ
a
,
where a
T is the semimajor axis of the transfer orbit O
T, which is given by
(8)
a
T
=
a
+
a
2
(
e
2
-
e
1
)
.The second impulse is applied at the periapsis of the nominal orbit and it is in the opposite direction with respect to the motion of the satellite, thus reducing the energy and the semimajor axis of the orbit to return it to its original value. Its magnitude is the difference between the velocities of the satellite at this point when it is in the orbitO
T and its velocity at the same point but when the satellite is at the orbit O2. The equation is
(9)
Δ
V
2
=
2
μ
a
(
1
-
e
1
2
)
-
μ
a
T
-
2
μ
a
(
1
-
e
1
2
)
-
μ
a
.The total variation of velocity required by the maneuver is given byΔ
V
T
=
Δ
V
1
+
Δ
V
2. The conversion to get the fuel mass consumed is obtained from the well-known Tsiolkovskii equation [21]:
(10)
m
0
=
m
f
e
(
Δ
V
i
/
g
0
I
sp
)
,
where g
0 is the acceleration of gravity at the surface of the Earth, assumed to be 9.8 m/s2; I
sp is the specific impulse of the engine (considered to be 340 s); the subscript i assumes the values 1 and 2 to represent both impulses; m
0 is the mass of the satellite before the impulse; and m
f is the mass of the satellite after the impulse, such that m
0
-
m
f is the fuel consumed by the maneuver. So, it is necessary to have a value for the initial mass of the satellite and the parameters of the combination engine/fuel in order to obtain the total consumption in terms of fuel mass. This equation is used for both impulses separately, updating the value of the mass after the first impulse, not to Δ
V
T, to get more accuracy.


### 2.4. The Impulsive Transfer to Correct the Inclination

In this situation, the simplest strategy is to make a single impulsive transfer to rotate the velocity of the satellite by the desired angle. This impulse has to be applied when the satellite passes by the intersection of the two planes: the one where the satellite is orbiting and the one where it should go after the maneuver. The magnitude of this impulse is given by [22]
(11)
Δ
V
i
=
2
v
sin
(
Δ
i
2
)
,
where v is the velocity of the satellite at the moment of the application of the impulse and Δ
i is the plane change required by the maneuver. Since the nominal orbit is not circular, the magnitude of this impulse depends on the true anomaly of the point of the application of the impulse, which depends on the geometry of the problem (values of the argument of the periapsis and the right ascension of the ascending node) and how the perturbation changes those variables. To avoid specific cases and considering that the eccentricity is small so that the magnitude of the velocity does not oscillate too much, it is possible to consider the worst case, which is when making the maneuver at the periapsis of the orbit, the point where the velocity is larger. In this case, the maximum velocity variation required by the maneuver that changes the inclination (Δ
V
max
⁡) is given by
(12)
Δ
V
i
max
⁡
=
2
sin
Δ
i
2
2
μ
a
(
1
-
e
1
2
)
-
μ
a
.This equation is obtained just by replacing the term in velocity by the velocity at periapsis. The conversion to fuel mass is obtained once again by using (10).


### 2.5. Low Thrust Maneuvers

The same transfers described above can be performed using low thrust. The reason to do that is the smaller fuel consumption required by the low thrust propulsion system, due to its higher specific impulse. It is assumed that the satellite can arrive and leave at any point in the final and initial orbits, respectively. The maneuver is assumed to be realized by an engine that delivers a thrust with variable direction and constant magnitude.The dynamical model assumes that the satellite does not suffer any perturbation during the maneuver. This assumption is valid because the durations of the maneuvers are very short when compared to the orbital period of the satellite. The thrust follows the rules of fixed magnitude, constant ejection velocity, and free angular motion. To solve this problem, first order necessary conditions for a local minimum are written to obtain the adjoint equations. Pontryagin’s Maximum Principle is used to give the control angles at each instant. The approach used here is the same hybrid approach showed in Biggs [8, 9] and later used by Gomes and Prado [10], where all the details of the method are shown, including the variables used and the equations of motion. Constraints can be added, if necessary. The solution of the problem is the time histories of the direction of the thrust and the instants to start and finish the application of the thrust. The fuel consumed and the times that the engine is active are also part of the solution.To allow more savings in the fuel consumption, two situations are simulated when using the low thrust propulsion system. The first one performs the maneuver in a single propulsive arc and the second one allows the possibility of two thrusting arcs, so the engine is turned on and off twice. This maneuver allows a better use of the fuel, by applying the thrust in the most favorable parts of the orbit, which are the apoapsis regions. Figure2 shows the results of a two-arc maneuver. The horizontal axis shows the true anomaly of the spacecraft where the engine should be fired and the vertical axis shows the direction of the thrust, both in degrees. One figure is made for each propulsion arc. This example is the maneuver to correct an error of 0.01 in the eccentricity, thus changing it from 0.02 to 0.01. The other situations are similar, so only the fuel consumption and the duration of the maneuvers are shown in the tables and the figures are omitted. The duration of the application of the thrust is short, so the plots of the trajectories are not interesting, since they are Keplerian orbits. After that, the errors of the propulsion system were considered, using a very simple approach. The solution is obtained as a set of points, which represents the whole trajectory. Then, we inserted a fixed error in the magnitude (5%) and in the direction of the solution obtained (2 degrees) at every point obtained as a solution of the maneuver. The constraints of reaching the final orbit specified were the same of the maneuvers without errors. So, using this simple approach, it was possible to see that the errors in the fuel consumed were in the order of 5–7%. This result is expected because the maneuvers have small magnitudes and so there is an almost linear dependence between errors and extra consumption.Low thrust to correct the eccentricity for the maneuver from 0.02 to 0.01 using two arcs of propulsion.
(a)
(b)


## 2.1. Single-Averaged Model

In Domingos et al. [16], Section 2, a complete review and more details about the development of the disturbing potential and the average process are shown.The equations of motion are used as in Domingos et al. [17] to compute the effects of the perturbing body on the satellite. Using an expansion in Legendre’s polynomials (considering r
′
≫
r), the second order part of the disturbing potential, considering an average over the eccentric anomaly of the satellite, is given by
(1)
〈
R
2
〉
=
μ
′
a
2
n
′
2
2
(
a
′
r
′
)
3
×
{
(
1
+
3
2
e
2
)
[
3
2
(
α
2
+
β
2
)
-
1
]
iiiiii
+
15
4
(
α
2
-
β
2
)
e
2
}
,
where r
′
=
a
′
(
1
-
e
′
2
)
/
(
1
+
e
′
cos
⁡
f
′
), cos
⁡
f
′
=
cos
⁡
M
′ + e
′
(
cos
⁡
2
M
′
-
1
) + (
9
/
8
)
e
′
2
(
cos
⁡
3
M
′
-
cos
⁡
M
′
) + ⋯, and sin
f
′
=
sin
M
′ + e
′
sin
2
M
′ + e
′
2
(
(
9
/
8
)
sin
3
M
′
-
(
7
/
8
)
sin
M
′
) + ⋯.When considering elliptic orbits for the perturbing body, the values ofα and β are given by α
=
cos
⁡
ω
cos
⁡
D
-
cos
⁡
i
sin
ω
sin
D, β
=
-
sin
ω
cos
⁡
D
-
cos
⁡
i
cos
⁡
ω
sin
D, and D
=
Ω
-
f
′
-
ω
′. The mean anomaly of the perturbing body M
′ is given by M
′
=
M
o
′
+
n
′
t. Thus, the variations in the orbital elements of the perturbed body are obtained. To do this, we derived Lagrange’s planetary equations that describe the variations of the mean orbital elements of the spacecraft. They are given by
(2)
d
e
d
t
=
-
(
1
-
e
2
)
1
/
2
n
a
2
e
∂
〈
R
〉
∂
ω
,
d
i
d
t
=
cos
⁡
i
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
ω
-
1
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
Ω
,
d
Ω
d
t
=
1
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
i
,
d
ω
d
t
=
-
cos
⁡
i
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
∂
〈
R
〉
∂
i
+
(
1
-
e
2
)
1
/
2
n
a
2
e
∂
〈
R
〉
∂
e
.The semimajor axis is constant since the mean anomalyM was eliminated in the perturbing function. Following this idea, we numerically integrate the equation of motion (3), obtained from Lagrange’s planetary equations. The results are shown below:
(3)
d
e
d
t
=
K
15
4
μ
′
n
′
2
e
1
-
e
2
n
×
[
sin
2
ω
(
cos
⁡
2
D
-
cos
⁡
2
i
sin
2
D
)
iiiii
(
cos
⁡
2
D
)
-
cos
⁡
i
cos
⁡
2
ω
sin
2
D
]
,
d
i
d
t
=
K
1
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
3
4
μ
′
n
′
2
a
2
×
{
e
2
[
-
5
cos
⁡
2
ω
sin
2
D
cos
⁡
2
i
3
2
iiiiiiiii
+
3
2
(
-
sin
2
D
+
sin
2
D
cos
⁡
2
i
)
iiiiiiiii
-
5
2
sin
D
cos
⁡
2
ω
(
1
+
cos
⁡
2
i
)
]
iiiiii
+
(
sin
D
-
sin
2
D
cos
⁡
2
i
)
iiiiii
3
2
-
5
cos
⁡
2
D
sin
2
ω
cos
⁡
i
}
,
d
Ω
d
t
=
K
3
4
μ
′
n
′
2
a
2
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
×
{
(
1
+
3
2
e
2
)
(
-
sin
2
D
sin
2
i
)
iiiii
+
[
cos
⁡
2
i
cos
⁡
2
ω
sin
2
D
+
sin
i
sin
2
D
sin
2
ω
]
iiiii
×
(
5
2
e
2
)
}
,
d
ω
d
t
=
K
3
4
μ
′
n
′
2
×
{
-
a
2
sin
i
[
μ
a
(
1
-
e
2
)
]
1
/
2
iiiii
×
[
(
1
+
3
2
e
2
)
(
-
sin
2
D
sin
2
i
)
iiiiiiiiii
+
(
5
2
e
2
)
(
cos
⁡
2
i
sin
2
D
cos
⁡
2
ω
)
iiiiiiiiii
+
1
2
sin
2
i
sin
2
D
sin
2
ω
]
iiiii
+
1
-
e
2
n
[
5
(
cos
⁡
2
D
-
cos
⁡
2
i
sin
2
D
)
cos
⁡
2
ω
6
2
iiiiiiiiiiiiiiiiiiiiiii
-
5
cos
⁡
i
sin
2
D
sin
2
ω
iiiiiiiiiiiiiiiiiiiiiii
+
a
2
[
(
e
2
)
]
1
/
2
6
2
(
cos
⁡
2
D
+
cos
⁡
2
i
sin
2
D
)
]
}
,
where K
=
(
1
+
3
e
′
cos
⁡
f
′
+
3
e
′
2
cos
⁡
2
f
′
)
/
(
1
-
6
e
′
+
15
e
′
2
).This set of equations has singularities for equatorial orbits, because the inclination appears in the denominator of some of them. There is no singularity for small eccentricities, because it does not appear in the denominator. So, it is necessary to check the results very carefully against the other models, which have no singularities, when near equatorial orbits are considered. In any case, exactly equatorial orbits were not used in the present research.


## 2.2. Double-Averaged Model

Using (1) and considering an average over the eccentric anomaly of the perturbing body (Domingos et al. [16]), the equation obtained for the double-averaged disturbing function is
(4)
〈
〈
R
2
〉
〉
=
K
′
[
2
(
3
cos
⁡
2
i
-
1
)
+
3
e
2
(
3
cos
⁡
2
i
-
1
)
iiiii
+
15
e
2
sin
2
i
cos
⁡
2
ω
]
.
Now, it is possible to quantify the resulting variations in the orbital elements of the satellite. To obtain this result, again Lagrange’s planetary equations are derived using the perturbing potential shown by (4). The results are given by
(5)
d
e
d
t
=
K
′
15
μ
′
n
′
2
e
1
-
e
2
8
n
sin
2
(
i
)
sin
(
2
ω
)
,
d
i
d
t
=
K
′
-
15
μ
′
n
′
2
e
2
16
n
1
-
e
2
sin
(
2
i
)
sin
(
ω
)
,
d
ω
d
t
=
K
′
3
μ
′
n
′
2
8
n
1
-
e
2
[
(
5
cos
⁡
2
(
i
)
-
1
+
e
2
)
iiiiiiiiiiiiiiiiiiiiii
+
5
(
1
-
e
2
-
cos
⁡
2
(
i
)
)
cos
⁡
(
2
ω
)
]
,
d
Ω
d
t
=
K
′
3
μ
′
n
′
2
cos
⁡
(
i
)
8
n
1
-
e
2
[
5
e
2
cos
⁡
(
2
ω
)
-
3
e
2
-
2
]
,
d
M
d
t
=
-
K
′
μ
′
n
′
2
8
n
[
(
3
e
2
+
7
)
(
3
cos
⁡
2
(
i
)
-
1
)
iiiiiiiiiiiiiiiii
+
15
(
1
+
e
2
)
sin
2
(
i
)
cos
⁡
2
(
ω
)
]
,
where K
′
=
[
1
+
(
3
/
2
)
e
′
2
+
(
15
/
8
)
e
′
4
].This set of equations has no singularities, so it is important to use it to verify the results obtained by the previous method, in particular when near equatorial orbits are considered.


## 2.3. The Impulsive Transfers to Correct the Eccentricity

In present study, we have considered the time evolution and maneuvers to correct only inclination and eccentricity, because they were considered the most important elements that impact the orbit of the satellite (see [13, 14]). As shown by the mathematical models used here, the semimajor axis is not affected by third-body perturbations, except for short periodic oscillations of small amplitudes. The argument of periapsis and the right ascension of the ascending node do not affect very much the mission under analysis, since they do not change the regions of the celestial body that is covered by the satellite. Those variables have a secular perturbation that causes a circular motion for both angles (Solórzano and Prado [20]).To return the eccentricity of an orbit to its original value and keeping the same semimajor axis, one possibility is to make a bi-impulsive transfer to send the satellite from an orbit O2, with a given semi-major axis “a” and eccentricity “e
2”, to an orbit O1, with the same semimajor axis “a” but with eccentricity “e
1”. Another possibility would be to perform a single impulsive maneuver at the point where the orbits intersect each other. The two-impulsive maneuver will be used, because it uses the apsis of the orbits, which are better in terms of fuel consumption, and also because splitting the maneuver in two parts makes the task of correcting errors in the impulses applied easier. Since the semimajor axis is not changed by the third-body perturbations, in terms of average, it is also desired to perform orbital maneuvers that do not change this parameter. So, when correcting the eccentricity, this is the only orbital parameter that needs to be changed by the maneuver. Figure 1 shows this situation. A good solution is a Hohmann type transfer that uses the most distant apoapsis. So, it is necessary to apply the first impulse, with magnitude Δ
V
1, to the satellite when it is at the apoapsis of the orbit O2 that is a point where the distance from the center of the Earth is r
a
2. This impulse puts the satellite in a transfer orbit O
T that is an ellipse tangent to the apoapsis of O2 and to the perigee of O1, where the distance from the center is r
p
1. When the satellite reaches the periapsis of the orbit O1, the second impulse, with magnitude Δ
V
2, is applied to the satellite and this impulse puts the satellite in the final desired orbit O2. This second impulse is applied opposite to the direction of the motion of the satellite, in order to remove energy from the satellite, thus returning its semimajor axis to the original value.Figure 1
Geometry of the bi-impulsive transfer.To obtain the magnitude of the impulses, it is possible to start with the vis viva equation (6), where v is the velocity of the satellite, r is the distance between the satellite and the center of the Earth, and μ is the gravitational parameter of the Earth, whose numerical value used here is 398600 km3/s2:
(6)
v
=
2
μ
r
-
μ
a
.Based on that equation, the magnitude of the first impulse is the difference between the velocity of the satellite when it is at the pointr
a
2 in the orbit O2 and the velocity when the satellite is at the same point but in the orbit O
T. The equation is
(7)
Δ
V
1
=
2
μ
a
(
1
-
e
2
)
-
μ
a
T
-
2
μ
a
(
1
-
e
2
)
-
μ
a
,
where a
T is the semimajor axis of the transfer orbit O
T, which is given by
(8)
a
T
=
a
+
a
2
(
e
2
-
e
1
)
.The second impulse is applied at the periapsis of the nominal orbit and it is in the opposite direction with respect to the motion of the satellite, thus reducing the energy and the semimajor axis of the orbit to return it to its original value. Its magnitude is the difference between the velocities of the satellite at this point when it is in the orbitO
T and its velocity at the same point but when the satellite is at the orbit O2. The equation is
(9)
Δ
V
2
=
2
μ
a
(
1
-
e
1
2
)
-
μ
a
T
-
2
μ
a
(
1
-
e
1
2
)
-
μ
a
.The total variation of velocity required by the maneuver is given byΔ
V
T
=
Δ
V
1
+
Δ
V
2. The conversion to get the fuel mass consumed is obtained from the well-known Tsiolkovskii equation [21]:
(10)
m
0
=
m
f
e
(
Δ
V
i
/
g
0
I
sp
)
,
where g
0 is the acceleration of gravity at the surface of the Earth, assumed to be 9.8 m/s2; I
sp is the specific impulse of the engine (considered to be 340 s); the subscript i assumes the values 1 and 2 to represent both impulses; m
0 is the mass of the satellite before the impulse; and m
f is the mass of the satellite after the impulse, such that m
0
-
m
f is the fuel consumed by the maneuver. So, it is necessary to have a value for the initial mass of the satellite and the parameters of the combination engine/fuel in order to obtain the total consumption in terms of fuel mass. This equation is used for both impulses separately, updating the value of the mass after the first impulse, not to Δ
V
T, to get more accuracy.


## 2.4. The Impulsive Transfer to Correct the Inclination

In this situation, the simplest strategy is to make a single impulsive transfer to rotate the velocity of the satellite by the desired angle. This impulse has to be applied when the satellite passes by the intersection of the two planes: the one where the satellite is orbiting and the one where it should go after the maneuver. The magnitude of this impulse is given by [22]
(11)
Δ
V
i
=
2
v
sin
(
Δ
i
2
)
,
where v is the velocity of the satellite at the moment of the application of the impulse and Δ
i is the plane change required by the maneuver. Since the nominal orbit is not circular, the magnitude of this impulse depends on the true anomaly of the point of the application of the impulse, which depends on the geometry of the problem (values of the argument of the periapsis and the right ascension of the ascending node) and how the perturbation changes those variables. To avoid specific cases and considering that the eccentricity is small so that the magnitude of the velocity does not oscillate too much, it is possible to consider the worst case, which is when making the maneuver at the periapsis of the orbit, the point where the velocity is larger. In this case, the maximum velocity variation required by the maneuver that changes the inclination (Δ
V
max
⁡) is given by
(12)
Δ
V
i
max
⁡
=
2
sin
Δ
i
2
2
μ
a
(
1
-
e
1
2
)
-
μ
a
.This equation is obtained just by replacing the term in velocity by the velocity at periapsis. The conversion to fuel mass is obtained once again by using (10).


## 2.5. Low Thrust Maneuvers

The same transfers described above can be performed using low thrust. The reason to do that is the smaller fuel consumption required by the low thrust propulsion system, due to its higher specific impulse. It is assumed that the satellite can arrive and leave at any point in the final and initial orbits, respectively. The maneuver is assumed to be realized by an engine that delivers a thrust with variable direction and constant magnitude.The dynamical model assumes that the satellite does not suffer any perturbation during the maneuver. This assumption is valid because the durations of the maneuvers are very short when compared to the orbital period of the satellite. The thrust follows the rules of fixed magnitude, constant ejection velocity, and free angular motion. To solve this problem, first order necessary conditions for a local minimum are written to obtain the adjoint equations. Pontryagin’s Maximum Principle is used to give the control angles at each instant. The approach used here is the same hybrid approach showed in Biggs [8, 9] and later used by Gomes and Prado [10], where all the details of the method are shown, including the variables used and the equations of motion. Constraints can be added, if necessary. The solution of the problem is the time histories of the direction of the thrust and the instants to start and finish the application of the thrust. The fuel consumed and the times that the engine is active are also part of the solution.To allow more savings in the fuel consumption, two situations are simulated when using the low thrust propulsion system. The first one performs the maneuver in a single propulsive arc and the second one allows the possibility of two thrusting arcs, so the engine is turned on and off twice. This maneuver allows a better use of the fuel, by applying the thrust in the most favorable parts of the orbit, which are the apoapsis regions. Figure2 shows the results of a two-arc maneuver. The horizontal axis shows the true anomaly of the spacecraft where the engine should be fired and the vertical axis shows the direction of the thrust, both in degrees. One figure is made for each propulsion arc. This example is the maneuver to correct an error of 0.01 in the eccentricity, thus changing it from 0.02 to 0.01. The other situations are similar, so only the fuel consumption and the duration of the maneuvers are shown in the tables and the figures are omitted. The duration of the application of the thrust is short, so the plots of the trajectories are not interesting, since they are Keplerian orbits. After that, the errors of the propulsion system were considered, using a very simple approach. The solution is obtained as a set of points, which represents the whole trajectory. Then, we inserted a fixed error in the magnitude (5%) and in the direction of the solution obtained (2 degrees) at every point obtained as a solution of the maneuver. The constraints of reaching the final orbit specified were the same of the maneuvers without errors. So, using this simple approach, it was possible to see that the errors in the fuel consumed were in the order of 5–7%. This result is expected because the maneuvers have small magnitudes and so there is an almost linear dependence between errors and extra consumption.Low thrust to correct the eccentricity for the maneuver from 0.02 to 0.01 using two arcs of propulsion.
(a)
(b)


## 3. Results

Based on the mathematical models described in the previous sections, performed some numerical simulations are now performed using (3), for the single-averaged model, and (5), for the double-averaged model, to verify the evolution of the eccentricity and inclination due to the third-body perturbations for high-altitude circular orbits (semimajor axis of 42284 km). These orbits have the altitude of the geostationary orbits around the Earth, and they are chosen because they have strong effects from the third-body perturbations. Thus, the choice of high-altitude orbits in addition to the system with mass parameters equivalent to the Earth-Moon system with variable eccentricity can better show the validity of the results presented. Near equatorial (i
0
=
10
-
3 degrees), near the critical inclination of the third-body perturbation (i
0
=
39
°~0.6803 radians), and high inclined orbits (i
0
=
80
°~1.3955 radians) are used for the simulations. The initial orbit of the satellite is assumed to have eccentricity e
0
=
0.01 and right ascension of the ascending node and argument of periapsis equal to zero.


### 3.1. Finding the Consumptions and the Durations of the Maneuvers

After verifying the evolution of the eccentricity and inclination for ~35 years (near 478 orbits of the perturbing body), it is possible to obtain the exact times where those variables reach the limits of variation accepted by the mission from the original orbit of the satellite. Several values are used for those limits of variation for the eccentricity,Δ
e, and for the inclination, Δ
i, thus studying the influence of this parameter. Those parameters represent different requirements that can be made by different missions. The values used are 0.0005, 0.001, 0.005, 0.01, 0.02, and 0.05 for Δ
e and 0.0001, 0.0005, 0.001, and 0.005 radians for Δ
i. When those limits are reached (e
0
±
Δ
e and/or i
0
±
Δ
i), a maneuver is used to return the satellite to its original orbit. The evolutions of the eccentricity and the inclination are shown in Figures 3 to 5.These figures show the evolution of the eccentricity as a function of time for near equatorial orbits. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). In (a), the results are for the full elliptic restricted three-body problem. (b) is for single-averaged model and (c) is for double-averaged model. The dashed curve represents the (t, e) pair to perform orbital maneuvers, depending on the limit of the eccentricity variation allowed by the mission.
(a)
(b)
(c)As shown in Section2, considering the instantaneous propulsive system, the maneuvers used to correct the deviations in eccentricity are the bi-impulsive planar version and the maneuver used to correct the inclination is the single impulsive maneuver. When the low trust system is considered, the maneuver to correct the errors in eccentricity is performed under two different hypotheses: with one or two propulsion arcs. The same system corrects the errors in inclination, but now using only one propulsion arc. Tables 1, 2, 3, and 4 show the consumptions, in terms of kg of fuel consumed, for each maneuver, which have fixed values for each limit chosen from the list of possibilities considered and the durations of the propulsion arcs for the low thrust maneuvers, expressed in minutes. These results are shown in Tables 5, 6, 7, 8, and 9. In the information that follows, “Impulse 1” is obtained from (7), “Impulse 2” is obtained from (9), the “total Impulse” is the sum of both impulses, and the values for the mass of propellant are obtained from (10). The low thrust maneuvers are obtained from the algorithm described in [10].Table 1
Fuel consumption for the maneuvers to correct the eccentricity.
Δ
e
Impulse 1 (m/s)
Impulse 2 (m/s)
Total impulse forΔ
e (m/s)
Fuel consumed (kg)
LT with 1 arc (kg)
Arc 1 LT (kg)
Arc 2 LT (kg)
Total LT with 2 arcs (kg)
0.0005
0.388
0.380
0.769
0.231
0.067
0.030
0.030
0.060
0.0010
0.777
0.761
1.537
0.461
0.135
0.060
0.060
0.120
0.0050
3.889
3.793
7.682
2.308
0.687
0.301
0.300
0.601
0.0100
7.793
7.563
15.356
4.619
1.350
0.605
0.605
1.210
0.0200
15.644
15.033
30.677
9.249
2.714
1.210
1.210
2.421
0.0500
39.547
36.904
76.451
23.210
6.779
3.040
3.040
6.080Table 2
Time required by the low thrust maneuvers to correct the eccentricity.
Limit reached forΔ
e
Thrusting time Arc 1 (min)
Thrusting time for Arc 2 (min)
Thrusting time with 2 arcs (min)
Thrusting time with 1 arc (min)
0.0005
0.094
0.094
0.188
0.210
0.0010
0.188
0.188
0.376
0.419
0.0050
0.940
0.938
1.878
2.091
0.0100
1.89
1.89
3.78
4.22
0.0200
3.781
3.781
7.562
8.439
0.0500
9.497
9.497
18.994
21.285Table 3
Fuel consumption for the maneuvers to correct the inclination.
Δ
i (rad)
Impulse forΔ
i (m/s)
Fuel consumed (kg)
LT with 1 arc (kg)
0.0001
0.310
0.0932
0.0244
0.0005
1.553
0.4660
0.1221
0.0010
3.106
0.9330
0.2443
0.0050
15.523
4.6700
1.2235Table 4
Time required by the low thrust maneuvers to correct the inclination.
Δ
i (rad)
Thrusting time (min)
0.0001
0.085
0.0005
0.426
0.0010
0.852
0.0050
4.270Table 5
Times when the limits are reached and fuel consumption per year of mission.
Eccentricity of the perturbing body
Time to reach theΔ
e and Δ
i limits (year)
Yearly consumption for impulsive maneuvers (kg)
e
′
=
0.0
Δ
e
=
0.0005
full model = 22.82
0.0101
Δ
e
=
0.001   full model = 34.78
0.0133
e
′
=
0.1
Δ
e
=
0.0005
full model = 19.93
0.0116
Δ
e
=
0.001   full model = 29.95
0.0154
e
′
=
0.2
Δ
e
=
0.0005 full model = 4.76
0.0485
Δ
e
=
0.001   full model = 6.71
0.0687
Δ
e
=
0.005 full model = 16.62
0.1389
Δ
e
=
0.01  full model = 26.08
0.1771Table 6
Times when the limits are reached and average yearly fuel consumption for the mission.
Eccentricity of the perturbing body
Time to reach the Δe and Δi limits (year)
Yearly consumption for impulsive maneuvers (kg)
Fore
′
=
0.0
full model = 5.50
Δ
e
=
0.0005   single averaged model = 5.75      double averaged model = 5.76      full model = 7.63
Δ
e
=
0.001 single averaged model = 8.23     double averaged model = 8.24     full model = 18.11
Δ
e
=
0.005
single averaged model = 20.12     double averaged model = 20.13      full model = 28.0
Δ
e
=
0.01
single averaged model = 31.30     double averaged model = 31.30      full model = 0.012
Δ
i
=
0.0001
single averaged model = 0.10     double averaged model = 22.92
0.04200.04020.0401 0.06040.05600.0559  0.12740.11470.1146 0.16500.14760.1476 7.76670.93200.0041
Fore
′
=
0.1
full model = 8.77
Δ
e
=
0.0005
single averaged model = 5.68     double averaged model = 5.67      full model =  12.84
Δ
e
=
0.001
single averaged model = 8.11      double averaged model = 8.11       full model = 30.98
Δ
e
=
0.005
single averaged model = 19.82     double averaged model = 19.83
Δ
e
=
0.01  single averaged model = 30.82     double averaged model = 30.84      full model = 0.012
Δ
i
=
0.0001  single averaged model = 0.083     double averaged model = 25.58
0.02630.04060.0407 0.03590.05680.0568  0.07450.11640.11630.01500.0150 7.76671.12290.0036
Fore
′
=
0.2
full model = 13.64
Δ
e
=
0.0005 single averaged model = 5.31     double averaged model = 5.41       full model = 21.10
Δ
e
=
0.001
single averaged model = 7.62     double averaged model = 7.75 Δ
e
=
0.005
single averaged model = 18.59     double averaged model = 18.94 Δ
e
=
0.01 single averaged model = 28.98     double averaged model = 29.45      full model = 0.012
Δ
i
=
0.0001     single averaged model = 0.071     double averaged model =  21.57 Δ
i
=
0.0005
single averaged model  = 13.03 Δ
i
=
0.001 single averaged model = 22.62
0.01690.04350.0427  0.02180.06050.05950.12420.1219 0.15940.15687.76671.31270.0043 0.03580.0412Table 7
Times when the limits are reached and the average yearly fuel consumption for the mission assuming a circular orbit for the perturbing body.
Limit reached
Time to reach the limit (year)
Yearly consumption for impulsive maneuvers (kg)
Δ
e
=
0.0005
Full model = 3.69Single-averaged model = 3.65Double-averaged model = 3.64
0.06260.06330.0635
Δ
e
=
0.001
Full model = 5.13Single-averaged model = 5.16Double-averaged model = 5.16
0.08990.08930.0893
Δ
e
=
0.005
Full model = 11.44Single-averaged model = 11.63 Double-averaged model = 11.63
0.20170.19850.1985
Δ
e
=
0.01
Full model = 16.01Single-averaged model = 16.40 Double-averaged model = 16.41
0.28850.28160.2815
Δ
e
=
0.02
Full model = 22.19Single-averaged model = 22.70 Double-averaged model = 22.71
0.41680.040740.4073
Δ
e
=
0.05
Full model = 32.36Single-averaged model = 33.20Double-averaged model = 33.21
0.71720.69910.6989
Δ
i
=
0.0001
Full model = 0.012Single-averaged model = 0.012Double-averaged model = 30.013
7.76677.76670.0031
Δ
i
=
0.0005
Full model = 29.98Single-averaged model = 32.15
0.01550.0145Table 8
Times when the limits are reached and the fuel consumption per year of mission using the value of 0.1 for the eccentricity of the orbit of the perturbing body.
Limit reached
Time to reach the limit (year)
Yearly consumption for impulsive maneuvers (kg)
Δ
e
=
0.0005
Full model = 3.29Single-averaged model = 3.58 Double-averaged model = 3.58
0.07020.06450.0645
Δ
e
=
0.001
Full model = 4.61Single-averaged model = 5.08  Double-averaged model = 5.08
0.10000.09070.0907
Δ
e
=
0.005
Full model =  10.38Single-averaged model = 11.44 Double-averaged model = 11.46
0.22240.20170.2014
Δ
e
=
0.01
Full model = 14.68Single-averaged model = 16.14double-averaged model = 16.15
0.31460.28620.2860
Δ
e
=
0.02
Full model = 20.46Single-averaged model = 22.35  Double-averaged model = 22.37
0.45210.41380.4135
Δ
e
=
0.05
Full model = 30.22Single-averaged model = 32.68Double-averaged model = 32.71
0.76800.71020.7096
Δ
i
=
0.0001
Full model = 0.012Single-averaged model = 0.009Double-averaged model = 29.56
7.766710.35560.0032
Δ
i
=
0.0005
Single-averaged model = 29.85
0.0156Table 9
Times when the limits are reached and the fuel consumption per year of mission assuming the value of 0.2 for the eccentricity of the orbit of the perturbing body.
Limit reached
Time to reach the limit (year)
Yearly consumption for impulsive maneuvers (kg)
Δ
e
=
0.0005
Full model = 3.05Single-averaged model = 3.37Double-averaged model = 3.42
0.07570.06860.0675
Δ
e
=
0.001
Full model = 4.0Single-averaged model = 4.78Double-averaged model = 4.85
0.11530.09640.0951
Δ
e
=
0.005
Full model = 9.1Single-averaged model = 10.77   Double-averaged model = 10.94
0.25360.21430.2110
Δ
e
=
0.01
Full model = 12.95Single-averaged model = 15.18 Double-averaged model = 15.42
0.35670.30430.2995
Δ
e
=
0.02
Full model = 18.28Single-averaged model = 21.01 Double-averaged model = 21.36
0.50600.44020.4330
Δ
e
=
0.05
Full model = 27.38Single-averaged model = 30.74 Double-averaged model = 31.24
0.84770.75500.7430
Δ
i
=
0.0001
Full model = 0.012Single-averaged model = 0.008 Double-averaged model = 28.24
7.766711.6500.0033
Δ
i
=
0.0005
Full model = 27.5Single-averaged model = 17.91Double-averaged model = 35.71
0.01690.02600.0130
Δ
i
=
0.001
Full model = 27.73Single-averaged model = 30.0
0.03360.0311To obtain the mass of fuel required by each maneuver, as an example, a satellite with an initial mass of 1000 kg is used. The hydrazine is assumed to be the fuel for the impulsive maneuvers, withI
sp
=
340 seconds. For the low thrust maneuvers, a force of 10 newtons and a specific impulse of the 1300 s were used for the propulsive system.The durations of the propulsion arcs were not too large, going from about 0.2 minutes to nearly 21 minutes, depending on the amplitude of the maneuver. Since all the corrections are small, the consumptions are nearly proportional to the amplitudes of the maneuvers. The maneuvers using two propulsion arcs divide the fuel consumption near equally between two arcs, sometimes showing a little larger consumption in the first thrusting arc, because this arc has distances to the Earth that are a little smaller, so it spends more fuel for the maneuver. Anyway, these differences in distances from the Earth are too small and not always implied in difference in consumption between both arcs considering the number of digits shown in the tables. The maneuvers with only one arc take more time and require more fuel, since the thrust is not applied in the optimal parts of the orbit that are the regions near the apoapsis. Figure2 shows the direction of the thrust to correct the eccentricity from 0.02 to 0.01, using two arcs of propulsion, as an example of the calculations made in the present research. Note that the direction has near linear variations for both arcs, close to the direction of motion for the first arc and close to the direction opposite to the motion for the second arc. It means that the low thrust maneuvers follow the main geometry of the impulsive maneuvers, but replacing the instantaneous variation of velocity by a continuous force acting in the satellite during some time. The main advantage of this system is the fuel economy, which is about four times smaller using the propulsion system showed here, due to the higher value of the specific impulse. It can be even larger using more efficient propulsive system, which means fuel with higher specific impulses. When using two arcs of propulsion, there is a fuel economy of around 10% with respect to the use of a single arc. The reason for these extra savings is the application of the thrust in the regions where they are more efficient, that is, the apoapsis of the trajectories involved. At those points, the gravitational force of the Earth is smaller, thus allowing larger effects due to the application of the thrust.Tables3 and 4 show the equivalent results of the maneuvers used to correct the inclination. Those maneuvers require less fuel than the ones to correct eccentricities, because the limits to perform an orbital correction are very small for the inclination. The fuel consumption goes from 0.0244 to 1.2235 kg and the duration of the propulsion goes from 0.08 to 4.27 minutes. There is an almost linear relation between the desired angle change and the fuel consumption, due to the small values involved in the variations of the inclination. In this situation, only one propulsion arc was used for each maneuver, since it is the best solution, even for the impulsive maneuvers.


### 3.2. Near Equatorial Orbits

Figure3 shows the evolution of the eccentricity as a function of time for the near equatorial orbits (i
0
=
10
-
3 degrees, since real equatorial orbits show singularities in the single-averaged equations of motion). The times are shown in years. Figures 3(a), 3(b), and 3(c) show the results of the cases of the full elliptic restricted problem and single- and double-averaged models, respectively. On the top of each figure, the corresponding eccentricity of the perturbing body (e
′) is shown. This parameter is used in order to study its effects on this problem and it represents the same system used before. It is clear that none of the specified limits is reached when using the single-averaged and the double-averaged techniques (see Figures 3(b) and 3(c)). The similar results, in terms of reaching the limits to make a maneuver, which comes from both averaged models, are expected because both averaged models have results that are very close to each other when the perturbing body is in a circular orbit (Domingos et al. [17]). This effect decreases with the increase in the eccentricity of the primaries. In general, both averaged models have results that are much closer to each other than closer to those of the full model. It is visible from the full model that, during the integration time, the evolution of the eccentricity has oscillations that increase in amplitude with the increase in the eccentricity of the primaries. The averaged models are not able to predict this property. This fact changes the times when the maneuvers are required, so it has a strong influence on the estimation of the fuel consumption for the station keeping maneuvers.On the other hand, the results show that, for the full elliptic restricted three-body problem, in Figure3(a), the limits of variation of the eccentricity 0.0005, 0.001, 0.005, and 0.01 are reached, so it is necessary to perform orbital maneuvers. As an example, in Figure 3(a), we can see that when e
′
=
0.2 only the limit of variation for the eccentricity Δ
e
=
0.05 is not reached during the timescale of the numerical integrations. As a result, in this latter case, we do evidence that there are four points where it is required to perform orbital maneuvers to correct the orbit. The inclinations have zero variations during this period, so it is not required to apply any propulsion system to correct the orbit. For cases with e
′
<
0.2, the results shown in Figure 3 evidence that only two orbital maneuvers are required to correct the orbit of the satellite. So, the general conclusion is that averaged models are not good enough to predict the evolution of the eccentricity and the orbital maneuvers required by the mission for near equatorial orbits, for those values of the accuracy required in eccentricity for the mission, because there is a variation of up to 0.015 in eccentricity and an average yearly consumption of up to 0.1771 kg that are not predicted by both averaged models. This average is obtained by dividing the total fuel expenditures during the mission by the total time considered for the evolution of the orbits. So, even if the spacecraft spends more than one year without the need of a maneuver, it is still possible to compute a “yearly average consumption.”Table5 shows the times, in years, when the limits of eccentricity are reached for each situation, as well as the estimated fuel consumption per year of mission for the several limits considered for the eccentricity and for the values of the eccentricity of the primaries.The fuel consumptions are not too large, going from 0.0101 to 0.1771 kg per year. It is also noted that there is a strong dependence of the fuel consumption on the eccentricity of the primaries, in the order of five times, when comparing the circular case to the case with eccentricity 0.2, for the same limits in the eccentricity variations. In the worst case considered here, with eccentricity of primaries of 0.2, the times for performing an orbital maneuver goes from 4.76 years, if the limit used for the variation of eccentricity isΔ
e
=
0.0005, to 26.08 years, if the limit used for the variation of eccentricity is Δ
e
=
0.01.


### 3.3. Near Critical Inclination Orbits

Figure4 shows the evolution of the eccentricity and inclination, respectively, as a function of time when using the three models for the case of near critical orbits (i
0
=
39
°~0.6803 radians). In these figures, Figures 4(a), 4(b), and 4(c) also show the results of the cases of the full elliptic restricted problem and the single- and double-averaged models, respectively. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). Regarding the full elliptic restricted three-body problem, Figure 4(a), the results show that the inclination reaches the first limit of 0.0001 rad at 0.012 years, while the eccentricity reached several limits up to 0.02 in 28 years for the circular case for the perturbing body. It is also interesting to note that, when the eccentricity increases, there is a decrease in the inclination. By focusing on the value of the eccentricity of the perturbing body, we can see that the inclination remains close to the initial value all the time.Evolution of the eccentricity and inclination as a function of time for the casei
0
=
39
° (~0.6803 radians). In (a), the results are for the full elliptic restricted three-body problem. (b) is for the single-averaged model and (c) is for the double-averaged model. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). The dashed curves represent the (t, e) and (t, i) pairs for eccentricity and inclination, respectively, to perform orbital maneuvers, depending on the limit of the eccentricity and inclination variations allowed by the mission. In the plots for the inclination, it showed the dashed lines for the limits of variation for times higher than 1 year for better graphic visibility.
(a)
(b)
(c)Figure 5
Evolution of the eccentricity as a function of time for the casei
0
=
80
° (~1.3955 radians). The results are for the full elliptic restricted three-body problem. Similar results have been obtained for single- and double-averaged models. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). The dashed curves represent the (t, e) pair to perform orbital maneuvers for each limit assumed for the variation in eccentricity.It can also be seen, in Figures4(b) and 4(c), that the evolution of the eccentricity has a similar behavior for single- and double-averaged models and that they do not suffer significant changes with the increase in e
′. In the two models, the times for reaching each eccentricity limit are very close. As an example, the eccentricity 0.02 was reached in approximately 30 years in both cases. The same does not occur when the study of the inclination is made, because it reaches its limits in different times for each model.The evolutions of the orbits also show another interesting feature of the averaged models. Opposite to what happens for the eccentricity, when the inclination is studied, the averaged models predict variations that are stronger than the ones observed in the real restricted three-body problem. It means that averaged models will predict orbital maneuvers and so will fuel consumption, in situations where they are really not required. Still, regarding the orbital evolution predicted by the three models, the results shown here evidence some differences. The decrease in the eccentricity shown by the full model for the case where the eccentricity of the primaries is 0.2 is not predicted by both averaged models. The dependence of the results on the eccentricity of the primaries is much stronger when considering the elliptical restricted three-body case than when it is predicted by the averaged models, where the evolutions are similar to each other in all cases studied. Both averaged models are very similar to each other regarding predictions for the eccentricity, but they have differences from each other for the evolution of the inclination.Considering the limits of eccentricity and inclinations reached during the integration period, it is clear that it is necessary to perform orbital maneuvers in most of the cases. Table6 shows the times, in years, when these limits are reached for each situation, as well as the estimated fuel consumption per year of the mission for the several limits considered for the eccentricity and inclination and for different values of the eccentricity of the primaries.Limits up to 0.01 in eccentricity are reached for all the models. The value 0.02 was never reached. Differences in the prediction of the times for the maneuvers and the average yearly fuel consumption to correct eccentricities are of the order of 10%, so not too large. Regarding the inclination, for the circular case only, the limitΔ
i
=
0.0001 was reached by the three models but the differences in the times of the maneuver and the average yearly fuel consumption are very large among the models, with the elliptical restricted problem requiring much more fuel. When the perturbing body is not in a circular orbit, the full model takes more time to reach the limits of the eccentricity, so the averaged models predict earlier times for the maneuvers, so more yearly fuel consumption occurs when compared to the real cases. Some larger limits for the eccentricity are not reached by the full model. For the inclination, the full model predicted again higher fuel consumptions and faster orbital maneuvers than did the averaged models, so those models under evaluate the fuel consumption to correct the inclination. It is cheaper to keep lower values of the limits for the eccentricity, because the velocity of the increase in the eccentricity gets higher with the eccentricity. So, it is better to perform the maneuvers earlier, not allowing the eccentricity to grow too much. For the inclination, those costs alternate from one limit to the other.So, the general conclusion is that averaged models are also not good enough to predict the evolution of the eccentricity and inclination also in this region of near critical inclination orbits because there are considerable differences in the prediction of the full model and both averaged models.


### 3.4. High Inclined Orbits

Figures5 and 6 show the temporal evolutions of the eccentricity and inclination, respectively, for the case of high inclined orbits, considered here to have i
0
=
80
° (~1.3955 radians). The same representations for the models given in Figure 4 are used here.Evolution of the inclination as a function of time for the casei
0
=
80
° (~1.3955 radians). In (a), the results are for the full elliptic restricted three-body problem. (b) is for the single-averaged model and (c) is for the double-averaged model. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). The dashed curves represent the (t, i) pair to perform orbital maneuvers, for each limit assumed for the inclination.
(a)
(b)
(c)In Figure5, only the results of the evolution of the eccentricity when considering the full elliptic restricted three-body problem are showed, because similar results have been obtained for the single- and double-averaged models. This similarity is an expected result of initial inclinations above the critical value (see Domingos et al. [17]). The eccentricity oscillates with large amplitude, reaching the maximum limit of 0.06 (a variation of 0.05 from the initial value 0.01) near the time 32.5 years, while the inclination reaches only the smaller values of the limit of variation (<10−3 rad) during the integration time. We note that, for the present case of i
0
=
80
°, all the limits of eccentricity are reached when using the full problem and both averaged models, so it is necessary to perform orbital maneuvers in all the situations. Considering the inclinations, see Figure 6, the three limits of inclination are only reached for e
′
=
0.2. The single-averaged model predicts earlier maneuvers than the full model does, so it overestimates the fuel consumption. The double-averaged model follows the full model better in this topic. The eccentricity of the primaries always increases the perturbations, so it also increases the average yearly fuel consumption. Tables 5 to 9 show the times when the limits are reached for each situation, as well as the estimated fuel consumption per year of mission for the several limits considered for the eccentricity and inclination for the values of 0.0, 0.1, and 0.2 for the eccentricity of the primaries.To correct the eccentricity, the times predicted by the three models are very similar, which generates very close fuel consumption estimations. There is a strong relation between the average yearly consumption and the limits imposed on the variation in eccentricity. The cost using the limitΔ
e
=
0.05 is ten times more than the cost when using the limit Δ
e
=
0.0005. So, the more frequent the maneuvers are performed the lower the annual cost is. Again, this is a consequence of the increase in the speed of the variation of the eccentricity with respect to the initial eccentricity. When considering the corrections in the inclination, the opposite occurs, so allowing larger limits reduces the total cost of the maneuvers per year of mission.In this situation, the corrections in eccentricity are made in times that still have similar predictions by the three models, but differences around 10% start to occur in some cases, generating fuel consumption estimates with differences of the same order. The relation between the yearly consumption and the limits imposed on the variation in eccentricity is kept in similar values, being more economical by an order of ten when using the strategy of making frequent maneuvers. The same previous pattern observed for the circular case occurs for the corrections in the inclination, and allowing larger limits still reduces the total cost of the maneuvers per year of mission.The same analysis made before applies here. The only difference is that there is an increase in the differences of the times to perform maneuvers and in the fuel consumptions, as a consequence, between the full and the averaged models, going to differences close to 15% now. It is a consequence of the fact that the averaged models loss quality when the eccentricity of the primaries increases, as shown in detail by Domingos et al. [17].


## 3.1. Finding the Consumptions and the Durations of the Maneuvers

After verifying the evolution of the eccentricity and inclination for ~35 years (near 478 orbits of the perturbing body), it is possible to obtain the exact times where those variables reach the limits of variation accepted by the mission from the original orbit of the satellite. Several values are used for those limits of variation for the eccentricity,Δ
e, and for the inclination, Δ
i, thus studying the influence of this parameter. Those parameters represent different requirements that can be made by different missions. The values used are 0.0005, 0.001, 0.005, 0.01, 0.02, and 0.05 for Δ
e and 0.0001, 0.0005, 0.001, and 0.005 radians for Δ
i. When those limits are reached (e
0
±
Δ
e and/or i
0
±
Δ
i), a maneuver is used to return the satellite to its original orbit. The evolutions of the eccentricity and the inclination are shown in Figures 3 to 5.These figures show the evolution of the eccentricity as a function of time for near equatorial orbits. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). In (a), the results are for the full elliptic restricted three-body problem. (b) is for single-averaged model and (c) is for double-averaged model. The dashed curve represents the (t, e) pair to perform orbital maneuvers, depending on the limit of the eccentricity variation allowed by the mission.
(a)
(b)
(c)As shown in Section2, considering the instantaneous propulsive system, the maneuvers used to correct the deviations in eccentricity are the bi-impulsive planar version and the maneuver used to correct the inclination is the single impulsive maneuver. When the low trust system is considered, the maneuver to correct the errors in eccentricity is performed under two different hypotheses: with one or two propulsion arcs. The same system corrects the errors in inclination, but now using only one propulsion arc. Tables 1, 2, 3, and 4 show the consumptions, in terms of kg of fuel consumed, for each maneuver, which have fixed values for each limit chosen from the list of possibilities considered and the durations of the propulsion arcs for the low thrust maneuvers, expressed in minutes. These results are shown in Tables 5, 6, 7, 8, and 9. In the information that follows, “Impulse 1” is obtained from (7), “Impulse 2” is obtained from (9), the “total Impulse” is the sum of both impulses, and the values for the mass of propellant are obtained from (10). The low thrust maneuvers are obtained from the algorithm described in [10].Table 1
Fuel consumption for the maneuvers to correct the eccentricity.
Δ
e
Impulse 1 (m/s)
Impulse 2 (m/s)
Total impulse forΔ
e (m/s)
Fuel consumed (kg)
LT with 1 arc (kg)
Arc 1 LT (kg)
Arc 2 LT (kg)
Total LT with 2 arcs (kg)
0.0005
0.388
0.380
0.769
0.231
0.067
0.030
0.030
0.060
0.0010
0.777
0.761
1.537
0.461
0.135
0.060
0.060
0.120
0.0050
3.889
3.793
7.682
2.308
0.687
0.301
0.300
0.601
0.0100
7.793
7.563
15.356
4.619
1.350
0.605
0.605
1.210
0.0200
15.644
15.033
30.677
9.249
2.714
1.210
1.210
2.421
0.0500
39.547
36.904
76.451
23.210
6.779
3.040
3.040
6.080Table 2
Time required by the low thrust maneuvers to correct the eccentricity.
Limit reached forΔ
e
Thrusting time Arc 1 (min)
Thrusting time for Arc 2 (min)
Thrusting time with 2 arcs (min)
Thrusting time with 1 arc (min)
0.0005
0.094
0.094
0.188
0.210
0.0010
0.188
0.188
0.376
0.419
0.0050
0.940
0.938
1.878
2.091
0.0100
1.89
1.89
3.78
4.22
0.0200
3.781
3.781
7.562
8.439
0.0500
9.497
9.497
18.994
21.285Table 3
Fuel consumption for the maneuvers to correct the inclination.
Δ
i (rad)
Impulse forΔ
i (m/s)
Fuel consumed (kg)
LT with 1 arc (kg)
0.0001
0.310
0.0932
0.0244
0.0005
1.553
0.4660
0.1221
0.0010
3.106
0.9330
0.2443
0.0050
15.523
4.6700
1.2235Table 4
Time required by the low thrust maneuvers to correct the inclination.
Δ
i (rad)
Thrusting time (min)
0.0001
0.085
0.0005
0.426
0.0010
0.852
0.0050
4.270Table 5
Times when the limits are reached and fuel consumption per year of mission.
Eccentricity of the perturbing body
Time to reach theΔ
e and Δ
i limits (year)
Yearly consumption for impulsive maneuvers (kg)
e
′
=
0.0
Δ
e
=
0.0005
full model = 22.82
0.0101
Δ
e
=
0.001   full model = 34.78
0.0133
e
′
=
0.1
Δ
e
=
0.0005
full model = 19.93
0.0116
Δ
e
=
0.001   full model = 29.95
0.0154
e
′
=
0.2
Δ
e
=
0.0005 full model = 4.76
0.0485
Δ
e
=
0.001   full model = 6.71
0.0687
Δ
e
=
0.005 full model = 16.62
0.1389
Δ
e
=
0.01  full model = 26.08
0.1771Table 6
Times when the limits are reached and average yearly fuel consumption for the mission.
Eccentricity of the perturbing body
Time to reach the Δe and Δi limits (year)
Yearly consumption for impulsive maneuvers (kg)
Fore
′
=
0.0
full model = 5.50
Δ
e
=
0.0005   single averaged model = 5.75      double averaged model = 5.76      full model = 7.63
Δ
e
=
0.001 single averaged model = 8.23     double averaged model = 8.24     full model = 18.11
Δ
e
=
0.005
single averaged model = 20.12     double averaged model = 20.13      full model = 28.0
Δ
e
=
0.01
single averaged model = 31.30     double averaged model = 31.30      full model = 0.012
Δ
i
=
0.0001
single averaged model = 0.10     double averaged model = 22.92
0.04200.04020.0401 0.06040.05600.0559  0.12740.11470.1146 0.16500.14760.1476 7.76670.93200.0041
Fore
′
=
0.1
full model = 8.77
Δ
e
=
0.0005
single averaged model = 5.68     double averaged model = 5.67      full model =  12.84
Δ
e
=
0.001
single averaged model = 8.11      double averaged model = 8.11       full model = 30.98
Δ
e
=
0.005
single averaged model = 19.82     double averaged model = 19.83
Δ
e
=
0.01  single averaged model = 30.82     double averaged model = 30.84      full model = 0.012
Δ
i
=
0.0001  single averaged model = 0.083     double averaged model = 25.58
0.02630.04060.0407 0.03590.05680.0568  0.07450.11640.11630.01500.0150 7.76671.12290.0036
Fore
′
=
0.2
full model = 13.64
Δ
e
=
0.0005 single averaged model = 5.31     double averaged model = 5.41       full model = 21.10
Δ
e
=
0.001
single averaged model = 7.62     double averaged model = 7.75 Δ
e
=
0.005
single averaged model = 18.59     double averaged model = 18.94 Δ
e
=
0.01 single averaged model = 28.98     double averaged model = 29.45      full model = 0.012
Δ
i
=
0.0001     single averaged model = 0.071     double averaged model =  21.57 Δ
i
=
0.0005
single averaged model  = 13.03 Δ
i
=
0.001 single averaged model = 22.62
0.01690.04350.0427  0.02180.06050.05950.12420.1219 0.15940.15687.76671.31270.0043 0.03580.0412Table 7
Times when the limits are reached and the average yearly fuel consumption for the mission assuming a circular orbit for the perturbing body.
Limit reached
Time to reach the limit (year)
Yearly consumption for impulsive maneuvers (kg)
Δ
e
=
0.0005
Full model = 3.69Single-averaged model = 3.65Double-averaged model = 3.64
0.06260.06330.0635
Δ
e
=
0.001
Full model = 5.13Single-averaged model = 5.16Double-averaged model = 5.16
0.08990.08930.0893
Δ
e
=
0.005
Full model = 11.44Single-averaged model = 11.63 Double-averaged model = 11.63
0.20170.19850.1985
Δ
e
=
0.01
Full model = 16.01Single-averaged model = 16.40 Double-averaged model = 16.41
0.28850.28160.2815
Δ
e
=
0.02
Full model = 22.19Single-averaged model = 22.70 Double-averaged model = 22.71
0.41680.040740.4073
Δ
e
=
0.05
Full model = 32.36Single-averaged model = 33.20Double-averaged model = 33.21
0.71720.69910.6989
Δ
i
=
0.0001
Full model = 0.012Single-averaged model = 0.012Double-averaged model = 30.013
7.76677.76670.0031
Δ
i
=
0.0005
Full model = 29.98Single-averaged model = 32.15
0.01550.0145Table 8
Times when the limits are reached and the fuel consumption per year of mission using the value of 0.1 for the eccentricity of the orbit of the perturbing body.
Limit reached
Time to reach the limit (year)
Yearly consumption for impulsive maneuvers (kg)
Δ
e
=
0.0005
Full model = 3.29Single-averaged model = 3.58 Double-averaged model = 3.58
0.07020.06450.0645
Δ
e
=
0.001
Full model = 4.61Single-averaged model = 5.08  Double-averaged model = 5.08
0.10000.09070.0907
Δ
e
=
0.005
Full model =  10.38Single-averaged model = 11.44 Double-averaged model = 11.46
0.22240.20170.2014
Δ
e
=
0.01
Full model = 14.68Single-averaged model = 16.14double-averaged model = 16.15
0.31460.28620.2860
Δ
e
=
0.02
Full model = 20.46Single-averaged model = 22.35  Double-averaged model = 22.37
0.45210.41380.4135
Δ
e
=
0.05
Full model = 30.22Single-averaged model = 32.68Double-averaged model = 32.71
0.76800.71020.7096
Δ
i
=
0.0001
Full model = 0.012Single-averaged model = 0.009Double-averaged model = 29.56
7.766710.35560.0032
Δ
i
=
0.0005
Single-averaged model = 29.85
0.0156Table 9
Times when the limits are reached and the fuel consumption per year of mission assuming the value of 0.2 for the eccentricity of the orbit of the perturbing body.
Limit reached
Time to reach the limit (year)
Yearly consumption for impulsive maneuvers (kg)
Δ
e
=
0.0005
Full model = 3.05Single-averaged model = 3.37Double-averaged model = 3.42
0.07570.06860.0675
Δ
e
=
0.001
Full model = 4.0Single-averaged model = 4.78Double-averaged model = 4.85
0.11530.09640.0951
Δ
e
=
0.005
Full model = 9.1Single-averaged model = 10.77   Double-averaged model = 10.94
0.25360.21430.2110
Δ
e
=
0.01
Full model = 12.95Single-averaged model = 15.18 Double-averaged model = 15.42
0.35670.30430.2995
Δ
e
=
0.02
Full model = 18.28Single-averaged model = 21.01 Double-averaged model = 21.36
0.50600.44020.4330
Δ
e
=
0.05
Full model = 27.38Single-averaged model = 30.74 Double-averaged model = 31.24
0.84770.75500.7430
Δ
i
=
0.0001
Full model = 0.012Single-averaged model = 0.008 Double-averaged model = 28.24
7.766711.6500.0033
Δ
i
=
0.0005
Full model = 27.5Single-averaged model = 17.91Double-averaged model = 35.71
0.01690.02600.0130
Δ
i
=
0.001
Full model = 27.73Single-averaged model = 30.0
0.03360.0311To obtain the mass of fuel required by each maneuver, as an example, a satellite with an initial mass of 1000 kg is used. The hydrazine is assumed to be the fuel for the impulsive maneuvers, withI
sp
=
340 seconds. For the low thrust maneuvers, a force of 10 newtons and a specific impulse of the 1300 s were used for the propulsive system.The durations of the propulsion arcs were not too large, going from about 0.2 minutes to nearly 21 minutes, depending on the amplitude of the maneuver. Since all the corrections are small, the consumptions are nearly proportional to the amplitudes of the maneuvers. The maneuvers using two propulsion arcs divide the fuel consumption near equally between two arcs, sometimes showing a little larger consumption in the first thrusting arc, because this arc has distances to the Earth that are a little smaller, so it spends more fuel for the maneuver. Anyway, these differences in distances from the Earth are too small and not always implied in difference in consumption between both arcs considering the number of digits shown in the tables. The maneuvers with only one arc take more time and require more fuel, since the thrust is not applied in the optimal parts of the orbit that are the regions near the apoapsis. Figure2 shows the direction of the thrust to correct the eccentricity from 0.02 to 0.01, using two arcs of propulsion, as an example of the calculations made in the present research. Note that the direction has near linear variations for both arcs, close to the direction of motion for the first arc and close to the direction opposite to the motion for the second arc. It means that the low thrust maneuvers follow the main geometry of the impulsive maneuvers, but replacing the instantaneous variation of velocity by a continuous force acting in the satellite during some time. The main advantage of this system is the fuel economy, which is about four times smaller using the propulsion system showed here, due to the higher value of the specific impulse. It can be even larger using more efficient propulsive system, which means fuel with higher specific impulses. When using two arcs of propulsion, there is a fuel economy of around 10% with respect to the use of a single arc. The reason for these extra savings is the application of the thrust in the regions where they are more efficient, that is, the apoapsis of the trajectories involved. At those points, the gravitational force of the Earth is smaller, thus allowing larger effects due to the application of the thrust.Tables3 and 4 show the equivalent results of the maneuvers used to correct the inclination. Those maneuvers require less fuel than the ones to correct eccentricities, because the limits to perform an orbital correction are very small for the inclination. The fuel consumption goes from 0.0244 to 1.2235 kg and the duration of the propulsion goes from 0.08 to 4.27 minutes. There is an almost linear relation between the desired angle change and the fuel consumption, due to the small values involved in the variations of the inclination. In this situation, only one propulsion arc was used for each maneuver, since it is the best solution, even for the impulsive maneuvers.


## 3.2. Near Equatorial Orbits

Figure3 shows the evolution of the eccentricity as a function of time for the near equatorial orbits (i
0
=
10
-
3 degrees, since real equatorial orbits show singularities in the single-averaged equations of motion). The times are shown in years. Figures 3(a), 3(b), and 3(c) show the results of the cases of the full elliptic restricted problem and single- and double-averaged models, respectively. On the top of each figure, the corresponding eccentricity of the perturbing body (e
′) is shown. This parameter is used in order to study its effects on this problem and it represents the same system used before. It is clear that none of the specified limits is reached when using the single-averaged and the double-averaged techniques (see Figures 3(b) and 3(c)). The similar results, in terms of reaching the limits to make a maneuver, which comes from both averaged models, are expected because both averaged models have results that are very close to each other when the perturbing body is in a circular orbit (Domingos et al. [17]). This effect decreases with the increase in the eccentricity of the primaries. In general, both averaged models have results that are much closer to each other than closer to those of the full model. It is visible from the full model that, during the integration time, the evolution of the eccentricity has oscillations that increase in amplitude with the increase in the eccentricity of the primaries. The averaged models are not able to predict this property. This fact changes the times when the maneuvers are required, so it has a strong influence on the estimation of the fuel consumption for the station keeping maneuvers.On the other hand, the results show that, for the full elliptic restricted three-body problem, in Figure3(a), the limits of variation of the eccentricity 0.0005, 0.001, 0.005, and 0.01 are reached, so it is necessary to perform orbital maneuvers. As an example, in Figure 3(a), we can see that when e
′
=
0.2 only the limit of variation for the eccentricity Δ
e
=
0.05 is not reached during the timescale of the numerical integrations. As a result, in this latter case, we do evidence that there are four points where it is required to perform orbital maneuvers to correct the orbit. The inclinations have zero variations during this period, so it is not required to apply any propulsion system to correct the orbit. For cases with e
′
<
0.2, the results shown in Figure 3 evidence that only two orbital maneuvers are required to correct the orbit of the satellite. So, the general conclusion is that averaged models are not good enough to predict the evolution of the eccentricity and the orbital maneuvers required by the mission for near equatorial orbits, for those values of the accuracy required in eccentricity for the mission, because there is a variation of up to 0.015 in eccentricity and an average yearly consumption of up to 0.1771 kg that are not predicted by both averaged models. This average is obtained by dividing the total fuel expenditures during the mission by the total time considered for the evolution of the orbits. So, even if the spacecraft spends more than one year without the need of a maneuver, it is still possible to compute a “yearly average consumption.”Table5 shows the times, in years, when the limits of eccentricity are reached for each situation, as well as the estimated fuel consumption per year of mission for the several limits considered for the eccentricity and for the values of the eccentricity of the primaries.The fuel consumptions are not too large, going from 0.0101 to 0.1771 kg per year. It is also noted that there is a strong dependence of the fuel consumption on the eccentricity of the primaries, in the order of five times, when comparing the circular case to the case with eccentricity 0.2, for the same limits in the eccentricity variations. In the worst case considered here, with eccentricity of primaries of 0.2, the times for performing an orbital maneuver goes from 4.76 years, if the limit used for the variation of eccentricity isΔ
e
=
0.0005, to 26.08 years, if the limit used for the variation of eccentricity is Δ
e
=
0.01.


## 3.3. Near Critical Inclination Orbits

Figure4 shows the evolution of the eccentricity and inclination, respectively, as a function of time when using the three models for the case of near critical orbits (i
0
=
39
°~0.6803 radians). In these figures, Figures 4(a), 4(b), and 4(c) also show the results of the cases of the full elliptic restricted problem and the single- and double-averaged models, respectively. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). Regarding the full elliptic restricted three-body problem, Figure 4(a), the results show that the inclination reaches the first limit of 0.0001 rad at 0.012 years, while the eccentricity reached several limits up to 0.02 in 28 years for the circular case for the perturbing body. It is also interesting to note that, when the eccentricity increases, there is a decrease in the inclination. By focusing on the value of the eccentricity of the perturbing body, we can see that the inclination remains close to the initial value all the time.Evolution of the eccentricity and inclination as a function of time for the casei
0
=
39
° (~0.6803 radians). In (a), the results are for the full elliptic restricted three-body problem. (b) is for the single-averaged model and (c) is for the double-averaged model. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). The dashed curves represent the (t, e) and (t, i) pairs for eccentricity and inclination, respectively, to perform orbital maneuvers, depending on the limit of the eccentricity and inclination variations allowed by the mission. In the plots for the inclination, it showed the dashed lines for the limits of variation for times higher than 1 year for better graphic visibility.
(a)
(b)
(c)Figure 5
Evolution of the eccentricity as a function of time for the casei
0
=
80
° (~1.3955 radians). The results are for the full elliptic restricted three-body problem. Similar results have been obtained for single- and double-averaged models. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). The dashed curves represent the (t, e) pair to perform orbital maneuvers for each limit assumed for the variation in eccentricity.It can also be seen, in Figures4(b) and 4(c), that the evolution of the eccentricity has a similar behavior for single- and double-averaged models and that they do not suffer significant changes with the increase in e
′. In the two models, the times for reaching each eccentricity limit are very close. As an example, the eccentricity 0.02 was reached in approximately 30 years in both cases. The same does not occur when the study of the inclination is made, because it reaches its limits in different times for each model.The evolutions of the orbits also show another interesting feature of the averaged models. Opposite to what happens for the eccentricity, when the inclination is studied, the averaged models predict variations that are stronger than the ones observed in the real restricted three-body problem. It means that averaged models will predict orbital maneuvers and so will fuel consumption, in situations where they are really not required. Still, regarding the orbital evolution predicted by the three models, the results shown here evidence some differences. The decrease in the eccentricity shown by the full model for the case where the eccentricity of the primaries is 0.2 is not predicted by both averaged models. The dependence of the results on the eccentricity of the primaries is much stronger when considering the elliptical restricted three-body case than when it is predicted by the averaged models, where the evolutions are similar to each other in all cases studied. Both averaged models are very similar to each other regarding predictions for the eccentricity, but they have differences from each other for the evolution of the inclination.Considering the limits of eccentricity and inclinations reached during the integration period, it is clear that it is necessary to perform orbital maneuvers in most of the cases. Table6 shows the times, in years, when these limits are reached for each situation, as well as the estimated fuel consumption per year of the mission for the several limits considered for the eccentricity and inclination and for different values of the eccentricity of the primaries.Limits up to 0.01 in eccentricity are reached for all the models. The value 0.02 was never reached. Differences in the prediction of the times for the maneuvers and the average yearly fuel consumption to correct eccentricities are of the order of 10%, so not too large. Regarding the inclination, for the circular case only, the limitΔ
i
=
0.0001 was reached by the three models but the differences in the times of the maneuver and the average yearly fuel consumption are very large among the models, with the elliptical restricted problem requiring much more fuel. When the perturbing body is not in a circular orbit, the full model takes more time to reach the limits of the eccentricity, so the averaged models predict earlier times for the maneuvers, so more yearly fuel consumption occurs when compared to the real cases. Some larger limits for the eccentricity are not reached by the full model. For the inclination, the full model predicted again higher fuel consumptions and faster orbital maneuvers than did the averaged models, so those models under evaluate the fuel consumption to correct the inclination. It is cheaper to keep lower values of the limits for the eccentricity, because the velocity of the increase in the eccentricity gets higher with the eccentricity. So, it is better to perform the maneuvers earlier, not allowing the eccentricity to grow too much. For the inclination, those costs alternate from one limit to the other.So, the general conclusion is that averaged models are also not good enough to predict the evolution of the eccentricity and inclination also in this region of near critical inclination orbits because there are considerable differences in the prediction of the full model and both averaged models.


## 3.4. High Inclined Orbits

Figures5 and 6 show the temporal evolutions of the eccentricity and inclination, respectively, for the case of high inclined orbits, considered here to have i
0
=
80
° (~1.3955 radians). The same representations for the models given in Figure 4 are used here.Evolution of the inclination as a function of time for the casei
0
=
80
° (~1.3955 radians). In (a), the results are for the full elliptic restricted three-body problem. (b) is for the single-averaged model and (c) is for the double-averaged model. The eccentricity of the perturbing body (e
′) is shown in color codes: 0.0 (red), 0.1 (green), and 0.2 (blue). The dashed curves represent the (t, i) pair to perform orbital maneuvers, for each limit assumed for the inclination.
(a)
(b)
(c)In Figure5, only the results of the evolution of the eccentricity when considering the full elliptic restricted three-body problem are showed, because similar results have been obtained for the single- and double-averaged models. This similarity is an expected result of initial inclinations above the critical value (see Domingos et al. [17]). The eccentricity oscillates with large amplitude, reaching the maximum limit of 0.06 (a variation of 0.05 from the initial value 0.01) near the time 32.5 years, while the inclination reaches only the smaller values of the limit of variation (<10−3 rad) during the integration time. We note that, for the present case of i
0
=
80
°, all the limits of eccentricity are reached when using the full problem and both averaged models, so it is necessary to perform orbital maneuvers in all the situations. Considering the inclinations, see Figure 6, the three limits of inclination are only reached for e
′
=
0.2. The single-averaged model predicts earlier maneuvers than the full model does, so it overestimates the fuel consumption. The double-averaged model follows the full model better in this topic. The eccentricity of the primaries always increases the perturbations, so it also increases the average yearly fuel consumption. Tables 5 to 9 show the times when the limits are reached for each situation, as well as the estimated fuel consumption per year of mission for the several limits considered for the eccentricity and inclination for the values of 0.0, 0.1, and 0.2 for the eccentricity of the primaries.To correct the eccentricity, the times predicted by the three models are very similar, which generates very close fuel consumption estimations. There is a strong relation between the average yearly consumption and the limits imposed on the variation in eccentricity. The cost using the limitΔ
e
=
0.05 is ten times more than the cost when using the limit Δ
e
=
0.0005. So, the more frequent the maneuvers are performed the lower the annual cost is. Again, this is a consequence of the increase in the speed of the variation of the eccentricity with respect to the initial eccentricity. When considering the corrections in the inclination, the opposite occurs, so allowing larger limits reduces the total cost of the maneuvers per year of mission.In this situation, the corrections in eccentricity are made in times that still have similar predictions by the three models, but differences around 10% start to occur in some cases, generating fuel consumption estimates with differences of the same order. The relation between the yearly consumption and the limits imposed on the variation in eccentricity is kept in similar values, being more economical by an order of ten when using the strategy of making frequent maneuvers. The same previous pattern observed for the circular case occurs for the corrections in the inclination, and allowing larger limits still reduces the total cost of the maneuvers per year of mission.The same analysis made before applies here. The only difference is that there is an increase in the differences of the times to perform maneuvers and in the fuel consumptions, as a consequence, between the full and the averaged models, going to differences close to 15% now. It is a consequence of the fact that the averaged models loss quality when the eccentricity of the primaries increases, as shown in detail by Domingos et al. [17].


## 4. Conclusions

The present paper studied the effects of the eccentricity of the orbit of a third-body that is disturbing a spacecraft orbiting the same central body. A hypothetical system considered here has the same masses and distances of the Earth-Moon system but has varied values for the eccentricity, in order to study its influence in the station keeping consumption. The idea is studying potential systems of star-planet or planet-moon that exist outside the solar system in terms of learning the necessity and the cost of keeping a permanent basis orbiting those systems. In the solar system, the dwarf planet Haumea is also a good candidate to use the ideas shown here. Another contribution of this paper is to study the differences of the single- and double-averaged models, both compared with the full elliptical restricted problem, in the times to require an orbital maneuver as well in the cost required. Emphasis is given on the estimation of the fuel consumption for station keeping maneuvers. Impulsive, with one or two impulses, and low thrust maneuvers, with one or two propulsion arcs, were used for a satellite with an initial mass of 1000 kg.The main advantage of the low thrust system is the fuel economy, which is near four times smaller for the propulsion systems used here, due to the higher specific impulse. When using two propulsion arcs, there is an extra saving near 10%.For the near equatorial orbits, the averaged models are not good in predicting the evolution of the eccentricity and the orbital maneuvers. There is a variation of up to 0.015 in the eccentricity and a yearly consumption of up to 0.1771 kg that occur using the full model but are not predicted by the averaged models. It is also shown that there is a strong dependence of the results on the eccentricity of the primaries, in the order of five, in terms of the fuel consumption required by the circular case when compared to the case where the eccentricity is 0.2. This estimation was one of the goals of the paper.For the near critical inclination orbits, an interesting feature of the averaged models is found. Opposite to what happens for the eccentricity, the averaged models predict variations that are stronger than the ones observed in the real restricted three-body problem. The differences in the prediction of the times for the maneuvers and the yearly fuel consumption to correct eccentricities are in the order of 10%.For the high inclined orbits, the times predicted by the three models are very similar when studying the eccentricity, which generates very close fuel consumption estimates. So, the averaged models are good to make predictions in this range of orbits.In general, it is noticed that there is a strong relation between the yearly consumption and the limits imposed on the variation in eccentricity. The cost using the limitΔ
e
=
0.05 is ten times more than the cost when using the limit Δ
e
=
0.0005, which indicates that the use of more frequent maneuvers reduces the annual cost of station keeping. For the corrections in inclination, the opposite occurs.



---

*Source: 359845-2014-07-03.xml*