How Rockets Work

The fuel in a rocket can be solid fuel, liquid fuel, or gaseous fuel. When the rocket burns fuel, it creates very hot gas which is allowed to escape through a nozzle at the bottom of the rocket. Each kind of fuel burns to give a particular average velocity to the gas and therefore a particular escape velocity when it leaves the nozzle. This fuel was moving with the rocket, but it is now moving away from the rocket with a speed $v_e$, that is, the escape velocity. If the escaping gas in some time interval has a mass $m_e$, then when it leaves the nozzle it will have experienced a change in momentum of $\Delta p_{exhaust} = -m_e v_e$. This must also equal the impulse exerted on the gas by the rocket. By Newton's third law, the rocket must have experienced an equal and opposite impulse and also an equal and opposite change in momentum.

If the rocket plus fuel had a mass $m$ and was moving at a velocity of $v$ before this time interval, then it's change in velocity $\Delta v$ can be found using momentum conservation such that $-\Delta p_{exhaust} = \Delta p_{rocket}$ and $m_e \times v_e = m \times \Delta v$.

Discovering the Rocket Equation

In a very short time interval, $dt$, then the amount of fuel burnt and exhausted will be $dm_e$. The change in the rocket's velocity will be $dv$. So momentum conservation in this short time interval becomes $dm_e \times v_e = m \times dv$.

But if we burn and exhaust an amount of fuel $dm_e$, the mass of the rocket must be reduced by the same amount, $dm = -dm_e$.

Making this substitution gives $-dm \times v_e = m \times dv$. We can rearrange this into $$dv=-v_e\frac{dm}{m}.$$

If we integrate this from the start where $v=0$ and $m=m_0$ to some later point we can find what the velocity $v$ must be if the remaining mass of the rocket and fuel is $m$, $$v=-v_e \ln \frac{m}{m_0}.$$

This result is called the rocket equation.

From Momentum to Force

For our simulation, we want to use the Euler method as we did previously, so we need to calculate the thrust force from the rocket. We know that the impulse on the rocket in a short time interval, $dt$, which we already calculated must be equal to the thrust force times the time interval, $J = dm_e \times v_e = F_{thrust} \times dt$. So the force must equal $$F_{thrust} = v_e \frac{dm_e}{dt},$$ where $\frac{dm_e}{dt}$ is the burn rate of the fuel.

Making a Simulation of a Rocket Launch

We'll be simulating the motion of a V2 rocket (https://en.wikipedia.org/wiki/V-2_rocket), the world's first long-range guided ballistic missile. It was an incredible technical feat using liquid oxygen and alcohol propellents to eventually achieve a range of 220 miles, a ceiling of 50-60 miles and a speed of 3400 mph.

The original V2 rocket had an empty mass of $4000 kg$ and was loaded with $8800 kg$ of fuel. It burned fuel at $129.4 kg/s$ and had an exhaust velocity of $2050 m/s$.

We'll start by assigning all the rocket variables including mass of empty rocket, starting fuel amount, burn rate of the fuel, and exhaust velocity. We also need to give initial values for position and velocity of the rocket and time as well as the size of the time steps in the simulation and the maximum time.

The simulation loop runs until time runs out. We'll assume the the rocket has a constant burn rate for it's fuel, so we need to run as long as there is fuel, but something has to change when the fuel runs out. Notice the if and else statements inside the while loop.

Your job is to put in the correct expressions for the thrust, net force on the rocket, and the updated mass of the rocket into the simulation below. You may have different expressions before and after fuel runs out.

At the end we display the final velocity of the rocket and compare to the prediction of the rocket equation.

from pylab import * # Downloads libraries

#Variables
Mr = 4000   #mass of empty rocket (kg)
Mf = 8800   #mass of initial fuel (kg)
Mtot = Mr+Mf    #total mass of the rocket and fuel (kg)
burn = 129.4    #change in mass as rocket uses fuel (dm/s)
Ve = 2050   #exhaust velocity (m/s)
t = 0
dt = 0.1    #time step
tmax = 90   #max time that simulation will run
v = 0   #V2 intial velocity
d = 0   #initial distance traveled (starting altitude at launch)

