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.
Last week, PSS Mike Paluszek visited ITER, the international fusion research experiment under construction in France. In light of Mike’s recent visit to ITER, we wanted to showcase an application of our tokamak Fusion Reactor Design function to the design of ITER. This function is part of the Fusion Energy Toolbox for MATLAB, a toolbox that includes a variety of physics and engineering tools for designing fusion reactors and studying plasma physics. We will also compute design parameters for ITER’s successor, the DEMOnstration power plant (DEMO), a fusion reactor currently in the design phase which is planned to achieve net electricity output.
We first apply the Fusion Reactor Design function to ITER. Note that ITER is expected to produce 500 Megawatts (500 MW) of fusion power, but this will not be converted into electric power, the power that goes into the electrical grid. DEMO, on the other hand, is planned to produce 500 MW of electric power from 2000 MW of fusion power. The Fusion Reactor Design function asks for the net electric power output of the reactor, P_E, as an input, so we generate a value for P_E for ITER by using the same ratio of electric-to-fusion power as in DEMO, giving us a P_E of 125 MW for ITER. The inputs used for the ITER design are shown below (see references [1,2]), where we use a data structure “d_ITER”:
d_ITER.a = 2; % plasma minor radius (m)
d_ITER.B_max = 13; % maximum magnetic field at the coils (T)
d_ITER.P_E = 125; % electric power output of the reactor (MW)
d_ITER.P_W = 0.57; % neutron wall loading (MW/m^2)
d_ITER.H = 1; % H-mode enhancement factor
d_ITER.consts.eta_T = 0.25; % thermal conversion efficiency
d_ITER.consts.T_bar = 8; % average ion temperature (keV)
d_ITER.consts.k = 1.7; % plasma elongation
d_ITER.consts.f_RP = 0.25; % recirculating power fraction
The first five inputs were described in our original post on the Fusion Reactor Design function. The function can be called to perform a parameter sweep over any of these inputs. We also specify values for some constants: the thermal conversion efficiency ‘eta_T’, the average ion temperature ‘T_bar’, the plasma elongation ‘k’, which is a measure of how elliptical the plasma cross-section is, and the recirculating power fraction ‘f_RP’. We can perform a parameter sweep over the minor radius (from a = 1.8 meters to a = 2.2 meters, with 100 points in between) and display a table of results simply with two lines of code:
d_ITER = FusionReactorDesign(d_ITER,'a',1.8,2.2,100); % run function
d_ITER.parameters % show table of resulting parameters
Looking at the results table from d_ITER.parameters, we see overall agreement with parameters for ITER [1,2]. The plasma major radius (essentially the tokamak radius) R_0 output is about 5 m, which is in the ballpark of the 6.2 m radius of ITER design, and the magnetic field at R_0 (on plasma axis) output is 4.8 Tesla, close to the ITER design value of 5.3 Tesla. The plasma current output is 17.5 MegaAmps, which is also close to ITER’s design of 15 MegaAmps.
The Fusion Reactor Design function also outputs plots that show whether or not the reactor satisfies key operational constraints for tokamaks, see the figure below. The first three curves check various constraints to ensure the plasma is stable, which we see are met as they are located in the unshaded region (though the green curve is marginally close to the constraint boundary). The blue curve’s position deep into the shaded region indicates that the reactor is far from producing enough electric current to sustain itself. The designers of ITER anticipated this, which is why ITER will additionally use a pulsed inductive current and test a combination of other techniques to drive the plasma current.
We now consider DEMO, which is in the design phase with the goal of net electrical power output. Similarly to running the ITER case, we set up a data structure (now called ‘d_DEMO’) with known DEMO input parameters [3] and perform a parameter sweep over the minor radius ranging from a = 2.7 meters to a = 3.1 meters:
d_DEMO.a = 2.9; % plasma minor radius (m)
d_DEMO.B_max = 13; % maximum magnetic field at the coils (T)
d_DEMO.P_E = 500; % electric power output of the reactor (MW)
d_DEMO.P_W = 1.04; % neutron wall loading (MW/m^2)
d_DEMO.H = 0.98; % H-mode enhancement factor
d_DEMO.consts.eta_T = 0.25; % thermal conversion efficiency
d_DEMO.consts.T_bar = 12.5; % average ion temperature (keV)
d_DEMO.consts.k = 1.65; % plasma elongation
d_DEMO.consts.f_RP = 0.25; % recirculating power fraction
d_DEMO = FusionReactorDesign(d_DEMO,'a',2.7,3.1,100); % run function
d_DEMO.parameters % show table of resulting parameters
The outputs for the DEMO case also show overall agreement with DEMO parameters [3]. The plasma major radius R_0 output is 7.8 m, which is not far from the 9 m design radius for DEMO. The resulting on-axis magnetic field output is 6.2 T, close to the 5.9 T of the DEMO design. The plasma current output is now 21 MegaAmps, which is less than 20% away from the design value of 18 MegaAmps. It is important to note that in each of these parameters, we see an increase going from ITER to DEMO, which is consistent both in our model’s output and the actual design parameters in the papers [1-3].
The operational constraints plot for DEMO is shown in the figure below. DEMO is a larger reactor than ITER, and given the favorable scaling of tokamak operation with size, we expect improved results for operational constraints in DEMO. The three curves which check plasma stability are all satisfied. Unlike in the case of ITER which had the green curve close to the shaded region, the green curve in the case of DEMO stays safely in the unshaded region. The blue curve is still in the unshaded region, but much closer to the boundary of the unshaded region than ITER (now ~1.8, much closer to 1 than in the case of ITER which was ~4). This shows an improvement for DEMO compared to ITER as it is closer to producing enough self-sustaining plasma current, though it will still need some help from other current-generating techniques which will be tested on ITER.
This function is part of release 2022.1 of the Fusion Energy Toolbox. Contact us at info@psatellite.com or call us at +01 609 275-9606 for more information.
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
Now, your ship mass includes your payload, your engine, your fuel tanks and your fuel. Assume we want to move a payload of 50,000 kg, somewhat larger than the NASA Deep Space Habitat. The engine mass is computed using a parameter called the specific power, in units of W/kg. The fuel tank mass is scaled from the fuel mass, typically adding another 10%. When we run the numbers, we find that the engine needs to have a specific power of about 1×108 W/kg, and an exhaust velocity of about 5000 km/s results in the maximum payload fraction. We can compute the fuel mass and trajectory using our MassFuelElectricConstantUE and StraightLineConstantAccel toolbox functions:
Distance, Velocity and Fuel MassThrust, Power, and Exhaust Velocity
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;
mEngine = mLeft - mPayload;
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('Payload mass is %g kg\n',mPayload);
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
Payload fraction is 0.416667
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.
Distance, velocity, and linear accelerationFuel mass, exhaust velocity, and thrust
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:
Earth to Mars Trajectory, 2.1 days, 0.5 AU traversed
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!
The Fusion Energy Toolbox for MATLAB is a toolbox for designing fusion reactors and for studying plasma physics. It includes a wide variety of physics and engineering tools. The latest addition to this toolbox is a new function for designing tokamaks, based on the paper in reference [1]. Tokamaks have been the leading magnetic confinement devices investigated in the pursuit of fusion net energy gain. Well-known tokamaks that either have ongoing experiments or are under development include JET, ITER, DIII-D, KSTAR, EAST, and Commonwealth Fusion Systems’ SPARC. The new capability of our toolboxes to conduct trade studies on tokamaks allows our customers to take part in this exciting field of fusion reactor design and development.
The Fusion Reactor Design function checks that the reactor satisfies key operational constraints for tokamaks. These operational constraints result from the plasma physics of the fusion reactor, where there are requirements for the plasma to remain stable (e.g., not crash into the walls) and to maintain enough electric current to help sustain itself. The tunable parameters include: the plasma minor radius ‘a’ (see figure below), the H-mode enhancement factor ‘H’, the maximum magnetic field at the coils ‘B_max’, the electric power output of the reactor ‘P_E’, and the neutron wall loading ‘P_W’, which are all essential variables to tokamak design and operation. H-mode is the high confinement mode used in many machines.
Illustration of the toroidal plasma of a tokamak. R is the major radius while a is the minor radius of the plasma. The red line represents a magnetic field line which helically winds along the torus. Image from [2].
This function captures all figure and table results in the original paper. We implemented a numerical solver which allows the user to choose a variable over which to perform a parameter sweep. A ‘mode’ option has been incorporated which allows one to select a desired parameter sweep variable (‘a’, ‘H’, ‘B_max’, ‘P_E’, or ‘P_W’) when calling the function. Some example outputs of the function are described below.
As an example, we will consider the case of tuning the maximum magnetic field at the coils ‘B_max’. The figure below plots the normalized operation constraint parameters for a tokamak as functions of B_max from 10 Tesla to 25 Tesla. The unshaded region, where the vertical axis is below the value of 1, is the region where operational constraints are met. We see that for magnetic fields below about 17.5 Tesla there is at least one operation constraint that is not met, while for higher magnetic fields all operation constraints are satisfied, thus meeting the conditions for successful operation. This high magnetic field approach is the design approach of Commonwealth Fusion Systems for the reactor they are developing [3].
Operational constraint curves as a function of B_max. Successful operation occurs if all of the curves are in the unshaded region. Note, f_B/f_NC, a ratio of the achievable to required bootstrap current, is set equal to 1. In this case P_E = 1000 MW, P_W = 4 MW/m2, and H = 1. For more details on the plotted parameters and how they function as operational plasma constraints, see reference [1].
Note, however, that there is a material cost associated with achieving higher magnetic fields, as described in reference [1]. This is illustrated in the figure below, which plots the cost parameter (the ratio of engineering components volume V_I to electric power output P_E) against B_max. There is a considerable increase in cost at high magnetic fields due to the need to add material volume that can structurally handle the higher current loads required.
Cost parameter (units of volume in cubic meters per megawatt of power, m3/MW) as a function of B_max.
In this post we illustrated the case of a tunable maximum magnetic field at the coils, though as mentioned earlier, there are other parameters you can tune. This function is part of release 2022.1 of the Fusion Energy Toolbox. Contact us at info@psatellite.com or call us at +01 609 275-9606 for more information.
Thank you to interns Emma Suh and Paige Cromley for their contributions to the development of this function.
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.
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.
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.
Solar Sail and Earth paths in the heliocentric frame.Navigation camera view.
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.
We’ve added some new tools to the Aircraft Control Toolbox for our upcoming 2019 release. The first is a new GUI for creating aircraft models. You import a Wavefront OBJ files and then you point and click to define leading edges, wing areas, engine locations and so forth. This makes it easier to import the geometric data. The GUI is shown below. It illuminates the view that you need to use for a given geometric element in red. The inertia matrix is generated from the mass and the surface geometry.
The new Model Creation GUI
A new simulation function was added to use the data from this GUI. It has a flat Earth aircraft model with a plugins architecture. You can add your own lift, drag and thrust models or use the simple built-in models. It is much simpler than AC.m which is designed to be a comprehensive high-fidelity simulation. We’ve added a new animation GUI to show you the results of your simulations.
We expect 2019.1 to be available in June. You can get a demo with previews of the new functions now.
New functions in the Lunar Cube module in 2016.1 allow you to easily plan lunar insertion and orbit change maneuvers. In the following pictures you can see a lunar orbit insertion from a hyperbolic orbit. In all figures the lunar terrain is exaggerated by a factor of 10.
The same maneuver looking down on the orbit plane. The green arrows are the force vectors.
The following figure shows a two maneuver sequence. The first puts the spacecraft into an elliptical orbit. The second circularizes the orbit.
We are adding the Lunar Cube Module in 2016.1 to our CubeSat Toolbox for MATLAB! It allows users to analyze and simulateCubeSats in lunar transfer and lunar orbit. It includes a new dynamical model for CubeSats that includes:
Earth, Moon and Sun gravity based on the JPL ephemerides
Spherical harmonic lunar gravity model
Reaction wheels
Thrusters
Power generation from solar panels
Battery energy storage
Variable mass due to fuel consumption
Solar pressure disturbances
Lunar topographic model
New graphics functions for lunar orbit operations
Lunar targeting function
Lunar mission control function for attitude control and orbit control
The module includes a script with a simulation of a 6U Cubesat leaving Earth orbit and reaching the moon. The following figure shows the Earth to Moon trajectory.
This figure shows the transfer orbit near the moon. The lunar topography is exaggerated by a factor of 10 to make it visible. It is based on Clementine measurements.
Here are results from the new LunarTargeting function. It finds optimal transfers to lunar orbits. The first shows the transfer path to the Moon’s sphere of influence.
The next shows the lunar hyperbolic orbit. In this case the transfer is into a high inclination lunar orbit.
Here is a picture of our DFD transfer vehicle. You can see the lander on the front and two Deep Space Optical Communication System (DSOC) assemblies mounted on trusses. There are 2 DFD engines.
A picture of the Pluto Lander. The solar panels are illuminated by a laser from the orbiter. The lander has a dry mass of 150 kg.
