# Rankine Cycle MATLAB Software

Rankine cycles are used to produce power from coal and nuclear fission power plants. They also serve as the bottoming cycle for combined cycle power plants. In the future, the Rankine cycle may be used to produce power from nuclear fusion power plants.

Princeton Satellite Systems has added Rankine cycle functions to release 2023.1 of its MATLAB toolboxes, which is coming soon. The following figures show the output from the function RankineCycle.m. The first diagram is the temperature entropy (T-s) diagram. Entropy can be viewed as the state of disorder in the system.

The red lines indicate the simple Rankine cycle while the bell-shaped curve is the T-s curve for water and steam. The second output is the cycle diagram, shown below.

Entropies are displayed on this diagram. You only need two pressures to define the cycle: the base and the peak. The pump compresses water which is fed to the boiler. The boiler produces steam that is converted to mechanical energy in the steam turbine. This drives a generator to produce electric power. The fluid output of the turbine passes through the condenser to become water.

This is the simplest Rankine cycle. The numbers on the diagram correspond to the numbers on the T-s diagram. The cycle is defined per unit mass of water flowing in the cycle so it can be scaled to any size. The MATLAB code has functions for handling steam tables, which are essential to studying or designing steam engines.

The Rankine cycle can be improved by adding reheat or regeneration. Organic liquids can also be used in place of water.

# Aerospace Engineering Winter Internship

I worked on two projects during my winter internship at Princeton Satellite Systems: a two-stage-to-orbit (TSTO) launch vehicle design proposal related to the NASA Space Launch System (SLS) and a satellite conjunction maneuver demo. These both used the Spacecraft Control Toolbox for MATLAB.

One of the main ideas behind the TSTO launch vehicle project is to propose an all-liquid variant of the SLS. Currently, the SLS first stage is mostly powered by two solid rocket boosters (SRB) upgraded from the Space Shuttle SRBs. However, our proposal is to replace the two SRBs with five liquid boosters (LB), each mated with an RS-25 engine. The second stage would remain the same. Using MATLAB, I analyzed the launch and trajectory performance of both variants and found similar performance. Additionally, the total mass of the all-liquid SLS variant would be approximately two-thirds the mass of the SRB-powered spacecraft. An approximate CAD model of the all-liquid SLS version is shown below.

In addition, the LBs can be used independently to power smaller high-performance TSTO launch vehicles that carry around 8,000 kg of payload to low earth orbit. Trajectory plots and a preliminary CAD model are shown below.

My other project this internship was to help out with a satellite conjunction avoidance demo with Ms. Stephanie Thomas. The goal was to create a solution in MATLAB to identify potential satellite-debris conjunctions and develop a method/algorithm to avoid the conjunctions. I mainly worked on testing the code and relevant functions and providing feedback about the solution’s comprehensiveness.

Overall, I greatly enjoyed this internship and the opportunity to work at PSS. I saw firsthand how even a small company can make significant contributions to aerospace and engineering through diverse interests yet specific, impressive skill sets.

# Doing the Mars run with fusion propulsion at 1 G

We received a comment on LinkedIn about how fast the “Mars run” could be achieved with a sustained 1 G acceleration. The reader suggested this could be done in 40 hours. What engine parameters would be required to make that happen?

Using a simple constant-acceleration, straight-line analysis, you can indeed compute that the trip should take only a couple of days. Assuming a Mars conjunction, the straight distance is about 0.5 AU. At this speed you can ignore the gravitational effects of the sun and so the distance is a simple integral of the acceleration: d = 1/2 at2. The ship accelerates for half the time then decelerates, and the change in velocity is ΔV = at. Combining the two halves of the trip, at an acceleration of 9.8 m/s2, the trip takes about 2.1 days.

% straight line: distance s = 0.5*at^2
acc = 9.8;             % accel, m/s^2
aU  = Constant('au');  % km
dF  = 0.5*aU*1000;     % distance, m
t   = sqrt(4*dF/acc);  % time for dF, s
dV  = t*acc/1000; % km/s
fprintf('\nAccel: %g m/s^2\n',acc)
fprintf('Time: %g days\n',t/86400)
fprintf('Delta-V: %g km/s\n',dV)

