Controlling the Eccentricity of Polar Lunar Orbits with Low-Thrust Propulsion

1 Grupo de Dinâmica Orbital & Planetologia, UNESP, Universidade Estadual Paulista, Avenida Ariberto Pereira da Cunha, 333, CEP, 12516-410, Guaratinguetá, SP, Brazil 2 Instituto Nacional de Pesquisas Espaciais—INPE, Avenida dos Astronautas 1752, CEP, 12227-010, São José dos Campos-SP, Brazil 3 Universidade Federal do ABC, UFAB, Santo André, 09210-170, SP, Brazil 4 Instituto de Fı́sica, Universidade de Brası́lia—UnB, CP 04455, CEP, 70919-970, Brası́lia, DF, Brazil


Introduction
Recently, several nations presented plans to reach the Moon.Satellites have been launched and many more are planned for following years see, e.g., 1 .The expectations are that in the near future there will be a lunar base.The lunar poles are particularly of interest since it seems to be where water can be found.Therefore, long living satellites in polar lunar orbits will be needed.It is well known that lunar satellites in polar orbits suffer a strong gravitational perturbation from the Earth.That effect is a natural consequence of the Lidov-Kozai resonance.
It is well known that the Lidov-Kozai resonance introduces equilibrium configurations.In the case of lunar polar orbits disturbed by the Earth's gravitational field, this can be used as an advantage to implement constellations of satellites with elliptic highly inclined orbits 2, 3 .On the other hand it causes instability for near circular highly inclined orbits.Wytrzyszczak et al. 4 studied the regular and chaotic motion of geosynchronous satellites disturbed by the Moon's gravitational field.They found that the chaotic nature of high inclination satellites is caused due to the significant eccentricity growth caused by the Lidov-Kozai resonance.
Similarly, the final fate of polar lunar near circular satellites is the collision with the Moon.Therefore, the control of the orbital eccentricity leads to the control of the satellite's lifetime.
In this paper we propose the control of the eccentricity using an electrical thruster, similar to the one that is in development at the University of Brasília.Electric propulsion is basically a technique of space propulsion which involves the conversion of electrical power into the kinetic power or thrust of the exhaust beam of ionized particles.The ability to obtain high exhaust velocities with ionized particles enables plasma thrusters to perform high specific impulse mission in space 5 .The main innovation of the thruster that is being developed at the University of Brasília is the use of a permanent magnet, saving energy during the mission.Preliminary results in the laboratory show that it is possible to obtain more than 100 mN with this technology.Inspired on this thruster project, in this work, we assume a constant exhaust velocity, and we can control the switch of the thruster during the mission.
In the present work we introduce an approach in order to keep the orbital eccentricity of lunar polar satellites at low values.The approach is based on the use of low-thrust propulsion in order to introduce a correction of the eccentricity.
In the next section we introduce the Lidov-Kozai resonance.In Section 3 we show the evolution of the eccentricities form our numerical integrations.The approach proposed to control the eccentricity and its application is presented in Section 3. In the final section we present our final comments.

The Lidov-Kozai Resonance
Lidov 6 , studying the dynamics of artificial satellites, and Kozai 7 , studying the dynamics of asteroids, independently discovered what is now called the Lidov-Kozai resonance.Following, we introduce the basic features of such resonance.
In this section we adopted a simple model see, e.g., 2 for the orbital evolution of an artificial satellite disturbed by a third body in circular equatorial orbit around the primary.It was obtained by double averaging the system 8 taking into account the disturbing function expanded in Legendre polynomials up to second order and the eccentricity of the disturbing body also up to the second order.The disturbing function of the problem is averaged independently over the mean longitudes of the satellite and the third body.The standard definition for average used in this work is where M is the mean anomaly, which is proportional to the time.

Mathematical Problems in Engineering 3
Following such approach one can find the double averaged disturbing function given by see, e.g., 9 where a, e, ω, and i are, respectively, the semimajor axis, eccentricity, argument of pericenter and inclination, G is the gravitational constant, a E is the semimajor axis of the Earth with respect to the Moon, and m 1 and m 2 are the masses of the Earth and Moon respectively.Substituting R in Lagrange's planetary equations see, e.g., 10 , we find where n is the mean motion and γ m 1 / m 1 m 12 .
Considering the case when de/dt 0 and dω/dt 0 one can find three first integrals: where a, e, i, and w are the semi-major axis, eccentricity, inclination, and argument of pericenter of the satellite, and a 0 , k 1 , and k 2 are the constants of motion.This system has a set of fixed points given by Therefore, for a system with e, i, and ωsatisfying conditions 2.10 , the satellite would be in what can be called a frozen orbit; that is, apart from short-period oscillations, the orbit would be kept fixed in size and location.
A simple analysis of 2.2 and 2.3 shows that 9 i for k 2 > 0 and any value of k 1 , w circulates, ii for k 2 < 0 and k 1 < 3/5, w librates around 90 So, that is the Lidov-Kozai resonance.When i > i * , the system behaves like a pendulum, with stable fixed points, librations around such points and circulation.
Following another approach let us consider a satellite orbiting the Moon and disturbed by the Earth in an elliptical orbit with respect to the Moon.Taking into account only the term P 2 of the Legendre polynomial, the disturbing potential is given by where e E and n E are the eccentricity, and mean motion of the Earth with respect to the Moon, and μ E m E / m E m M , where m E and m M are the masses of the Earth and Moon, respectively.
One can identify the Lidov-Kozai resonant features by numerically integrating Lagrange's planetary equations for the temporal variation of the argument of pericentre, w, and the eccentricity, e, with the disturbing potential given by 2.10 .In Figure 1 we present a sample of these numerical integrations in a diagram e versus w.There are two sets of initial values of eccentricity and inclination e o , i o .One for low inclination i o 20 • and other for high inclination i o 56.2 • .This figure shows a clear dependence of the eccentricity on the argument of pericentre for an orbit with high inclination.All satellites inclined to the orbital plane of the third body the Earth by more than 39.2 • , the critical inclination, experience a considerable growth of eccentricity.The Earth causes the Lidov-Kozai resonance driving the eccentricity growth.For all trajectories from the set with initial inclination higher than the critical value, the argument of pericentre librates, while it circulates for trajectories from the set with initial inclination lower than the critical value.

