Propulsion Matlab Code
This code implemented the numerical simulation of the trajectory of the spacecraft.
tic
clear all
close all
%Speed profile with Nuclear thermal Rocket : Phase 1
%Explicit Euler Method
propulsion=true;
Mdry1=100e3;
Mdry2=50e3;
Payloadmass=25e3;
Mdry=Mdry1+Mdry2+Payloadmass;
% halfmass=1-0.5/(1+Mm/Mfuel_0);% this is a percentage
%%FUEL MASS FOR DIFFERENT PHASES
Mfuel_01=700e3;
Mfuel_02=300e3;
Mfuel_03=250e3;
%%TOTAL MASS
M0=Mdry+Mfuel_01+Mfuel_02+Mfuel_03;
AU=150e9;
%% Number of iterations
kmax=1e8;
%% Initial speed
v=42e3;
t=0;
r=150e9;
%%Step time
h=800;
%%SPECIFIC IMPULSE
IS1=550000; %realistic value 5000s (110 times the reality...)
IS2=450000;
IS3=330000;
g0=9.81;
%%FLOW RATE
alpha1=-100;
alpha2=-100;
alpha3=-100;
%%Thrust in different phases
Fthrust1=g0*(-alpha1)*IS1;
Fthrust2=g0*(-alpha2)*IS2;
Fthrust3=g0*(-alpha3)*IS3;
%%THE EFFECTIVE EXHAUST VELOCITY
Vp1=-Fthrust1/alpha1;
Vp2=-Fthrust2/alpha2;
Vp3=+Fthrust3/alpha2;
R_eridani=9.9e16-150e9;
Msun=1.989e30;
Meridani=0.85*Msun;
G=6.67e-11;
%% Init loop paremeters
k=1;
tp=0;
M=M0;
deltaV1=(Vp1*alpha1/M+G*Msun/(r^2)-G*Meridani/((R_eridani-r)^2));
print=1500;
kprint=0;
Mfuel1=Mfuel_01;
Mfuel2=Mfuel_02;
Mfuel3=Mfuel_03;
M=Mdry+Mfuel1+Mfuel2+Mfuel3;
while (k<=kmax) &&r<R_eridani
% tp=[tp,t];
%v(k+1)=v(k)*(1/h+G*Msun*h/r(k)^3)/(1/h-G*Msun*h/r(k)^3)-(G*Msun*h/r(k)^2+Vp*alpha/(M0-alpha*t))/(1/h-G*Msun*h/r(k)^3);
%%to avoid out of memory
if print==1500
kprint=kprint+1;
vv(kprint)=v;
rr(kprint)=r;
MM(kprint)=M;
% MMfuel(kprint)=Mfuel;
tt(kprint)=t;
print=0;
end
print=print+1;
if Mfuel1>=0&&propulsion
deltaV1=(Vp1*alpha1/M+G*Msun/(r^2)-G*Meridani/(R_eridani-r)^2);
v1=v;
v=v-h*deltaV1;
r=r+(v+v1)*h/2;
t=t+h;
k=k+1;
Mfuel1=Mfuel1+alpha1*h;
M=Mdry+max(0,Mfuel1)+max(0,Mfuel2)+max(0,Mfuel3);
else
if Mfuel2>=0
deltaV2=(Vp2*alpha2/M+G*Msun/(r^2)-G*Meridani/((R_eridani-r)^2));
v1=v;
v=v-h*deltaV2;
r=r+(v+v1)*h/2;
t=t+h;
k=k+1;
Mfuel2=Mfuel2+alpha2*h;
M=Mdry2+Payloadmass+max(0,Mfuel2)+max(0,Mfuel3);
else
if r<R_eridani-10000*AU
% M=Mdry2+Mfuel3+Payloadmass;
deltaV=G*Msun/(r^2)-G*Meridani/((R_eridani-r)^2);
v1=v;
v=v-h*deltaV;
r=r+(v+v1)*h/2;
t=t+h;
k=k+1;
else
if Mfuel3>0 &&v>100
deltaV3=(Vp3*alpha2/M+G*Msun/(r^2)-G*Meridani/((R_eridani-r)^2));
v1=v;
v=v-h*deltaV3;
r=r+(v+v1)*h/2;
t=t+h;
k=k+1;
Mfuel3=Mfuel3+alpha3*h;
M=Mdry2+Payloadmass+max(0,Mfuel3);
else
% M=Mdry2+Mfuel3+Payloadmass;
deltaV=G*Msun/(r^2)-G*Meridani/((R_eridani-r)^2);
v1=v;
v=v-h*deltaV;
r=r+(v+v1)*h/2;
t=t+h;
k=k+1;
end
end
end
end
end
% tp1=tp(1:end-1);
tmonth=tt/(3600*24*30);
vkm=vv/1000;
rAU=rr/150e9;
Mtonnes=MM/1000;
figure(1)
plot(tmonth,vkm,'LineWidth',2);
Xlabel('time (months)','FontSize',20)
Ylabel('Spacecraft velocity (km/s)','FontSize',20)
title('Spacecraft Velocity versus time','FontSize',20)
set(gca,'FontSize',20)
grid on
figure(2)
plot(tmonth,Mtonnes,'LineWidth',2);
Xlabel('time (months)','FontSize',20)
Ylabel('Mass (tonnes)','FontSize',20)
title('Spacecraft Mass versus time','FontSize',20)
set(gca,'FontSize',20)
grid on
figure(3)
plot(tmonth,rAU,'LineWidth',2);
Xlabel('time (months)','FontSize',20)
Ylabel('Spacecraft distance traveled (AU)','FontSize',20)
title('Spacecraft distance traveled versus time','FontSize',20)
set(gca,'FontSize',20)
grid on
figure(4)
plot(rAU,vkm,'LineWidth',2);
Xlabel('distance traveled (AU)','FontSize',20)
Ylabel('Spacecraft velocity (km/s)','FontSize',20)
title('Spacecraft Velocity versus distance traveled','FontSize',20)
set(gca,'FontSize',20)
grid ons
toc