#Creation of graph
figure(figsize=(18,6), dpi=160)

#Calculation loop
while t < tmax:
  if Mtot > Mr: #parameter set to a half step furter to account for last remaining fuel under a full step

    #Forces acting on the system
    thrust = Ve*burn
    Fnet=thrust

    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot-burn*dt

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'ro', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'ro', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'ro', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

  else:
    #Forces acting on the system
    thrust = 0
    Fnet=0

    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'bo', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'bo', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'bo', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

print('Final velocity was', '{:.2E}'.format(v), 'm/s')

#We'll compare to the velocity predicted by the rocket equation
Rv=-Ve*log(Mr/(Mr+Mf))

print('Final velocity predicted by rocket equation is', '{:.2E}'.format(Rv), 'm/s')

Final velocity was 2.39E+03 m/s
Final velocity predicted by rocket equation is 2.38E+03 m/s

Exercise 2. Next we will work to make our simulation more realistic by first adding the force of gravity. First we'll assume a constant force of gravity, $F_G = mg$. We should also change the simulation to run until the rocket returns to ground.

from pylab import * # Downloads libraries

#Variables
Mr = 4000   #mass of empty rocket (kg)
Mf = 8800   #mass of initial fuel (kg)
Mtot = Mr+Mf    #total mass of the rocket and fuel (kg)
burn = 129.4    #change in mass as rocket uses fuel (dm/s)
Ve = 2050   #exhaust velocity (m/s)
t = 0
dt = 0.1    #time step
tmax = 450   #max time that simulation will run
v = 0   #V2 intial velocity
d = 0   #initial distance traveled (starting altitude at launch)

#Creation of graph
figure(figsize=(18,6), dpi=160)

#Calculation loop
while t < tmax:
  if Mtot > Mr: #parameter set to a half step furter to account for last remaining fuel under a full step

    #Forces acting on the system
    thrust = Ve*burn
    gravity=-Mtot*9.81
    Fnet=thrust+gravity


    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot-burn*dt

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'ro', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'ro', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'ro', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

  else:
    #Forces acting on the system
    thrust = 0
    Fnet=Mr*-9.81

    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'bo', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'bo', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'bo', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

print('Final velocity was', '{:.2E}'.format(v), 'm/s')

#We'll compare to the velocity predicted by the rocket equation
Rv=-Ve*log(Mr/(Mr+Mf))

print('Final velocity predicted by rocket equation is', '{:.2E}'.format(Rv), 'm/s')

Final velocity was -2.04E+03 m/s
Final velocity predicted by rocket equation is 2.38E+03 m/s

Exercise 3. We know that the gravitational force changes if the distance from the center of the earth changes, so we will modify to include the more accurate force of gravity equation, $F_G = GMm/r^2$.

from pylab import * # Downloads libraries

#Variables
Mr = 4000   #mass of empty rocket (kg)
Mf = 8800   #mass of initial fuel (kg)
Mtot = Mr+Mf    #total mass of the rocket and fuel (kg)
burn = 129.4    #change in mass as rocket uses fuel (dm/s)
Ve = 2050   #exhaust velocity (m/s)
t = 0
dt = 0.1    #time step
tmax = 450   #max time that simulation will run
v = 0   #V2 intial velocity
d = 0   #initial distance traveled (starting altitude at launch)

#Creation of graph
figure(figsize=(18,6), dpi=160)

