plot graph using Tsiolkovsky’s rocket equations.

Implement a program that can be used to calculate the speed and position of a rocket whose motion is expressed through the differential equation. A plot that shows a test of your implementation compared to the solution of Tsiolkovsky’s rocket equations. Use data from the tables below but set G and CD to zero.A plot of how the velocity changes in the first 1000 s of the rocket’s flight according to the solution of (1) using the parameters stated below.

m*dV/dt = ve*dm/dt -GMm/r^2 -0.5*rho*A*V^2*cd(1)

Initial values (At t=0)

V speed =0 ms-1

r radius =R

h altitude =0 km

this is my code so far.

clear;
clf;

%TIME CONDITION
End_time = 1000;
Time_step = 1;

%INTIAL CONDITION
Intial_speed = 0;
Intial_altitude = 0;

%FIXED PARAMETERS
G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = 75;                 % Area of the rocket 
Cd = 0.4;               % Drag Coefficient for the rocket
m_e = 54000;            % Empty mass of the rocket
m_0 = 894000;           % Intial mass of the rocket
v_e = 4500;             % Exhaust gas speed
dmdt_f = 5000;          % Rate of change of mass with fuel left
dmdt_o = 0;             % Rate of change of mas with all fuel being used

%STTING ARRAYS
t = 0:Time_step:End_time;
Speed = zeros(1,length(t));
Altitude = zeros(1,length(t));
Speed(1) = Intial_speed;
Altitude(1) = Intial_altitude;
r(1) = Intial_altitude+R;
m(1)= m_0;
r(1) = Intial_altitude+R;
rho(1) = 1.225;

%FOR LOOP
for i=2:length(t)
  %CALCULATION FOR dmdt BEING NOT CONSTANT THROUGHT  
    if m(i-1)> m_e
        dmdt = dmdt_f;
    else
       dmdt = dmdt_o;
    end
    
    m(i) = m(i-1) - Time_step*dmdt;             %Calculate the current speed
    
    rho(i)=1.225*10.^((-3*Altitude(i-1)/50000));      %Calculate the current air density
  
    Speed(1) = v_e + Speed(i-1) - v_e*m(i)/m(i-1) - Time_step*G*M/(r(i-1))^2 - Time_step*0.5*rho(i)*A*Cd*Speed(i-1)^2/m(i-1);  %Calculate the current Velocity
    
    r(i) = r(i-1) + Time_step*0.5*(Speed(i)+Speed(i-1));    %Calculate the current radius
    
    Altitude(i) = Altitude(i-1) + Time_step*0.5*(Speed(i)+Speed(i-1));  %Calculate the current altitude
    
end
  
%DISPLAY THE RESULTS
subplot(2,1,1)
plot (t,Speed,'m')
xlabel('Time(s)')
ylabel('Speed(m/s)')
subplot(2,1,2)
plot(t,Altitude)
xlabel('Time(s)')
ylabel('Altitude(m)')

This are my graphs

As a result for this question i need to get the following graphs but i don't know how to change my code so that i can get these graphs

I need help to fix my code so that i can get the following graph as a result

[Xavier MICHEL]

Hi, did you fix your code after all and how did you plot the theoritical height with the rocket equation ? (Orange curve)

Answer 1

[James Tursa]

This

Speed(1) = ...

needs to be this

Speed(i) = ...

and this

... v_e*m(i)/m(i-1) ...

should be this

... Time_step*v_e*dmdt/m(i-1) ...

Also, that V^2 term only works if the rocket is going up. If the rocket ever starts falling down that term will not work as is. You should change V^2 to abs(V)*V for this to work in that situation.

[Amparo Dias]

I still don't get the follwoing graphs for my code and i don't know how to fix it

[James Tursa]

Can you point me to the source of your differential equation? Also can you post your current code?

[Amparo Dias]

This is the equation

This is my code

clear;
clf;

%TIME CONDITION
End_time = 1000;
Time_step = 1;

%INTIAL CONDITION
Intial_speed = 0;
Intial_altitude = 0;

%FIXED PARAMETERS
G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = 75;                 % Area of the rocket 
Cd = 0.4;               % Drag Coefficient for the rocket
m_e = 54000;            % Empty mass of the rocket
m_0 = 894000;           % Intial mass of the rocket
v_e = 4500;             % Exhaust gas speed
dmdt_f = 5000;          % Rate of change of mass with fuel left
dmdt_o = 0;             % Rate of change of mas with all fuel being used

%SETTING ARRAYS
t = 0:Time_step:End_time;
Speed = zeros(1,length(t));
Altitude = zeros(1,length(t));
Speed(1) = Intial_speed;
Altitude(1) = Intial_altitude;
r(1) = Intial_altitude+R;
m(1)= m_0;
r(1) = Intial_altitude+R;
rho(1) = 1.225;

%FOR LOOP
for i=2:length(t)
  %CALCULATION FOR dmdt BEING NOT CONSTANT THROUGHT  
    if m(i-1)> m_e
        dmdt = dmdt_f;
    else
       dmdt = dmdt_o;
    end 
    m(i) = m(i-1) - Time_step*dmdt;             %Calculate the current speed
    rho(i)=1.225*10.^((-3*Altitude(i-1)/50000));      %Calculate the current air density
    Speed(i) = v_e + Speed(i-1) - Time_step*v_e*dmdt/m(i-1) - Time_step*G*M/(r(i-1))^2 - Time_step*0.5*rho(i)*A*Cd*Speed(i-1)^2/m(i-1);  %Calculate the current Velocity 
    r(i) = r(i-1) + Time_step*0.5*(Speed(i)+Speed(i-1));    %Calculate the current radius 
    Altitude(i) = Altitude(i-1) + Time_step*0.5*(Speed(i)+Speed(i-1));  %Calculate the current altitude
end
    
%DISPLAY THE RESULTS
subplot(2,1,1)
plot (t,Speed,'m')
xlabel('Time(s)')
ylabel('Speed(m/s)')
subplot(2,1,2)
plot(t,Altitude)
xlabel('Time(s)')
ylabel('Altitude(m)')

In addition, can you help me plot the compared graphs of veclocity vs time and altitude vs time with discretised graphs and Tsiolkovsky’s graphs just like the above graphs. I tried but it just won't work

[James Tursa]

So, that differential equation is a bit confusing because that dmdt term isn't being used as the derivative of mass with respect to time (which would be negative), it is being used as a mas flow rate (a positive number and the opposite of the mass rate of change). If we go with that, and I modify your code to calculate dvdt directly and then apply it to the Speed, I get something that looks reasonable.

Change this

Speed(i) = etc.

to this

    dvdt = v_e * dmdt / m(i-1) - G*M/r(i-1)^2 - 0.5*rho(i)*A*Cd*Speed(i-1)*abs(Speed(i-1))/m(i-1);
    Speed(i) = Speed(i-1) + Time_step*dvdt;

and I get this:

You might want to put in a check to see if you are on the ground with not enough thrust to lift off, in which case Speed and Altitude should simply be set to 0 and not propagate.

[Amparo Dias]

Thank you

[Harshita]

Hey James, I am Harshita from India. I am analysing some graphs for my uni project. I want to make the trajectroy and velocity graphs for a liquid propellent rocket, but I am not able to correctly figure out the code for that. I wanted to ask if you could please help me with that. I'll be really grateful if you see this message and revert back.

Regrads, Harshita