Accel: 9.8 m/s^2
Time: 2.02232 days
Delta-V: 1712.34 km/s

d = StraightLineDataStructure;
d.dF = 0.5*aU;   % 0.5 AU
d.tF = 2*86400;  % 2.1 days
d.uE = 5000;     % km/s
d.eta = 1;       % jet power sigma
d.sigma = 1e8;   % W/kg
d.mP = 50000;    % kg
dOut = MassFuelElectricConstantUE( d )
StraightLineConstantAccel(dOut)
fprintf('Acceleration: %g m/s^2\n',dOut.a*1e3)
fprintf('Exhaust velocity: %g km/s\n',d.uE)
fprintf('Initial power: %g GW\n',dOut.p(1)*1e-9)
fprintf('Engine mass: %g kg\n',dOut.mE)

Acceleration: 10.02 m/s^2
Exhaust velocity: 5000 km/s
Specific Power: 1e+08 W/kg
Initial power: 2832.61 GW
Engine mass: 28326.1 kg
Payload fraction: 0.442172

This produces the following plots:

The power needed is… over 2.8 terawatts! That’s about equal to the total power output of the entire Earth, which had an installed power capacity of 2.8 terawatts in 2020. And the engine would need to weigh less than 30 tons, about the size of a loaded tractor-trailer truck. For comparison, we estimate a Direct Fusion Drive would produce about 1 MW per ton, which is a specific power of 1×103 W/kg. So, this is why you see us trying to design an engine that can do the Mars transfer in 90 days and not 3 days!

Now, there is another consideration here. Namely, constant acceleration at 1 G is not the optimal solution by any means. The optimal solution for a fast, light transfer is actually a linear acceleration profile. This knowledge goes way back: 1961! Here’s a reference:

Leitmann, George. "Minimum Transfer Time for a Power-Limited Rocket." Journal of Applied Mechanics 28, no. 2 (June 1, 1961): 171-78. https://doi.org/10.1115/1.3641648.

This would mean that the engine changes its exhaust velocity during trip, passing through infinity at the switch point. We compute this in our “straight-line, power-limited” or SLPL function series. While this can’t be done physically, even an approximation of this with a variable impulse thruster will one day be more efficient than constant acceleration or thrust. How much better? The power needed is nearly 1/2 the constant acceleration solution, 1.5 TW, and the specific power needed is reduced by half, to 5.6×107 W/kg. However, those are still insane numbers!

mD = 80000;  % dry mass: engine, tanks, payload
m0 = 1.5*mD; % wet mass: with fuel
tF = 3*86400;
vF = 0;

[Pj,A,tau] = SLPLFindPower( aU, tF, vF, mD, m0 );

mTank = 0.05*(m0-mD); % tanks, scale with fuel
mLeft = mD-mTank;

disp('Straight-line Power-limited (linear accel)')
fprintf('Engine power is %g GW\n',Pj*1e-9);
fprintf('Engine mass is %g kg\n',mEngine);
fprintf('sigma is %g W/kg\n',Pj/mEngine);

SLPLTrajectory( A, tau, Pj, m0, tF )

Straight-line Power-limited (linear accel)
Engine power is 1573.26 GW
Engine mass is 28000 kg
sigma is 5.6188e+07 W/kg

The trajectory and engine output are plotted below. The linear acceleration results in a curved velocity plot, while in the constant acceleration case, we saw a linear velocity plot. You can see the spike in exhaust velocity at the switch point, which occurs exactly at the halfway point.

After all, who needs 1G gravity when the trip only takes 2 days?

For even more fun though, we computed a planar trajectory to Mars using the parameters we found – just to confirm the straight-line analysis is in fact a good approximation. This figure shows the paths the optimization takes:

It is in fact approximately a straight line!

In reality though, these power system numbers are not even remotely plausible with any technology we are aware of today. That’s why we are designing engines to reduce the Mars trip time to 90 days from 8 or 9 months – still a big improvement!

# Hohmann Transfer Simulation with the Spacecraft Control Toolbox

