1D Rocket Launch Simulation

The Physics: Combining Thrust, Gravity, and Mass

As you know from your physics background, the motion of an object is dictated by the net force acting on it, according to Newton's Second Law, $F_{net}=ma$. For our rocket launching straight up, we have two primary forces to consider during the engine burn:

1. Thrust (T): A powerful upward force generated by the motor. We will consider this positive.
2. Weight (W): The force of gravity pulling the rocket down. $W=mg$, where $g$ is the acceleration due to gravity (approximately 9.81 m/s²). We will consider this negative.

Therefore, the net force on the rocket is:
$F_{net}(t)=T-W(t)=T-m(t)g$

Notice that because the mass $m(t)$ is changing, the weight of the rocket is also changing. To find the rocket's acceleration, we rearrange Newton's second law:
$a(t)=\frac{F_{net}(t)}{m(t)}=\frac{T}{m(t)}-g$
This simple but powerful equation is the heart of our simulation. It tells us the rocket's acceleration at any instant, given its current mass. As mass decreases, the same thrust produces a greater acceleration.

The following resource provides an excellent overview of these governing equations. While it also includes drag, which we'll cover later, its derivation of the core acceleration equation is directly relevant.

From Physics to Code: The Simulation Loop

With the acceleration equation defined, we can now use Euler's method, which we studied in Module 2, to simulate the flight step-by-step. The process involves iterating through time in small increments (dt) and updating the rocket's state (its position, velocity, and mass) at each step.

Here is the logic for our simulation loop during the engine burn phase:

1. Calculate Acceleration: At the current time step, use the current mass m to find the acceleration: a = (T / m) - g.
2. Update Velocity: The new velocity is the old velocity plus the change caused by the acceleration over the time step: v_new = v_old + a * dt.
3. Update Altitude: Similarly, the new altitude is the old altitude plus the distance covered: y_new = y_old + v_old * dt.
4. Update Mass: The rocket's mass decreases as it burns fuel: m_new = m_old - mass_flow_rate * dt.
5. Advance Time: Move to the next time step: t_new = t_old + dt.
6. Repeat: Continue this loop until the motor burns out (t >= burn_time).

After the motor burns out, thrust T becomes zero, and the mass becomes constant at its "dry" value. The rocket will continue to coast upwards, but it will be decelerating due to gravity. We'll simulate this "coast phase" with a second loop that runs until the rocket's velocity drops to zero, which marks its highest point, or apogee.

The Illinois tutorial you just looked at implements this exact structure.

A Complete Python Implementation

Now, let's write the code. The script below combines all these elements into a functional 1D simulation. It uses rocket parameters similar to those in the "Making a simple flight simulation" article you reviewed in the last lesson. As a hands-on learner, the best way to understand this is to run it, study it, and see the output for yourself.

import numpy as np
import matplotlib.pyplot as plt

# --- 1. Simulation Parameters ---
# Rocket properties
m_wet = 1.0       # Wet mass (kg) - total mass with fuel
m_dry = 0.906     # Dry mass (kg) - mass without fuel
m_prop = m_wet - m_dry # Propellant mass (kg)

# Motor properties (approximating an Estes F15 motor)
total_impulse = 49.6 # N*s
burn_time = 3.4      # s
T = total_impulse / burn_time # Average Thrust (N), assumed constant
mass_flow_rate = m_prop / burn_time # Mass flow rate (kg/s)

# Environment and Simulation settings
g = 9.81          # Gravity (m/s^2)
dt = 0.1          # Time step (s)

# --- 2. Initial Conditions & Data Storage ---
t = 0.0
y = 0.0           # Altitude
v = 0.0           # Velocity
m = m_wet         # Mass

# Lists to store simulation data, initialized with starting values
time_data = [t]
altitude_data = [y]
velocity_data = [v]
mass_data = [m]

# --- 3. Simulation Loop (Burn Phase) ---
print("Simulating Burn Phase...")
while t < burn_time:
    # Calculate net force and acceleration at the start of the timestep
    F_net = T - m * g
    a = F_net / m

    # Update state using Euler's method
    v_new = v + a * dt
    y_new = y + v * dt  # Use velocity from start of step for position update
    m_new = m - mass_flow_rate * dt
    t_new = t + dt

    # Update variables for the next iteration
    v, y, m, t = v_new, y_new, m_new, t_new

    # Store the new values
    velocity_data.append(v)
    altitude_data.append(y)
    mass_data.append(m)
    time_data.append(t)

print(f"Burnout at t={t:.2f}s, altitude={y:.2f}m, velocity={v:.2f}m/s.")

# --- 4. Simulation Loop (Coast Phase) ---
print("Simulating Coast Phase to Apogee...")
# Thrust is now 0, mass is constant at dry mass
T_coast = 0
m = m_dry 

while v > 0: # Loop until apogee (velocity becomes non-positive)
    # Calculate acceleration (only gravity acts now)
    a = (T_coast / m) - g
    
    # Update state using Euler's method
    v_new = v + a * dt
    y_new = y + v * dt
    t_new = t + dt

    # Update variables
    v, y, t = v_new, y_new, t_new
    
    # Store new values
    velocity_data.append(v)
    altitude_data.append(y)
    mass_data.append(m) # Mass is constant now
    time_data.append(t)

print(f"Apogee at t={t:.2f}s, max altitude={y:.2f}m.")

# --- 5. Plotting Results ---
print("Generating plots...")
plt.style.use('default')

# Plot 1: Altitude vs. Time
plt.figure(figsize=(10, 8))
plt.subplot(3, 1, 1)
plt.plot(time_data, altitude_data, label='Altitude')
plt.axvline(x=burn_time, color='r', linestyle='--', label='Motor Burnout')
plt.title("1D Rocket Simulation")
plt.ylabel("Altitude (m)")
plt.grid(True)
plt.legend()

# Plot 2: Velocity vs. Time
plt.subplot(3, 1, 2)
plt.plot(time_data, velocity_data, label='Velocity')
plt.axvline(x=burn_time, color='r', linestyle='--', label='Motor Burnout')
plt.ylabel("Velocity (m/s)")
plt.grid(True)
plt.legend()

# Plot 3: Mass vs. Time
plt.subplot(3, 1, 3)
plt.plot(time_data, mass_data, label='Mass')
plt.axvline(x=burn_time, color='r', linestyle='--', label='Motor Burnout')
plt.xlabel("Time (s)")
plt.ylabel("Rocket Mass (kg)")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

When you run this code, pay close attention to the shape of each graph:

• Mass: Decreases linearly during the burn and then flatlines at the dry mass.
• Velocity: Increases non-linearly during the burn (why non-linear?). After burnout, it decreases linearly due to constant gravitational deceleration.
• Altitude: Rises, with its slope (velocity) first increasing and then decreasing, peaking at apogee.

Conclusion

You have now successfully built and executed a 1D vertical rocket launch simulation from first principles. This is a significant milestone that integrates physics concepts with your Python programming skills.

Here are the key takeaways from this lesson:

• The acceleration of a rocket is determined by the net force (thrust minus weight) divided by its instantaneous mass: $a(t)=T/m(t)-g$.
• A simulation can be broken into distinct phases, such as the burn phase (thrust on, mass changing) and the coast phase (thrust off, mass constant).
• Euler's method provides a straightforward, iterative way to solve the equations of motion numerically by stepping through time and updating the rocket's state at each increment.

In our next lesson, we will explore the famous Tsiolkovsky rocket equation. This is an analytical equation that provides an ideal solution for the change in velocity. We will compare the results of our numerical simulation to this equation and discuss why they differ, highlighting the limitations of the ideal equation for atmospheric flight.