Javier Chivite Sierra
Space Engineering, master's level (120 credits)
Lule˚ a tekniska universitet
Internship / Master Thesis
MSc in Space Technique & Instrumentation MSc in Space Science and Technology
Multiple fly-by for interplanetary missions
Author:
Javier Chivite Sierra
Supervisors:
Roberto Armellin
Francesco Sanfedino
2 R´esum´e 3
3 Introduction 4
4 Fly-by 8
5 Mission Development 12
5.1 Launch and Arrival to the Moon . . . 12
5.2 Fly-By . . . 15
5.3 Second Fly-By . . . 20
5.4 Third Fly-by and trajectory to reach the target . . . 23
6 Optimisation 26
7 Future Work 28
8 Conclusions 29
References 30
1 Abstract
Current state-of-the-art of propulsion system for space vehicles does not allow to deliver the required payload for the mission to all bodies in the Solar System. Therefore, alternatives have been developed to reach those bodies without having the necessary technology. Gravity Assist Manoeuvres take advantage of the encounter of the spacecraft with one or more celestial bodies to modify the velocity vector of the spacecraft. These manoeuvres have already been used previously to reach high ∆v targets with a very low propellant consumption.
This thesis models a Gravity Assist Manoeuvre to later apply the model to a space mission to reach a target with multiple gravity assist manoeuvres around the Moon to reduce the fuel consumption.
In the first part, the gravity assist manoeuvre is designed based on the angle that the normal vector of the plane of the fly-by forms with the perpendicular vector to the velocity of the spacecraft on the Moon reference frame and the velocity of the Moon. The second designed parameter for the fly-by is the angle that the velocity of the spacecraft relative to the Moon rotates about the plane of the fly-by modifying the direction of the spacecraft’s velocity.
The second part of the project applies the previous concept of gravity assist manoeuvres to a space mission. A spacecraft orbiting on a Geostationary Transfer Orbit is injected into an orbit to the Moon. Once the spacecraft reaches the Moon, it flies by the Moon modifying the direction and magnitude of the velocity of the spacecraft in the Earth reference frame. The orbit obtained after the fly-by is then propagated for a given period of time before injecting the spacecraft again into an orbit to the Moon. After arrival to the Moon, the direction and magnitude of the velocity of the spacecraft in the Earth reference frame is modified through a second fly-by. Afterwards, the previous process is repeated again for a third fly-by before transforming the velocity of the spacecraft into the Heliocentric reference frame. The orbit is propagated for a period of time before getting injected into an orbit to the target of the mission.
The third part of the project aims to optimise the mission developed in the second part, though only two fly-bys will be considered in this part to simplify the optimisation process. The previous mission has been developed in several sections and the minimum fuel consumption has been determined individually for each section, obtaining a local minimum. Unfortunately, the global minimum fuel consumption determined as the addition of those local minimum fuel consumption might be different, i.e. different sections might influence other sections leading to a lower global minimum fuel consumption that has not been considered before.
The fourth and last part shows future work that might be done to the project. It includes the modification and application of the code for the optimisation, the study of powered fly-bys to modify the previous parts, and the addition of perturbations and space interactions to develop a more realistic mission.
2 R´ esum´ e
L’´etat de l’art actuel des syst`emes de propulsion pour v´ehicules spatiaux ne permet pas d’envoyer les diff´erentes charges utiles vers tous les objets qui composent le syst`eme solaire. Des alternatives ont ´et´e d´evelopp´ees afin d’atteindre ces objets sans la technologie n´ecessaire. Les manœuvres d’assistance par gravit´e prennent avantage de la rencontre entre le v´ehicule spatial avec un ou plusieurs objets c´elestes pour modifier la v´elocit´e du v´ehicule spatial. Ces manœuvres ont d´ej`a ´et´e utilis´ees auparavant afin d’atteindre des ∆V tr`es ´elev´es pour une consommation de combustible tr`es faible.
Ce rapport mod`ele une manœuvre d’assistance par gravit´e pour plus tard l’appliqu´e `a une mission spatiale qui a pour but d’atteindre une cible `a l’aide de plusieurs manœuvres d’assistance par gravit´e autour de la Lune afin de r´eduire la consommation de propergol. Dans la premi`ere partie, la manœuvre d’assistance gravitationnelle est con¸cue en fonction de l’angle que forme le vecteur normal du plan du survol avec le vecteur perpendiculaire `a la vitesse de l’engin spatial sur le cadre de r´ef´erence de la Lune et la vitesse de la Lune. Le deuxi`eme param`etre con¸cu pour le survol est l’angle de rotation de la vitesse de l’engin spatial par rapport `a la Lune autour du plan du survol en modifiant sa direction.
La deuxi`eme partie du projet applique le dit concept des manœuvres d’assistance par gravit´e `a une mission spatiale. Un engin spatial en orbite sur une orbite de transfert g´eostationnaire est inject´e sur une orbite de la Lune. Une fois que l’engin spatial atteint la Lune, il vole vers la Lune en modifiant la direction et la magnitude de sa vitesse en utilisant la Terre comme r´ef´erence.
L’orbite obtenue apr`es le survol est alors propag´ee pendant une dur´ee donn´ee avant de r´einjecter l’engin spatial sur une orbite vers la Lune. Apr`es l’arriv´ee sur la Lune, la direction et la magnitude de la vitesse de l’engin spatial dans la r´ef´erence terrestre sont modifi´ees par un deuxi`eme survol.
Ensuite, le processus pr´ec´edent est r´ep´et´e pour un troisi`eme survol avant de modifier la vitesse de l’engin spatial en utilisant le rep`ere h´eliocentrique. L’orbite se propage pendant un certain temps avant d’ˆetre inject´ee sur une orbite vers la cible de la mission.
La troisi`eme partie a pour objectif d’optimiser la mission d´evelopp´ee dans la deuxi`eme partie, mais seulement deux survols seront pris en consid´eration dans cette partie pour simplifier le processus d’optimisation. La mission pr´ec´edente a ´et´e d´evelopp´ee en plusieurs sections et la consommation minimale de carburant a ´et´e d´etermin´ee individuellement pour chaque section, en obtenant un minimum local. Malheureusement, la consommation globale d´etermin´ee comme l’ajout de cette consommation de carburant minimale locale pourrait ˆetre diff´erente, c’est-`a-dire que diff´erentes sections pourraient influencer d’autres sections, conduisant `a une consommation de carburant minimale globale plus faible qui n’avait pas ´et´e prise en compte auparavant.
La quatri`eme et derni`ere partie montre le travail futur qui pourrait ˆetre r´ealis´e. Il comprend la modification et l’application du code pour l’optimisation, l’´etude des survols motoris´es pour mod- ifier les parties pr´ec´edentes, et l’ajout de perturbations et d’interactions spatiales pour d´evelopper
3 Introduction
Space exploration has always fascinated mankind but, it was not until the first half of the 20th century that many studies were carried out to develop solid and liquid-propellant rockets. In 1903, Konstantin Tsiolkovski presented his work ”In Exploring Space with Reactive Devices” with severa mathematical calculations and design concepts that laid out the foundations for liquid- fueled rocketry [1]. A few decades later, the German V-2 missile, powered by a liquid-propellant rocket engine and developed as a Vengeance Weapon against the Allies, became the first human artefact to achieve sub-orbital spaceflight, crossing the Karman line on June 20, 1944 [2].
Over the next decades, the space race between USA and USSR made a breakthrough in the history of space exploration. On April 12, 1961, Soviet Yuri Gagarin reached the outer space and orbited the Earth [3], and eight years later, on July 21, 1969, American Neil Armstrong and Edwin E.
Aldrin set foot on the Moon on a mission aboard Apollo 11 [4].
Since then, several missions for the exploration of space and universe have been carried out.
However, propulsion current state-of-the-art is not enough to deliver the required payload of the mission to all bodies within the Solar System. Reaching Mercury, Jupiter or the outer part of the Solar System is not possible just using current propulsion capabilities. Thus, gravity assist manoeuvres have been developed to gain the required change in velocity to reach those bodies without fuel consumption.
A gravity assist manoeuvre (GAM), also called swing-by, fly-by or gravitational slingshot, is the use of the relative movement and gravity of a planet or other massive celestial body to modify the velocity of a spacecraft. As the spacecraft gets close to the celestial body, the gravity of the celestial body produces a change in the velocity vector of the spacecraft [5]. This concept has been extensively studied and applied over the last fifty years and, since 1970, several fly-bys and orbital flights have been conducted around planets such as Venus, Jupiter, Mercury, Saturn and Mars, as well as different small bodies such as comets.
Mariner 10 probe was the first spacecraft that used a fly by to reach another planet. In this case, the spacecraft swung by Venus on February 5, 1974 on its way to Mercury, becoming the first spacecraft that explored this planet. In addition, one of the main goals of the mission was to prove the possibility of using gravity assist manoeuvres [6].
Pioneer 11 was the second mission of the Pioneer program to study Jupiter and the outer solar system, but the first one that investigated Saturn and its main rings. Pioneer 11 probe did a fly by around Jupiter in order to obtain the necessary energy to reach Saturn. Subsequently, the spacecraft swung by Saturn following an escape trajectory from the Solar system [7].
Figure 1: Trajectories of Pioneer 10 and 11 [7]
NASA continued its interplanetary exploration program with the Voyager 1 and Voyager 2 probes.
Even though these spacecraft were going to investigate Jupiter and Saturn, both were able to continue their mission into the outer Solar system. Both probes were injected into a direct transfer to Jupiter in 1977. Voyager 2 encountered Jupiter, Saturn, Uranus and Neptune, whereas Voyager 1 was launched on a faster trajectory that enabled the probe to reach Jupiter and Saturn sooner, avoiding flying by Uranus and Neptune. Both probes have currently exited the Solar system. A particular fact to highlight is that both missions encountered a favourable alignment of the outer planets (Jupiter, Saturn, Uranus and Neptune), which happened from 1970 to 1990 [8].
Figure 2: Trajectories of Voyager 1 and Voyager 2 (credit: NASA)
Galileo spacecraft aimed to reach Jupiter injected into a direct transfer. Due to some modifications in the safety protocol as a result of the Challenger accident in 1986, the trajectory of the probe was
complex MGA trajectory, including two swing-bys around Venus, a swing-by around the Earth, and one around Jupiter. Once the Cassini probe was orbiting around Saturn, it carried out several swing-bys around Saturn’s moons to modify its orbit in order to study the Saturn system [10].
NASA launched another mission to Mercury in 2004, the Messenger mission. In order to reach Mercury, several swing-bys were carried out. The target would be reached after a swing-by around the Earth and two additional swing-bys around Venus. Messenger’s performed three swing-bys around Mercury to decelerate its velocity relative to Mercury and ease the orbital insertion ma- noeuvre. The gravity assist manoeuvres were used to overcome the problem of massive acceleration that accompanies flight toward the Sun, changing the in-plane velocity components of Messenger as well as the orbital inclination. In this case, the inclination of Mercury is about 7 degrees over the ecliptic.
Figure 3: MESSENGER trajectory
The Rosetta mission developed by the European Space Agency (ESA) was launched in 2004 and aimed to rendezvous the comet 67P/Churyumov–Gerasimenko and release the lander, Philae, on it. The spacecraft went through several swing-bys to perform the rendezvous with minimum propellant mass. In order to reach the comet, the spacecraft carried out a fly-by around the Earth, a fly-by around Mars and two more swing-bys around the Earth. The rendezvous with the comet took place on August 6, 2014, reducing Rosetta’s relative velocity. Rosetta entered the orbit on September 10 at an altitude of 30 km from the nucleus [11].
A mission to highlight is the New Horizons, developed by NASA and launched in 2006. This mission aimed to study the dwarf planet Pluto, its moons and other objects contained in the Kuiper Belt. In this case, the spacecraft flew by Jupiter in order to increase its velocity by about 14.000 km/h, shortening its trip to Pluto by three years [12].
Figure 4: New Horizon trajectory [13]
BepiColombo mission developed by ESA in partnership with Japanese Aerospace Exploration Agency (JAXA) aims to investigate Mercury. In order to reach Mercury, the spacecraft launched in 2018 has performed a swing-by around the Earth on April 10, 2020, a swing-by around Venus on October 15, 2020, and is expected to do another fly-by around Venus in 2021 and six swing-bys around Mercury before getting inserted into Mercury orbit in 2025.
Previous missions shown above have demonstrated the feasibility of gravity assist manoeuvres for space missions. The goal of the current project is to propose low-energy interplanetary trajectories using multiple fly-bys around the Moon for a scientific observation mission of a celestial body such as an asteroid.
4 Fly-by
As stated in ”Section 3”, a fly-by or gravity assist manoeuvre modifies the velocity of the spacecraft when the space vehicle swings by close to a celestial body, so its gravity changes the velocity vector of the spacecraft. In this passage, there is a momentum exchange between the spacecraft and the celestial body. The spacecraft can either increase, decrease or simply rotate its inertial velocity whereas the celestial body will barely lose orbital momentum due to its massive mass compared to the one of the spacecraft.
Figure 5: Gravity assist trajectories. Decreasing spacecraft’s energy relative to the main body (Left) Increasing spacecraft’s energy relative to the main body (Right) [14]
To model the Fly-by, the procedure stated at [11] has been followed. It will be assumed that the spacecraft is located at the same position as the celestial body, in the case of the current project, the moon. The gravitational parameter of the moon and its mean radius will be given respectively as µM and RM. The velocity of the spacecraft on the arrival to the moon will be v− and the velocity of the celestial body, i.e. the Moon, will be given as vm. Both velocities will be considered to be in the Geocentric Reference Frame. In the Moon reference system (parallel to the Geocentric Reference Frame), the spacecraft approaching the Moon from a great distance (r−→ ∞) with a velocity relative to the Moon will be calculated as:
v−∞= v−− vm (1)
As the spacecraft approaches from infinity, the orbit of the spacecraft about the moon is hyperbolic.
The objective is to find out the outgoing velocity of the spacecraft at the end of the fly by, v+∞. In the case of a hyperbolic orbit, the eccentricity will be greater than one (e>1), so the semimajor axis ”a” is negative. The Vis-viva (energy) equation can be expressed as:
E=v 2−µm
r = (v∞−)2
2 =(v+∞)2
2 = −µm 2a
From the previous equation, it can be noted that for a hyperbolic orbit, the magnitude of the incoming and outgoing velocity of the legs of the hyperbola must be the same.
v∞− = v∞+ = v∞
Thus, the magnitude of the velocity relative to the moon, before and after the swing-by, remains constant. The fly-by will only change the direction of the relative velocity vector, but not its magnitude.
During the fly-by, the spacecraft follows a hyperbola, so the fly-by is planar. The hyperbola lies on a plane that contains the incoming relative velocity vector (v−∞) and the centre of mass of the moon. To determine the plane attitude, an additional angle ”γ” is added as a parameter of the fly-by. It is used to calculate the point where the spacecraft pierces the sphere of influence of the moon. Otherwise, infinite different planes could be considered.
Figure 6: (a) Hyperbola plane, containing the incoming velocity of the leg of the hyperbola, the body and angle γ of rotation about v−∞ with respect to a reference. (b) Geometry to compute γ (credit: [11])
In order to model the fly-by, the first step is to determine the plane of the hyperbola. In this case, the normal vector of the plane taken as reference (shown in Figure 6(a)) will be assessed as the cross product of the incoming velocity of the leg of the hyperbola and the velocity of the moon.
The Reference will be considered an auxiliary plane. Therefore, the auxiliary reference system is determined as:
n1= v+∞
|v+∞| (2)
n3= v+∞×vm
|v+∞×vm| (3)
n2= n3× n1 (4)
where n3is the third component of the auxiliary reference frame, and also the reference to calculate the required plane for the fly-by.
In order to get the normal vector of the plane of the fly-by, a rotation angle ”γ” of n3 about the n1 vector (direction of the incoming velocity of the leg of the hyperbola (v−∞)) is required.
Figure 7: Rotation about the unit vector with the direction of the incoming velocity of the leg of the hyperbola
The new reference system is given through the rotation matrix:
[R(γ)] =
1 0 0
0 cos(γ) sin(γ) 0 −sin(γ) cos(γ)
Thus, the new reference system can be expressed as:
u1 u2
u3
=
1 0 0
0 cos(γ) sin(γ) 0 −sin(γ) cos(γ)
n1 n2
n3
= [UN] (5)
v+∞ is found by rotating v−∞ an angle φ about the axis u3, also known as nΠ in Figure 6, but first v−∞ has to be expressed in the new reference system. Denoting the coordinates in the new reference system as ”2”, v−∞ is expressed as:
v−∞2= [UN]v−∞ (6)
The rotation of v∞− about the axis u3 is determined through the rotation matrix:
[R(φ)] =
cos(φ) sin(φ) 0
−sin(φ) cos(φ) 0
0 0 1
Therefore, v+∞ can be calculated as:
v+∞2= [R(φ)]v−∞2 (7)
The value of v+∞ in the Moon reference system can be assessed through the inverse of the [UN]
matrix. As [UN] matrix is an orthogonal matrix, in this case:
[UN]−1= [UN]T
Thus, v+∞ in the Moon reference system is evaluated as:
v+∞= [UN]−1v+∞2= [UN]Tv+∞2 (8) The outgoing velocity of the spacecraft v+ in the Earth reference system is calculated as the addition of the velocity of the moon plus the outgoing velocity relative to the moon (v+∞).
v+= v+∞+ vm (9)
The increment of velocity of the spacecraft in the Earth reference system due to the fly by is given by:
∆V= v+− v− (10)
In addition, from the previous value of the angle φ, the eccentricity of the hyperbolic orbit of the fly-by is determined as:
sin(φ 2) = 1
e −→ e= 1
sin(φ2) (11)
And the radius of the perigee of the hyperbolic orbit can be obtained as:
e= 1 +rpV∞2
µM −→ rp= (e − 1)µM
V∞2 (12)
5 Mission Development
The aim of the mission is to reach a target, an asteroid, with the minimum fuel possible, saving weight and reducing the cost of the mission. This is achieved through several fly-bys around the Moon. The procedure followed will be explained in the following subsections.
All code written and provided for the current project has been developed using MATLAB software.
The ephemeris of the Moon with respect to the Earth has been obtained using CSPICE toolkit developed by NASA for different programming languages and software.
5.1 Launch and Arrival to the Moon
The vehicle proposed for the mission would be the ARIANE 6 launched from Kourou, which is located in the French Guiana (longitude: −52.77◦), to a Geostationary Transfer Orbit (GTO) with the following orbit parameters:
Parameters of the GTO Orbit
a 24467 km
e 0.73
i 5.13◦
ω 178◦
θ 0◦
where ”a” is the semimajor axis, ”e” is the eccentricity, ”i” is the inclination, ”ω” is the argument of the perigee and ”θ” is the true anomaly.
The Right Ascension of the Ascending Node (RAAN) is missing in the previous table. That is because this orbit parameter depends on the insertion time into orbit. Therefore, the local sidereal time (degrees) for Kourou’s longitude has been determined adding the longitude of Kourou to the Greenwich sidereal time (degrees) at the specified UT, which depends on the Julian Date. The Julian Date (JD) for periods between 1900 and 2100 is obtained as:
J0 = 367(year) − INT (7 · year + INT (month+912 )
4 ) + INT (275 · month
9 ) + day + 1721013.5
J D= J0 +
5
60+minute
60 + hour
24 (13)
The number of centuries since the epoch J2000 can be expressed as:
j= jd0 − 2451545
36525 (14)
where ”jd0” depends on the difference between the Julian Date (JD) and the floor value of the Julian Date. The Greenwich sidereal time (degrees) at 0 hr UT (g0) can be determined as:
g0 = 100.4606184 + 36000.77004 · j + 0.000387933 · j2−2.583 · 10−8· j3; (15) The Greenwich sidereal time ”gst” (degrees) at the specified UT is given by:
gst= g0 + 360.98564724 · (JD − jd0); (16)
And the Local Sidereal time can be calculated adding the longitude of Kourou to the Greenwich sidereal time:
LST = gst + (Longitudeof Kourou) (17)
which must be reduced to a range between 0◦ and 360◦ to obtain the value of the Right Ascension of the Ascending Node (RAAN) of the orbit.
Figure 8: GTO Orbit
For the previous orbit parameters, the state vector of the spacecraft (position and velocity) can be assessed first in the Perifocal coordinates and transformed into the Geocentric Equatorial coordinates afterwards. First, the Specific angular momentum can be obtained as:
h=qµE· a ·(1 − e2) (18)
In the perifocal frame, the position vector can be written as:
r= xp + yq
where x=r · cos(θ) and y=r · sin(θ), and r, is the magnitude of r, given by the orbit equation:
r= h2 µE
1 1 + e · cos(θ)
Therefore, the position of the spacecraft in the perifocal frame is calculated as:
h2 1
Figure 9: Perifocal frame pqw (credit: Curtis)
In order to transform the coordinates from perifocal to geocentric equatorial frame, a 3-1-3 rotation is carried out. The rotation matrix [Q]=[R3(ω)][R1(i)][R3(Ω)]
[Q] =
−sin(Ω)cos(i)sin(ω) + cos(Ω)cos(ω) −sin(Ω)cos(i)cos(ω) − cos(Ω)sin(ω) sin(Ω)sin(i) cos(Ω)cos(i)sin(ω) + sin(Ω)cos(ω) cos(Ω)cos(i)cos(ω) − sin(Ω)sin(ω) −cos(Ω)sin(i)
sin(i)sin(ω) sin(i)cos(ω) cos(i)
where Ω is the Right Ascension of the Ascending Node, ”i” is the inclination and ω is the argument(21)
of the perigee. The position and velocity in the geocentric equatorial frame can be expressed as:
rGEF= [Q]rPF (22)
vGEF= [Q]vPF (23)
After the initial position and the velocity are calculated in the Geocentric Equatorial frame, the state vector has been propagated for an increment of time ∆t=10.000 seconds before injecting the spacecraft on an orbit to the Moon. This is done through the Kepler’s problem, using Lagrange coefficients. Given the initial position and velocity, the time of flight between the two states and the gravitational parameter, the position and velocity after such time of flight can be determined.
For the current project, the code written by R. Armellin, P. Di Lizia, and F. Topputo for the Kepler’s problem was used. The position and velocity of the spacecraft after propagating the orbit are r1 and v1 respectively.
For the arrival to the Moon, the position and velocity of the Moon are computed using the CSPICE toolkit. The arrival date time and the observer, which is the reference system for the position and velocity of the Moon, i.e., the Earth, are the inputs.
In order to link the spacecraft that has been propagated on the orbit throughout a ∆t time, and the Moon, the Lambert arc is used. The Lambert Arc is a coast arc that is calculated solving the Lambert’s Problem [15]. The Lambert’s problem is the problem to find out the orbital parameters of the conic arc connecting two points r1 and r2 in space in a given time ∆t and number of full revolutions. After determining the orbital parameters, it is calculated the velocity vectors v1 and v2 at the boundaries of the arc. The code used in the project to solve the Lambert’s problem has been developed by Izzo, and Lancaster, Blanchard & Gooding.
The necessary increment of velocity for the Lambert arc is calculated as the difference between the velocity at the boundary of the GTO for the Lambert arc (vL1) and the velocity of the spacecraft
propagated on the GTO orbit (v1).
∆VL1= vL1− v1 (24)
The previous increment of velocity ∆VL1 has been calculated for a specific departure and arrival date time. However, changing the departure and arrival dates, it would result in a new trajectory with new values for the orbit parameters and the velocities at the boundaries of the arc. The plot of ∆VL1 depending on the departure and arrival date times is known as Porkchop Plot. In the project, it was created an array with different departure times and another array with different arrival times. The Lambert’s problem was solved for each trajectory and the ∆VL1 determined.
The results are shown in the figure below:
DELTA-V CONTOUR (km/s)
11-May-2020 00:00 12-May-2020 00:00
DEPARTURE DATE 15-May-2020
16-May-2020 17-May-2020 18-May-2020 19-May-2020 20-May-2020 21-May-2020 22-May-2020 23-May-2020
ARRIVAL DATE
2.5 3 3.5 4 4.5 5 5.5 6 6.5
Figure 10: Porkchop Plot
In order to continue with the fly-by once the spacecraft reaches the other boundary of the Lambert arc, the solution selected was the departure and arrival date times for which the value of ∆VL1
was the lowest.
5.2 Fly-By
As the aim of the project is to use the fly-bys to modify the velocity of the spacecraft in order to decrease the total ∆V required for the mission, i.e., saving fuel, the model of the fly-by was explained in Section 4. In the current section, the previous model is applied for the mission.
of the auxiliary reference frame about the direction of the incoming velocity relative to the moon vector v−∞, i.e., n1. The incoming velocity relative to the moon vector is later turned an angle ”φ”
about the normal vector of the plane (u3), obtaining the outgoing velocity vector relative to the moon v+∞. The outgoing velocity relative to the moon vector v+∞ in the Moon reference system plus the velocity of the moon vmgive the outgoing velocity vector v3 in the Earth reference frame.
v3= v+∞+ vm (26)
The increment of velocity obtained for the fly-by is determined as:
∆VFB= v3− v2 (27)
A similar procedure to the one carried out for the Porkchop has been done in this section. A different value of ”γ” and ”φ” define a different plane for the hyperbola and different rotation of incoming velocity relative to the moon, modifying the eccentricity, the radius of the perigee and the direction of the outgoing velocity relative to the moon. Therefore, the ∆VFB, eccentricity and radius of the perigee have been plotted depending on both ”γ” and ”φ” parameters:
(a) ∆VFB of the Fly By (b) Eccentricity of the orbit after the Fly By Figure 11
Figure 12: Radius of the perigee of the Hyperbolic orbit
However, some constraints must be added to the model to make it feasible:
• Solutions where the radius of the perigee is lower than the radius of the Moon (rper< rM) are not possible as the spacecraft would crash into the Moon.
• The sphere of influence (SOI) of the Moon, which is the spheroid region where the Moon is the primary gravitational influence on an orbiting object. The radius of the sphere of influence is calculated as:
rSOIM E = d · (mM
mE)2/5 (28)
where d is the semimajor axis of the Moon orbit around the Earth, mM is the mass of the Moon and mE is the mass of the Earth.
If the radius of the perigee of the hyperbolic orbit is greater than the radius of the sphere of influence of the Moon (rper > rSOIM E), then, the Moon has no influence on the spacecraft.
• As the objective of the project is to do several fly-bys around the Moon, if the solutions after the fly-by imply a Hyperbolic orbit or an orbit that is not under the gravitational influence of the Earth, another fly-by around the Moon will not be feasible. Thus, in order to do another fly-by, orbits around the Earth with an eccentricity greater than one (e > 1) will be discarded. In addition, if the new orbit of the spacecraft around the Earth after the fly-by has an apogee which is greater than the radius of the sphere of influence of the Earth (rap> rSOIES, the Earth will not have primary influence upon the spacecraft. In that case the body with primary gravitational influence would be the Sun. As the aim is to do another fly-by, these solutions will also be discarded for the current fly-by.
The radius of the sphere of influence of the Earth has been determined as:
rSOIES = dES·(mE
mS)2/5 (29)
where dES is the semimajor axis of the orbit of the Earth about the Sun, and mS is the mass of the Sun.
After the previous restraints are applied to the the Fly-By, the solutions for ∆VFBand eccentricity depending on ”γ” and ”φ” that meet the previous requirements are plotted below:
Assuming that the position of the spacecraft during the fly-by is the same as the one for the Moon with respect to the Earth, rM, and obtaining the velocity of the spacecraft after the fly-by, v3, the parameters of the resulting orbit can be obtained from the state vector [16].
The specific angular momentum can be obtained as:
h= r × v (30)
The inclination of the orbit can be calculated as:
i= hz
h (31)
The node line is defined as:
N= K × h (32)
where K is the third unit vector of the Geocentric reference system. From the previous equation, the right ascension of the ascending node (RAAN or Ω) can be determined:
Ω =
cos−1(NNX) (NY ≥0)
360◦− cos−1(NNX) (NY <0) (33) The eccentricity vector is assessed as:
e= 1
µE[(v2−µE
r )r − r · vrv] (34)
where vr is the radial velocity determined as: vr = r · v/r The argument of the perigee is calculated as:
ω=
cos−1(N·eN e) (eZ ≥0)
360◦− cos−1(N·eN e) (eZ <0) (35) And the true anomaly as:
θ=
cos−1(e·rer) (vr≥0)
360◦− cos−1(e·rer) (vr<0) (36) The semimajor axis can be determined from Vis-viva equation:
a= µE
2(µrE−v22) (37)
The orbit parameters depending on the angles ”γ” and ”φ” applying the previous constraints are plotted in the following figures:
(a) Semimajor axis of the orbit after the Fly-By (b) Inclination of the orbit after the Fly-By Figure 14: Semimajor axis and inclination of the orbit of the spacecraft around the Earth after the Fly-By
(a) Specific angular momentum of the orbit after the Fly-By
(b) True anomaly of the orbit after the Fly-By
Figure 15: Specific angular momentum and true anomaly of the orbit of the spacecraft around the Earth after the Fly-By
(a) RAAN of the orbit after the Fly-By (b) Argument of the perigee of the orbit after the Fly-By
Figure 16: RAAN and Argument of the perigee of the orbit of the spacecraft around the Earth after the Fly-By
5.3 Second Fly-By
To continue with the mission, the selected values of the angles ”γ” and ”φ” for the fly-by have been the ones that provide the highest value of ∆V once the constraints for the fly-by are met.
However, as it will be discussed in the next section, the author acknowledges that these values might not be the optimum values to reach the target afterwards.
As the aim was to do another fly-by, different options were proposed:
1. Propagation of both orbits of the spacecraft and the moon throughout time solving Kepler’s problem. The aim was to determine the time were the spacecraft and the moon crossed at the same position. Neglecting all interactions with the environment, this option would not include fuel consumption. However, the crossing of the Moon and the spacecraft might take years or decades, resulting on the invalidity of this option for the mission.
2. Propagation of the orbit of the spacecraft for a period of time and use of the Lambert arc to link the position of the spacecraft at that time and the position of the moon after a given time of flight. Unlike the first option, this implies fuel consumption but save a large amount of time compared to the previous option. Therefore, the selected option to continue was the second option.
For the second option, a similar procedure to the one carried out on the way to the Moon for the first fly-by was followed. Assuming the spacecraft was located at the same position as the Moon after the first fly-by and having the velocity vector of the spacecraft after the fly-by (v3), Kepler’s problem was used. Given the position and velocity of the spacecraft (state vector) and a given period of time ∆t, the state vector of the spacecraft after ∆t can be obtained. The velocity of the spacecraft after propagating the orbit for a ∆t time is denoted as (v4)
Afterwards, given a time of flight to go to the Moon after propagating the orbit , the position of the Moon was determined for the arrival time (addition of the previous propagation time, ∆t, and time of flight to the date time of the first fly by). Once the initial and final positions of the
required arc as well as the time of flight between both positions are known, Lambert’s problem is solved to calculate the Lambert arc that links both positions.
For the Lambert arc, as stated in previous subsections, orbit parameters are determined and so are the velocity vectors at the boundaries of the arc. The necessary increment of velocity for the Lambert arc is calculated as the difference between the velocity at the boundary of the orbit of the spacecraft for the Lambert arc (vL4) and the velocity of the spacecraft propagated on the orbit after the first fly-by (v4).
∆VL2= vL4− v4 (38)
As it was done for the calculation of the first Lambert arc, modifying the propagation time of the orbit after the first fly-by, ∆t, and the time of flight for the Lambert arc will lead to new trajectories and new values for the orbit parameters as well as different velocity vectors at the boundaries of the arc. Again, depending on the departure and arrival date times on the way to the moon, ∆VL2 can be plotted.
Figure 17: Porkchop Plot
For the next fly-by, the selected ∆t and Time of Flight were the ones for which the value of ∆VL2
was the lowest.
The second Fly-By was calculated following the same procedure done in section 5.2 and applying the constraint where the radius of the perigee must be greater than the radius of the Moon
(a) ∆VFB for the second Fly-By (b) Eccentricity of the orbit after the second Fly- By
Figure 18
For Figure 18, only the constraint related to the radius of the Moon has been applied. If only two fly-bys about the Moon are required, as the target (the asteroid) orbits about the Sun, its velocity is given in the Heliocentric reference frame. Therefore, in order to determine the outgoing velocity of the spacecraft after the fly-by in the Heliocentric reference system, the velocity of the Earth relative to the Sun must be added to the previous outgoing velocity of the spacecraft in the Earth reference frame. It must be highlighted cases where the orbit about the Earth after the fly-by is elliptical but the apogee of the orbit is larger than the radius of the sphere of influence of the Earth. In the apogee of the orbit, the spacecraft will not be primarily under the gravitational influence of the Earth anymore. The velocity of the spacecraft at such position in the Sun reference frame is calculated adding the velocity of the Earth relative to the Sun plus the velocity of the spacecraft relative to the Earth .
(a) ∆VFB for the second fly-by if more fly-bys are required
(b) Eccentricity of the orbit after the second fly- by if more fly-bys are required
Figure 19
In case that other fly-bys about the moon are required, solutions where either the orbit is hyperbolic
around the Earth (e > 1) or the apogee of the orbit is greater than the radius of the sphere of influence of the Earth, have been discarded. Otherwise, another fly-by would not be feasible. In order to do a third fly-by, it has been selected the values of ”γ” and ”φ” that allow the highest increment of velocity for the second fly-by (∆VFB) in Figure 19. However, the author must acknowledge that those values might not be the optimum values to reach the target afterwards.
5.4 Third Fly-by and trajectory to reach the target
For the third fly-by, the procedure carried out was the same as the one done in the previous subsection 5.3. The orbit of the spacecraft was propagated after the second fly-by for a period of time (∆t2). After that, the Lambert arc that links such position of the spacecraft on its orbit and the position of the Moon at the arrival time (addition of the propagation time of the spacecraft after the second fly-by (∆t2) plus the time of flight (TOF) from the position of the spacecraft to the moon) is determined. Different values of ∆t and TOF give different trajectories for the Lambert arc. Thus, different velocities at the boundaries of the Lambert arc. The necessary increment of velocity on the boundary of the Lambert arc where the spacecraft is located depending on different departure and arrival times has been plotted in a Porkchop plot below.
Figure 20: Porkchop Plot
From figure 20, the selected departure and arrival times for the Lambert arc to continue with
(a) ∆VFBfor the third Fly-By (b) Eccentricity of the orbit after the third Fly- By
Figure 21
After the fly-by and assuming that no more fly-bys around the Moon will be done (otherwise, the same procedure would have to be followed), the velocity of the Earth relative to the Sun is added to the velocity of the spacecraft relative to the Earth after the fly-by. The final velocity of the spacecraft will be expressed in a Heliocentric reference frame.
vs/cOrbit/Sun= vs/cOrbit/Earth+ vs/cEarth/Sun (39) where
vEarth/Suns/c = vEarthEarth/Sun+ ωEarth/Sun× dEarth/Spacecraft (40) and ωEarth/Sun is the angular velocity vector of the Earth relative to the Sun.
As the distance from the Earth to the Sun is much larger than the distance from the spacecraft to the Earth, the last term of the previous equation has been neglected leading the previous equation to vs/cEarth/Sun≈ vEarthEarth/Sun.
In order to continue with the mission and reach the target, a similar procedure to the one carried out previously to go to the Moon would have to be followed. The orbit of the spacecraft in the Sun reference frame would have to be propagated using Kepler’s problem throughout a period of time (∆tSun), obtaining the new state vector of the spacecraft for a given time.
The state vector of the target should be provided for an arrival time. The arrival time is calculated as the addition of a time of flight (T OFSun) to go from the spacecraft to the target once the orbit of the spacecraft has been propagated a period of time (∆tSun). Given both positions of the spacecraft and the target as well as the time of flight between them (T OFSun), the Lambert arc that links both boundaries can be determined. The Lambert arc provides the velocity at the boundaries of the arc. The necessary increment of velocity at the boundary of the Lambert arc is calculated as the difference between the spacecraft’s velocity vector at the boundary of the Lambert arc and the velocity vector of the spacecraft propagated previously using Kepler’s problem.
In addition, on the other boundary of the Lambert arc, the one belonging to the target, another manoeuvre would be necessary to meet the requirements for injection into the orbit around the target.
Unfortunately, the ephemeris of the asteroid or another target was not provided, so calculations could not be carried out for this part of the thesis. However, the same process could be followed to reach any other celestial body.
6 Optimisation
In the previous Section 5, it was selected the minimum value of the increment of velocity vector when solving Lambert’s problem or the maximum increment of velocity vector in the case of the fly-by. However, the problem was analysed individually and not globally. Thus, the solution selected for each part might not be the most optimum one globally to reach the target.
In order to set the optimisation problem, the first question arisen was how to model it? A necessary requirement that came up to simplify the problem was to define the number of fly-bys. Therefore, to include several fly-bys around the Moon but to avoid having too many variables, two fly-bys were considered for the optimisation problem. For the case of three or more fly-bys, the same variables used for a fly-by in the problem should be used but, as the number of fly-bys would increase, the number of variables of the decision vector would be greater too.
The aim of the mission is to reach the target with the minimum fuel possible saving weight and money. Therefore, fuel is only consumed to modify the velocity vector for the Lambert arc.
The considered fly-by is an unpowered fly-by, so no fuel is used. Thus, the objective would be to minimise the different increments of velocity when the Lambert arc is calculated. Considering only two fly-bys, the objetive function or cost function tries to get the minimum increment of velocity for the two Lambert arcs before both fly-bys plus the increment of velocity for the Lambert arc to reach the target after the second fly-by.
to minimise: J(x) =X2
i=1
∆VLi(x) + ∆VT(x) (41)
where x is the decision vector that will be explained later, ∆VLi is the increment of velocity for each Lambert arc and ∆VT is the increment of velocity after the second fly-by and after having propagated the orbit throughout time in order to determine the Lambert arc to reach the target.
A huge issue for the optimisation problem is the dimensionality, affecting the computational cost.
[17] states that a grid of N points in one dimension would have ND points in D dimensions.
Nevertheless, the number of points can differ for each dimension. For this project, for the two fly-by case, there are 10 dimensions considered which are gathered in the decision vector shown below:
x= [Tdep, T OFL1, γ1, φ1, P T1, T OFL2, γ2, φ2, P T2, T OFL3] (42) where Tdep is the departure date time from the GTO, T OFL1 and T OFL2 are respectively the Time Of Flight for the first Lambert arc that links the GTO and the Moon and the Time Of Flight for the second Lambert arc, γ and φ are the angles that define the fly-by. ”γ” defines the angles that the normal vector of the plane forms with the reference and φ defines the rotation of the incoming velocity in the Moon reference system about the normal vector to the plane. As there are two fly-bys, it is considered twice. P T1 and P T2 are the period of time for which the outgoing orbit is propagated after the first and second fly-by respectively. T OFL3 is the Time Of Flight for the Lambert arc that links the propagated orbit after the second fly-by and the target in the Heliocentric reference frame.
At this point, the author must acknowledge that all work developed in Section 5 had been done for date times in ”Day-Month-Year” format. That complicates the iteration process for the opti- misation. Therefore, the author suggest modifying the format to Julian date, where all values are numerical, so easing the addition of different period of times to a specific date time.
Several constraints must be added to the optimisation problem:
• In the previous section, thrust limitations have not been considered and it is likely that it will not affect to the current mission. Either way, the maximum amount of thrust that the propulsion system of the spacecraft can provide to modify the velocity must be taken into account. In case the cost function exceeds the maximum increment of velocity that the propulsion system can provide, that trajectory will be unfeasible and thus, pruned out.
• Minimum Altitude: As stated previously for the Fly-By, the radius of the perigee of the hyperbolic orbit on the Moon reference frame must be greater than the radius of the Moon.
In case the radius of the perigee is lower than the one of the Moon, the trajectory will be considered unfeasible and will be discarded.
• A constraint only applied to the first fly-by is that if the outgoing orbit on the Earth reference frame is hyperbolic or the apogee of the orbit exceeds the radius of the Sphere of Influence of the Earth, then, it is discarded as it makes another fly-by around the Moon impossible.
The solution step summary for the optimisation is:
1. Compute the state vector of the spacecraft on the GTO orbit for a departure time (Tdep).
2. Adding the Time of Flight (T OFL1) to the departure time (Tdep), the arrival time to the Moon is computed. For the arrival time, the state vector (position and velocity) of the Moon can be computed using the ephemeris routine from CSPICE software.
3. To solve the Lambert’s problem determining ∆VL1 on the boundary of the Lambert arc located at the GTO orbit.
4. To solve the fly-by depending on ”γ1” and ”φ1” angles and propagate the outgoing orbit in the Earth reference frame for a period of time (P T1).
5. To solve the Lambert’s problem determining ∆VL2 on the boundary of the Lambert arc located at the previous propagated orbit. The objective is to find the combination of ”γ1” and ”φ1” angles for the fly-by and the propagation time (P T1) that minimise the increment of velocity vector for the Lambert arc.
6. To solve the fly-by depending on ”γ2” and ”φ2” angles.
7. Compute the velocity of the spacecraft in the Heliocentric reference frame. Therefore, com- pute the velocity of the Earth relative to the Sun for the given time using CSPICE toolkit and add it to the velocity obtained for the second fly-by for the spacecraft relative to the Earth.
8. To propagate the outgoing orbit in the Heliocentric reference frame for a period of time (P T2).
9. To solve the Lambert’s problem determining ∆VT on the boundary of the Lambert arc located at the propagated orbit of the spacecraft. The objective is to find the combination
7 Future Work
• Even though an approach to the optimisation process has been done, due to lack of time, the whole optimisation problem has not been solved. First, the code structure should be adapted for the optimisation process. Initially, all the code programmed aimed to model the fly-by and check the development of the different elements of the mission. However, it was not organised in a way that it could be used for optimisation. Therefore, the structure of the code must be changed to allow optimisation. Second, the mission must be optimised to determine the decision vector that minimises the cost function.
• For the current project, only unpowered fly-bys have been used. However, powered fly-bys combining gravity and propulsive action should also be studied as well as the combination of powered and unpowered fly-bys [18].
• The current model does not take into account any interaction between the spacecraft and the environment or any perturbation, simplifying computation. All previous studies assumed that the Earth was spherical and it was the only body exerting a force on the spacecraft (apart of the fly-by). However, the Earth is not spherical, nor homogeneous. All other solar system bodies exert a force on the spacecraft and the spacecraft is moving in interplanetary space which is not really empty (plasma, solar wind, electromagnetic field and radiation [19].
In order to obtain a more realistic mission, previous perturbations should be consideration when existed.
8 Conclusions
In this thesis, a method to reduce fuel consumption for space missions has been investigated. As current state-of-the-art of propulsion systems does not allow space vehicles to reach every body in the Solar System, Multiple Gravity Assist Manoeuvres (MGA) are used to modify the velocity of the spacecraft reaching high-energy targets with a limited amount of propellant. Therefore, MGA manoeuvres have been studied, modelled and applied for a mission to a target (e.g. an asteroid).
The Fly-By has been modelled in order to allow the user choose the position (based on the angle that the normal vector of the plane forms with respect to a reference) where the incoming velocity of the spacecraft relative to the Moon pierces the sphere of influence of the Moon. It also allows the user to select the angle that the previous velocity vector rotates about the normal vector of the plane. Designing the Fly-By based on the previous two parameters, which are two angles, allows the user to plot different outputs such us eccentricity, radius of the perigee and increment of velocity depending on the the previous angles. That eases the selection of the previous parameters to meet any given requirement.
Given a few initial scripts to solve Kepler’s Problem, Lambert’s Problem and some parameters for the development of the mission, the author had to deal with concepts of Astrodynamics that he had not studied before the beginning of the internship/research. Thus, the development of the project was delayed. In addition, as stated in Section 7, the structure of the code was oriented to calculate each section individually (e.g. selecting a solution from a section of the code and applying it to the next section without considering the influence that a different solution might have on other sections of the code which might lead to better results). On the other side, even though the mission was not optimised to get the minimum fuel consumption, the results plotted on Section 5 would be feasible. Therefore, the procedure proposed for the development of the mission would be accepted in the absence of the optimised solution or in case no optimal solution is required.
Due to lack of time, the optimisation part has not been completed, but an initial approach has been carried out on how to model it. The objective function to optimise has been defined and all variables taking part in the process have been gathered in a decision vector. In addition, perturbations and interactions between the spacecraft and environment should be added to the previous mission to obtain a more realistic model.
Overall, the project has been an opportunity to take a deep look at the development of space missions designing and combining several aspects of the orbit such as propagation of the orbit using Kepler’s problem, linking different positions in space through Lambert’s problem or modifying both direction and magnitude of the velocity of a spacecraft flying by the Moon through Gravity assist manoeuvres. In addition, the project has shown current technological limitations of space missions and how to deal with them to reach the target.
References
[1] D. W. Reynolds. Apollo, The epic journey to the Moon, 1963-1972. Zenith Press, 2013. isbn:
978-0-7603-4452-1.
[2] G. P. Kennedy. Germany’s V-2 Rocket. Schiffer Publishing Ltd, 2006. isbn: 978-07643245297.
[3] J. Doran P. Bizony. Starman: The truth behind the legend of Yuri Gagarin. Bloomsbury Pub- lishing PLC, 1999. isbn: 978-0747542674.
[4] A. Chaikin. A man on the Moon: The voyages of the Apollo astronauts. Penguin, 1998. isbn:
978-0143112358.
[5] M. A. Minovitch. The invention that opened the solar system to exploration. 2010. doi: 10.
1016/j.pss.2010.01.008.
[6] D. L. Shirley. Mariner 10 mission analysis - Application of a black art. In Proceedings of American Institute of Aeronautics and Astronautics, Aerospace Sciences Meeting California (USA), 1975.
[7] V. T. Toth S. G. Turyshev. The Pioneer Anomaly. 2010. doi: 10.12942/lrr-2010-4.
[8] P. A. Penzo E. C. Kohlhase. Voyager mission description. Vol. 21. 1977, pp. 77–101.
[9] P. A. Penzo R. E. Diehl D. I. Kaplan. Satellite tour design for the Galileo mission. 1983.
[10] S. Flanagan F. Peralta. The Pioneer Anomaly. Vol. 3. 1995, pp. 1603–10. doi: 10.1016/
0967-0661(95)00171-P.
[11] M. Ceriotti. Global optimisation of multiple gravity assist trajectories. PhD thesis, Univer- sity of Glasgow, 2010.
[12] NewHorizons. http://solarsystem.nasa.gov/missions/new-horizons/in-depth/#:˜:
text=New%20Horizons%20is%20a%20NASA, 50%20AU%20from%20the%20Sun.&text=The%
20spacecraft%20was%20now%20set,reaches%20of%20the%20solar%20system.. [13] C.T. Russell. New Horizons. Springer, 2001. isbn: 978-0-387-89517-8.
[14] D. A. Vallado. Fundamentals of Astrodynamics and Applications. Space Technology Library, 2013. isbn: 978-1881883180.
[15] R. H. Battin. An Introduction to the Mathematics and Methods of Astrodynamics. Ameri- can Institute of Aeronautics and Astronautics, 1999. isbn: 1-56347-342-9.
[16] H. D. Curtis. Orbital Mechanics for Engineering Students. Elvesier, 2005. isbn: 0 7506 6169 [17] K. Qadir. Multi-gravity assist design tool for interplanetary trajectory optimisation. MSc the-0.
sis, Lulea University of Technology- Cranfield University, 2008.
[18] D.R. Myatt D. Izzo V.M. Becerra. “Search space pruning and global optimisation of multi- ple gravity assist spacecraft trajectories”. In: J Glob Optim 38:283–296 (). doi: 10.1007/
s10898-006-9106-0.
[19] H. Garrett D. Hastings. Spacecraft-Environment Interactions. Cambridge University Press, 2004. isbn: 0 521 47128 1.