Eccentricity Growth
In this section we present numerical integrations considering two dynamical systems: the 3body problem, Moon-Earth-satellite, and the 4-body problem, Moon-Earth-Sun-satellite.In all simulations the satellite is initially in polar orbit i 90 • .
In the case of the 3-body problem, considering a coordinate system centered in the barycenter of the Earth-Moon system X, Y, Z , the equations of motion of the satellite are given by where m is mass and the index i 1 refers to the Earth and i 2 refers to the Moon.In this system, the equations of motion for the moon and the Earth are given by

3.2
First, we simulated the system, integrating numerically 3.1 and 3.2 , considering a satellite with initial eccentricity equals to 0.0001 and a range of initial altitudes between 100 km and 5000 km. Figure 2 shows the evolution of the satellite's eccentricity for the 3-body simulations, considering altitudes h 100, 200, 500, 1000, and 5000 km.The plots show an exponential evolution of the eccentricity.We computed the time needed in order to reach the eccentricity that corresponds to the collision of the satellite with the Moon.A fit of the collision time, T collision , in Earth days, as a function of the altitude, h, is given by the expression:

3.3
The same set of simulations was performed considering the 4-body problem, adding the perturbations of the Sun.Now we considered a new coordinate system X Y Z , centered at the barycenter of Sun-Earth-Moon system.In this case we integrate numerically the equations similar using prime in the variables to equations 3.1 and 3.2 , but we add the index i 3, where the fourth term refers to Sun.However, the results found for the 4-body model were not significantly different from those found for the 3-body model.The empirical expression for the length of time needed to occur the collision with the Moon as a function of the initial altitude is given by T collision 2494e 0.063h−3.94×10−4 h 2 . 3.4 A comparison of the two sets of simulations and 2.5 and 2.6 is shown in Figure 3.

Controlling the Eccentricity
In order to control the satellite's eccentricity we will use low-thrust propulsion.Following the work of Sukhanov 11 , we use the locally optimal thrust for each orbital element.This development is based on the performance index, through the minimization of a functional in the direction of the orbital element to be changed.In our case, the eccentricity is the parameter to be minimized.The result is a vector, called Lawden's primer vector, P, which gives the direction of the thruster to be turned on.The eccentricity of the satellite relative to the Moon is given by where c is the magnitude of the angular momentum, and h is the integral energy.The primer vector is given by where p a c 2 /μ, p r , and p n are the radial and the tangential components, respectively.The approach we are proposing is a very simple one.The idea is to introduce a correction on the eccentricity every time it reaches a certain limit.The procedure is as follows.Following, we present the results of some simulations assuming e o 0.04 and Δe 0.01.These simulations were made considering a set of different thrust values, from 0.1 N up to 0.4 N. In each run we measured the length of time T Thruster , needed to correct the eccentricity value from e 0.04 to e 0.05 .From these results we obtained empirical expressions of T Thruster as a function of the initial altitude and as a function of the thrust value.As an example, in Figure 4 is shown the temporal evolution of the eccentricity and of the orbital radius for a satellite with an initial altitude of 500 km and using a thruster of 0.2 N.
In Figure 5 we present the propellant consumption per year of lifetime for the whole set of simulations, that is, different initial altitudes and different values of the thruster.The time intervals that the thrusters are turned on and off are shown in Figures 6 and 7.

Final Comments
In the present work we have studied the problem of polar lunar satellites in near circular orbits under the gravitational perturbations of the Earth and the Sun.The problem is dominated by the Lidov-Kozai resonance, which forces satellites with near circular orbits to have an exponential growth of its eccentricity.In order to keep the satellite with low eccentricity we proposed to use low-thrust propulsion every time the eccentricity reaches a limiting eccentricity, acceptable by the requirements of the satellite mission.The results show   that the satellite's lifetime can be reasonably extended several years at a not so expensive cost.Therefore, it is shown that low-thrust propulsion is very adequate for this kind of purpose.

Figure 2 :
Figure 2: Time evolution of the eccentricity for the 3-body problem.The colour code indicates the initial altitude in kilometers.The time is in Earth days.

Figure 3 :Figure 4 :
Figure 3: Collision time as a function of the initial altitude.The red crosses are for the 3-body simulations and green are for the 4-body simulations.The blue curve corresponds to 2.5 and the purple curve corresponds to 2.6 .

Figure 5 :
Figure 5: The propellant consumption per year of lifetime for the whole set of simulations, that is, different initial altitudes and different values of the thruster.

Figure 6 :
Figure 6: Time interval that the thrusters are turned on and off, for the whole set of simulations.
There are two sets of initial values of eccentricity and inclination e o , i o : One for low inclination, i o 20 • thick lines and other for high inclination, i o 56.2 • , thin lines .The values of the initial eccentricities are represented by the colour code in the bottom of the figure.The argument of pericentre circulates in the case of initial low inclination, while librates in the case of initial high inclination.
• or 270 • , iii for k 2 0 and w 90 • or 270 • , i i * ∼ 39.2 • , where i * 39.2 • is the critical value, which corresponds to the frozen orbit.Figure 1: Sample of the satellite's orbital evolution in a diagram e versus w.