Hohmann transfers are a well-known maneuver used to change the semi-major axis of an orbit. The Spacecraft Control Toolbox allows you to compute the required velocity changes, and integrate them into a full simulation.

In this demonstration, we create a 6U CubeSat that has 3 orthogonal reaction wheels and a single hydrazine thruster. The thruster is aligned with the body x-axis and must be aligned with the velocity vector to do the maneuver. An ideal Hohmann maneuver is done with impulsive burns at two points in the orbit. In reality, with a thruster, we have to do a finite burn.

The Hohmann transfer is computed with the following Spacecraft Control Toolbox code:

rI = [-7000;0;0];
vI = [0;-sqrt(mu/Mag(rI));0];
OrbMnvrHohmann(Mag(rI),rF);
[dV,tOF] = OrbMnvrHohmann(Mag(rI),rF);

The first time OrbMnvrHohmann is called, it generates the plot below of the planned Hohmann transfer. The function computes the delta-V and also the time of flight, which will be used to determine the start time of the second thruster burn.

We create a short script with numerical integration to implement the maneuver using a thruster. The burn durations are computed based on the thrust and the mass of the spacecraft. In this case, they are about three minutes long. The maneuver is quite small, so the mass change is not important. The attitude control system uses the PID3Axis function which is a general-purpose attitude control algorithm. The simulation is a for loop, shown below. The ECI vector for the burn is passed to the attitude control system, which updates every step of the simulation.

% Simulation loop
for k = 1:n

% Update the controller
dC.eci_vector = uBurn(:,kMnvr);
[tRWA, dC]    = PID3Axis( x(7:10), dC );

% Start the first burn
inMnvr = false;
if( t(k) > tStart(1) && t(k) < tEnd(1) )
inMnvr = true;
end

% Switch orientation
if( t(k) > tEnd(1) )
kMnvr = 2;
end

% Start the second burn
if( t(k) > tStart(2) && t(k) < tEnd(2) )
inMnvr = true;
end

if( inMnvr )
dRHS.force = thrustE*QTForm(x(7:10),dC.body_vector)*nToKN; % kN
else
dRHS.force = [0;0;0];
end
el = RV2El(x(1:3),x(4:6));
xP(:,k) = [x;tRWA;Mag(dRHS.force)/nToKN;el(1);el(5)];

% Right hand side
dRHS.torqueRWA = -tRWA;
x = RK4(@RHSRWAOrbit,x,dT,0,dRHS);
end

The maneuver logic just waits a quarter orbit then performs the first burn, by applying the thrust along the body vector. It then waits for the time of flight and then starts the next burn. The start and stop times are pre-computed. RK4 is Fourth Order Runge-Kutta, a popular numerical algorithm included with the toolbox.

At the final orbit radius an attitude maneuver is needed to reorient for the final burn.

The spacecraft body rates, in the body frame, during the maneuver are shown below.

The reaction wheel rates are shown below. The simulation does not model any particular wheel. Friction is not included in the simulation, although the right-hand-side function can include friction.

The wheel torques and rocket thrust are shown below. The thruster is a 0.2 lbf hydrazine thruster that is based on the Aerojet-Rocketdyne MR-103. The PID controller does not demand much torque.

The semi-major axis and eccentricity are shown below. The middle portion is during the transfer orbit.

The eccentricity is zero at the start and finish. Note the slope in both eccentricity and semi-major axis due to the finite acceleration. At the end of the simulation, we print the achieved orbital elements:

Final SMA        7099.72 km
SMA error         0.28 km
Final e          1.3e-05

The result is very close to the ideal solution!

This post shows how you can easily integrate attitude and orbit control. Email us for more information! We’d be happy to share the script. We can also offer a 30 day demo to let you explore the software.

# A Third Planet Discovered Orbiting Proxima Centauri

## Introduction

A third planet, as large as 26% of the mass of Earth, has been discovered orbiting our nearest stellar neighbor, Proxima Centauri .Astronomer João Faria and his collaborators detected Proxima Centauri d using the Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations.

It would be exciting to send a spacecraft to enter the Alpha-Centauri system and orbit this planet. At Princeton Satellite System we’ve looked at interstellar flight using the Direct Fusion Drive nuclear fusion propulsion system.

