connor/thesis
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

She did. Now I'm done!

Browse Source
This commit is contained in:
Connor
2022-03-15 22:20:45 -06:00
parent d00b977581
commit 298eb38ff1
7 changed files with 198 additions and 299 deletions
Download Patch File Download Diff File Expand all files Collapse all files

266
LaTeX/trajectory_design.tex
View File

@@ -16,24 +16,24 @@
very high-fidelity force models that account for aerodynamic pressure, solar radiation
pressure, multi-body effects, and other forces may be too time intensive for a
particular application. Initial surveys of the solution space often don't require such
complex models in order to gain valuable insight.
complex models in order to gain valuable preliminary insight.
Therefore, a common approach (and the one utilized in this implementation) is to first
use a lower-fidelity dynamical model that captures only the gravitational force due to
the primary body around which the spacecraft is orbiting. This approach can provide an
A common approach (and the one utilized in this implementation) is to first use a
lower-fidelity dynamical model that captures only the gravitational force due to the
primary body around which the spacecraft is orbiting. This approach can provide an
excellent low-to-medium fidelity model that is useful as an underlying model in an
algorithm for quickly categorizing a search space for initial mission feasibility
explorations.
In order to explore the Two Body Problem, we must first examine the full set of
assumptions associated with the force model\cite{vallado2001fundamentals}. Firstly, we
are only concerned with the nominative two bodies: the spacecraft and the planetary body
around which it is orbiting. Secondly, both of these bodies are modeled as point masses
with constant mass. This removes the need to account for non-uniform densities and
asymmetry. Finally, for convenience in notation at the end, we'll also assume that the
mass of the spacecraft ($m_2$) is much much smaller than the mass of the planetary body
($m_1$) and enough so as to be considered negligible. The only force acting on this
system is then the force of gravity that the primary body enacts upon the secondary.
are only concerned with the gravitational influence between the nominative two bodies:
the spacecraft and the planetary body around which it is orbiting. Secondly, both of
these bodies are modeled as point masses with constant mass. This removes the need to
account for non-uniform densities and asymmetry. Finally, for convenience in notation at
the end, we'll also assume that the mass of the spacecraft ($m_2$) is much much smaller
than the mass of the planetary body ($m_1$) and enough so as to be considered
negligible.
\begin{figure}[H]
\centering
@@ -45,7 +45,6 @@
Under these assumptions, the force acting on the body due to the law of universal
gravitation is:
\begin{align}
F_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\
F_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|}
@@ -53,7 +52,6 @@
And by Newton's second law (force is the product of mass and acceleration), we can
derive the following differential equations for $r_1$ and $r_2$:
\begin{align}
m_2 \ddot{\vec{r}}_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\
m_1 \ddot{\vec{r}}_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|}
@@ -65,7 +63,6 @@
inertial frame. $G$ is the universal gravitational parameter, $m_1$ is the mass of the
planetary body, and $m_2$ is the mass of the spacecraft. From these equations, we can
then determine the acceleration of the spacecraft relative to the planet:
\begin{equation}
\ddot{\vec{r}} = \ddot{\vec{r}}_2 - \ddot{\vec{r}}_1 =
- \frac{G \left( m_1 + m_2 \right)}{r^2} \frac{\vec{r}}{\left| r \right|}
@@ -76,27 +73,19 @@
negligible $m_2$ term. We can also introduce, for convenience, a gravitational parameter
$\mu$ which represents the gravity constant for the system about the center of motion
($\mu = G (m_1 + m_2) \approx G m_1$). Doing so and simplifying produces:
\begin{equation}
\ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r}
\end{equation}
We may also wish to utilize the total orbital energy for a spacecraft within this model.
Since the spacecraft is acting only under the gravitational influence of the planet and
no other forces, we can define the total specific mechanical energy as:
no other forces, we can define the total specific mechanical energy as
\cite{vallado2001fundamentals}:
\begin{equation} \label{energy}
\xi = \frac{v^2}{2} - \frac{\mu}{r}
\end{equation}
\noindent
Where the first term represents the kinetic energy of the spacecraft and the second term
represents the gravitational potential energy.
\subsection{Kepler's Laws}
Now that we've fully qualified the forces acting within the Two Body Problem, we can concern
ourselves with more practical applications of it as a force model. It should be noted,
firstly, that the spacecraft's position and velocity (given an initial position and velocity
@@ -105,6 +94,8 @@
one-dimensional equations (one for each component of the three-dimensional space) and
three unknowns (the three components of the second derivative of the position).
\subsection{Kepler's Laws}
In the early 1600s, Johannes Kepler produced just such a solution, by taking advantages of
what is also known as ``Kepler's Laws'' which are\cite{murray1999solar}:
@@ -113,68 +104,61 @@
expanded to any orbit by re-wording as ``all orbital paths follow a conic section
(circle, ellipse, parabola, or hyperbola) with a primary mass at one of the foci''.
Specifically the path of the orbit follows the trajectory equation:
The conic trajectory equation explains this observation and offers a description
of the path as:
\begin{equation}
r = \frac{\sfrac{h^2}{\mu}}{1 + e \cos(\theta)}
\end{equation}
Where $h$ is the angular momentum of the satellite, $e$ is the
where $h$ is the angular momentum of the satellite, $e$ is the
eccentricity of the orbit, and $\theta$ is the true anomaly, or simply
the angular distance the satellite has traversed along the orbit path.
the angular distance the satellite has traversed along the orbit path from
periapsis.
\item The area swept out by the imaginary line connecting the primary and secondary
bodies increases linearly with respect to time. This implies that the magnitude of the
orbital speed is not constant. For the moment, we'll just take this
value to be a constant:
\begin{equation}\label{swept}
\frac{\Delta t}{T} = \frac{k}{\pi a b}
\end{equation}
Where $k$ is the constant value, $a$ and $b$ are the semi-major and
where $k$ is the constant value, $a$ and $b$ are the semi-major and
semi-minor axis of the conic section, and $T$ is the period. In the
following section, we'll derive the value for $k$.
\item The square of the orbital period is proportional to the cube of the semi-major
axis of the orbit, regardless of eccentricity. Specifically, the relationship is:
axis of the orbit, regardless of eccentricity. For an elliptical orbit this
observation connects to the following known expression for the orbit period:
\begin{equation}
T = 2 \pi \sqrt{\frac{a^3}{\mu}}
\end{equation}
Where $T$ is the period and $a$ is the semi-major axis.
where $T$ is the period and $a$ is the semi-major axis.
\end{enumerate}
\subsection{Kepler's Equation}
Kepler was able to produce an equation to represent the angular displacement of an
orbiting body around a primary body as a function of time, which we'll derive now for
the elliptical case\cite{vallado2001fundamentals}. Since the total area of an ellipse is
the product of $\pi$, the semi-major axis, and the semi-minor axis ($\pi a b$), we can
relate (by Kepler's second law) the area swept out by an orbit as a function of time, as
we did in Equation~\ref{swept}. This leaves just one unknown variable $k$, which we can
determine through use of the geometric auxiliary circle, which is a circle with radius
equal to the ellipse's semi-major axis and center directly between the two foci, as in
Figure~\ref{aux_circ}.
the elliptical case\cite{vallado2001fundamentals}. Because the total area of an ellipse
is the product of $\pi$, the semi-major axis, and the semi-minor axis ($\pi a b$), we
can relate (by Kepler's second law) the area swept out by an orbit as a function of
time, as we did in Equation~\ref{swept}. This leaves just one unknown variable $k$,
which we can determine through use of the geometric auxiliary circle, which is a circle
with radius equal to the ellipse's semi-major axis and center directly between the two
foci, as in Figure~\ref{aux_circ}.
\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{fig/kepler}
\caption{Geometric Representation of Auxiliary Circle}\label{aux_circ}
\caption{Geometric representation of auxiliary circle}\label{aux_circ}
\end{figure}
In order to find the area swept by the spacecraft\cite{vallado2001fundamentals}, $k$, we
can take advantage of the fact that that area is the triangle $k_1$ subtracted from the
elliptical segment $PCB$:
\begin{equation}\label{areas_eq}
k = area(seg_{PCB}) - area(k_1)
\end{equation}
\noindent
Where the area of the triangle $k_1$ can be found easily using geometric formulae:
\begin{align}
area(k_1) &= \frac{1}{2} \left( ae - a \cos E \right) \left( \frac{b}{a} a \sin E \right) \\
&= \frac{ab}{2} \left(e \sin E - \cos E \sin E \right)
@@ -186,7 +170,6 @@
can find the area for the elliptical segment $PCB$ by first finding the circular segment
$POB'$, subtracting the triangle $COB'$, then applying the fact that an ellipse is
merely a vertical scaling of a circle by the amount $\frac{b}{a}$.
\begin{align}
area(PCB) &= \frac{b}{a} \left( area(POB') - area(COB') \right) \\
&= \frac{b}{a} \left( \frac{a^2 E}{2} - \frac{1}{2} \left( a \cos E \right)
@@ -197,26 +180,20 @@
By substituting the two areas back into Equation~\ref{areas_eq} we can get the $k$ area
swept out by the spacecraft:
\begin{equation}
k = \frac{ab}{2} \left( E - e \sin E \right)
\end{equation}
Which we can then substitute back into the equation for the swept area as a function of
time (Equation~\ref{swept}) for period of time since the spacecraft left periapsis:
\begin{equation}
\frac{\Delta t}{T} = \frac{t_2 - t_{peri}}{T} = \frac{E - e \sin E}{2 \pi}
\end{equation}
Which is, effectively, Kepler's equation. It is commonly known by a different form:
\begin{equation}
M = \sqrt{\frac{\mu}{a^3}} \Delta t = E - e \sin E
\end{equation}
Where we've defined the mean anomaly as $M$ and used the fact that $T =
\sqrt{\frac{a^3}{\mu}}$. This provides us a useful relationship between Eccentric Anomaly
where we've defined the mean anomaly as $M$ and used the fact that $T =
\sqrt{\frac{a^3}{\mu}}$. This provides us a useful relationship between eccentric anomaly
($E$) which can be related to spacecraft position, and time, but we still need a useful
algorithm for solving this equation in order to use this equation to propagate a
spacecraft.
@@ -224,34 +201,25 @@
\subsection{LaGuerre-Conway Algorithm}\label{laguerre}
For this thesis, the algorithm used to solve Kepler's equation was the general numeric
root-finding scheme first developed by LaGuerre in the 1800s and first applied to
Kepler's equation by Bruce Conway in 1985\cite{laguerre_conway}. In his paper, Conway
makes a compelling argument for utilizing the less common LaGuerre method over higher
order Newton or Newton-Raphson methods.
The Newton-Raphson methods, while found to generally have quite impressive convergence
rates (generally successfully solving Kepler's equation correctly within 5 iterations),
were prone to failures in convergence given certain specific initial conditions.
Therefore LaGuerre's algorithm is proposed as an alternative.
The algorithm can be relatively easily derived by examining the polynomial equation with
$m$ roots:
root-finding scheme first developed by LaGuerre in the 1800s and first applied to Kepler's
equation by Bruce Conway in 1985\cite{laguerre_conway}. In his paper, Conway makes a
compelling argument for utilizing the less common LaGuerre method over higher order Newton
or Newton-Raphson methods. The Newton-Raphson methods, while found to generally have quite
impressive convergence rates (generally successfully solving Kepler's equation correctly
within 5 iterations), were prone to failures in convergence given certain specific initial
conditions. Therefore LaGuerre's algorithm is proposed as an alternative.
The algorithm can be derived by examining the polynomial equation with $m$ roots:
\begin{equation}
g(x) = (x - x_1) (x - x_2) ... ( x - x_m)
\end{equation}
\noindent
We can then generate some useful convenience functions as:
\begin{align}
\ln|g(x)| &= \ln|(x - x_1)| + \ln|(x - x_2)| + ... + \ln|( x - x_m)| \\
\frac{d\ln|g(x)|}{dx} &= \frac{1}{x - x_1} + \frac{1}{x - x_2} + ... + \frac{1}{x -
x_m} = G_1(x)
\end{align}
and
\begin{align}
\frac{-d^2\ln|g(x)|}{dx^2} &= \frac{1}{(x - x_1)^2} + \frac{1}{(x - x_2)^2} + ... +
\frac{1}{(x - x_m)^2} = G_2(x)
@@ -259,42 +227,32 @@
Now we define the targeted root as $x_1$ and make the approximation that all of the
other roots are equidistant from the targeted root, which means:
\begin{equation}
x - x_i = b, i=2,3,...,m
\end{equation}
\noindent
We can then rewrite $G_1$ and $G_2$ as:
\begin{align}
G_1 &= \frac{1}{a} + \frac{n-1}{b} \\
G_2 &= \frac{1}{a^2} + \frac{n-1}{b^2}
\end{align}
\noindent
Which may be solved for $a$ in terms of $G_1$, $G_2$:
\begin{equation}
a = \frac{n}{G_1 \pm \sqrt{(n-1)(nG_2 - G_1^2)}}
\end{equation}
\noindent
With corresponding iteration function:
\begin{equation}
x_{i+1} = x_i - \frac{n g(x_i)}{g'(x_i) \pm \sqrt{(n-1)^2 f'(x_i)^2 - n (n-1) f(x_i)
f''(x_i)}}
\end{equation}
This iteration scheme can be shown to be globally convergent, regardless of the initial
guess. More relevantly, Conway also showed that the application of this method to
Kepler's equation was shown to converge with similar speed to many of the best common
higher order Newton-Raphson solvers. However, LaGuerre's method was also found to be
incredibly robust, converging to the correct value for every one of Conway's 500,000
tests. Because of this robustness, it is very useful for propagating spacecraft states.
guess. Conway also showed that the application of this method to Kepler's equation was shown
to converge with similar speed to many of the best common higher order Newton-Raphson
solvers. However, LaGuerre's method was also found to be incredibly robust, converging to
the correct value for every one of Conway's 500,000 tests. Because of this robustness, it is
useful for solving Kepler's equation.
\section{Interplanetary Considerations}\label{interplanetary}
\section{Interplanetary Trajectories}\label{interplanetary}
In interplanetary travel, the primary body most responsible for gravitational forces might
be a number of different bodies, dependent on the phase of the mission. In fact, at some
@@ -346,14 +304,15 @@
This effectively breaks the trajectory into a series of arcs each governed by a distinct
Two-Body problem patched together by distinct transition points. These transition points
occur along the spheres of influence of the planets nearest to the spacecraft.
occur along the spheres of influence of the planets nearest to the spacecraft. A
conceptual example of this process, labeled the method of patched conics, appears in
Figure~\ref{patched_conics_fig}.
Therefore, we must understand how to convert our spacecraft's state from the Sun frame
to the planetary frame as it crosses this boundary. An elliptical orbit about the sun
will have enough orbital energy to represent a hyperbolic orbit around the planet. So we
first need to determine the velocity of the spacecraft relative to the planet as it
crosses the SOI, which we can determine by subtraction \cite{vallado2001fundamentals}:
\begin{equation}
\vec{v}_{sc/p} = \vec{v}_{sc/sun} - \vec{v}_{planet/sun}
\end{equation}
@@ -361,8 +320,8 @@
Since the orbit around the planet is hyperbolic, in order to characterize the hyperbola
we must determine the velocity of the spacecraft when it has infinite distance relative
to the planet. Since this never occurs, a further approximation is made that the
velocity that the spacecraft has (relative to the planet) as it crosses the SOI can be
modeled as the $\vec{v}_\infty$ of that hyperbolic arc.
velocity of the spacecraft (relative to the planet) as it crosses the SOI can be modeled
as the $\vec{v}_\infty$ of that hyperbolic arc.
As an example, we may wish to determine the velocity relative to the planet that the
spacecraft has at the periapsis of its hyperbolic trajectory during the flyby. This
@@ -371,14 +330,12 @@
around its target planet. For a given incoming hyperbolic $\vec{v}_\infty$, we can first
determine the specific mechanical energy of the hyperbola at infinite distance by using
Equation~\ref{energy}:
\begin{equation}
\xi = \frac{v^2}{2} - \frac{\mu}{r} = \frac{v_\infty^2}{2}
\end{equation}
We can then leverage the conservation of energy to determine the velocity at a
particular point, $r_{ins}$:
\begin{align}
\xi_{ins} &= \frac{v_{ins}^2}{2} - \frac{\mu}{r_{ins}} \\
\xi_{ins} &= \xi_\infty = \frac{v_\infty^2}{2} \\
@@ -387,14 +344,13 @@
\subsection{Launch Considerations}
Generally speaking, an interplanetary mission begins with launch. For a satellite of
given size, a certain amount of orbital energy can be imparted to the satellite by the
launch vehicle. In practice, this value, for a particular mission, is actually
determined as a parameter of the mission trajectory to be optimized. The excess velocity
at infinity of the hyperbolic orbit of the spacecraft that leaves the Earth can be used
to derive the launch energy. This is usually qualified as the quantity $C_3$, which is
actually double the kinetic orbital energy with respect to the Sun, or simply the square
of the excess hyperbolic velocity at infinity\cite{wie1998space}.
For a satellite of given size, a certain amount of orbital energy can be imparted to the
satellite by the launch vehicle. In practice, this value, for a particular mission, is
actually determined as a parameter of the mission trajectory to be optimized. The excess
velocity at infinity of the hyperbolic orbit of the spacecraft that leaves the Earth can
be used to derive the launch energy. This is usually qualified as the quantity $C_3$,
which is actually double the kinetic orbital energy with respect to the Sun, or simply
the square of the excess hyperbolic velocity at infinity\cite{wie1998space}.
This algorithm will assume that the initial trajectory at the beginning of the mission
will be some hyperbolic orbit with velocity enough to leave the Earth. That initial
@@ -405,12 +361,12 @@
what the maximum mass any launch provider is capable of imparting that specific $C_3$
to.
A similar approach is taken at the end of the mission. This algorithm doesn't attempt to
exactly match the velocity of the planet. Instead, the excess hyperbolic velocity is
also treated as a parameter that can be minimized by the cost function. If a mission is
to then end in insertion, a portion of the mass budget can then be used for an impulsive
thrust engine, which can provide a final insertion burn. This approach also allows
flexibility for missions that might end in a flyby rather than insertion.
A similar approach is taken at the end of the trajectory. This algorithm doesn't attempt
to exactly match the velocity of the planet. Instead, the excess hyperbolic velocity is
also treated as a parameter that can be minimized by the cost function. If a trajectory
is to then end in insertion, a portion of the mass budget can then be used for an
impulsive thrust engine, which can provide a final insertion burn. This approach also
allows flexibility for missions that might end in a flyby rather than insertion.
\subsection{Gravity Assist Maneuvers}
@@ -441,7 +397,7 @@
\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{fig/flyby}
\caption{Visualization of velocity changes during a gravity assist}
\caption{Velocity changes during a gravity assist}
\label{grav_assist_fig}
\end{figure}
@@ -451,7 +407,7 @@
turning angle of this bend. In doing so, one can effectively achieve a (restricted) free
impulsive thrust event.
\subsection{Flyby Periapsis}
\subsection{Flyby Periapsis Altitude}
Now that we understand gravity assists, the natural question is then how to leverage
them for achieving certain velocity changes\cite{cho2017b}. But first, we must consider
@@ -460,7 +416,6 @@
mentioned in the previous section, given an excess hyperbolic velocity entering the
planet's sphere of influence ($\vec{v}_{\infty, in}$) and a target excess hyperbolic
velocity as the spacecraft leaves the sphere of influence ($\vec{v}_{\infty, out}$):
\begin{equation}\label{turning_angle_eq}
\delta = \arccos \left( \frac{\vec{v}_{\infty,in} \cdot
\vec{v}_{\infty,out}}{|\vec{v}_{\infty,in}| |\vec{v}_{\infty,out}|} \right)
@@ -470,12 +425,10 @@
that we must target in order to achieve the required turning angle. The periapsis of the
flyby, however, can provide a useful check on what turning angles are possible for a
given flyby, since the periapsis:
\begin{equation}\label{periapsis_eq}
r_p = \frac{\mu}{v_\infty^2} \left[ \frac{1}{\sin\left(\frac{\delta}{2}\right)} - 1 \right]
\end{equation}
Cannot be lower than some safe value that accounts for the radius of the planet and
cannot be lower than some safe value that accounts for the radius of the planet and
perhaps its atmosphere if applicable.
\subsection{Multiple Gravity Assist Techniques}
@@ -511,7 +464,9 @@
less than 180 degrees, which we classify as a Type I trajectory, and the second will
have a $\Delta \theta$ of greater than 180 degrees, which we call a Type II
trajectory. They will also differ in their direction of motion (clockwise or
counter-clockwise about the focus). This can be seen in Figure~\ref{type1type2}.
counter-clockwise about the focus). This can be seen in Figure~\ref{type1type2},
where both of the Lambert's solutions are presented for sample points in an orbit
around the Sun.
\begin{figure}[H]
\centering
@@ -523,7 +478,6 @@
The iteration used in this thesis will start by first calculating the change in true
anomaly, $\Delta \theta$, as well as the cosine of this value, which can be found
by:
\begin{align}
\cos (\Delta \theta) &= \frac{\vec{r}_1 \cdot \vec{r}_2}{|\vec{r}_1| |\vec{r}_2|} \\
\Delta \theta &= \arctan(y_2/x_2) - \arctan(y_1/x_1)
@@ -532,7 +486,6 @@
The direction of motion is then chosen such that counter-clockwise orbits are
considered, as travelling in the same direction as the planets is generally more
efficient. Next, the variable $A$ is defined:
\begin{equation}
A = DM \sqrt{|r_1| |r_2| (1 - \cos(\Delta \theta))}
\end{equation}
@@ -547,7 +500,6 @@
time of flight matches the expected value to within a provided tolerance. In order
to calculate the time of flight at each step, we must first calculate some useful
coefficients:
\begin{equation}\label{loop_start}
c_2 = \begin{cases}
\frac{1-\cos(\sqrt{\psi})}{\psi} \quad &\text{if} \, \psi > 10^{-6} \\
@@ -555,7 +507,6 @@
1/2 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6}
\end{cases}
\end{equation}
\begin{equation}
c_3 = \begin{cases}
\frac{\sqrt{\psi} - \sin \sqrt{\psi}}{\psi^{3/2}} \quad &\text{if} \, \psi > 10^{-6} \\
@@ -563,23 +514,18 @@
1/6 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6}
\end{cases}
\end{equation}
\noindent
Where the conditions of this piecewise function represent the elliptical,
hyperbolic, and parabolic cases, respectively. Once we have these, we can calculate
another variable, $y$:
\begin{equation}
y = |r_1| + |r_2| + \frac{A (c_3 \psi - 1)}{\sqrt{c_2}}
\end{equation}
We can then finally calculate the variable $\chi$, and from that, the time of
flight:
\begin{equation}
\chi = \sqrt{\frac{y}{c_2}}
\end{equation}
\begin{equation}
\Delta t = \frac{c_3 \chi^3 + A \sqrt{y}}{\sqrt{c_2}}
\end{equation}
@@ -592,22 +538,17 @@
The resulting $f$ and $g$ functions (and the derivative of $g$, $\dot{g}$) can then
be calculated:
\begin{align}
f &= 1 - \frac{y}{|r_1|} \\
g &= A \sqrt{\frac{y}{\mu}} \\
\dot{g} &= 1 - \frac{y}{|r_2|}
\end{align}
And from these, we can calculate the velocities of the transfer points as:
\begin{align}
\vec{v}_1 &= \frac{\vec{r}_1 - f \vec{r}_2}{g} \\
\vec{v}_2 &= \frac{\dot{g} \vec{r}_2 - \vec{r}_1}{g}
\end{align}
\noindent
Fully constraining the connecting orbit.
Fully describing the connecting path with the specified flight time.
\subsubsection{Planetary Ephemeris}
@@ -620,8 +561,8 @@
The primary use of SPICE in this thesis, however, was to determine the planetary
ephemeris at a known epoch. Using the NAIF0012 and DE430 kernels, ephemeris in the
ecliptic plane J2000 frame (ICRF) could be easily determined for a given epoch, provided as
a decimal Julian Day since the J2000 epoch.
International Celestial Reference Frame could be easily determined for a given
epoch, provided as a decimal Julian Day since the J2000 epoch.
\subsubsection{Porkchop Plots}
@@ -641,10 +582,9 @@
Using porkchop plots such as the one in Figure~\ref{porkchop}, mission designers can
quickly visualize which natural trajectories are possible between planets. Using the
fact that incoming and outgoing $v_\infty$ magnitudes must be the same for a flyby,
a savvy mission designer can even begin to work out what combinations of flybys
might be possible for a given timeline, spacecraft state, and planet selection.
a mission designer can even begin to work out what combinations of flybys might be
possible for a given timeline, spacecraft state, and planet selection.
%TODO: Create my own porkchop plot
\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{fig/porkchop}
@@ -652,13 +592,7 @@
\label{porkchop}
\end{figure}
However, this is an impulsive thrust-centered approach. The solution to Lambert's
problem assumes a natural trajectory. A natural trajectory is unnecessary when the
trajectory can be modified by a continuous thrust profile along the arc. Therefore,
for the hybrid problem of optimizing both flyby selection and thrust profiles,
porkchop plots are less helpful, and an algorithmic approach is preferred.
\section{Low Thrust Considerations} \label{low_thrust}
\section{Modeling Low Thrust Control} \label{low_thrust}
In this section, we'll discuss the intricacies of continuous low-thrust trajectories in
particular. There are many methods for optimizing such profiles and we'll briefly discuss
@@ -666,7 +600,7 @@
as introduce the concept of a control law and the notation used in this thesis for modelling
low-thrust trajectories more simply.
\subsection{Specific Impulse}
\subsection{Engine Model}
The primary advantage of continuous thrust methods over their impulsive counterparts is
in their fuel-efficiency in generating changes in velocity. Put specifically, all
@@ -678,45 +612,34 @@
This efficiency is often captured in a single variable called specific impulse, often
denoted as $I_{sp}$. We can derive the specific impulse by starting with the rocket
thrust equation\cite{sutton2016rocket}:
\begin{equation}
F = \dot{m} v_e + \Delta p A_e
\end{equation}
\noindent
Where $F$ is the thrust imparted, $\dot{m}$ is the fuel mass rate, $v_e$ is the exhaust
velocity of the fuel, $\Delta p$ is the change in pressure across the exhaust opening,
and $A_e$ is the area of the exhaust opening. We can then define a new variable
$v_{eq}$, such that the thrust equation becomes:
\begin{align}
v_{eq} &= v_e + \frac{\Delta p A_e}{\dot{m}} \\
F &= \dot{m} v_{eq} \label{isp_1}
\end{align}
\noindent
And we can then take the integral of this value with respect to time to find the total
impulse, dividing by the weight of the fuel to derive the specific impulse:
\begin{align}
I &= \int F dt = \int \dot{m} v_{eq} dt = m_e v_{eq} \\
I_{sp} &= \frac{I}{m_e g_0} = \frac{m_e v_{eq}}{m_e g_0} = \frac{v_{eq}}{g_0}
\end{align}
Plugging Equation~\ref{isp_1} into the previous equation we can derive the following
formula for $I_{sp}$:
\begin{equation} \label{isp_real}
I_{sp} = \frac{F}{\dot{m} g_0}
\end{equation}
\noindent
Which is generally taken to be a value with units of seconds and effectively represents
the efficiency with which a thruster converts mass to thrust.
\subsection{Sims-Flanagan Transcription}
This thesis chose to use a model well suited for modeling low-thrust paths: the
In this thesis the following approach is used for modeling low-thrust paths: the
Sims-Flanagan transcription (SFT)\cite{sims1999preliminary}. The SFT allows for
flexibility in the trade-off between fidelity and performance, which makes it very
useful for this sort of preliminary analysis.
@@ -752,7 +675,7 @@
continuous low-thrust trajectory within the Two-Body Problem, with only
linearly-increasing computation time\cite{sims1999preliminary}.
\subsection{Low-Thrust Control Laws}
\subsection{Low-Thrust Control Vector Description}
In determining a low-thrust arc, a number of variables must be accounted for and,
ideally, optimized. Generally speaking, this means that a control law must be determined
@@ -765,24 +688,23 @@
The methods for determining this direction varies greatly depending on the particular
control law chosen for that mission. Often, this process involves first determining a
useful frame to think about the kinematics of the spacecraft. In this case, we'll use a
frame often used in these low-thrust control laws: the spacecraft $\hat{R} \hat{\theta}
\hat{H}$ frame. In this frame, the $\hat{R}$ direction is the radial direction from the
center of the primary to the center of the spacecraft. The $\hat{H}$ hat is
perpendicular to this, in the direction of orbital momentum (out-of-plane) and the
$\hat{\theta}$ direction completes the right-handed orthonormal frame.
frame often used in these low-thrust control laws: the spacecraft-centered $\hat{R}
\hat{\theta} \hat{H}$ frame. In this frame, the $\hat{R}$ direction is the radial
direction from the center of the primary to the center of the spacecraft. The $\hat{H}$
hat is perpendicular to this, in the direction of orbital momentum (out-of-plane) and
the $\hat{\theta}$ direction completes the right-handed orthonormal triad.
This frame is useful because, for a given orbit, especially a nearly circular one, the
$\hat{\theta}$ direction is nearly aligned with the velocity direction for that orbit at
that moment. This allows us to define a set of two angles, which we'll call $\alpha$ and
$\beta$, to represent the in and out of plane pointing direction of the thrusters. This
convention is useful because a $(0,0)$ set represents a thrust force more or less
directly in line with the direction of the velocity, a commonly useful thrusting
direction for most effectively increasing (or decreasing if negative) the angular
momentum and orbital energy of the trajectory.
convention is useful because, in a near-circular path, a $(0,0)$ set represents a thrust
force more or less directly in line with the direction of the velocity, a commonly
useful thrusting direction for most effectively increasing (or decreasing if negative)
the angular momentum and orbital energy of the trajectory.
Using these conventions, we can then redefine our thrust vector in terms of the $\alpha$
and $\beta$ angles in the chosen frame:
\begin{align}
F_r &= F \cos(\beta) \sin (\alpha) \\
F_\theta &= F \cos(\beta) \cos (\alpha) \\
@@ -791,12 +713,12 @@
\subsubsection{Thrust Magnitude}
However, there is actually another variable that can be varied by the majority of
electric thrusters. Either by controlling the input power of the thruster or the duty
cycle, the thrust magnitude can also be varied, limited by the maximum thrust available
to the thruster. Not all control laws allow for this fine-tuned control of the thruster.
There is another variable that can be varied by the majority of electric thrusters.
Either by controlling the input power of the thruster or the duty cycle, the thrust
magnitude can also be varied, limited by the maximum thrust available to the thruster.
Not all control laws allow for this fine-tuned control of the thruster.
The algorithm used in this thesis does vary the magnitude of the thrust control. In
The approach used in this thesis does vary the magnitude of the thrust control. In
certain cases it actually can be useful to have some fine-tuned control over the
magnitude of the thrust. Since the optimization in this algorithm is automatic, it is
relatively straightforward to consider the control thrust as a 3-dimensional vector in