Rocket Velocity Example

 

·       The physics of a rocket launched straight up (neglecting air resistance) tells us that the velocity of the rocket is given by the following mathematical relationship

 

 

·       The physical quantities represented by each of the variables are:

 

1.    g – acceleration due to gravity (feet/second2)  (g=32.17)

2.    m – Total initial rocket mass (slugs, slugs*g=weight)

3.    Note: m = mass of fuel + mass of empty rocket (m=mf+me)

4.    u – speed of burning fuel mass exiting the rocket (feet/second)

5.    q – rate of fuel mass loss (slugs/second)

6.    t – time since launch (seconds)

7.    b – the time at which we run out of fuel (stop burn time) (seconds).

 

·       Let’s assume we have been given this much (at this point, don’t worry about how this equation was derived).  What we want to do is investigate what this equation is “telling” us.  We want to understand what our rocket will be doing (in particular how fast will it be going and ultimately where is it over time).

 

·       Let’s start by trying to plot this.  Note that this equation is unlike any we have discussed so far because it is piece-wise defined.  It has two distinct parts!  One way to deal with this is to use a for loop over t, with an if statement to test which part of the function above we are in for a given time.

 

·       The break point in the equation above is at time b, when we stop the burn (shut down the engines).  It now rides on momentum and is at the mercy of gravity.  Because the situation has changed at time b, a new equation describes its behavior (that makes sense…).

 

·       If we burn until we lose all our fuel mass, then

 

and

 

 

    

·       Let’s try to write MATLAB code to evaluate this equation and plot velocity vs. time for a particular rocket.  Let’s use these values: mf=100 slugs, me=100 slugs, m=200 slugs, u=8000 ft/s, q=1 slug/s, g=32.17 ft/s2, b=100 s.

 

clear; close all; clc

m=200;

u=8000;

q=1;

g=32.17;

b=100;

 

count=1;

for t=0:3:300

   

    if t<=b

        v(count)=u*log(m/(m-q*t))-g*t;

    else  

        v(count)=u*log(m/(m-q*b))-g*t;

    end

   

    tt(count)=t;

    count=count+1;

   

end

 

figure(1)

plot(tt,v)

xlabel('Time (sec) ');

ylabel('Velocity (ft/sec) ');

title('Rocket Velocity Curve ');

axis tight

 

 

·       Here is a flowchart of this program.

 

·       What we probably really want to know is the height of the rocket at each time.  We can compute this from the velocity if we are clever.  The height will start at a reference of 0 feet.  The height at the next time “slice” (after a set time interval), can be estimated by the velocity multiplied by the length of the time interval.   More specifically,

 

new height = old height + (old velocity)x(length of time interval).

 

          We also know that the height should never go negative (it will hit the ground at height 0).

 

·       Let’s try to add this relationship to our code.

 

clear; close all; clc

m=200;

u=8000;

q=1;

g=32.17;

b=100;

 

count=1;

h(1)=0;

for t=0:3:300

   

    if t<=b

        v(count)=u*log(m/(m-q*t))-g*t;

    else  

        v(count)=u*log(m/(m-q*b))-g*t;

    end

   

 

    h(count+1)=h(count)+v(count)*3;

   

    if h(count+1)<0;

        h(count+1)=0;

    end

   

    tt(count)=t;

    count=count+1;

   

end

 

figure(1)

subplot(121)

plot(tt,v)

xlabel('Time (sec) ');

ylabel('Velocity (ft/sec) ');

title('Rocket Velocity Curve ');

grid

axis tight

 

subplot(122)

N=length(tt);

% h array is one longer than tt.

plot(tt,h(1:N))

xlabel('Time (sec) ');

ylabel('Height (ft) ');

title('Rocket Height Curve ');

grid

axis tight

 

     % display the max height and when that happens in command window

[maxh,nmax]=max(h);

fprintf('Rocket reaches %.1f feet after %.1f seconds \n',maxh,tt(nmax));

 

 

·       This problem is perfect for an animation!!  Try this…

 

figure(2)    % new figure window

for n=1:N % remember N is the length of our time array

   plot(0,h(n),'r^');                          % plot the height vertically, use a red triangle

   axis([-1,1,0,maxh]);                    % fix the axes to show the full height of the rocket

             ylabel('Height (ft) ')

   title('Rocket Flight Animation')

   S=sprintf('t=%.2f sec',tt(n));    % create a string with the current time

   text(-.5,10000,S);                        % display the current time

   M(n)=getframe;                          % store figure for movie later

end

movie(M,1,10);           % play the frames as a movie

 

 

·       Notice that tt(34)=99 (that’s when we run out of fuel).  Let’s change the color of our rocket when it is out of fuel.

 

figure(2)   

for n=1:N

 

     if n<34

          plot(0,h(n),'r^'); % red while it has fuel

          text(.5,10000, 'BURN ON! ')

              else

plot(0,h(n),'b^'); % blue when it is gliding

text(.5,10000, 'BURN OFF! ')

              end

 

         

   axis([-1,1,0,maxh]);

             ylabel('Height (ft) ')

   title('Rocket Flight Animation')                 

   S=sprintf('t=%.2f sec',tt(n));   

   text(-.5,10000,S);                       

   M(n)=getframe;

                            

end

movie(M,1,10);