## Interstellar Fusion Propulsion

At the 2021 Breakthrough Energy Conference we presented findings for both flyby and orbital missions. Flyby missions are easier, but orbit entry would allow detailed study of the planet. A flyby gets your spacecraft close, but it is moving really fast!

The following charts give an outline of our talk. The first shows the optimal exhaust velocity based on sigma, the ratio of power to mass. Our designs have a sigma from 0.75 to 2 kW/kg. With 2 kW/kg, the optimal exhaust velocity is 4000 km/s. The mission would take about 800 years. Our current designs can’t get exhaust velocities higher than 200 km/s. We’d need another method to produce thrust.

The next plot shows a point mission that reaches Alpha Centauri in 500 years. This requires a sigma of about 20. The spacecraft accelerates and decelerates continuously. The mission could be improved by staging, much like on a rocket that launches from the Earth into orbit.

The next figure shows how the starship would enter the Alpha Centauri system.

The final plot shows the orbital maneuvers that lower the orbit and rendezvous with the planet.

Even 500 years is a long time! This is over ten times the lifetime of Voyager, but much less than some engineering marvels built on the Earth.

We hope to someday be able to build fusion powered spacecraft that will head into interstellar space!

# Moonfall – The Movie

Moonfall is a movie coming out in 2022. It creates a scenario where the Moon’s orbit is changed and set on a collision course with the Earth. It is fun to work out the orbital mechanics.

Let us assume that the Moon is in a circular orbit around the Earth. It is actually more influenced by the Sun than the Earth, but the circular orbit approximation is sufficient for our purposes. A mysterious force changes the orbit from circular to elliptical so that at closest approach it hits the Earth. The transfer orbit has an eccentricity of 0.9673 and a semi-major axis of 195000 km. The new orbital period is 9.9 days so it will hit the Earth in 5 days!

What kind of force is needed? The required velocity change is 0.83 km/s so a force of 6 x 1016 N applied over 10 seconds is required. To get an idea of how large that force really is, the Space Launch System (SLS) Block 2 vehicle produces about 10 million pounds of thrust , which is approximately 50 x 106 N (50 MN). Hence it would take 1.2 billion SLS rockets firing for 10 seconds to perform such a re-direction of the Moon! An image of the SLS is shown below (image from ).

As the Moon approaches the Earth it is going to raise the tides. A simple formula (really only valid when the Moon is far from the Earth) is

where is the gravitational constant for the moon, is the gravitational constant for the Earth, r is the distance between the Earth and Moon and a is the radius of the Earth. The distance during the approach and the wave height are shown in the following plot.

By around 3 days the tides started getting really big! We’d expect the Moon’s gravitational force also to pull on the solid part of the Earth’s surface, causing all sorts of trouble.

References:

 https://www.nasa.gov/sites/default/files/atoms/files/0080_sls_fact_sheet_10092018.pdf

# PSS Toolboxes 2021.1 Now Available!

Version 2021.1 of Princeton Satellite Systems toolboxes for MATLAB is now available! Over 50 new functions and scripts are included. Many other existing functions have been improved.

One new function is AtmNRLMSISE.m, an atmosphere function based on the NRL MSISE model. It is uses extensive flight data and includes sun effects. It computes the overall density and the number density of all atmospheric constituents. Our function has an easy to use interface that automatically incorporates the sun information and lets you input your spacecrafts ECI coordinates. You can also choose to use the original interface. Here is a comparison with the existing scale height model.

We provide a complete set of functions for planning lunar missions in the Missions module. The software includes landing control systems and trajectory optimizaton tools. You can use our Optical Navigation system for your cis-lunar missions and explore our cutting-edge neural network terminal descent software.

Here are two images from an optical navigation simulation for a solar sail.

The Spacecraft Control Toolbox provides you with a lot of ways to do things, so you can use your own creativity to perform analyses or design a mission.

# The Optical Navigation Module for the Spacecraft Control Toolbox is Now Available

