Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Mission design study of an RTG powered, ion engine equipped interstellar spacecraft
(USC Thesis Other)
Mission design study of an RTG powered, ion engine equipped interstellar spacecraft
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
MISSION DESIGN STUDY OF AN RTG POWERED, ION ENGINE EQUIPPED
INTERSTELLAR SPACECRAFT
by
Joshua A. Fogel
______________________________________________________________________________
A Thesis Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
MASTER OF SCIENCE
(ASTRONAUTICAL ENGINEERING)
May 2014
Copyright 2014 Joshua A. Fogel
Table of Contents
List of Figures ................................................................................................................................ 4
List of Tables .................................................................................................................................. 6
Abbreviations ................................................................................................................................. 7
Abstract .......................................................................................................................................... 9
Chapter 1: Introduction ................................................................................................................ 11
1.1 The Importance of the Interstellar Precursor Mission ............................................... 12
1.2 Mission Simulator Requirements ............................................................................... 14
1.3 Scope of This Thesis Work ........................................................................................ 16
Objectives ............................................................................................................ 16
Assumptions ......................................................................................................... 17
Chapter 2: Simulating the IPM .................................................................................................... 18
2.1 Mission Physics ......................................................................................................... 18
The Solar Gravitational Environment .................................................................. 19
The Ion Engine ..................................................................................................... 20
2.2 Simulator Implementation ......................................................................................... 23
Simulator Structure and Operation ...................................................................... 23
Initial Conditions & Escape Trajectory Characterization .................................... 27
Orbit Propagation: The Fourth Order Runge-Kutta Method ............................... 33
Total Acceleration Calculation ............................................................................. 36
2.3 Software Validation .................................................................................................... 37
Chapter 3: Simulation Methods ................................................................................................... 39
3.1 Simulator Settings ...................................................................................................... 39
3.2 Input Parameter Sets and Scenarios ........................................................................... 40
Chapter 4: Simulation Results and Conclusions .......................................................................... 42
4.1 Example Simulation Outputs...................................................................................... 42
Output Plots from Scenario C4 ............................................................................ 42
Example Fringe Scenario: Plots from Scenario E3 .............................................. 48
4.2 Orbit and Mission Architecture Analysis ................................................................... 52
Total Flight Time to 200 AU ................................................................................ 52
Heliocentric Distance at ECO .............................................................................. 54
Orbit Characterization .......................................................................................... 57
4.3 Ion Engine Model Analysis ........................................................................................ 62
4.4 Initial Position Comparison: Circular Earth Orbit, Perihelion, Aphelion .................. 71
2
Perihelion Launch ................................................................................................ 71
Aphelion Launch .................................................................................................. 72
4.5 Conclusions ................................................................................................................ 73
Mission Architecture for Circular Earth Orbit Launches ..................................... 74
Ion Engine Performance ...................................................................................... 76
Initial Position Comparison ................................................................................. 77
4.6 Suggestions for Future Work ..................................................................................... 77
References .................................................................................................................................... 79
Appendix A: Raw Simulation Output Data: ................................................................................ 80
Initial Position: Earth Orbit, Circular ............................................................................... 80
A.1 Simulation Output Data: 5 year burn time .................................................... 80
A.2 Simulation Output Data: 10 year burn time .................................................. 81
A.3 Simulation Output Data: 15 year burn time .................................................. 82
Appendix B: Raw Simulation Output Data: ................................................................................ 84
Initial Position: Earth Orbit, Perihelion ........................................................................... 84
B.1 Simulation Output Data: 15 year burn time .................................................. 84
Initial Position: Earth Orbit, Aphelion ............................................................................. 85
B.2 Simulation Output Data: 15 year burn time .................................................. 85
Appendix C: ................................................................................................................................ 87
Matlab Code: .................................................................................................................... 87
C.1 Main Script - “Orbit_prop.m” ..................................................................... 87
C.2 Fourth Order Runge-Kutta Function - “rk4_Orb.m” ................................... 93
C.3 Total Acceleration Function - “A_tot.m” ...................................................... 95
C.4 Acceleration Function: Gravity - “A_grav.m” ............................................. 96
C.5 Acceleration Function: Thruster - “A_thrust.m” ......................................... 96
C.6 Instantaneous Thruster Parameters Function - “thruster_params.m” ........ 96
3
List of Figures
Figure 1: Solar System Structure (logarithmic scale in units of AU) [5] .................................... 11
Figure 2: Accelerations felt by the spacecraft during engine-burn, with depictions of
instantaneous orbits at various moments in time ......................................................................... 18
Figure 3: IPM Spacecraft Simulator operational block diagram ................................................ 23
Figure 4: Heliocentric Inertial Reference Frame and Initial Conditions .................................... 28
Figure 5: Hyperbolic Orbit Geometry for -[ey]unit ≤ 0 ................................................................ 32
Figure 6: Simulation hierarchy, illustrating simulation Set C ..................................................... 40
Figure 7: Scenario C4: Trajectory ............................................................................................... 42
Figure 8: Scenario C4: Heliocentric Speed vs. Position ............................................................. 43
Figure 9: Scenario C4: Orbit Semi-major Axis vs. Time ............................................................ 45
Figure 10: Scenario C4: Orbit Eccentricity vs. Time .................................................................. 46
Figure 11: Set C: Spacecraft Total Mass vs. Time ...................................................................... 48
Figure 12: Scenario E3: Trajectory ............................................................................................. 49
Figure 13: Scenario E3: Heliocentric Speed vs. Position ........................................................... 50
Figure 14: Scenario E3: Heliocentric Speed vs. Time ................................................................ 51
Figure 15: Scenario E3: Orbit Eccentricity vs. Time .................................................................. 51
Figure 16: Total Flight Time to 200 AU vs. Launch Velocity for burn time of 5 ....................... 52
Figure 17: Flight Time to 200 AU vs. Launch Velocity for burn time of 10 years ..................... 53
Figure 18: Flight Time to 200 AU vs. Launch Velocity for burn time of 15 years ..................... 53
Figure 19: Heliocentric Position @ ECO vs. Launch Velocity, 5 year burn time ....................... 55
Figure 20: Heliocentric Position @ ECO vs. Launch Velocity, 10 year burn time ..................... 56
4
Figure 21: Heliocentric Distance @ ECO vs. Launch Velocity, 15 year burn time .................... 56
Figure 22: Heliocentric Excess Velocity vs. Launch Velocity for burn time of 15 years ........... 58
Figure 23: Solar Escape Target Angle vs. Launch Velocity for burn time of 15 years ............... 59
Figure 24: Flight Time to Escape vs. Launch Velocity for burn time of 15 years ...................... 60
Figure 25: Final Orbit Semi-major Axis vs. Launch Velocity for burn time of 15 years ........... 61
Figure 26: Final Orbit Eccentricity vs. Launch Velocity for burn time of 15 years ................... 62
Figure 27: Ion Engine Beam Current vs. Burn Time (analytic) .................................................. 63
Figure 28: Ion Engine Accelerating V oltage vs. Burn Time (analytic) ....................................... 64
Figure 29: Propellant Mass Flow Rate vs. Burn Time (analytic) ............................................... 65
Figure 30: Ion Engine Thruster Force vs. Burn Time (analytic) ................................................. 66
Figure 31: Ion Engine Beam Exit Velocity vs. Burn Time (analytic) ......................................... 68
Figure 32: Ion Engine Isp vs. Burn Time (analytic) .................................................................... 69
Figure 33: Ion Engine ∆V vs. Burn Time (analytic) ................................................................... 70
Figure 34: Difference in Total Flight Times between circular and perihelion launches ............. 71
Figure 35: Difference in Total Flight Times between circular and aphelion launches ............... 72
5
List of Tables
Table 1: Set A Trajectory Results (tb = 5 years, P0 = 1000 W) .................................................... 80
Table 2: Set B Trajectory Results (tb = 5 years, P0 = 1500 W) ................................................... 80
Table 3: Set C Trajectory Results (tb = 5 years, P0 = 2000 W) .................................................... 80
Table 4: Engine Outputs for (tb = 5 years), Sets A, B & C .......................................................... 81
Table 5: Set D Trajectory Results (tb = 10 years, P0 = 1000 W) .................................................. 81
Table 6: Set E Trajectory Results (tb = 10 years, P0 = 1500 W) .................................................. 81
Table 7: Set E Trajectory Results (tb = 10 years, P0 = 2000 W) .................................................. 81
Table 8: Engine Outputs for (tb = 10 years), Sets D, E & F ......................................................... 82
Table 9: Set G Trajectory Results (tb = 15 years, P0 = 1000 W) ................................................. 82
Table 10: Set H Trajectory Results (tb = 15 years, P0 = 1500 W) ............................................... 82
Table 11: Set I Trajectory Results (tb = 15 years, P0 = 2000 W) ................................................. 82
Table 12: Engine Outputs for (tb = 15 years), Sets G, H & I ....................................................... 83
Table 13: Set Ap Trajectory Results (tb = 15 years, P0 = 1000 W) .............................................. 84
Table 14: Set Bp Trajectory Results (tb = 15 years, P0 = 1500 W) .............................................. 84
Table 15: Set Cp Trajectory Results (tb = 15 years, P0 = 2000 W) .............................................. 84
Table 16: Engine Outputs for (tb = 15 years), Sets Ap, Bp & Cp ................................................... 85
Table 17: Set Aa Trajectory Results (tb = 15 years, P0 = 1000 W) .............................................. 85
Table 18: Set Ba Trajectory Results (tb = 15 years, P0 = 1500 W) .............................................. 85
Table 19: Set Ca Trajectory Results (tb = 15 years, P0 = 2000 W) .............................................. 85
Table 20: Engine Outputs for (tb = 15 years), Sets Aa, Ba & Ca ................................................... 86
6
Abbreviations
Symbols:
UG gravitational potential
µ standard gravitational parameter (m
3
/s
2
)
C3 characteristic energy (m
2
/s
2
)
R heliocentric position (m)
V heliocentric velocity (m/s)
A linear acceleration (m/s
2
)
Rp radius of perihelion (m)
Ra radius of aphelion (m)
a orbital semi-major axis (m)
p orbital semi-parameter (m)
e orbital eccentricity
∆V “delta-V”, change in velocity (m/s)
Vesc heliocentric escape velocity
V∞ hyperbolic excess velocity (m/s)
ξ solar escape target angle (deg)
kn Runge-Kutta integration increment
h integration step size (sec/step)
tb engine burn time
P0 total input electric power (W)
Pb beam power (W)
Vb accelerating voltage (V)
Ib beam current (A)
FThrust ion engine thruster force (N)
Φ total impulse (N•sec)
Isp specific impulse (sec)
Ve ion engine exit velocity (m/s)
Mxe atomic mass of Xenon (kg)
Mo spacecraft structural dry mass (kg)
Mp propellant mass (kg)
m propellant mass flow rate (kg/s)
Acronyms:
IPM Interstellar Precursor Mission
RTG Radioisotope Thermoelectric Generator
LISM Local Interstellar Medium
ECO Engine Cut-Off
SMA Semi-Major Axis
7
EP Electric Propulsion
LEP Low-thrust Electric Propulsion
Constants:
AU astronomical unit: AU = 149597870700 (m)
echarge fundamental charge unit: echarge = 1.602176565*10
-19
(C)
AMU atomic mass unit: AMU = 1.66053886*10
-27
(kg)
Mxe atomic mass of Xenon: Mxe = 131 (AMU)
yr Earth year: yr = 3.15569088*10
7
(sec)
Planetary Values, Taken on July 29, 2013 [7]:
µsun solar standard gravitational constant: µsun =1.32712440018*10
20
(m
3
/s
2
)
RP-Earth Earth perihelion: RpEarth = 0.983123587067824 (AU)
RA-Earth Earth aphelion: RaEarth = 1.016517960841986 (AU)
a Earth Earth semi-major axis: aEarth = 0.9998207739549034 (AU)
8
Abstract
This research explores a variety of mission and system architectures for an unmanned
Interstellar Precursor Mission (IPM) spacecraft with a Radioisotope Thermoelectric Generator
(RTG) powered Ion Engine using Xenon propellant, traveling on a (direct) ballistic escape
trajectory to the undisturbed Interstellar Medium (~200 AU). The main goal of this work was to
determine the relationship between the propulsion system design parameters and the ensuing
escape trajectory. To do this, an orbit simulator was created in Matlab using a fourth order
Runge-Kutta numerical integration method to propagate the thrusting spacecraft’s trajectory
through time. The accelerations due to the Sun’s gravity and the Ion Engine thrust were modeled
separately and then combined into a single total acceleration vector at each time step, with the
thrust direction assumed to be in the direction of the spacecraft’s instantaneous velocity vector.
The propellant of the thruster was also designed to be completely consumed by the time of
engine cut-off (ECO), meaning a constant propellant mass flow rate. Simulations were run for
burn times of 5, 10 & 15 years, with heliocentric launch velocities of 0, 5, 7, 10 & 12 km/sec
from a circular 1 AU Earth orbit, and with RTG supplied engine input powers of 1000, 1500 &
2000 W. A total of 45 simulations were run for the circular 1 AU case, as well as additional
comparison simulations for launches from an elliptical Earth orbit at perihelion and aphelion.
The results of these simulations yielded many interesting results on the total fly-out times
to 200 AU, which ranged dramatically from ~35 to ~140 years depending on the propulsion
system settings and orbital initial conditions, as well as descriptions of the ECO distances from
the Sun for each mission. The simulations also revealed the inherent gravitational maneuver
9
inefficiency felt by all low thrust spacecraft, which becomes more apparent under certain
conditions. Relations between launch velocity and the various characteristics of the escape
trajectory were also shown. Finally, the propulsion system model showcased many details about
the performance of the Ion Engine, including the consequence that increasing the input power
level to the engine will result in a greater increase of thrust when operating at shorter engine burn
times.
10
Chapter 1: Introduction
Figure 1: Solar System Structure (logarithmic scale in units of AU) [5]
On August 25th, 2012, the deep space probe Voyager 1 officially crossed the heliopause,
the boundary shock separating the heliosphere, the region of space dominated by the solar wind,
and the Local Interstellar Medium (LISM), becoming the first manmade object to reach
interstellar space. This historic milestone of the Voyager 1 mission, launched on September 5th,
1977, marks a nearly 122 Astronomical Units (AU) journey, which included the capture of the
first detailed images of the Jovian and Saturnian systems as part of V oyager’s grand tour of the
outer planets. Now in its extended mission phase, Voyager 1 will continue its travel out of the
solar system, its suite of instruments continuing to study the nature and composition of the LISM
until its power generators can no longer support its scientific systems (expected ~2025). While
the V oyager mission represents one of the greatest achievements in deep space exploration to
11
date, current scientific models of the heliopause and LISM that describe the interaction of matter
and electromagnetic fields between the Sun and the surrounding galaxy are still rather primitive,
and require further investigation.
1.1 The importance of the Interstellar Precursor Mission
Even today, very little is known about the nature of the LISM. In order to better
understand the composition and dynamics of the LISM, which would allow for more confident
travel through interstellar space in the future, experimental data must be collected on these
regions. Unfortunately, the highest fidelity scientific measurement methods currently available,
such as spectroscopy, interferometry, and particle sampling, can only be accomplished with
instruments in situ. Significant questions relating to the composition of matter fundamental to
the formation of habitable planets and life could be answered by direct sampling of the LISM.
Due to the expense of such ambitious missions, as well as the technical difficulties of
maintaining spacecraft hardware during the decades-long travel time, investigations into a
reference mission called the “Interstellar Precursor Mission” (IPM) are being conducted to
determine the suitable mission architecture necessary to optimize such an endeavor. The mission
would require a high speed launched, low mass spacecraft to travel to roughly 200 AU from the
Sun in order to reach a region of interstellar space unperturbed by the solar wind and magnetic
field, and to do so within a reasonable time frame (a “reasonable” time frame can be defined as
about half the professional lifetime of a scientist or engineer, or on the order of ~20 years). Since
12
plasma interactions between the heliosphere and the LISM may occur as near as ~70-80 AU, all
propulsion systems with engine exhaust must be shut off well before reaching this boundary so
as to not interfere with the probe’s array of sensors. Because the heliophysics side of the IPM is
already well defined, meaning no significant studies on scientific instruments and their related
topics are necessary, further studies must be conducted on the capabilities of high energy launch
vehicles, and more importantly the spacecraft propulsion requirements.
The main engineering problem associated with any interstellar mission is how to design
the spacecraft’s propulsion system to adequately accelerate it to high solar escape velocities.
Voyager 1 is currently traveling at roughly ~3.6 AU/yr, thanks primarily to the gravitational assist
fly-by’s of Jupiter and Saturn, but a future IPM would require at least 2-3 times this fly-out
speed. This will likely be achieved through a combination of gravity assists and a variety
propulsion and launch methods, but a simple direct ballistic escape method could just as easily
be effective. While many propulsion architectures have been proposed, ranging from solar sails
to nuclear-pulse propulsion, the type gaining popularity is Electric Propulsion (EP) where an
external power source (typically solar arrays or nuclear reactors) supplies energy to an engine,
where it is converted into kinetic energy and thereby thrust (the methods of which vary
depending on EP type). The benefit of EP engines over traditional chemical engines (used
primarily in launch vehicles, which have high-thrust, low-efficiency) is their high-efficiencies,
measured in Specific Impulse, Isp (units of seconds), the impulse per unit of propellant. This
allows the system to achieve the same total change in velocity (∆V) as a comparable chemical
rocket (which have Isp’s on the order of ~400 sec, compared to EP systems which can have Isp’s
upwards of 1000 sec), but with far less propellant. The draw back is that EP systems provide far
13
lower thrust than chemical systems, which means they can only be used once placed in orbit. As
a result, the required engine burn times required to achieve a given ∆V for EP systems can be on
the order of years, compared to chemical engines which typically burn on the order of minutes.
For the Interstellar Precursor Mission, the typical low-thrust EP (LEP) system of choice
is the Ion Engine paired with a Radioisotope Thermoelectric Generator (RTG), a power source
that uses thermocouples to convert heat from a radioactive fuel into electricity by way of the
Seebeck effect. RTG’s are favored by deep space probes due to their relatively high power
output (~10
3
W) and long operational life spans, which with modern advancement (and the 80
year half-life of the standard plutonium fuel) can maintain a suitable output power threshold up
to 40 years. This system combination demands a low-mass spacecraft due to its low thrust
output, but boasts an Isp considerably higher than alternative propulsion methods (with values
upwards of 3000 seconds depending on propellant type and flow rate). While Ion Engines have
become an increasingly popular design choice for deep space probes and geosynchronous
satellites (for station keeping) over the past 30 years, the relationship between the propulsion
system’s design parameters and the resultant escape trajectories are not readily apparent.
1.2 Mission Simulator Requirements
Before a full-scale mission can be undertaken, thorough mission and system level
simulations must be conducted to determine the best combination of design parameters for the
spacecraft’s hardware and mission architectures. To do this, mathematical models of each
14
spacecraft system and environmental influence, such as the Sun’s gravitational acceleration, must
be constructed. It should be noted that while these models are designed to be as realistic as
possible, they are still approximations that contain simplifying assumptions, mostly due to the
inherent constraining nature of a computerized representation of the natural universe. Once the
modeling equations are selected, they must be combined into a set of differential equations, the
solutions to which represent the spacecraft’s motion through time. Since it is difficult to obtain
closed form solutions to the differential equations of motion for a thrusting spacecraft with
variable mass, numerical integration methods must be employed which range in complexity and
accuracy. Instead of producing a continuous analytical solution, numerical methods calculate an
approximate solution as a set of values at a discrete number of points (in the case of the thrusting
spacecraft under the influence of gravity, the numerical solution is an array of points representing
the spacecraft’s three-dimensional position at different moments in time, known collectively as
its trajectory). The user is free to choose the step size in time between these points, which allows
control over the accuracy of the numerical solution and the real-time duration of the integration
process. The most common of these methods is the fourth order Runge-Kutta process, which can
be tailored specifically to the thrusting spacecraft problem.
After the models and integrator are built, the simulator’s accuracy can be verified by
comparing the integration results of a spacecraft with zero thrust in a given orbit with the
analytical solution of the known orbit. Once confirmed, the simulator can then be confidently
used to demonstrate the effects of different mission and spacecraft designs on the ensuing
trajectory, as well as to characterize the performance of the modeled systems. These results can
then be applied to future mission design efforts, allowing a spacecrafts’ hardware and mission
15
architecture to be more appropriately designed. For the IPM, simulation results should include
heliocentric positions and velocities over time, flight times to the LISM, characterizations of the
pre-engine cut-off (ECO) and post-ECO trajectories, and descriptions of the performance of the
propulsion system (including the thrust, propellant flow rate, engine efficiency, and electrical
performance).
1.3 Scope of This Thesis Work
Objectives
The primary goal of this thesis work is to simulate a LEP interstellar probe mission to the
LISM via a ballistic, or direct, escape trajectory (no gravity assists) in order to develop a top-
level understanding of the relationship between the propulsion system’s design and the ensuing
escape trajectory. Additional objectives which enable this goal are as follows:
1. Construct a basic model of a spacecraft with an RTG powered Ion Engine.
2. Construct a basic model of the solar gravity field.
3. Construct a fourth order Runge-Kutta numerical integrator, tailored for
orbit propagation, with reasonable accuracy.
4. Run the simulator for different spacecraft designs and initial conditions
that enable the probe to reach the LISM.
5. Obtain useful output data to characterize the probe’s resultant trajectory
and engine performance for each set of design parameters.
16
Assumptions
This thesis work makes the following simplifying assumptions:
1. The simulation begins after the spacecraft has already left the Earth’s
gravitational sphere of influence, so the Earth’s gravitational influence is
not considered during the simulation.
2. Besides the acceleration provided by the spacecraft’s propulsion
system, the only perturbing environmental acceleration felt by the
spacecraft is the Sun’s gravitational acceleration.
3. The Sun’s gravity field is modeled by Newton’s law of universal
gravitation; relativistic effects are not considered.
4. In the solar gravity model, the spacecraft is considered to be a point-mass
with negligible mass. In the propulsion system model, the spacecraft is
considered to be a point-mass with positive mass; a real-(simulation)-time
combination of structural dry mass (burnout), Mb, and instantaneous
propellant mass, Mp.
5. The direction of the spacecraft’s thruster acceleration vector is designed to
always be along the spacecraft’s instantaneous heliocentric velocity
vector; it is assumed that the necessary attitude control system
(and self-contained power supply) to maintain this pointing is part of the
spacecraft’s structural dry mass, as defined in the spacecraft model. This
also applies to the launch velocity, which is assumed to be in the direction
of the initial orbital velocity vector.
17
Chapter 2: Simulating the IPM
2.1 Mission Physics
Before a simulator can be constructed, the underlying physics of the spacecraft’s motion
must be defined. For the IPM, the dynamics of the entire system is assumed to be determined
solely by the solar gravitational environment through which the spacecraft travels, and the RTG
powered Ion Engine which imparts an acceleration to the spacecraft and changes its mass during
the burn phase.
X
Y
Instantaneous Orbits
AGravity
Trajectory
Sun
AThrust
Figure 2: Accelerations felt by the spacecraft during engine-burn, with depictions of
instantaneous orbits at various moments in time.
18
The Solar Gravitational Environment
From the basic derivation of two-body motion, the gravitational acceleration of a
spacecraft whose mass, mSC, is significantly less than the mass of the primary gravitating body is
determined by the spacecraft’s three-dimensional inertial position relative to the primary at any
instant in time. This is shown in the gravitation potential, which can be expressed in terms of the
primary body’s standard gravitational parameter, µ, which is the product of the universal
gravitation constant G and the primary’s mass:
With the mass of the spacecraft, mSC , Newton’s Law of Gravitation can then be written as:
Using Newton’s second law (eqn 2.1.3), Equation 2.1.1, and Equation 2.1.2, the acceleration of
the spacecraft due to gravity (eqn 2.1.4), which points along the direction of the relative position
vector towards the primary, can be derived:
Equation 2.1.4 can be expressed using Equations 2.1.5a-c, which decompose the gravitation
acceleration into its inertial, or in this case heliocentric, x, y, and z scalar components:
19
The Ion Engine
At any moment in time, the performance of an RTG powered Ion Engine can be
characterized by a series of equations describing the engine’s design parameters. The rate at
which the engine expels the propellant from the thruster, the propellant mass flow rate, can be set
to any function of time, position, or velocity, etc, but for the purposes of this study is set to a
constant such that all of the propellant is consumed by the end of the defined burn time, tb, which
maintains a constant thrust over time:
The beam power, Pb, can be calculated from the total input power, P0, supplied by the RTG, and
the engine power efficiency, εpower:
20
If the ion exhaust beam contains multiply charged species, then the ratio of doubly charged ions
to singly charged ions can be represented by I
++
/I
+
with Ib = I
+
+ I
++
. The thrust correction
factor, α, is then defined by Equation 2.1.8:
Assuming uniform beam divergence from the thruster, along with a constant ion current density
profile accelerated by a uniform electric field, the total thrust correction factor can be written as
the product of the thrust correction factor, α, and the beam divergence correction factor for a
given average half-angle beam divergence, Θ:
For the purposes of this study, both I
++
/I
+
and Θ are assumed to be zero for all time, due to their
typically minimal impact on the engine performance. This keeps the total correction factor γ = 1,
which idealizes the beam conditions. The magnitude of the propellant exit velocity from the
thruster, Ve, is defined in terms of the beam power, Pb, as:
With the exit velocity obtained, the thruster’s instantaneous specific impulse, Isp, can be found:
In order to calculate the beam current, Ib, of the thruster, the ion flow rate of the beam must first
be computed from the propellant flow rate and the atomic mass of the propellant (in kg), which
21
in the case of this study is Xenon due to its high atomic mass:
The beam current and accelerating voltage can now be calculated:
The total efficiency of the engine, εtotal, is defined as:
which can also be written in terms of an expression for jet power, Pjet:
From this expression of total efficiency, the Ion Engine thrust can now be written as a function of
total efficiency, propellant flow rate, and total input power:
Since the propellant flow rate is constant, the total mass of the spacecraft, Mo, at any given time
(for t ≤ tb) can be calculated from the burnout mass, Mb, and the propellant flow rate:
Finally, the ∆V of the Ion Engine itself can be computed using the Tsiolkovsky Rocket Equation:
22
From Newton’s second law, the acceleration due to the engine thrust is then simply
AThrust = (FThrust / Mo) at any instant. Based on the definition of propellant flow rate, once
t = tb, all the propellant has been consumed leaving only the structural dry mass, Mb. After the
engine burn phase is complete, t > tb, the control parameters mdot, Ve, Isp, Fthrust, Ndot, Ib, and Vb
must be set to 0, which effectively “turns off” the thruster so it no longer imparts any
acceleration to the spacecraft.
2.2 Simulator Implementation
Simulator Structure and Operation
RK4
Orbit Propagator
Total
Acceleration
Gravitational
Acceleration
Thruster
Acceleration
Instantaneous
Thruster
Parameters
Main Program:
•Initial Conditions
•Escape Trajectory
Characterization
! Single Function Call
! Iterative Function Call
Figure 3: IPM Spacecraft Simulator operational block diagram
23
The simulator developed for this thesis simulates the three-dimensional motion of a
thrusting spacecraft with variable mass under the influence of a central gravitating body. Since
analytical solutions to the differential equations modeling a thrusting spacecraft with changing
total mass are difficult to obtain, numerical methods must be employed. This is done using a
fourth-order Runge-Kutta numerical integration method which solves the equations of motion to
propagate the trajectory through time. As shown in Figure 3, the simulator is comprised of six
separate routines: a Main Program that handles the initial conditions and post-propagation
analysis, a Numerical Orbit Propagator function, three Acceleration functions, and an Ion
Engine function that calculates the Instantaneous Thruster Parameters. To operate, the user
must enter the following input parameters into two of these subroutines:
Main Program User Inputs:
1. The standard gravitational parameter of the central body, µ
2. Ion Engine burn time, tb
3. Time span of orbit propagation, tspan
4. Estimate of the numerical integration time step, hestimate
5. Launch velocity relative to Earth in heliocentric coordinates, Vlaunch
6. Earth’s heliocentric orbit radius at launch, Rparking
7. Semi-major axis of Earth’s heliocentric orbit, aparking
Instantaneous Thruster Parameters User Inputs:
1. Atomic mass of propellant, Melement
2. Structural dry mass of spacecraft, Mb
3. Total mass of stored propellant, Mp
24
4. Total engine input power provided by RTG, P0
5. Power efficiency, εpower
6. Mass utilization efficiency, εmass util
7. Average half-angle divergence of the ion exhaust beam, Θ
8. Fraction of double ion current in the ion exhaust beam, I
++
/I
+
9. Propellant mass flow rate function
During the simulation, all calculations are done with values expressed in base units (m, kg, sec,
etc), and later converted to more appropriate units (km, AU, years, etc) during the post-
propagation data analysis. Upon completion, the simulator provides the following scalar outputs:
1. Did the spacecraft achieve escape velocity?
2. Flight time from engine ignition to moment of escape (yr)
3. Magnitude of heliocentric position at ECO (AU)
4. Magnitude of heliocentric velocity at ECO (AU/yr)
5. Total flight time to 200 AU (yr)
6. Magnitude of heliocentric velocity at 200 AU (AU/yr)
7. Hyperbolic excess velocity (AU/yr)
8. Final orbit eccentricity
9. Final orbit semi-major axis (AU)
10. Solar escape target angle (deg)
In addition to these values, the simulator also produces the following plots characterizing
different aspects of the simulated mission:
1. Three-dimensional heliocentric position (trajectory)
25
2. Magnitude of heliocentric velocity vs. time
3. Magnitude of heliocentric velocity vs. heliocentric distance
4. Orbit eccentricity vs. time
5. Orbit semi-major axis vs. time
6. Total spacecraft mass vs. time
7. Propellant mass flow rate vs. time
8. Engine Thrust vs. time
9. Thruster specific impulse vs. time
10. Thruster exit velocity vs. time
11. Accelerating voltage vs. time
12. Beam current vs. time
The Main Program organizes the initial conditions and sends them to the Numerical Orbit
Propagator through a single function call. The Orbit Propagator uses the fourth order Runge-
Kutta process which takes the spacecraft under the influence of a single total acceleration vector
at one state (position and velocity), and predicts the next state at the defined time step later. This
requires iterative function calls to the Total Acceleration function which in-turn require the
Gravitational and Thruster Acceleration functions to determine the instantaneous accelerations
felt by the spacecraft. The Total Acceleration function adds the two separate accelerations
together to create a single total acceleration vector, which is then passed back to the propagator.
Each call to the Thruster Acceleration function requires an additional call to the Instantaneous
Thruster Parameters function to calculate the thrust and total mass at the current state. These
values are then used by the Thruster Acceleration function to calculate the magnitude of thruster
26
acceleration vector. (It should be noted that Thruster Acceleration function assumes the thrust
vector to be along the direction of the instantaneous heliocentric velocity vector). Once the
numerical integration process is complete, the orbit propagator sends the state and parameter
matrices back to the Main Program, where the data is processed. The algorithms used by each
function are described in detail in subsequent subsections of this chapter.
The time it takes to run a complete simulation varies greatly depending on the desired
integration time span, but for a typical 7 year propagation with an integration time step of ~100
sec/step, the real elapsed time of the simulation is on the order of 13 minutes. A 17 year
propagation using the same step takes roughly 50 minutes to complete.
Initial Conditions & Escape Trajectory Characterization
While the user provides most of the initial conditions required to initiate the orbit
propagation, the velocity of the spacecraft in its initial orbit must be calculated as well. The
relation known as the Vis-Viva Equation (eqn 2.2.1) is used to calculate the magnitude of the
velocity of the Earth (eqn 2.2.2) based on the user specified Earth orbital parameters (launching
from a circular orbit, or an elliptical orbit from perihelion or aphelion):
27
The initial three-dimensional state in the heliocentric inertial reference frame combines Vlaunch
and VEarth (eqn 2.2.3) by the assuming that the Earth’s position is on the positive heliocentric
inertial X-axis, and that both VEarth and Vlaunch are in the heliocentric +Y direction (Figure 4).
X
Y
Sun
Earth Orbit
VEarth
Vlaunch
Escape
Trajectory
Target Angle, ξ
Direction
of Escape
Figure 4: The Heliocentric Inertial Reference Frame and Initial Conditions
28
The program then calculates the integer number of steps necessary for the Runge-Kutta
integration process based on the user supplied estimate to ensure the number of integration steps
best covers the desired propagation time:
The program then calls the Numerical Orbit Integrator, which requires inputs of tspan, Nsteps,
and the y0 state, and provides the following outputs:
1. t : (Nsteps + 1)x1 array of time values at each integration step in units of (sec)
2. y : (Nsteps + 1)x6 state matrix of position and velocity vectors in base units
3. tescape index : Integer variable which indicates at which integration step the
spacecraft achieved Solar escape velocity (if at all).
4. mdot : (Nsteps + 1)x1 array of thruster mass flow rate values in .
5. Ve : (Nsteps + 1)x1 array of thruster exit velocity values in units of .
6. Isp : (Nsteps + 1)x1 array of thruster specific impulse values in units of (sec)
7. Fthrust : (Nsteps + 1)x1 array of thruster force in units of (N).
8. Mo : (Nsteps + 1)x1 array of spacecraft total mass values in units of (kg).
9. Ib : (Nsteps + 1)x1 array of thruster beam current in units of (A).
10. Vb : (Nsteps + 1)x1 array of thruster accelerating voltage values in units of (V).
Once propagation is complete, if the spacecraft has achieved escape velocity during the
simulated engine burn, its escape trajectory is characterized and the total flight time to 200 AU
computed from just a few simple parameters. Plugging the magnitude values of R and V taken
just after ECO into the Vis-Via equation (eqn 2.2.1), the definition of orbital angular momentum
29
(eqn 2.2.5), and an expression for the hyperbolic semi-parameter (eqn 2.2.6), the final hyperbolic
eccentricity, efinal, can be found.
The hyperbolic excess velocity, or V∞, can also be computed from the Vis-Viva equation by
setting R = ∞:
The total flight time from engine ignition to 200 AU (FT200 AU) can be described by two separate
time segments: the time from ignition to burn-out (tb), and the time from burn-out to 200 AU
(tECO_200AU). However, since the instantaneous orbit changes between each time step during the
engine burn, the final hyperbolic orbit’s perihelion, and corresponding epoch time, does not
coincide with the spacecraft’s initial position and epoch. This means that tECO_200AU is comprised
of the time from the theoretical hyperbolic epoch to ECO, subtracted from the time from
theoretical hyperbolic epoch to 200 AU. This yields the following expression for the total flight
time from ignition to 200 AU:
30
To calculate th-epoch to 200 AU, the spacecraft’s true anomaly at 200 AU (θ200 AU) and
hyperbolic anomaly at 200 AU (H200 AU) must be calculated from the orbit equation (eqn 2.2.9),
the semi-parameter definition (eqn 2.2.6), and the true anomaly/hyperbolic anomaly relation (eqn
2.2.10). th-epoch to 200 AU can then be calculated using the Hyperbolic Kepler Equation (eqn 2.2.11).
Similarly, th-epoch to ECO can be found by first calculating the true anomaly at ECO with Equation
2.2.9 and converting to hyperbolic anomaly at ECO with Equation 2.2.10, and solving the
Hyperbolic Kepler Equation (eqn 2.2.11) once again for ∆t, which in this case is th-epoch to ECO.
The total flight time from ignition to 200 AU can now be computed using Equation 2.2.8.
31
To calculate the Target Angle, ξ, the angle measured counter-clockwise from the
heliocentric X-axis to the V∞ direction (Figure 5), the direction of the final hyperbolic orbit’s
eccentricity vector, which points towards the final hyperbolic perihelion, must first be computed
using the following heliocentric vector equation (eqn 2.2.12):
X
Y
[e]
unit
ψ
Direction
of Escape
ξ
Figure 5: Hyperbolic Orbit Geometry for -[ey]unit ≤ 0
32
The half-angle between the hyperbolic asymptotes, ø, is defined by equation 2.2.13:
The angle between the -[e]unit direction and the heliocentric X axis, ψ, is then calculated by:
The target angle, ξ, can now be calculated based on the condition of -[ey]unit:
To keep the domain of ξ between 0 ≤ ξ ≤ 2π, the following correction is made:
Orbit Propagator: The Fourth Order Runge-Kutta Method
To solve the equations of motion of the thrusting spacecraft with variable mass, a
numerical integration method is employed, known as the fourth order Runge-Kutta process,
tailored specifically to the orbit problem. This method takes an initial orbital state, a desired
integration time span, and the number of integration steps, and approximates the state at each
time step using the using four pairs of discrete vector equations. In the simulator, the various
thruster parameters are also recorded within this function at each state via a function call to the
Instantaneous Thruster Parameters function.
33
The function begins by calculating the actual integration step size, h, from the calculated
number of steps and total propagation time (eqn 2.2.16):
It should be noted that all state and thruster parameter matrices are defined to have (Nsteps + 1)
rows since Matlab counts the first array index as “1”, not “0”. The initial position and velocity is
parsed from the initial state vector, and is used with the initial time t = 0 in a function call to the
Total Acceleration function which calculate the initial total acceleration (1x3 vector). The
initial position, velocity, and acceleration vectors are then assigned to the variables “Rprev”,
“Vprev”, and “Aprev”, respectively. For the first iteration, these “previous” values are the initial
parameters in the Runge-Kutta process, which runs iteratively from the second time step (step
count j = 2) to the final step (j = Nsteps + 1) by steps of h seconds. Each step, the time array “t” is
updated by adding one step size “h” to the previous time value (eqn 2.2.17):
The Fourth Order Runge-Kutta method for numerical orbit propagation uses four pairs of
1x3 vector variables, kn (eqn 2.2.18a-h), which require the spacecraft’s total instantaneous
acceleration vector in order to numerically approximate the next state vector from the previous
state (eqn 2.2.19):
34
The newly calculated position and velocity vectors are then recorded in the state matrix y. If
Equation 2.2.20 is satisfied, and the spacecraft has not escaped in a previous iteration, then the
spacecraft has achieved escape velocity at the current step. The time of escape and
corresponding time array index value are recorded.
35
After Rnext and Vnext are computed, the Instantaneous Thruster Parameters function is
called with an input of the current time value to calculate the set of thruster parameters at the
current state. These values are recorded in the current index of the thruster parameter matrices.
Finally, the Rprev and Vprev variables are set to the newly calculated Rnext and Vnext, and Aprev is
recalculated using the Total Acceleration function using Rnext and Vnext, and the current time
value t(j). The Runge-Kutta process is then repeated at the next step, and so on; iterations
continue until the step counter j = (Nsteps + 1).
Total Acceleration Calculation
To determine the spacecraft’s total acceleration vector at any given time, the simulator’s
Total Acceleration function combines the individual accelerations governed by the solar
gravitational model and the Ion Engine model, as described in section 2.1. The total acceleration
felt by the spacecraft is calculated by taking a sum of the (n) number of separate acceleration
vectors acting on the spacecraft (eqn 2.2.21), which in the case of this study is just two:
As derived in section 2.1, the spacecraft’s gravitational acceleration is a function of heliocentric
position, R, whose scalar accelerations can be represented as the vector AGrav (eqn 2.5.6),
calculated in the simulator by an independent function:
36
Similarly, the thruster acceleration vector can be found by first determining the direction of the
thrust vector, which in the case of this study is assumed to be along the direction of the
heliocentric velocity unit vector:
Using Newton’s second law (eqn 2.1.3), the acceleration due to the thruster is then calculated as
the instantaneous thruster force over the instantaneous total mass, in the direction of V
unit
:
2.3 Software Validation
The orbit propagator was validated by analytically calculating two reference orbits with
the Ion Engine switched off, and comparing the spacecraft’s positions at various true anomalies
(θ) in the orbit with simulated trajectories based on the same initial conditions. A circular orbit
about the Sun at 1 AU and an elliptical 1x8 AU orbit were used in the test. To find the analytical
positions at the corresponding times from the numerical solutions, the Elliptical Kepler Equation
was solved using the Newton-Raphson iterative process. For the 1 AU circular case, the
difference in position between the analytical and numerical solutions near specific true anomaly
angles are as follows:
37
For the 1 AU circular case:
∆Rθ = 90º = 17.030 m
∆Rθ = 180º = 34.000 m
∆Rθ = 270º = 17.060 m
∆Rθ = 360º = 0.870 m
For the 1x8 AU elliptical case:
∆Rθ = 90º = 4.210 m
∆Rθ = 180º = 8.901 m
∆Rθ = 270º = 4.108 m
∆Rθ = 360º = 0.603 m
This implies that the total relative error is on the order of 10
-10
, which follows the expected result
that the fourth order Runge-Kutta total accumulated error should be on the order of O(h
4
).
38
Chapter 3: Simulation Methods
3.1 Simulator Settings
The trajectory simulator developed for this thesis gives the user the freedom to tailor the
Ion Engine model and initial conditions to their specific simulation needs. For the purposes of
this study, the simulations were conducted using the following settings:
1. hestimate = 100 (sec/step) for all simulations
2. Spacecraft structural dry-mass (burnout), Mb = 790 (kg)
3. Initial propellant mass, Mp = 440 (kg)
4. Xenon Propellant, Mxe
5. Propellant mass flow rate = a constant for all simulations, as determined
by Equation 2.1.6
6. εpower = 0.6
7. εmass util = 0.9
8. I
++
/I
+
= 0
9. Θ = 0º
While no perturbations are currently included in the simulator, the total acceleration function is
structured in such a way that additional accelerations could easily be incorporated into the
simulator architecture.
39
3.2 Input Parameter Sets and Scenarios
Each simulation run, or “scenario”, used a different combination of control parameters:
1. Rlaunch set to RCirc-Earth = 1 AU
2. Vlaunch set to either 0, 5, 7, 10, or 12 km/sec
3. P0 set to either 1000, 1500, or 2000 W
4. tb set to either 5, 10, or 15 years
In order to accommodate all possible combinations, each scenario was categorized first by the
burn time, tb, then by the simulation “set” determined by input power, P0, and then by the launch
velocity, Vlaunch (Figure 6):
Scenario: C2
Scenario: C3
Scenario: C4
Scenario: C1
Burn Time:
t b = 5 yr
Set A:
P 0 = 1000 W
Set B:
P 0 = 1500 W
Set C:
P 0 = 2000 W
V launch = 0 km/sec
V launch = 5 km/sec
V launch = 7 km/sec
V launch = 10 km/sec
V launch = 12 km/sec
Burn Time:
t b = 10 yr
Burn Time:
t b = 15 yr
Scenario: C5
Set D: P 0 = 1000 W
Set E: P 0 = 1500 W
Set F: P 0 = 2000 W
Set G: P 0 = 1000 W
Set H: P 0 = 1500 W
Set I: P 0 = 2000 W
Figure 6: Simulation hierarchy, illustrating simulation Set C.
40
Each simulation set is unique, meaning that sets A through C correspond exclusively with
tb = 5 yr, sets D through F correspond with tb = 10 yr, and sets G through I correspond with
tb = 15 yr. This architecture calls for a total of 45 individual simulations, whose results are
detailed in subsequent chapters and Appendix A. While these scenarios only represent missions
under the initial circular Earth orbit assumption, additional scenarios were simulated with
Rlaunch = RP-Earth & RA-Earth (with set designations of Ap - Cp for perihelion and Aa - Ba for
aphelion) in order to compare the perihelion and aphelion results with the circular case
(Appendix B). These comparison sets were run using only tb = 15 yr since longer burn times
were found to have a larger number of scenarios that achieved escape velocity; since only
scenarios for which the spacecraft achieved escape velocity were considered for post-simulation
analysis, using the highest burn time provided the most usable data points.
41
Chapter 4: Simulation Results and Analysis
4.1 Example Simulation Outputs
Output Plots from Scenario C4
To demonstrate a typical simulation run, let us examine the results from mission scenario
C4, simulated with a 5 year burn time, an RTG supplied power of P0 = 2000 W, and launched
with a velocity of Vlaunch = 10 km/sec from an initial circular Earth orbit:
Figure 7: Scenario C4: Trajectory
42
As seen in the output trajectory, the combined initial conditions and input power provided
sufficient energy into the system to cause the spacecraft to escape, and to do so in less than one
complete revolution around the Sun. The spacecraft achieved escape velocity in a little less than
half of the burn time, which allowed it to spend more of the burn phase increasing its hyperbolic
excess velocity rather than attempting to escape. In other words, the more time the spacecraft
can spend burning post-escape, the faster it can reach 200 AU (as it will be shown, the launch
velocity has a significant impact on when escape occurs). The calculated total flight time to 200
AU is just under 60 years, which while comparable to the total flight time of the Voyager 1 probe
to 200 AU (~56.4 years), it is triple the defined “reasonable time frame” of ~20 years for such a
fly-out mission.
Figure 8: Scenario C4: Heliocentric Speed vs. Position
43
As a companion plot to Figure 7, Figure 8 illustrates that during the pre-escape burn
phase, the gravitational dynamics of the system naturally decrease the spacecraft’s velocity as it
moves farther away from the Sun, while the Ion Engine simultaneously attempts increase its
velocity. This is seen in Figure 7 where the spacecraft’s actual velocity is initially very close to
the required escape velocity at that position, but then diverges from this boundary during the
burn. Eventually the acceleration due to the Ion Engine overtakes this gravitational effect, and
the spacecraft begins to approach, and finally break through, the escape barrier. This represents a
source of inefficiency in all low-thrust propulsion architectures, known as the maneuver
inefficiency, since a portion of the thrust is used to counter the slowing gravitational effect. This
phenomena will be re-examined in the eccentricity plot (Figure 10) and the scenario E3 results.
Figure 8 also demonstrates that after engine cut-off the spacecraft is solely under the influence of
gravitational acceleration, meaning its heliocentric velocity will continue to decrease as it
approaches the hyperbolic excess velocity asymptote characterized by its final orbit.
44
Figure 9: Scenario C4: Orbit Semi-major Axis vs. Time
As theory would predict, the spacecraft’s semi-major axis increases following an inverse
relation, with a discontinuity at the point of escape. This is because the boundary between a
captured and escaped orbit is a parabolic orbit, which is technically a closed orbit with
aparabolic = ∞. Once through this point, the semi-major axis of any hyperbolic orbit is negative,
with “larger”, or more accurately put “higher energy”, orbital semi-major axes closer to 0 (from
the -∞ direction).
45
Figure 10: Scenario C4: Orbit Eccentricity vs. Time
The spacecraft’s motion once again acts as expected, with its orbital eccentricity
increasing through the point of escape, and remaining constant after ECO. While it is difficult to
see in this specific scenario’s eccentricity plot, Figure 10 does provide another piece of evidence
to support the maneuver inefficiency demonstrated in Figure 8. In simple elliptical orbit
mechanics, theory predicts that if a ∆V maneuver is done at periapsis, the resulting effect is an
increase, or “raise”, of the radius of apoapsis, and vice-versa for a burn at apoapsis. If such a
maneuver is done near, but not on, either of these locations, then same effect is seen to a lesser
extent, but with the addition of an incurred rotation of the entire orbit in-plane.
46
In the case of the simulations, the moment the engine ignited at the start of the
simulation, the spacecraft’s orbit is immediately no longer circular, and while it’s orbit size and
shape continuously change, the orbit remains in the elliptical regime prior to escape, simply
becoming more and more elliptical. But similar to the simple theoretical case, during the pre-
escape burn phase, as the spacecraft travels from it’s initial position, or relative “perihelion”,
towards its relative “aphelion”, the thrust causes a continuous rotation of the trajectory in plane,
until it reaches its relative aphelion (if it does so at all). The implication of this is that if the
spacecraft does not achieve escape velocity prior to a approaching it’s relative aphelion, upon
reaching it the spacecraft’s thrust will cease raising aphelion, as it had been since ignition, and
begin to raise perihelion instead. This is the source of the inherent inefficiency of a low-thrust
mission, since ideally the spacecraft would perform an impulse burn at ignition, or relative
perihelion, in order to use all of the propellant mass to exclusively increase aphelion past the
point of escape instantaneously. While this effect is negligible for scenario C4 , as seen by the
slight concavity change of eccentricity during the pre-escape burn phase and the spacecraft never
reaching its relative aphelion, it has a prominent impact on scenarios with longer burn times and
low energy initial conditions where complete revolutions of the Sun occur prior to escape. This
will be explored in scenario E3.
47
Figure 11: Set C: Spacecraft Total Mass vs. Time
As expected, the total mass of the spacecraft decreases in a linear fashion due to the constant
propellant mass flow rate, and stops decreasing once all of the propellant has been depleted.
Example Fringe Scenario: Plots from Scenario E3
The results of this study have revealed that longer burn times allow scenarios with
previously unsuccessful lower-energy initial conditions to escape. This is because despite a
lower propellant flow rate (eqn 2.1.6), the Ion Engine has more time to accelerate the spacecraft
to escape velocity. The draw back to this is that these “fringe” scenarios that just barely manage
to escape often have total flight times to 200 AU significantly larger than the the ~20 year goal.
48
That being said, the results of such scenarios still merit further analysis, so the following are
highlights from mission scenario E3, simulated with a 10 year burn time, RTG supplied power of
P0 = 1500 W, launched with a velocity of Vlaunch = 7 km/sec from an initial circular Earth orbit:
Figure 12: Scenario E3: Trajectory
While the spacecraft in scenario E3 was able to escape, several key differences can be seen
between the scenario C4 and E3 trajectories. The most obvious difference is that the spacecraft
made almost one and a half revolutions around the Sun prior to ECO. This meant that less time
was spent thrusting in the general direction of the hyperbolic excess velocity, which resulted in a
49
lower ECO velocity and therefore a longer flight time to 200 AU. More importantly, the
complete revolution about the Sun meant that the maneuver inefficiency mentioned in the
scenario C4 discussion was more prominent. As seen in the scenario plots, the spacecraft
reached its “relative aphelion” Rrel-ap = ~4 AU at roughly 225º measured counter-clockwise from
the X-axis. The orbital radius then decreased towards a new “relative perihelion” Rrel-peri = ~2.2
AU, before increasing through the point of escape. Over the course of this relative orbit phase,
the engine’s thrust alternated between raising aphelion and raising perihelion, resulting in distinct
increases and decreases of eccentricity (Figure 15). So while the size of the orbit continued to
increase throughout, it was not dynamically efficient since only aphelion needed to be raised.
Figure 13: Scenario E3: Heliocentric Speed vs. Position
50
Figure 14 (top): Scenario E3: Heliocentric Speed vs. Time
Figure 15 (bottom): Scenario E3: Orbit Eccentricity vs. Time
51
4.2 Orbit and Mission Architecture Analysis
Once all 45 initial circular Earth orbit scenarios were simulated, the results were
organized and plotted to illustrate the relationships between the input control parameters and the
characteristics of the resulting mission (NOTE: only scenarios where the spacecraft achieved
escape velocity were included in this analysis). Design characteristics of the Voyager 1 mission,
such as its total flight time to 200 AU from launch, are included for reference wherever relevant:
Total Flight Time to 200 AU
To examine the primary mission design parameter, the total flight times to 200 AU were
plotted based on burn time with all corresponding power sets included in each:
Figure 16: Total Flight Time to 200 AU vs. Launch Velocity for burn time of 5 years
52
Figure 17 (top): Flight Time to 200 AU vs. Launch Velocity for burn time of 10 years
Figure 18 (bottom): Flight Time to 200 AU vs. Launch Velocity for burn time of 15 years
53
As shown in Figure 17, a 5 year burn time only allowed scenarios with high energy
launch velocities, Vlaunch, greater than 10 km/sec to escape. Among those, only those launched
with Vlaunch = 12 km/sec reached 200 AU faster than V oyager 1. However, even these scenarios
are more than double the desired ~20 year fly-out time. For burn times of 10 and 15 years
(Figures 17 & 18), many more scenarios were able to escape, but only a fraction of those
reached 200 AU before Voyager 1. The highest energy scenarios for both 10 and 15 year burn
times were closer to desired fly-out time, with the fastest being ~35.5 years with a 15 burn time,
but this could still be considered out of the desired range. Extrapolating from these figures, it
may be possible to achieve the ~20 year fly-out with either a larger power input P0 or launch
velocity Vlaunch, or more likely a combination of both. While increasing the burn time would also
reduce the fly-out time, the resulting ECO distances will likely exceed the desired RECO limit of
RECO < ~80 AU. It is apparent that while increasing the input power certainly decreases the fly-
out time, as the launch velocity increases the benefit of higher power levels becomes less
significant. On average, there is roughly a ~30 year fly-out time difference between power levels
at low launch velocities, but less than ~10 year difference between power levels at high launch
velocities. This trade off can be used as a powerful design factor; a launch vehicle with a high
C3 value could be paired with a lower power engine, or a launch vehicle with a lower C3 could be
paired with a higher powered engine to achieve comparable fly-out times.
Heliocentric Distance at ECO
The second main design parameter for direct solar escape trajectories is the heliocentric
position of engine cut-off. Since the exhaust from the Ion Engine can interfere with the
54
spacecraft’s suite of scientific instruments and sensors, it is imperative that ECO occur well
before the science phase of the mission is planned to begin in order to allow the exhaust
particulates time to dissipate. Since plasma interactions between the Sun and the Interstellar
Medium may occur as near as 70 - 80 AU from the Sun, depending on solar activity and the
direction of escape, this range can be used as a secondary restriction to assess the simulated
missions.
Figure 19: Heliocentric Position @ ECO vs. Launch Velocity, 5 year burn time
55
Figure 20 (top): Heliocentric Position @ ECO vs. Launch Velocity, 10 year burn time
Figure 21 (bottom): Heliocentric Distance @ ECO vs. Launch Velocity, 15 year burn time
56
As shown in Figures 19, 20 & 21, all simulated trajectories achieved an ECO distance
within the required boundary range. However, by examining the ECO distance for missions with
P0 = 2000 W launched with Vlaunch = 12 km/sec, it appears that every 5 year increase of burn time
increases the ECO distance by about ~22 AU. It can therefore be surmised that if the burn time
were increased to 20 years for the highest energy initial conditions case, the ECO distance would
be ~92 AU, which is well above the desired limit. So while increasing the burn time would
likely decrease the total flight time to 200 AU, the corresponding ECO distances for the high
energy initial conditions cases will likely be unacceptable.
Orbit Characterization
In order to characterize the typical final orbit of direct solar escape trajectories, the
following plots were produced exclusively for the 15 year burn time cases since they had the
most number of successful escape trajectories:
57
Figure 22: Heliocentric Excess Velocity vs. Launch Velocity for burn time of 15 years
As expected, scenarios with higher input power levels and launch velocities resulted in larger
excess velocities upon escape. Figure 22 also shows that the scenarios with excess velocities
well above Voyager 1’s were also the scenarios that reached 200 AU well before it (Figure 18).
58
Figure 23: Solar Escape Target Angle vs. Launch Velocity for burn time of 15 years
The comparison of solar escape target angles yields two interesting results: first, it reveals that as
the launch velocity increases, it becomes the dominant factor in determining the target angle;
variations in the input power have a lessening effect as the energy provided by launch increases.
Second, the two scenarios that appear as outliers, with angles far below the higher input power
cases for the same given launch velocity, represent the “fringe” scenarios discussed earlier.
Because these cases just barely managed to achieve escape velocity, they took longer to escape
compared to the similar higher input power cases and thus revolved further about the Sun. Since
the domain of the target angle is 0 < 360, this resulted in comparatively low angles.
59
Figure 24: Flight Time to Escape vs. Launch Velocity for burn time of 15 years
Figure 24 demonstrates that the launch velocity is the dominant factor in the time it takes for the
spacecraft to achieve escape velocity. Since the required escape velocity from a 1 AU circular
orbit is roughly Vreq-esc = 42.1 km/sec, the minimum launch velocity required to launch directly
into an escape trajectory from this orbit is Vlaunch-min = ~12.337 km/sec. For launch vehicle
selection, this corresponds to a minimum C3 value of:
Such launches would allow for shorter burn times since the engine would be thrusting during the
more efficient post-escape phase, which would help maintain acceptable ECO distances.
60
Figure 25: Final Orbit Semi-major Axis vs. Launch Velocity for burn time of 15 years
Figure 25 confirms that higher energy initial conditions result in “larger”, or more energetic,
hyperbolic escape orbits. (NOTE: While V oyage 1’s semi-major axis appears to be zero, it is
actually -0.0032 AU, signifying a high energy hyperbolic trajectory. This was achieved through
its various planetary gravity assists during its primary mission phase and transition into its
extended mission for Solar System fly-out.)
61
Figure 26: Final Orbit Eccentricity vs. Launch Velocity for burn time of 15 years
As the initial conditions become more energetic, the resulting escape trajectory becomes more
eccentric. The slight decrease in eccentricities for Vlaunch = 12 km/sec can likely be attributed to
the faster escape times, which mean more of the thrust time is spent in the direction of final
escape, locking the spacecraft into a final eccentricity sooner than in the lower energy cases.
4.3 Ion Engine Model Analysis
To showcase the engine model for design parameters corresponding to a wide range of
burn times, the model equations were analytically plotted against burn time:
62
Figure 27: Ion Engine Beam Current vs. Burn Time (analytic)
Equation:
For a given type and amount of propellant, assuming constant propellant flow rate, the
corresponding beam current is solely determined by the burn time. Figure 27 shows that to keep
the beam current below 1 A, burn times of roughly 10 years or more are desirable.
63
Figure 28: Ion Engine Accelerating Voltage vs. Burn Time (analytic)
Equation:
As burn time increases, the difference in accelerating voltage between adjacent input
power levels diverges linearly. This means lower burn times correspond to lower voltages for
any given input power level.
64
Figure 29: Propellant Mass Flow Rate vs. Burn Time (analytic)
Equation:
The propellant mass flow rate function is designed to ensure that all propellant is
completely exhausted at ECO. This means the mass flow rate is only determined by the total
amount of propellant and the burn time. Advanced functions could easily be implemented in
place of this, such as making the propellant flow rate a function of heliocentric position to
balance the gravitational and thruster accelerations over time in order to minimize the maneuver
inefficiency.
65
Figure 30: Ion Engine Thruster Force vs. Burn Time (analytic)
Equation:
Since the thruster force drops off as a square root of an inverse, larger burn times
correspond with a smaller difference between adjacent input power levels. This means
implementing a higher input power level will result in a greater thrust increase at lower burn
times. Using the equation above and the definition of impulse, a useful expression can be
derived that relates the the thruster force at one burn time, tb, with the thruster force at another
66
burn time, (n·tb), given the same input power. Since FThrust is constant over the course of any
given burn time, the total impulse due to the thruster is:
So for two different burn times:
Applying the equation for FThrust and dividing the above equations, we get:
Finally, plugging in the original impulse definitions yields the following simple expression
relating the thruster force for one burn time to the thruster force for another burn time:
67
Figure 31: Ion Engine Beam Exit Velocity vs. Burn Time (analytic)
Equation:
68
Figure 32: Ion Engine Isp vs. Burn Time (analytic)
Equation:
For low-thrust electric propulsion systems, the Isp for a given propellant mass is only
limited by the thruster’s input power and propellant flow rate, which in the case of this study is
only dependent on burn time. Longer burn times grant a greater increase in Isp between adjacent
power levels.
69
Figure 33: Ion Engine ∆V vs. Burn Time (analytic)
Equation:
The above ∆V’s were calculated using the same masses which were held constant
throughout all of the simulations:
Mp = 440 [kg]
Mb = 790 [kg]
Mo = 1230 [kg]
70
4.4 Initial Position Comparison: Circular Earth Orbit, Perihelion, Aphelion
While all previous results were based on launches from an assumed circular heliocentric
Earth orbit, scenarios were also run for launches from Earth perihelion and aphelion with 15 year
burn times to asses the flight time advantage, or disadvantage, of launching from these locations:
Perihelion Launch
Figure 34: Difference in Total Flight Times between circular and perihelion launches
In Figure 34, a ∆FT > 0 implies that the circular version of that scenario took longer to
reach 200 AU than the perihelion case, and conversely a ∆FT < 0 implies that the circular version
took less time to reach 200 AU than the perihelion case. Instinctually all perihelion cases should
71
take less time to reach 200 AU than the circular cases since the initial Earth aphelion distance is
naturally larger than 1 AU, and while this holds true for the majority of the simulated scenarios,
two pairs of scenarios stand contrary to this with the circular cases having faster fly-out times
than the perihelion cases. These outliers turn out to once again be the “fringe” scenarios
discussed in previous sections; because they only barely managed to escape, scenarios launched
from perihelion escaped later, or closer to the ECO time, than those from the circular cases. This
means the circular fringe cases had longer post-escape burn phases, resulting in larger velocities
at ECO and shorter fly-out times than the perihelion fringe cases.
Aphelion Launch
Figure 35: Difference in Total Flight Time between circular and aphelion launches
72
Similar to the perihelion comparisons, in Figure 34 a ∆FT > 0 implies that the circular
version of that scenario took longer to reach 200 AU than the aphelion case, and conversely a
∆FT < 0 implies that the circular version took less time to reach 200 AU than the aphelion case.
Instinctually all aphelion cases should take more time to reach 200 AU than the circular cases
since the initial position is greater than 1 AU (meaning the initial burn phase raises perihelion,
which is less than 1 AU), and while this holds true for the majority of the simulated scenarios,
the same two fringe scenarios stand out with the circular cases having longer fly-out times than
the aphelion cases. This is caused by the same mechanism discussed in the perihelion
comparisons; fringe scenarios launched from aphelion escape faster than in the circular cases,
resulting in more time spent burning in the post-escape burn phase and faster fly-outs. When
combined with the results from the perihelion comparisons, this implies that for all scenarios that
just manage to escape, initial positions farther from the Sun will result in faster fly-out times than
those closer to the Sun, a seemingly backwards conclusion.
4.5 Conclusions
The following is a summery of the conclusions found by this study for the given
assumptions and simulator settings described in Chapter 3. The results have been broken down
into the categories of the mission architecture, Ion Engine performance, and launch position
comparisons:
73
Mission Architecture for Circular Earth Orbit Launches
1. The faster the spacecraft can escape, and therefore spend more time burning
post-escape, the faster it can reach 200 AU.
2. The gravitational dynamics of the system naturally decrease the spacecraft’s
velocity as it moves farther away from the Sun, which causes an inefficiency for
all low-thrust spacecraft. During the pre-escape burn phase, if the spacecraft does
not achieve escape velocity prior to approaching it’s relative aphelion, the
spacecraft’s thrust will cease raising aphelion, and instead begin to raise
perihelion. This means scenarios with longer burn times and low energy initial
conditions which reach their relative aphelion at least once (“fringe” scenarios)
are more susceptible to this inefficiency.
3. Longer burn times allow scenarios with previously unsuccessful lower-energy
initial conditions to escape. The draw back to this is that these “fringe” scenarios
just barely manage to escape, resulting in total flight times to 200 AU
significantly longer than the the ~20 year goal.
4. Burn times of 5 years only resulted scenarios with high energy launch velocities,
Vlaunch ≥ 10 km/sec, to escape. Only those launched with Vlaunch = 12 km/sec
reached 200 AU faster than Voyager 1. However, these scenarios are more than
double the desired ~20 year fly-out time. For burn times of 10 and 15 years,
many more scenarios were able to escape, but only a fraction of those reached
200 AU before Voyager 1. The highest energy scenarios for both 10 and 15 year
74
burn times were closer to desired fly-out time, with the fastest being ~35.5 years
with a 15 burn time, which could still be considered too long.
5. All simulated trajectories have an ECO distance within the required ~70 - 80 AU
boundary. However, by examining the ECO distances for trajectories with
P0 = 2000 W launched with Vlaunch = 12 km/sec, it appears that for every 5 year
increase of burn time, the ECO distance increases by about ~22 AU. It can
therefore be surmised that if the burn time were increased to 20 years for the
highest energy initial condition case, the ECO distance would be ~92 AU, which
is well above the desired limit. So while increasing the burn time would likely
decrease the total fly-out time to within the desired time frame, the corresponding
ECO distances for the high energy initial conditions cases would likely be
unacceptable.
6. As the launch velocity increases, it becomes the dominant factor in determining
the escape target angle; variations in the input power have less of an impact on
the target angle as the launch velocity increases.
7. The launch velocity is the dominant contributor to the time it takes to achieve
escape velocity. The minimum launch velocity required to launch directly into an
escape trajectory from a circular 1 AU orbit is Vlaunch-min = ~12.337 km/sec. For
launch vehicle selection, this corresponds to a minimum C3 value of:
75
Such launches would allow for shorter burn times since the engine would be
thrusting during the more efficient post-escape burn phase, which would aid in
maintaining acceptable ECO distances.
8. Higher energy initial conditions result in “larger”, or more energetic, hyperbolic
escape trajectories.
9. The final orbit eccentricity decreases from Vlaunch = 10 km/sec to
Vlaunch = 12 km/sec for high energy launches. This can likely be attributed to the
faster escape times, which mean more time is spent thrusting in the direction of
final escape, which locks the spacecraft into a final eccentricity sooner than in the
lower energy cases.
Ion Engine Performance
1. To minimize the beam current, increase the engine burn time
(to roughly ~10 years or more for the greatest effect)
2. To minimize the accelerating voltage for any given input power level, decrease
the burn time.
3. Higher input power levels will result in a greater increase of thrust for lower
engine burn times.
4. For low-thrust electric propulsion systems, the Isp for a given propellant mass is
only limited by the thruster’s input power and propellant flow rate (which in the
case of this study is only dependent on burn time). Longer burn times cause a
greater increase in Isp between adjacent power levels.
76
Initial Position Comparison
1. Most perihelion launches take less time to reach 200 AU than circular launches,
but circular “fringe” scenarios that just manage to escape have longer post-escape
burn phases, resulting in larger velocities at ECO and shorter fly-out times than
the corresponding perihelion fringe cases. The opposite is true for fringe launches
from aphelion. This means that fringe scenarios with initial positions farther from
the Sun in Earth’s eccentric orbit will result in faster fly-out times than those
closer to the Sun, a seemingly backwards conclusion.
4.6 Suggestions for Future Work
While this study has explored in detail a great number of direct ballistic escape trajectory
mission architectures, there are many scenarios that must still be run to fully simulate all possible
Interstellar Precursor Mission designs. Without having to develop additional code or simulation
tools, simulations can be run for many more input power levels and launch velocities to develop
a higher resolution of cross-scenario comparisons. Other propellant types and initial orbital
conditions can be tested, including any realistic inclination changes required to escape in the
desired direction. The most significant work that can be done with the existing simulator is to
adjust the thrust profile via the propellant mass flow rate function so it is no longer a constant.
The benefit of a control function based on heliocentric position (or possibly flight time or
instantaneous total mass) is that it would allow FThrust to be more efficiently tailored to the
77
instantaneous gravitational acceleration, meaning a higher thrust closer to the Sun, and lower
thrust farther away, the effect of which could significantly reduce flight time to 200 AU.
Since the current simulator is a simplification of the true dynamics of the IPM, namely
the lack of any perturbations, higher fidelity simulations could be run simply by adding
additional acceleration functions which more realistically model the interactions between the
spacecraft and the space environment. Accelerations due to solar radiation pressure, spacecraft
heating, plasma interactions, or even planetary gravity perturbations could be included in the
Total Acceleration function when the instantaneous acceleration vector of the spacecraft is
calculated. A more robust RTG model could also be constructed, enabling a more realistic
engine input power, P0, that decreases over time as the radioactive RTG fuel decays. These
relatively minor additions to the simulator could have a significant impact on the simulation
results, and would better represent the true dynamics such a mission.
The most significant work that must be conducted following this study is the inclusion of
gravity assist maneuvers. While these are by no means required for a spacecraft to reach 200 AU
in a reasonable time, as this research suggests, their contribution to the mission can dramatically
increase the spacecraft’s escape speed and decrease the fly-out time to 200 AU, likely saving
significant technical and financial cost since lower energy launch vehicles and spacecraft power
systems can be employed. Following the trail blazed by Voyager 1, these gravity assist
maneuvers would likely involve single or multiple fly-by’s of Jupiter, Saturn, or a combination
of the outer planets. This would require significant effort and alterations to the simulator
developed for this thesis, but such research would yield incredible results of great significance to
the Interstellar Precursor Mission design effort.
78
References
1. D. M. Goebel, I. Katz, “Fundamentals of Electric Propulsion: Ion and Hall Thrusters”, 1 -
34, Wiley & Sons, 2008.
2. G. L. Matloff, “Deep Space Probes: To the Outer Solar System and Beyond (Second
Edition)”, 30 - 57, Praxis Publishing Ltd, 2005.
3. M. Gruntman, R. McNutt, R. E. Gold, S. M. Krimigis, E. C. Roelof, J. C. Leary, G.
Gloeckler, P. L. Koehn, W. S. Kurth, S. R. Oleson, D. Fiehler, “Innovative Explorer
Mission To Interstellar Space”, International Academy of Astronautics, Italy (2005), 1 - 6.
4. R. L. McNutt, R. F. Wimmer-Schweingruber, The International Interstellar Probe Team,
“Enabling Interstellar Probe”, Acta Astronautica (2010), 790 - 801.
5. NASA/JPL-Caltech, “PIA 17046: V oyager Goes Interstellar (Artist Concept)”, September
12, 2013.
6. D.I. Fiehler, R. L. McNutt, “Mission Design for the Innovative Interstellar Explorer
Vision Mission”, Journal of Spacecraft and Rockets, V ol. 43, No. 6m November-
December 2006, 1239 - 1246.
7. Giorgini, J.D., Yeomans, D.K., Chamberlin, A.B., Chodas, P.W., Jacobson, R.A., Keesey,
M.S., Lieske, J.H., Ostro, S.J., Standish, E.M., Wimberly, R.N., "JPL's On-Line Solar
System Data Service", Bulletin of the American Astronomical Society, V ol 28, No. 3, p.
1158, 1996.
8. C.J. V oesenek, “Implementing a Fourth Order Runge-Kutta Method for Orbit
Simulation”, June 14, 2008.
79
Appendix A:
Raw Simulation Output Data:
Initial Position: Earth Orbit, Circular
A.1 Simulation Output Data: 5 year burn time:
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
A1 0 N/A N/A N/A N/A N/A N/A N/A N/A
A2 5 N/A N/A N/A N/A N/A N/A N/A N/A
A3 7 N/A N/A N/A N/A N/A N/A N/A N/A
A4 10 3.528 92.0855 14.5005 3.0047 1.893 1.2379 -11.0148 198.239
A5 12 0.27043 51.14037 19.9354 4.2716 3.7797 1.65588 -2.76307 161.426
Table 1: Set A Trajectory Results (tb = 5 years, P0 = 1000 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
B1 0 N/A N/A N/A N/A N/A N/A N/A N/A
B2 5 N/A N/A N/A N/A N/A N/A N/A N/A
B3 7 N/A N/A N/A N/A N/A 0.5635 9.747 N/A
B4 10 2.67281 68.4354 15.8102 3.53039 2.7331 1.508773 -5.28446 190.958
B5 12 0.21357 45.189 21.1008 4.757 4.346 1.90717 -2.0899 159.94
Table 2: Set B Trajectory Results (tb = 5 years, P0 = 1500 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
C1 0 N/A N/A N/A N/A N/A N/A N/A N/A
C2 5 N/A N/A N/A N/A N/A N/A N/A N/A
C3 7 N/A N/A N/A N/A N/A 0.67725 16.415 N/A
C4 10 2.187 57.8051 16.89441 3.9682 3.32768 1.760989 -3.56484 186.907
C5 12 0.18149 41.3269 22.076 5.1626 4.8037 2.13906 -1.71063 158.864
Table 3: Set C Trajectory Results (tb = 5 years, P0 = 2000 W)
80
Simulation
Set
Specific Impulse
(sec)
Accelerating
Voltage
(V)
Current
(A)
Mass Flow Rate
(kg/sec)
Thruster Force
(N)
Thruster
Exit Velocity
(m/sec)
A 2115.318 292.1278 2.053895 0.0000027886128 0.054878 20744.19156
B 2590.7259 438.1918 2.053895 0.0000027886128 0.067212 25406.342
C 2991.5126 584.2557 2.053895 0.0000027886128 0.07761 29336.717
Table 4: Engine Outputs for (tb = 5 years), Sets A, B & C
A.2 Simulation Output Data: 10 year burn time:
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final SMA
(AU)
Target
Angle
(deg)
D1 0 N/A N/A N/A N/A N/A N/A N/A N/A
D2 5 N/A N/A N/A N/A N/A N/A N/A N/A
D3 7 N/A N/A N/A N/A N/A 0.628 15.77932 N/A
D4 10 6.413 79.835 23.8292 2.9768 2.3555 1.7119 -7.11462 212.996
D5 12 0.416 46.134 37.9383 4.62927 4.39876 2.1683 -2.04015 163.808
Table 5: Set D Trajectory Results (tb = 10 years, P0 = 1000 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess
Velocity (AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
E1 0 N/A N/A N/A N/A N/A N/A N/A N/A
E2 5 N/A N/A N/A N/A N/A 0.887386 37.1803 N/A
E3 7 8.869 111.718 12.9966 2.925 1.5762 1.3146 -15.888 212.199
E4 10 4.765 59.8855 27.7483 3.739 3.3374 2.32144 -3.544 202.77
E5 12 0.323488 40.6077 41.1864 5.3009 5.1169 2.63529 -1.5077 162.442
Table 6: Set E Trajectory Results (tb = 10 years, P0 = 1500 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position @
ECO (AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
F1 0 N/A N/A N/A N/A N/A 0.448146 6.701 N/A
F2 5 9.417 126.732 8.25598 3.36742 1.3329 1.2232 -22.218 93.6641
F3 7 7.565 76.475 14.5499 3.50888 2.6241 1.844894 -5.7326 163.855
F4 10 3.8516 50.6792 30.9234 4.36709 4.0642 2.85642 -2.3898 196.996
F5 12 0.27192 37.0632 43.89719 5.86091 5.7054 3.06686 -1.2127 161.4
Table 7: Set E Trajectory Results (tb = 10 years, P0 = 2000 W)
81
Simulation
Set
Specific Impulse
(sec)
Accelerating
Voltage
(V)
Current
(A)
Mass Flow Rate
(kg/sec)
Thruster Force
(N)
Thruster
Exit Velocity
(m/sec)
D 2991.512 584.2557 1.0269 0.0000013943064 0.03880529 29336.717
E 3663.839 876.3836 1.0269 0.0000013943064 0.04752658 35929.993
F 4230.637 1168.511 1.0269 0.0000013943064 0.05487897 41488.383
Table 8: Engine Outputs for (tb = 10 years), Sets D, E & F
A.3 Simulation Output Data: 15 year burn time:
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
G1 0 N/A N/A N/A N/A N/A N/A N/A N/A
G2 5 N/A N/A N/A N/A N/A 0.745497 19.659 N/A
G3 7 14.2766 148.7209 9.6966 3.0506 1.079 1.1812 -33.8939 107.917
G4 10 9.209 75.159 31.7618 3.1027 2.6723 2.478336 -5.5278 228.538
G5 12 0.54294 43.8958 56.86198 5.0269 4.8869 2.698334 -1.6529 165.174
Table 9: Set G Trajectory Results (tb = 15 years, P0 = 1000 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
H1 0 N/A N/A N/A N/A N/A N/A N/A N/A
H2 5 14.351 147.081 10.44966 2.96348 1.10767 1.2233 -32.174 40.7016
H3 7 12.7714 101.155 16.2242 2.9482 1.9559 1.9229 -10.318 314.575
H4 10 6.762 56.4888 39.33329 4.04468 3.78843 3.46688 -2.7505 213.3
H5 12 0.41764 38.739 62.747 5.8359 5.72713 3.38042 -1.2035 163.855
Table 10: Set H Trajectory Results (tb = 15 years, P0 = 1500 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
I1 0 N/A N/A N/A N/A N/A 0.825187 30.1161 N/A
I2 5 12.426 88.80568 16.7044 3.1799 2.320747 2.025461 -7.3294 257.18
I3 7 10.5215 70.6534 23.483 3.569 3.06222 2.9385 -4.2097 240.61
I4 10 5.42092 47.9284 45.316 4.815 4.6315 4.30031 -1.8402 205.81
I5 12 0.3486 35.465 67.654 6.51029 6.42 4.01068 -0.9577 162.89
Table 11: Set I Trajectory Results (tb = 15 years, P0 = 2000 W)
82
Simulation
Set
Specific Impulse
(sec)
Accelerating
Voltage
(V)
Current
(A)
Mass Flow Rate
(kg/sec)
Thruster Force
(N)
Thruster
Exit Velocity
(m/sec)
G 3663.8397 876.383 0.684631 0.0000009295376 0.03168 35929.993
H 4487.2689 1314.575 0.684631 0.0000009295376 0.0388052 44005.075
I 5181.4518 1752.767 0.684631 0.0000009295376 0.044808 50812.684
Table 12: Engine Outputs for (tb = 15 years), Sets G, H & I
83
Appendix B:
Raw Simulation Output Data:
Initial Position: Earth Orbit, Perihelion
B.1 Simulation Output Data: 15 year burn time:
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
AP1 0 N/A N/A N/A N/A N/A N/A N/A N/A
AP2 5 N/A N/A N/A N/A N/A 0.77764 21.9201 N/A
AP3 7 14.463 161.322 9.3958 3.0437 0.92817 1.14714 -45.821 89.2525
AP4 10 8.449 70.682 33.916 3.2455 2.86459 2.492388 -4.8105 218.6211
AP5 12 0.2709 42.9977 58.2528 5.1317 4.997 2.69444 -1.5803 163.4953
Table 13: Set Ap Trajectory Results (tb = 15 years, P0 = 1000 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
BP1 0 N/A N/A N/A N/A N/A N/A N/A N/A
BP2 5 14.388 150.228 10.7748 2.909 1.066 1.2155 -34.734 23.9648
BP3 7 12.5831 98.357 17.4057 2.934 2.01868 1.97125 -9.687 294.359
BP4 10 6.2017 54.4254 41.2271 4.181 3.945 3.42469 -2.536 206.848
BP5 12 0.2134 38.0949 64.0958 5.936 5.831 3.36269 -1.1609 162.315
Table 14: Set Bp Trajectory Results (tb = 15 years, P0 = 1500 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
CP1 0 N/A N/A N/A N/A N/A 0.8205 29.47 N/A
CP2 5 12.291 86.954 17.1397 3.206 2.3819 2.053177 -6.958 244.686
CP3 7 10.2664 68.762 24.336 3.637 3.1597 2.96267 -3.954 228.368
CP4 10 4.9747 46.6249 47.065 4.9437 4.771 4.2269 -1.734 200.738
CP5 12 0.1811 34.956 68.971 6.607 6.5198 3.9806 -0.929 161.447
Table 15: Set Cp Trajectory Results (tb = 15 years, P0 = 2000 W)
84
Simulation
Set
Specific Impulse
(sec)
Accelerating
Voltage
(V)
Current
(A)
Mass Flow Rate
(kg/sec)
Thruster Force
(N)
Thruster
Exit Velocity
(m/sec)
Ap 3663.8397 876.383 0.684631 0.0000009295376 0.03168 35929.993
Bp 4487.2689 1314.575 0.684631 0.0000009295376 0.0388052 44005.075
Cp 5181.4518 1752.767 0.684631 0.0000009295376 0.044808 50812.684
Table 16: Engine Outputs for (tb = 15 years), Sets Ap, Bp & Cp
Initial Position: Earth Orbit, Aphelion
B.2 Simulation Output Data: 15 year burn time:
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
Aa1 0 N/A N/A N/A N/A N/A N/A N/A N/A
Aa2 5 N/A N/A N/A N/A N/A 0.706 17.603 N/A
Aa3 7 14.1412 142.0047 10.324 3.00069 1.1649 1.19726 -29.087 126.858
Aa4 10 9.9799 80.5242 29.426 2.965 2.4725 2.449484 -6.457 241.136
Aa5 12 0.8762 44.8652 55.434 4.9188 4.7719 2.70134 -1.7335 167.016
Table 17: Set Aa Trajectory Results (tb = 15 years, P0 = 1000 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU
(yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
Ba1 0 N/A N/A N/A N/A N/A N/A N/A N/A
Ba2 5 14.299 143.188 10.246 3.00859 1.1604 1.2342 -29.316 57.6566
Ba3 7 12.8902 102.666 15.0781 2.99225 1.928 1.883 -10.618 336.452
Ba4 10 7.34454 58.8347 37.324 3.9058 3.6249 3.50795 -3.0042 221.017
Ba5 12 0.6648 39.4268 61.366 5.733 5.619 3.397 -1.2499 165.536
Table 18: Set Ba Trajectory Results (tb = 15 years, P0 = 1500 W)
Scenario
ID
Launch
Velocity
(km/sec)
Time to
Escape
(yr)
Time to
200 AU (yr)
Position
@ ECO
(AU)
Velocity
@ ECO
(AU/yr)
Hyperbolic
Excess Velocity
(AU/yr)
Final
Eccentricity
Final
SMA
(AU)
Target
Angle
(deg)
Ca1 0 N/A N/A N/A N/A N/A 0.828 30.585 N/A
Ca2 5 12.557 90.678 16.217 3.16004 2.26222 1.9958 -7.7135 270.898
Ca3 7 10.7779 72.658 22.5202 3.5082 2.9667 2.9086 -4.4849 254.266
Ca4 10 5.886 49.373 43.479 4.6844 4.486 4.3775 -1.9611 211.697
Ca5 12 0.549 36.0034 66.307 6.411 6.317 4.0408 -0.9891 164.469
Table 19: Set Ca Trajectory Results (tb = 15 years, P0 = 2000 W)
85
Simulation
Set
Specific Impulse
(sec)
Accelerating
Voltage
(V)
Current
(A)
Mass Flow Rate
(kg/sec)
Thruster Force
(N)
Thruster
Exit Velocity
(m/sec)
Aa 3663.389 876.383 0.684631 0.0000009295376 0.03168439 35929.9937
Ba 4487.268 1314.575 0.684631 0.0000009295376 0.03880529 44005.0755
Ca 5181.45 1752.767 0.684631 0.0000009295376 0.044808 50812.684
Table 20: Engine Outputs for (tb = 15 years), Sets Aa, Ba & Ca
86
Appendix C:
Matlab Code
C.1 Main Script - “Orbit_Prop.m”
% Main Simulator Script:
% This program organizes the simulation initial conditions, and calls the
% Runge-Kutta numerical integrator, and conducts analysis of the simulation
% results (NOTE: use must alter the titles of the output plots to reflect
% the name and number of the scenario).
%
% Inputs:
% mu: standard gravitational parameter (m^3/s^2)
% tb: engine burn time (s)
% tspan: time span of orbit propagation
% h_estimate: estimate of integration time step (s)
% V_launch: launch velocity relative to the earth (m/s)
% r_parking: initial distance from the sun (m)
% a_parking: initial heliocentric orbital semi-major axis (m)
%
% Output Values:
% a_final_AU: final orbit semi-major axis (AU)
% eccentricity_final: final orbit eccentricity
% R_escape_ECO_AU: heliocentric radius at ECO (AU)
% V_escape_ECO_km_sec: heliocentric velocity at ECO (km/s)
% V_escape_ECO_AU_yr: heliocentric velocity at ECO (AU/yr)
% V_200_AU_in_AU_yr: heliocentric velocity at 200 AU (AU/yr)
% V_inf_AU_yr: hyperbolic excess velocity (AU/yr)
% Total_Flight_Time_to_200_AU: Total flight time to 200 AU (yr)
% Target_Angle: solar escape target angle (deg)
%
% Output Plots:
% 3D heliocentric trajectory, helio velocity vs. time, helio velocity
% vs. helio distance, eccentricity vs. time, SMA vs. time, total S/C
% mass vs. time, propellant mass flow rate vs. time, thruster force
% vs. time, thruster Isp vs. time, thruster exit velocity vs. time,
% accelerating voltage vs. time, and beam current vs. time.
clear all
close all
format long
clc
global tb mu h yr t1
fprintf(datestr(now))
fprintf('\n')
fprintf('Start the Clock! \n \n')
tic;
87
t1 = tic; % Starts "stop watch"
mu = 1.32712440018*10^(20); % Solar gravitational constant (m^3/s^2)
yr = (60*60*24*365.242); % 1 year in (sec)
AU = 149597870700; % 1 AU in (m)
rp_Earth = 0.9831235870678204; % in (AU) Taken on July 29, 2013
ra_Earth = 1.016517960841986; % in (AU) Taken on July 29, 2013
a_Earth = 0.9998207739549034; % in (AU) Taken on July 29, 2013
%%% User Defined: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tb = 15*yr; % Burn time (sec)
tspan = 17*yr; % integration time (seconds)
h_estimate = 100; % # of (sec/step)
V_launch = 12000; % Launch velocity relative to Earth (m/sec)
r_parking = rp_Earth*AU; % Initial parking orbit (m)
a_parking = a_Earth*AU; % (m) a = r_parking for circular orbit about sun
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V_Earth = sqrt(2*mu/r_parking - mu/a_parking); % (m/sec)
y0 = [r_parking 0 0 0 (V_Earth+V_launch) 0]; % initial (inertial) state:
% [rx ry rz vx vy vz] in [m m/s] (1AU circ) V = 29784.69183 m/s
N_steps = round(tspan/h_estimate); % estimate # of integration steps (steps)
[t y t_escape_index eccentricity a_inst mdot Ve Isp Force_thrust M_tot I Voltage] =
rk4_Orb(tspan, N_steps, y0); % Calls RK-4 Integrator
fprintf('Stop the Clock! \n \n')
Elapsed_Time_min = (toc(t1))/60 % Stops "stop watch", outputs elapsed time (min)
% Calculate final orbit right after ECO:
a_final = mu/((2*mu/norm(y((round(tb/h)+1),1:3)))-norm(y((round(tb/h)+1),4:6))^2);
a_final_AU = a_final/AU % final SMA
H_ang_final = norm(cross(y((round(tb/h)+1),1:3),y((round(tb/h)+1),4:6))); % magnitude
% of Angular Momentum (km^2/sec)
p_semi = (H_ang_final^2)/mu; % semi parameter (km)
eccentricity_final = sqrt(1 - p_semi/a_final) % final orbit eccentricity
if t_escape_index > 0 % if spacecraft has escaped...
R_escape_ECO_AU = (1/AU)*norm(y(round(tb/h),1:3)) % Position @ ECO (AU)
V_escape_ECO_km_sec = (10^-3)*norm(y(round(tb/h),4:6)) % Velocity @ ECO (km/sec)
V_escape_ECO_AU_yr = (yr/AU)*norm(y(round(tb/h),4:6)) % Velocity @ ECO (AU/yr)
V_200_AU_in_km_sec =(10^-3)*sqrt(2*mu/(200*AU)-mu/a_final) % Velocity@200AU (km/s)
V_200_AU_in_AU_yr=(yr/AU)*sqrt(2*mu/(200*AU)-mu/a_final) % Velocity@200 AU (AU/yr)
V_inf_AU_yr = (yr/AU)*sqrt(mu/(-a_final)) % Hyperbolic Excess Velocity (AU/yr)
True_Anom_200_AU = acos(((a_final*(1-eccentricity_final^2))/(200*AU) - 1)/
eccentricity_final); % True Anomaly @ 200 AU on final orbit (rad)
H_200_AU = 2*atanh(sqrt((eccentricity_final-1)/(eccentricity_final+1))*tan
(True_Anom_200_AU/2)); % Hyperboic Eccentric Anomaly @ 200 AU (rad)
t_hyper_200_AU = (eccentricity_final*sinh(H_200_AU)-H_200_AU)/sqrt(mu/(-
(a_final^3))); % Time from hyperbolic epoch to 200 AU (sec)
True_Anom_burnout = acos(((a_final*(1-eccentricity_final^2))/(norm(y((round(tb/h)
+1),1:3))) - 1)/eccentricity_final); % True Anomaly @ burnout on final orbit (rad)
H_burnout = 2*atanh(sqrt((eccentricity_final-1)/(eccentricity_final+1))*tan
(True_Anom_burnout/2)); % Hyperboic Eccentric Anomaly @ burnout (rad)
t_hyper_burnout = (eccentricity_final*sinh(H_burnout)-H_burnout)/sqrt(mu/(-
(a_final^3))); % Time from hyperbolic epoch to burnout (sec)
88
Total_Flight_Time_to_200_AU = (tb + (t_hyper_200_AU - t_hyper_burnout))/yr
% Total flight time to 200 AU (yr)
% Calculate Target Angle:
e_vec = (1/mu)*cross(y(end, 4:6), cross(y(end, 1:3), y(end, 4:6))) - (1/norm(y
(end, 1:3)))*y(end, 1:3); % eccentricity vector
e_hat = (1/norm(e_vec))*e_vec; % eccentricity unit vector
phi = acos(1/eccentricity_final); % 1/2 of hyper. approach/departure angle
Angle = acos(dot([1 0 0], (-e_hat))); % angle between -e_hat direction and X-axis
if (-e_hat(1, 2)) > 0
Target = Angle - phi; % angle between half hyper approach and Angle
else
Angle = 2*pi - Angle;
Target = Angle - phi;
end
if Target < 0! ! % Keeps domain: 0 <= Target angle <= 360
Target = 2*pi - abs(Target);
end
Target_Angle_deg = (180/pi)*Target
end
% Plots the Trajectory
figure
hold on
title('Scenario Bp5: Trajectory: 15 year burn, 17 year propagation', 'fontsize', 16,
'fontweight','b');
plot3((1/AU)*y(1,1), (1/AU)*y(1,2), (1/AU)*y(1,3), 'r.', 'LineWidth', 1.001,
'MarkerSize', 6); % plots ignition
plot3((1/AU)*y(1:round(tb/h),1), (1/AU)*y(1:round(tb/h),2), (1/AU)*y(1:round(tb/h),3),
'r'); % trajectory, converted from [m] to [AU]
if t_escape_index > 0
plot3((1/AU)*y(t_escape_index,1), (1/AU)*y(t_escape_index,2), (1/AU)*y
(t_escape_index,3), 'g.', 'LineWidth', 1.01, 'MarkerSize', 9); % if S/C has escaped,
% plot '.' to indicate when
plot3((1/AU)*y((round(tb/h)+1),1), (1/AU)*y((round(tb/h)+1),2), (1/AU)*y((round
(tb/h)+1),3), 'b.', 'LineWidth', 1.01, 'MarkerSize', 9);
plot3((1/AU)*y((round(tb/h)+1):end,1), (1/AU)*y((round(tb/h)+1):end,2), (1/AU)*y
((round(tb/h)+1):end,3), 'b');
plot3(0, 0, 0, 'y*', 'LineWidth', 1.001, 'MarkerSize', 12)
legend('Engine Ignition', 'Engine Burn','Escape Velocity Achieved', 'Engine Cut-
off (ECO)', 'Final Orbit' );
else
plot3((1/AU)*y((round(tb/h)+1),1), (1/AU)*y((round(tb/h)+1),2), (1/AU)*y((round
(tb/h)+1),3), 'b.', 'LineWidth', 1.01, 'MarkerSize', 9);
plot3((1/AU)*y((round(tb/h)+1):end,1), (1/AU)*y((round(tb/h)+1):end,2), (1/AU)*y
((round(tb/h)+1):end,3), 'b');
plot3(0, 0, 0, 'y*', 'LineWidth', 1.001, 'MarkerSize', 12)
legend('Engine Ignition', 'Engine Burn', 'Engine Cut-off (ECO)', 'Final Orbit' );
end
xlabel('X (AU)', 'fontsize', 14);
ylabel('Y (AU)', 'fontsize', 14);
zlabel('Z (AU)', 'fontsize', 14);
if t_escape_index > 0
text(-6, 1.8, sprintf(' Heliocentric Position @ ECO = %0.3f (AU)',
R_escape_ECO_AU), 'fontsize', 12, 'fontweight','b', 'BackgroundColor', 'w',
'EdgeColor', 'k')
89
text(-6, 0.9, sprintf(' Heliocentric Velocity @ ECO = %0.3f (AU / yr)',
V_escape_ECO_AU_yr), 'fontsize', 12, 'fontweight','b', 'BackgroundColor', 'w',
'EdgeColor', 'k')
text(-6, 0, sprintf(' Solar System Escape Target Angle = %0.3f (deg)',
Target_Angle_deg), 'fontsize', 12, 'fontweight','b', 'BackgroundColor', 'w',
'EdgeColor', 'k')
text(-6, -0.9, sprintf(' Total Flight Time to 200 AU = %0.3f (years)',
Total_Flight_Time_to_200_AU), 'fontsize', 12, 'fontweight','b', 'BackgroundColor',
'w', 'EdgeColor', 'k')
end
axis equal
grid on
hold off
% Plots Magnitude of Heliocentric Velocity vs. Time
figure
hold on
if t_escape_index > 0
plot((1/yr)*t(t_escape_index), (yr/AU)*norm(y(t_escape_index, 4:6)), 'g.',
'LineWidth', 1.01, 'MarkerSize', 10); % time, norm(V_escape)
end
plot((1/yr)*t(round(tb/h)), (yr/AU)*norm(y(round(tb/h), 4:6)), 'b.', 'LineWidth',
1.01, 'MarkerSize', 10); % plot(ECO_time, norm(V_burnout))
plot(((1/yr)*t), (yr/AU)*sqrt(y(:, 4).^2+y(:, 5).^2+y(:, 6).^2), 'r');
% plot(time, norm(velocity))
grid on
if t_escape_index > 0
legend('Escape Velocity Achieved', 'Engine Cut-off', 'Location', 'NorthEast')
text(1, 1.8, sprintf(' Heliocentric Excess Velocity = %0.3f (AU/yr)',
V_inf_AU_yr), 'fontsize', 12, 'fontweight','b', 'BackgroundColor', 'w', 'EdgeColor',
'k')
else
legend('Engine Cut-off', 'Location', 'NorthEast')
end
xlabel('Time (years)', 'fontsize', 14);
ylabel('Heliocentric Speed (AU/yr)', 'fontsize', 14);
title('Scenario Bp5: Heliocentric Speed vs. Time', 'fontsize', 16, 'fontweight','b');
axis([0 17 2 8])
hold off
% Plots Magnitude of Heliocentric Velocity vs. Position
figure
hold on
plot((1/AU)*sqrt(y(:, 1).^2+y(:, 2).^2+y(:, 3).^2), (yr/AU)*sqrt(2*mu*(y(:, 1).^2+y(:,
2).^2+y(:, 3).^2).^(-0.5)), 'r--'); % norm(position), norm(V_escape_required)
plot((1/AU)*sqrt(y(:, 1).^2+y(:, 2).^2+y(:, 3).^2), (yr/AU)*sqrt(y(:, 4).^2+y(:,
5).^2+y(:, 6).^2), 'm'); % norm(position), norm(V_actual)
if t_escape_index > 0
plot((1/AU)*norm(y(t_escape_index, 1:3)), (yr/AU)*norm(y(t_escape_index, 4:6)),
'g.', 'LineWidth', 1.01, 'MarkerSize', 10); % norm(escape_position), norm
(escape_velocity)
end
plot((1/AU)*norm(y(round(tb/h), 1:3)), (yr/AU)*norm(y(round(tb/h), 4:6)), 'b.',
'LineWidth', 1.01, 'MarkerSize', 10); % norm(burnout_position), norm
(burnout_velocity)
grid on
90
if t_escape_index > 0
legend('Required Escape Velocity', 'Actual Velocity', 'Escape Velocity Achieved',
'Engine Cut-off', 'Location', 'NorthEast')
text(1, 1.8, sprintf(' Hyperbolic Excess Velocity = %0.3f (AU / yr)',
V_inf_AU_yr), 'fontsize', 12, 'fontweight','b', 'BackgroundColor', 'w', 'EdgeColor',
'k')
else
legend('Required Escape Velocity', 'Actual Velocity', 'Engine Cut-off',
'Location', 'NorthEast')
end
xlabel('Radial Position (AU)', 'fontsize', 14);
ylabel('Heliocentric Speed (AU / yr)', 'fontsize', 14);
title('Scenario Bp5: Heliocentric Speed vs. Position', 'fontsize', 16,
'fontweight','b');
hold off
% Plots Instantaneous Orbit Eccentricity
figure
hold on
if t_escape_index > 0
plot(((1/yr)*t(1:t_escape_index)), eccentricity(1:t_escape_index), 'g');
one = ones(size(t,1),1);
plot(((1/yr)*t), one, 'r--');
plot((1/yr)*t(t_escape_index), eccentricity(t_escape_index), 'r.', 'LineWidth',
1.01, 'MarkerSize', 10); % if S/C has escaped, plot '.' to indicate when
plot(((1/yr)*t(t_escape_index:end)), eccentricity(t_escape_index:end), 'b');
plot((1/yr)*t(round(tb/h)), eccentricity(round(tb/h)), 'b.', 'LineWidth', 1.01,
'MarkerSize', 9); % plots burnout
legend('Elliptical', 'Parabolic', 'Escape Eccentricity Achieved', 'Hyperbolic',
'Engine Cut-off', 'Location', 'Southeast');
else
plot(((1/yr)*t), eccentricity, 'g');
one = ones(size(t,1),1);
plot(((1/yr)*t), one, 'r--');
plot((1/yr)*t(round(tb/h)), eccentricity(round(tb/h)), 'b.', 'LineWidth', 1.01,
'MarkerSize', 9); % plots burnout
legend('Elliptical', 'Parabolic', 'Engine Cut-off', 'Location', 'Southeast');
end
xlabel('Time (years)', 'fontsize', 14);
ylabel('Eccentricity', 'fontsize', 14);
title('Scenario Bp5: Instantaneous Eccentricity vs. Time', 'fontsize', 16,
'fontweight','b');
grid on
axis([0 17 0 1.4])
hold off
% Plots Instantaneous Semi-Major Axis
figure
hold on
plot((1/yr)*t(round(tb/h)), (1/AU)*a_inst(round(tb/h)), 'b.', 'LineWidth', 1.01,
'MarkerSize', 9);
plot((1/yr)*t(1:(t_escape_index-1)), (1/AU)*a_inst(1:(t_escape_index-1)), 'r');
plot((1/yr)*t((t_escape_index+1):end), (1/AU)*a_inst((t_escape_index+1):end), 'r');
legend('Engine Cut-off', 'Location', 'Northeast');
xlabel('Time (years)', 'fontsize', 14);
ylabel('Semi-major Axis (AU)', 'fontsize', 14);
91
title('Scenario Bp5: Instantaneous Semi-major Axis vs. Time', 'fontsize', 16,
'fontweight','b');
grid on
axis([0 16 -100 100])
hold off
% Plots total Mass
figure
hold on
plot(((1/yr)*t), M_tot, 'm');
xlabel('Time (years)', 'fontsize', 14);
ylabel('Total Mass (kg)', 'fontsize', 14);
title('Set Bp: Spacecraft Total Mass vs. Time', 'fontsize', 16, 'fontweight','b');
axis([0 17 750 1250])
grid on
hold off
% Plots Thruster Force
figure
hold on
plot(((1/yr)*t), Force_thrust, 'r');
xlabel('Time (years)', 'fontsize', 12);
ylabel('Force (N)', 'fontsize', 12);
title('Set Bp: Ion Engine Thruster Force vs. Time', 'fontsize', 14, 'fontweight','b');
grid on
axis([0 17 0 0.08])
hold off
% Plots Isp
figure
hold on
plot(((1/yr)*t), Isp, 'b');
xlabel('Time (years)', 'fontsize', 12);
ylabel('Isp (sec)', 'fontsize', 12);
title('Set Bp: Ion Engine Specific Impulse vs. Time', 'fontsize', 14,
'fontweight','b');
grid on
axis([0 17 0 5500])
hold off
% Plots Thruster Exit Velocity
figure
hold on
plot(((1/yr)*t), Ve, 'r');
xlabel('Time (years)', 'fontsize', 12);
ylabel('Exit Velocity (m/sec)', 'fontsize', 12);
title('Set Bp: Ion Engine Exit Velocity vs. Time', 'fontsize', 14, 'fontweight','b');
grid on
hold off
% Plots Mass Flow Rate
figure
hold on
plot(((1/yr)*t), mdot, 'm');
xlabel('Time (years)', 'fontsize', 12);
ylabel('Mass Flow Rate (kg/sec)', 'fontsize', 12);
title('Set Bp: Ion Engine Mass Flow Rate vs. Time', 'fontsize', 14, 'fontweight','b');
grid on
92
axis([0 17 0 1.8*10^-6])
hold off
% Plots Accelerating Voltage
figure
hold on
plot(((1/yr)*t), Voltage, 'r');
xlabel('Time (years)', 'fontsize', 12);
ylabel('Accelerating Voltage (V)', 'fontsize', 12);
title('Set Bp: Ion Engine Accelerating Voltage vs. Time', 'fontsize', 14,
'fontweight','b');
grid on
hold off
% Plots Current
figure
hold on
plot(((1/yr)*t), I, 'r');
xlabel('Time (years)', 'fontsize', 12);
ylabel('Current (A)', 'fontsize', 12);
title('Set Bp: Ion Engine Electrical Current vs. Time', 'fontsize', 14,
'fontweight','b');
grid on
hold off
C.2 Fourth Order Runge-Kutta Function - “rk4_Orb.m”
% 4th Order Runge Kutta Orbit Integrator:
%
% INPUTS:
% tspan = integration time interval (sec)
% N_steps = # of integration steps (integer)
% y0 = initial state [x y z xdot ydot zdot] in inertial cord, (m & m/sec)
%
% OUTPUTS:
% t: array of time values for each integration step (sec)
% y: matrix of states for each time value [x y z xdot ydot zdot] in
inertial cord, (km & km/sec)
% t_escape_yrs: time of escape (yr)
% t_escape_index: time array index value of time of escape
% Individual arrays of Ion Engine parameters
function [t y t_escape_index eccentricity a_inst mdot Ve Isp Force_thrust M_tot I
Voltage] = rk4_Orb(tspan, N_steps, y0)
global h mu t1;
h = tspan/N_steps; % # of seconds per step (seconds/step)
y = zeros(N_steps+1, 6); % establishes output state matrix
t = zeros(N_steps+1, 1); % establishes output array of time values
test = 1; % on/off switch for “check Escape Velocity” if loop
t_escape_index = 0; % Initializes index of time array "escape index"
% initialize engine parameter arrays:
mdot = zeros(N_steps+1,1); % array of mass flow rate values (kg/sec)
Ve = zeros(N_steps+1,1); % array of exit velocity values (m/sec)
93
Isp = zeros(N_steps+1,1); % array of specific Impulse values (sec)
Force_thrust = zeros(N_steps+1,1); % array of thrust force values (N)
M_tot = zeros(N_steps+1,1); % array of total mass values (kg)
I = zeros(N_steps+1,1); % array of current values (A)
Voltage = zeros(N_steps+1,1); % array of voltage values (V)
eccentricity = zeros(N_steps+1,1); % array if instantaneous orbit eccentricities
a_inst = zeros(N_steps+1,1); % array of instantaneous semi-major axis (m)
R0 = [y0(1) y0(2) y0(3)]; % initial position in inertial cord, m
V0 = [y0(4) y0(5) y0(6)]; % initial velocity in inertial cord, m/sec
A0 = A_tot(R0, V0, t(1)); % combined accel’s from functions A_Grav() & A_thrust()
y(1,1:3) = R0; % initial state
y(1,4:6) = V0;
R_prev = R0; % establishes initial conditions in RK-4 loop
V_prev = V0;
A_prev = A0;
[mdot(1) Ve(1) Isp(1) Force_thrust(1) M_tot(1) I(1) Voltage(1)] = thruster_params(0);
% initial thruster params
for j = 2:(N_steps + 1) % RK-4 Loop for N_steps, starting with second time value
(IC's count as first)
t(j) = t(j-1) + h; % next time value = previous time value + 1 time step in (sec)
% Runge-Kutta Coefficient (1x3) Vectors:
k1_v_next = A_prev;
k1_r_next = V_prev;
k2_v_next = A_tot((R_prev + k1_r_next * (h/2)), V_prev, t(j));
k2_r_next = V_prev + k1_v_next*(h/2);
k3_v_next = A_tot((R_prev + k2_r_next * (h/2)), V_prev, t(j));
k3_r_next = V_prev + k2_v_next*(h/2);
k4_v_next = A_tot((R_prev + k3_r_next * h), V_prev, t(j));
k4_r_next = V_prev + k3_v_next * h;
% Next calculated velocity & position:
V_next = V_prev + (h/6) * (k1_v_next + 2*k2_v_next + 2*k3_v_next + k4_v_next);
R_next = R_prev + (h/6) * (k1_r_next + 2*k2_r_next + 2*k3_r_next + k4_r_next);
y(j,1:3) = R_next; % recorded state matrix (output)
y(j,4:6) = V_next;
% Check to see if S/C has already reached escape velocity
if norm(V_next) >= sqrt(2*mu/norm(R_next)) && test == 1
fprintf(datestr(now))
fprintf('\n')
fprintf('S/C has escaped! \n \n')
t_escape_yrs = t(j)/(60*60*24*365.242) % time of escape (years)
t_escape_index = j; % time array "t" index of escape t_escape
test = 0; % switches off escape loop once S/C has escaped
end
94
% Simulation Time Status Loop:
if j == round((N_steps+1)/4)
fprintf(datestr(now))
fprintf('\n')
fprintf('%0.1f min: 25%% of Integration Complete... \n \n', toc(t1)/60)
elseif j == round((N_steps+1)/2)
fprintf(datestr(now))
fprintf('\n')
fprintf('%0.1f min: 50%% of Integration Complete... \n \n', toc(t1)/60)
elseif j == round(3*(N_steps+1)/4)
fprintf(datestr(now))
fprintf('\n')
fprintf('%0.1f min: 75%% of Integration Complete... \n \n', toc(t1)/60)
elseif j == (N_steps+1)
fprintf(datestr(now))
fprintf('\n')
fprintf('%0.1f min: 100%% of Integration Complete... \n \n', toc(t1)/60)
end
% Calculate instantaneous orbit:
a_inst(j) = mu/((2*mu/norm(R_next)) - norm(V_next)^2);
H_inst = norm(cross(R_next,V_next));
p_inst = (H_inst^2)/mu;
eccentricity(j) = sqrt(1 - p_inst/a_inst(j));
[mdot(j) Ve(j) Isp(j) Force_thrust(j) M_tot(j) I(j) Voltage(j)] = thruster_params
(t(j)); % updates output arrays with next set of thruster params
R_prev = R_next;
V_prev = V_next;
A_prev = A_tot(R_next, V_next, t(j)); % calculates next total acceleration
% from newly calculated state
end
C.3 Total Acceleration Function - “A_tot.m”
% Total Acceleration Function:
% Takes in R & V 1x3 vectors and time value, outputs 1x3 acceleration vector
function A_tot = A_tot(R, V, t)
A_tot = A_Grav(R) + A_thrust(V, t); % Combines accelerations from functions A_Grav()
and A_thrust(), in (m/sec^2)
95
C.4 Acceleration Function: Gravity - “A_grav.m”
% Gravity Acceleration Function:
% Takes in position vector, R, outputs 3D acceleration vector
function A_Grav = A_Grav(R)
global mu % Defined in Orbit_Prop.m
xddot = -(R(1)*mu)/(((R(1))^2 + (R(2))^2 + (R(3))^2)^(3/2)); % x acceleration
yddot = -(R(2)*mu)/(((R(1))^2 + (R(2))^2 + (R(3))^2)^(3/2)); % y acceleration
zddot = -(R(3)*mu)/(((R(1))^2 + (R(2))^2 + (R(3))^2)^(3/2)); % z acceleration
A_Grav = [xddot yddot zddot]; % acceleration in intertial cord, (m/sec^2)
C.5 Acceleration Function: Thruster - “A_thrust.m”
% Thruster Acceleration Function:
% Takes in velocity vector and time value, outputs 3D acceleration vector
function A_thrust = A_thrust(V, t)
[mdot Ve Isp Force_thrust M_tot I Voltage] = thruster_params(t); % Retrieve
% instantaneous thruster params
V_unit = (1/norm(V))*V; % Velocity unit vector (direction of thrust)
A_thrust = ((Force_thrust / M_tot)) * V_unit; % Acceleration due to thruster force
along V_unit direction, (m/sec^2)
C.6 Instantaneous Thruster Parameters Function - “thruster_params.m”
% Instantaneous Thruster Parameters Function:
% The Ion Engine Model that calculates the various engine parameters.
%
% User Inputs:
% Mxe: Atomic mass of xenon (or desired propellant) (kg)
% Mo: Structural dry mass of spacecraft (kg)
% Mp: Initial total mass of propellant (kg)
% P_total_in: Total input electric power from RTG (W)
% Eff_power: Power (electrical) efficiency
% Eff_mass_util: Mass utilization efficiency
% theta: Average half-angle divergence of the beam exhaust (rad)
% IPP_IP: Fraction of doubly chaged ion current in the beam
function [mdot Ve Isp Force_thrust M_tot I Voltage] = thruster_params(t)
global h tb yr;
e_charge = 1.602176565*10^(-19); % Fundemental charge unit (C)
AMU = 1.66053886*(10^(-27)); % 1 Atomic Mass Unit in (kg)
96
%%%% User Defined: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Mxe = 131*AMU; % Atomic mass of Xenon in kg
Mo = 790; % S/C structure dry mass (burnout) (kg)
Mp = 440; % Propellant mass (kg)
P_total_in = 1500; % Total input electric power (W)
Eff_power = 0.6; % Power (electric) efficiency
Eff_mass_util = 0.9; !! % Mass utilization efficiency
theta = 0*pi/180; % Average half-angle divergence of beam (rad)
IPP_IP = 0; % Fraction of double ion current in the beam
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pb = P_total_in * Eff_power; % Beam Power for thruster (W)
alpha = (1+ 0.707*IPP_IP)/(1+IPP_IP); % Trust correction factor
gamma = alpha*cos(theta); % total thrust correction
Force_thrust = 0; % Initialize Thrust Force
M_tot = Mo + Mp; % Initialize Total Mass
if t <= tb ! ! ! ! % If current time is within the burn phase...
mdot=(Mp/tb); ! ! ! ! % Mass flow rate function (kg/sec)
Ve = sqrt((2*Pb)/mdot); % Exit Velocity (m/sec)
Isp = Ve/9.80665; % Specific Impulse (sec)
Ndot = mdot/Mxe; % # of Xe+ ions per sec (#/sec)
I = Ndot*e_charge; % Beam current (A)
Voltage = Pb/I; % Accelerating voltage (V)
Eff_total = (gamma^2)*Eff_power*Eff_mass_util; % Total efficiency
Force_thrust = sqrt(2*Eff_total*mdot*P_total_in); % Force of thruster (N)
M_tot = Mo + (Mp - mdot*t); % Total S/C and Mp mass (kg)
if mdot*t > Mp % Loop that ensures that if Mp is used up before completed burn
time, M_tot = dry mass only and F_thrust = 0
M_tot = Mo;
Force_thrust = 0;
Ve = 0;
Isp = 0;
end
% Engine "off" loop for after the burn phase is complete:
elseif t > tb
mdot = 0; % Mass flow rate function (kg/sec)
Ve = 0; % Exit Velocity (m/sec)
Isp = 0; % Specific Impulse (sec)
Force_thrust = 0; % Force of thruster (N)
M_tot = Mo; % Total S/C mass without Mp (kg), assumes all Mp has been used up
Ndot = 0; % # of Xe+ ions per sec (#/sec)
I = 0; % Accelerating current (A)
Voltage = 0; % Accelerating voltage (V)
end
97
Abstract (if available)
Abstract
This research explores a variety of mission and system architectures for an unmanned Interstellar Precursor Mission (IPM) spacecraft with a Radioisotope Thermoelectric Generator (RTG) powered Ion Engine using Xenon propellant, traveling on a (direct) ballistic escape trajectory to the undisturbed Interstellar Medium (~200 AU). The main goal of this work was to determine the relationship between the propulsion system design parameters and the ensuing escape trajectory. To do this, an orbit simulator was created in Matlab using a fourth order Runge-Kutta numerical integration method to propagate the thrusting spacecraft's trajectory through time. The accelerations due to the Sun's gravity and the Ion Engine thrust were modeled separately and then combined into a single total acceleration vector at each time step, with the thrust direction assumed to be in the direction of the spacecraft's instantaneous velocity vector. The propellant of the thruster was also designed to be completely consumed by the time of engine cut-off (ECO), meaning a constant propellant mass flow rate. Simulations were run for burn times of 5, 10 & 15 years, with heliocentric launch velocities of 0, 5, 7, 10 & 12 km/sec from a circular 1 AU Earth orbit, and with RTG supplied engine input powers of 1000, 1500 & 2000 W. A total of 45 simulations were run for the circular 1 AU case, as well as additional comparison simulations for launches from an elliptical Earth orbit at perihelion and aphelion. ❧ The results of these simulations yielded many interesting results on the total fly-out times to 200 AU, which ranged dramatically from ~35 to ~140 years depending on the propulsion system settings and orbital initial conditions, as well as descriptions of the ECO distances from the Sun for each mission. The simulations also revealed the inherent gravitational maneuver inefficiency felt by all low thrust spacecraft, which becomes more apparent under certain conditions. Relations between launch velocity and the various characteristics of the escape trajectory were also shown. Finally, the propulsion system model showcased many details about the performance of the Ion Engine, including the consequence that increasing the input power level to the engine will result in a greater increase of thrust when operating at shorter engine burn times.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion
PDF
Trajectory mission design and navigation for a space weather forecast
PDF
Incorporation of mission scenarios in deep space spacecraft design trades
PDF
The space environment near the Sun and its effects on Parker Solar Probe spacecraft and FIELDS instrumentation
PDF
Systems engineering and mission design of a lunar South Pole rover mission: a novel approach to the multidisciplinary design problem within a spacecraft systems engineering paradigm
PDF
Techniques for analysis and design of temporary capture and resonant motion in astrodynamics
PDF
Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model
PDF
The development of an autonomous subsystem reconfiguration algorithm for the guidance, navigation, and control of aggregated multi-satellite systems
PDF
Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers
PDF
Numerical and experimental investigations of dust-plasma-asteroid interactions
PDF
An investigation into magnetically induced extractor-less electrospray propulsion devices
PDF
Experimental and numerical investigations of charging interactions of a dusty surface in space plasma
PDF
Ice: an impedance spectroscopy and atomistic simulation study
PDF
An experimental investigation of cathode erosion in high current magnetoplasmadynamic arc discharges
PDF
Study of rarefaction effects in gas flows with particle approaches
PDF
Advanced nuclear technologies for deep space exploration
PDF
Kinetic studies of collisionless mesothermal plasma flow dynamics
PDF
Numerical and experimental investigations of ionic electrospray thruster plume
PDF
Modeling human regulation of momentum while interacting with the environment
PDF
Modification of electrode materials for lithium ion batteries
Asset Metadata
Creator
Fogel, Joshua A.
(author)
Core Title
Mission design study of an RTG powered, ion engine equipped interstellar spacecraft
School
Viterbi School of Engineering
Degree
Master of Science
Degree Program
Astronautical Engineering
Publication Date
01/13/2014
Defense Date
10/16/2013
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
fourth order Runge-Kutta,interstellar,ion engine,mission simulation,OAI-PMH Harvest,RTG,trajectory design
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Gruntman, Michael (
committee chair
), Madni, Azad M. (
committee member
), Wang, Joseph (
committee member
)
Creator Email
joshuafo@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-360006
Unique identifier
UC11296210
Identifier
etd-FogelJoshu-2228.pdf (filename),usctheses-c3-360006 (legacy record id)
Legacy Identifier
etd-FogelJoshu-2228.pdf
Dmrecord
360006
Document Type
Thesis
Format
application/pdf (imt)
Rights
Fogel, Joshua A.
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
fourth order Runge-Kutta
interstellar
ion engine
mission simulation
RTG
trajectory design