#Calculation loop
while t < tmax:
  if Mtot > Mr: #parameter set to a half step furter to account for last remaining fuel under a full step

    #Forces acting on the system
    thrust = Ve*burn
    G=6.6743E-11
    M=5.972E24
    Re=6371000
    gravity=-(G*M*Mtot)/(Re+d)**2
    Fnet=thrust+gravity


    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot-burn*dt

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'ro', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'ro', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'ro', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

  else:
    #Forces acting on the system
    thrust = 0
    G=6.6743E-11
    M=5.972E24
    Re=6371000
    gravity=(G*M*Mtot)/(Re+d)**2
    Fnet=-gravity

    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'bo', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'bo', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'bo', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

print('Final velocity was', '{:.2E}'.format(v), 'm/s')

#We'll compare to the velocity predicted by the rocket equation
Rv=-Ve*log(Mr/(Mr+Mf))

print('Final velocity predicted by rocket equation is', '{:.2E}'.format(Rv), 'm/s')

Final velocity was -1.87E+03 m/s
Final velocity predicted by rocket equation is 2.38E+03 m/s

Exercise 4. We'll add air resistance to make it even more realistic. The formula for drag is $F_{drag} = \frac{1}{2} \rho v^2 C_d A$, where $\rho$ is the air density, $C_d$ is the drag coefficient, and $A$ is the cross-sectional area.

from pylab import * # Downloads libraries

#Variables
Mr = 4000   #mass of empty rocket (kg)
Mf = 8800   #mass of initial fuel (kg)
Mtot = Mr+Mf    #total mass of the rocket and fuel (kg)
burn = 129.4    #change in mass as rocket uses fuel (dm/s)
Ve = 2050   #exhaust velocity (m/s)
t = 0
dt = 0.1    #time step
tmax = 450   #max time that simulation will run
v = 0   #V2 intial velocity
d = 0   #initial distance traveled (starting altitude at launch)
G=6.6743E-11
M=5.972E24
Re=6371000
g=9.81
A=3.14*(1.65/2)**2
#Creation of graph
figure(figsize=(18,6), dpi=160)

#Calculation loop
while d>=0:
  if Mtot > Mr: #parameter set to a half step furter to account for last remaining fuel under a full step

    #Forces acting on the system
    Mtot==Mtot-(burn*dt)
    thrust = Ve*burn
    gravity=-(G*M*Mtot)/(Re+d)**2
    Fd=0.5*1.22*(-v)*abs(v)*0.125*A
    Fnet= thrust + gravity + Fd


    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot-burn*dt

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'ro', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'ro', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'ro', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

  else:
    #Forces acting on the system
    thrust = 0
    gravity=-(G*M*Mtot)/(Re+d)**2
    Fd=0.5*1.22*(-v)*abs(v)*0.125*A
    Fnet=Mr*-9.81+Fd

    #Acceleration
    a=Fnet/(Mtot)
    Mtot=Mtot

    #Graphs being plotted
    subplot(131)
    plot(t, v, 'bo', markersize=0.1)
    ylabel('Velocity [m/s]')
    xlabel("time [s]")

    subplot(132)
    plot(t, d, 'bo', markersize=0.1)
    ylabel('Distance [m]')
    xlabel("time [s]")

    subplot(133)
    plot(t, a, 'bo', markersize=0.1)
    ylabel('Acceleration [m/s^2]')
    xlabel("time [s]")

    #Velocity using euler approximation
    v=v+a*dt

    #Distance traveled using euler approximation
    d=d+v*dt+a*dt*dt/2

    #time step
    t= t+dt

print('Final velocity was', '{:.2E}'.format(v), 'm/s')

#We'll compare to the velocity predicted by the rocket equation
Rv=-Ve*log(Mr/(Mr+Mf))

print('Final velocity predicted by rocket equation is', '{:.2E}'.format(Rv), 'm/s')

Final velocity was -4.88E+02 m/s
Final velocity predicted by rocket equation is 2.38E+03 m/s

Exercise 5. In reality, the air density is not constant as the altitude changes. We'll need to somehow incorporate the changes in air density into our simulation. How can we do that?

Now we want to learn how to store the data from the simulation and write an output file that we can read into Excel or another program.