Space optical navigation employs a camera for attitude determination and a second high dynamic range camera on a pan/track mount for terrain and celestial body tracking. Navigation and attitude determination are performed in a Bayesian framework using anUnscented Kalman Filter with an IMU as the navigation and attitude base. The Optical Navigation Module provides MATLAB code for implementing optical navigation. Additional measurements can be added including a sun sensor for sun distance measurements in interplanetary space, Global Positioning System (GPS) measurements near the Earth, and range and range rate from ground stations or other spacecraft in deep space. The system is suitable for both lunar and Mars landing missions and icy moon and asteroid orbital missions such as Artemis, Lunar Orbital Platform Gateway, Orion Multi-Purpose Crew Vehicle, Europa Clipper, Lucy, Psyche. It is also applicable to any situation where GPS is not available.

The Optical Navigation Module allows you to implement an optical navigation system for any of these applications. It includes dynamical models for cis-lunar and deep space missions along with measurement models for all of these sensors. Several scripts provide examples to get you going quickly.

This picture shows the camera aimed at the horizon and the stars that it can see during Earth reentry. The step counter gives the integration step. The star numbers are sequential from the file of stars but the stars come from the Hipparcos catalog.

This pictures shows the laboratory hardware for an optical navigation camera on a pan/tilt mount. Flexible cables eliminate the need for slip rings simplifying the design. The platform is driven by orthogonal stepping motors with harmonic drives.

Note the size. As with all of our toolboxes, full source code is provided.

# What Makes a Reaction Wheel a Reaction Wheel?

A DC motor is the core of all momentum and reaction wheels. If you apply a voltage a, current will be produced which will cause the wheel to change speed. At the same time, the back electromotive force (EMF) will build up, eventually driving the motor torque to zero.

This is evident from the dynamical equation for a DC motor. $J\dot{\omega} = \frac{K_T}{R}\left(V - K_T\omega\right) + T_F$ $J$ is the inertia, $K_T$ is the torque constant, $V$ the voltage, $T_F$ the friction torque, $R$ the motor impedance and $\omega$ is the angular rate of the shaft.

You can turn this into a reaction wheel by adding current feedback as shown in the following block diagram. $G$ is the forward gain. The input is the desired torque. This is divided by the torque constant to get the desired current. The difference between the motor current and the desired current is integrated. How do you pick the gain? If you work through the equations you will get this equation for the voltage, $V$ $\dot{V} + \frac{G}{R}V = G\frac{T_C}{K_T} + \frac{G}{R}K_T\omega$ $R/G$ is the time constant. The response is shown in the following plot. Even as the speed increases, the difference between the desired torque and motor torque is nearly zero.

# Visiting Planet 9

In 2015, astronomers from Caltech determined that a giant ninth planet may be orbiting the Sun. It was called Planet X and then Planet 9. The discovery was based on perturbations in the orbits of TNOs, trans Neptunian Objects. The planet has about the mass of Neptune and is in a 10,000 to 20,000 year solar orbit. Jakub Scholtz of Durham University and James Unwin of University of Illinois at Chicago hypothesize that Planet 9 might be a black hole. The orbit of Planet 9 looks something like this.

We used a semi-major axis of 700 AU, an inclination of 30 degrees and an eccentricity of 0.6. The plot shows the full orbit of Planet 9, but the simulation only shows 150 years of the other planets.

It would be very interesting to visit Planet 9. One way is to use a solar sail. The sail would start on a trajectory aiming at perigee very close to the sun and then accelerate at high speed. Another approach is to use a spacecraft propelled by Direct Fusion Drive, a fusion propulsion system we’ve been working on for several years. A 26000 kg spacecraft with a 12 MW engine and 2000 kg of payload could rendezvous with Planet 9 (based on the above orbit) in just 11 years. This is the spacecraft trajectory

Direct Fusion Drive is based on the Princeton Field Reversed Configuration reactor invented by Dr. Samuel Cohen of the Princeton Plasma Physics Laboratory (PPPL). We have a experiment running at PPPL, funded by an ARPA-E OPEN grant, to perform critical ion-heating tests. Earlier work was funded by the NASA NIAC program. Hopefully we will be in a position to send a mission to Planet 9 in the not too distant future!