Notre Dame Rocket Team Rocket 2020

Launched at 19045-18879 Avery Rd, Three Oaks, MI 49128. Permission to use flight data given by Brooke Mumma, 2020.

Import Results (23rd feb)

1. Measured Stability Margin 2.875 cal

2. Official Target Altitude 4,444 ft

3. Measured Altitude 4,320 ft or 1316.736 m

4. Drift: 2275 ft

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

# Importing libraries
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import savgol_filter

from rocketpy import Environment, Flight, Function, Rocket, SolidMotor

plt.style.use("seaborn-v0_8-dark-palette")

RocketPy Simulation

Define a dictionary with the inputs for the simulation

parameters = {
    # Mass Details
    "rocket_mass": (18.998, 0.010),  # Rocket dry mass: 20.846 kg
    # propulsion details
    "motor_structure_mass": (1.848, 0.1),
    "burn_time": (3.433, 0.1),
    "nozzle_radius": (0.02475, 0.001),
    "throat_radius": (0.01075, 0.001),
    "grain_separation": (0.003, 0.001),
    "grain_density": (1519.708, 30),
    "grain_outer_radius": (0.033, 0.001),
    "grain_initial_inner_radius": (0.015, 0.002),
    "grain_initial_height": (0.12, 0.001),
    "grains_center_of_mass_position": (-0.35, 0.100),
    "nozzle_position": (0, 0.100),
    "motor_position": (3.391, 0.100),
    # aerodynamic details
    "center_of_mass_without_motor": (1.3, 0.100),
    "drag_coefficient": (0.44, 0.1),
    "inertia_i": (73.316, 0.3 * 73.316),
    "inertia_z": (0.15982, 0.3 * 0.15982),
    "radius": (0.1015, 0.001),
    "power_off_drag": (1, 0.033),
    "power_on_drag": (1, 0.033),
    ## nose cone
    "nose_length": (0.610, 0.001),
    "nose_radius": (0.1015, 0.001),
    "nose_position": (0, 0.100),
    ## fins
    "fin_span": (0.165, 0.001),
    "fin_root_chord": (0.152, 0.001),
    "fin_tip_chord": (0.0762, 0.001),
    "fin_sweep_angle": (13, 0.5),
    "fin_position": (3.050, 0.100),
    ## transitions
    "transition_top_radius": (0.1015, 0.010),
    "transition_bottom_radius": (0.0775, 0.010),
    "transition_length": (0.127, 0.010),
    "transition_position": (1.2, 0.010),
    # launch and environment details
    "wind_direction": (0, 3),
    "wind_speed": (1, 0.30),
    "inclination": (90, 1),
    "heading": (181, 3),
    "rail_length": (3.353, 0.001),
    # parachute details
    "cd_s_drogue": (1.5 * np.pi * (24 * 25.4 / 1000) * (24 * 25.4 / 1000) / 4, 0.1),
    "cd_s_main": (2.2 * np.pi * (120 * 25.4 / 1000) * (120 * 25.4 / 1000) / 4, 0.1),
    "lag_rec": (1, 0.5),
}

# rocket: nose_to_tail

Environment

Define the Environment object

# Environment conditions
env = Environment(
    gravity=9.81,
    latitude=41.775447,
    longitude=-86.572467,
    date=(2020, 2, 23, 16),
    elevation=206,
)

env.set_atmospheric_model(
    type="Reanalysis",
    file="../../tests/fixtures/acceptance/NDRT_2020/ndrt_2020_weather_data_ERA5.nc",
    dictionary="ECMWF",
)

Visualize the Environment object

# env.info()

Motor

Define the SolidMotor object

motor_l1395 = SolidMotor(
    thrust_source="../../tests/fixtures/acceptance/NDRT_2020/ndrt_2020_motor_Cesaroni_4895L1395-P.eng",
    burn_time=parameters.get("burn_time")[0],
    dry_mass=parameters.get("motor_structure_mass")[0],
    dry_inertia=(0, 0, 0),
    center_of_dry_mass_position=parameters.get("grains_center_of_mass_position")[0],
    grains_center_of_mass_position=parameters.get("grains_center_of_mass_position")[0],
    grain_number=5,
    grain_separation=parameters.get("grain_separation")[0],
    grain_density=parameters.get("grain_density")[0],
    grain_outer_radius=parameters.get("grain_outer_radius")[0],
    grain_initial_inner_radius=parameters.get("grain_initial_inner_radius")[0],
    grain_initial_height=parameters.get("grain_initial_height")[0],
    nozzle_radius=parameters.get("nozzle_radius")[0],
    throat_radius=parameters.get("throat_radius")[0],
    interpolation_method="linear",
    nozzle_position=parameters.get("nozzle_position")[0],
    coordinate_system_orientation="combustion_chamber_to_nozzle",  # combustion_chamber_to_nozzle"
)

motor_l1395.plots.draw()

motor_l1395.info()

Nozzle Details
Nozzle Radius: 0.02475 m
Nozzle Throat Radius: 0.01075 m

Grain Details
Number of Grains: 5
Grain Spacing: 0.003 m
Grain Density: 1519.708 kg/m3
Grain Outer Radius: 0.033 m
Grain Inner Radius: 0.015 m
Grain Height: 0.12 m
Grain Volume: 0.000 m3
Grain Mass: 0.495 kg

Motor Details
Total Burning Time: 3.433 s
Total Propellant Mass: 2.475 kg
Average Propellant Exhaust Velocity: 1977.740 m/s
Average Thrust: 1425.839 N
Maximum Thrust: 1800.0 N at 0.1 s after ignition.
Total Impulse: 4894.905 Ns

Rocket

Create the Rocket object

ndrt2020 = Rocket(
    radius=parameters.get("radius")[0],
    mass=parameters.get("rocket_mass")[0],
    inertia=(
        parameters.get("inertia_i")[0],
        parameters.get("inertia_i")[0],
        parameters.get("inertia_z")[0],
    ),
    power_off_drag=parameters.get("drag_coefficient")[0],
    power_on_drag=parameters.get("drag_coefficient")[0],
    center_of_mass_without_motor=parameters.get("center_of_mass_without_motor")[0],
    coordinate_system_orientation="nose_to_tail",
)
ndrt2020.set_rail_buttons(1.5, 2, 45)

ndrt2020.add_motor(motor=motor_l1395, position=parameters.get("motor_position")[0])

Adding aerodynamic surfaces

nose_cone = ndrt2020.add_nose(
    length=parameters.get("nose_length")[0],
    kind="tangent",
    position=parameters.get("nose_position")[0],
)
fin_set = ndrt2020.add_trapezoidal_fins(
    4,
    span=parameters.get("fin_span")[0],
    root_chord=parameters.get("fin_root_chord")[0],
    tip_chord=parameters.get("fin_tip_chord")[0],
    position=parameters.get("fin_position")[0],
    sweep_angle=parameters.get("fin_sweep_angle")[0],
    radius=parameters.get("transition_bottom_radius")[0],
)
transition = ndrt2020.add_tail(
    top_radius=parameters.get("transition_top_radius")[0],
    bottom_radius=parameters.get("transition_bottom_radius")[0],
    length=parameters.get("transition_length")[0],
    position=parameters.get("transition_position")[0],
)

Due to the chosen bluffness ratio, the nose cone length was reduced to 0.610 m.

Adding Parachute

drogue = ndrt2020.add_parachute(
    "Drogue",
    cd_s=parameters.get("cd_s_drogue")[0],
    trigger="apogee",
    sampling_rate=105,
    lag=parameters.get("lag_rec")[0],
    noise=(0, 8.3, 0.5),
)
main = ndrt2020.add_parachute(
    "Main",
    cd_s=parameters.get("cd_s_main")[0],
    trigger=167.64,
    sampling_rate=105,
    lag=parameters.get("lag_rec")[0],
    noise=(0, 8.3, 0.5),
)

ndrt2020.draw()

ndrt2020.info()

Inertia Details

Rocket Mass: 18.998 kg
Rocket Dry Mass: 20.846 kg (With Motor)
Rocket Mass: 23.321 kg (With Propellant)
Rocket Inertia (with motor, but without propellant) 11: 83.351 kg*m2
Rocket Inertia (with motor, but without propellant) 22: 83.351 kg*m2
Rocket Inertia (with motor, but without propellant) 33: 0.160 kg*m2
Rocket Inertia (with motor, but without propellant) 12: 0.000 kg*m2
Rocket Inertia (with motor, but without propellant) 13: 0.000 kg*m2
Rocket Inertia (with motor, but without propellant) 23: 0.000 kg*m2


Geometrical Parameters

Rocket Maximum Radius: 0.1015 m
Rocket Frontal Area: 0.032365 m2

Rocket Distances
Rocket Center of Dry Mass - Center of Mass withour Motor: 0.216 m
Rocket Center of Dry Mass - Nozzle Exit Distance: 1.875 m
Rocket Center of Dry Mass - Center of Propellant Mass: 2.225 m
Rocket Center of Mass - Rocket Loaded Center of Mass: 0.236 m


Aerodynamics Lift Coefficient Derivatives

Nose Cone Lift Coefficient Derivative: 2.000/rad
Tail Lift Coefficient Derivative: -0.834/rad
Fins Lift Coefficient Derivative: 8.674/rad

Aerodynamics Center of Pressure

Nose Cone Center of Pressure to CM: 0.282 m
Tail Center of Pressure to CM: 1.261 m
Fins Center of Pressure to CM: 3.097 m
Distance - Center of Pressure to Center of Dry Mass: -0.928 m
Initial Static Margin: 4.569 c
Final Static Margin: 5.733 c


Parachute Details

Parachute Drogue with a cd_s of 0.4378 m2
Ejection signal trigger: At Apogee
Ejection system refresh rate: 105.000 Hz
Time between ejection signal is triggered and the parachute is fully opened: 1.0 s


Parachute Details

Parachute Main with a cd_s of 16.0525 m2
Ejection signal trigger: 167.64 m (AGL)
Ejection system refresh rate: 105.000 Hz
Time between ejection signal is triggered and the parachute is fully opened: 1.0 s

Flight

# Flight
flight = Flight(
    rocket=ndrt2020,
    environment=env,
    rail_length=parameters.get("rail_length")[0],
    inclination=parameters.get("inclination")[0],
    heading=parameters.get("heading")[0],
)

flight.info()
flight.plots.trajectory_3d()

Initial Conditions

Position - x: 0.00 m | y: 0.00 m | z: 206.00 m
Velocity - Vx: 0.00 m/s | Vy: 0.00 m/s | Vz: 0.00 m/s
Attitude - e0: -0.009 | e1: -0.000 | e2: -0.000 | e3: -1.000
Euler Angles - Spin φ : -90.50° | Nutation θ: -0.00° | Precession ψ: -90.50°
Angular Velocity - ω1: 0.00 rad/s | ω2: 0.00 rad/s| ω3: 0.00 rad/s


Surface Wind Conditions

Frontal Surface Wind Speed: -5.04 m/s
Lateral Surface Wind Speed: 2.57 m/s


Launch Rail

Launch Rail Length: 3.353  m
Launch Rail Inclination: 90.00°
Launch Rail Heading: 181.00°


Rail Departure State

Rail Departure Time: 0.217 s
Rail Departure Velocity: 11.023 m/s
Rail Departure Static Margin: 4.634 c
Rail Departure Angle of Attack: 27.248°
Rail Departure Thrust-Weight Ratio: 6.617
Rail Departure Reynolds Number: 1.796e+05


Burn out State

Burn out time: 3.433 s
Altitude at burn out: 327.778 m (AGL)
Rocket velocity at burn out: 170.783 m/s
Freestream velocity at burn out: 172.406 m/s
Mach Number at burn out: 0.518
Kinetic energy at burn out: 3.040e+05 J


Apogee State

Apogee Altitude: 1502.771 m (ASL) | 1296.771 m (AGL)
Apogee Time: 16.681 s
Apogee Freestream Speed: 22.609 m/s


Parachute Events

Drogue Ejection Triggered at: 16.686 s
Drogue Parachute Inflated at: 17.686 s
Drogue Parachute Inflated with Freestream Speed of: 24.250 m/s
Drogue Parachute Inflated at Height of: 1292.102 m (AGL)
Main Ejection Triggered at: 61.276 s
Main Parachute Inflated at: 62.276 s
Main Parachute Inflated with Freestream Speed of: 26.515 m/s
Main Parachute Inflated at Height of: 142.148 m (AGL)


Impact Conditions

X Impact: 461.725 m
Y Impact: 188.331 m
Time of Impact: 93.981 s
Velocity at Impact: -4.322 m/s


Maximum Values

Maximum Speed: 173.100 m/s at 3.27 s
Maximum Mach Number: 0.524 Mach at 3.27 s
Maximum Reynolds Number: 2.467e+06 at 3.25 s
Maximum Dynamic Pressure: 1.837e+04 Pa at 3.25 s
Maximum Acceleration During Motor Burn: 67.083 m/s² at 0.10 s
Maximum Gs During Motor Burn: 6.841 g at 0.10 s
Maximum Acceleration After Motor Burn: 11.691 m/s² at 10.02 s
Maximum Gs After Motor Burn: 1.192 g at 10.02 s
Maximum Upper Rail Button Normal Force: 1.235 N
Maximum Upper Rail Button Shear Force: 0.404 N
Maximum Lower Rail Button Normal Force: 9.904 N
Maximum Lower Rail Button Shear Force: 3.214 N


Numerical Integration Settings

Maximum Allowed Flight Time: 600.000000 s
Maximum Allowed Time Step: inf s
Minimum Allowed Time Step: 0.000000e+00 s
Relative Error Tolerance:  1e-06
Absolute Error Tolerance:  [0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 1e-06, 1e-06, 1e-06, 1e-06, 0.001, 0.001, 0.001]
Allow Event Overshoot:  True
Terminate Simulation on Apogee:  False
Number of Time Steps Used:  456
Number of Derivative Functions Evaluation:  1294
Average Function Evaluations per Time Step: 2.837719

Comparison with the real flight data

Load the available flight data: altitude above ground level (m), vertical velocity (m/s), time (s)

flight_data = np.loadtxt(
    "../../tests/fixtures/acceptance/NDRT_2020/ndrt_2020_flight_data.csv",
    skiprows=1,
    delimiter=",",
    usecols=(3, 4),  # 3: Time (s); 4: Altitude (Ft-AGL)
)

# TODO: In the future, we could also use the axial acceleration data to compare against the simulation

Convert to Function objects

actual_z = Function(
    source=np.column_stack((flight_data[:, 0], flight_data[:, 1] / 3.281)),
    inputs="Time (s)",
    outputs="Altitude above ground level (m)",
    interpolation="linear",
    extrapolation="zero",
)  # the division by 3.281 is to convert from feet to meters

# Calculate the actual vertical velocity as the derivative of the altitude
actual_vz = actual_z.derivative_function()
actual_vz_filtered = Function(
    source=np.column_stack(
        (actual_vz.source[:, 0], savgol_filter(actual_vz.source[:, 1], 51, 3))
    ),
    inputs="Time (s)",
    outputs="Vertical velocity (m/s)",
    interpolation="linear",
    extrapolation="zero",
)

Function.compare_plots(
    plot_list=[(actual_vz, "Actual"), (actual_vz_filtered, "Filtered")],
    ylabel="Vertical velocity (m/s)",
    title="Vertical velocity",
    xlabel="Time (s)",
)

# the actual vertical velocity will be calculated as the derivative of the altitude
actual_az = actual_vz_filtered.derivative_function()

actual_az_filtered = Function(
    source=np.column_stack(
        (actual_az.source[:, 0], savgol_filter(actual_az.source[:, 1], 51, 3))
    ),
    inputs="Time (s)",
    outputs="Vertical acceleration (m/s^2)",
    interpolation="linear",
    extrapolation="zero",
)

Function.compare_plots(
    plot_list=[(actual_az, "Actual"), (actual_az_filtered, "Filtered")],
    ylabel="Vertical acceleration (m/s^2)",
    title="Vertical acceleration",
    xlabel="Time (s)",
)

Get the simulated results

simulated_z = flight.z - env.elevation
simulated_vz = flight.vz
simulated_az = flight.az
simulated_t_final = flight.t_final
simulated_apogee = flight.apogee - env.elevation

Plots comparison

plt.plot(actual_z[:, 0], actual_z[:, 1], label="Flight data")
plt.plot(simulated_z[:, 0], simulated_z[:, 1], label="RocketPy")
plt.xlabel("Time (s)")
plt.ylabel("Altitude (m)")
plt.ylim(0, 1390)
plt.xlim(0, round(simulated_t_final, -1))
plt.legend()
plt.grid()
plt.show()

plt.plot(actual_vz_filtered[:, 0], actual_vz_filtered[:, 1], label="Flight data")
plt.plot(simulated_vz[:, 0], simulated_vz[:, 1], label="RocketPy")
plt.xlabel("Time (s)")
plt.ylabel("Vertical velocity (m/s)")
plt.xlim(0, round(simulated_t_final, -1))
plt.legend()
plt.grid()
plt.show()

plt.plot(actual_az_filtered[:, 0], actual_az_filtered[:, 1], label="Flight data")
plt.plot(simulated_az[:, 0], simulated_az[:, 1], label="RocketPy")
plt.xlabel("Time (s)")
plt.ylabel("Vertical acceleration (m/s^2)")
plt.xlim(0, round(simulated_t_final, -1))
plt.legend()
plt.grid()
plt.show()

Numerical comparison

print("Apogee (AGL)")
print(f"RocketPy: {simulated_apogee:.2f} m")
print(f"Real data: {actual_z.max:.2f} m")
diff = abs(actual_z.max - simulated_apogee)
print(f"Absolute error: {diff:.2f} m")
print(f"Relative error: {diff / actual_z.max * 100:.2f} %")

Apogee (AGL)
RocketPy: 1296.77 m
Real data: 1320.29 m
Absolute error: 23.52 m
Relative error: 1.78 %

print("Max Velocity")
print(f"RocketPy:  {simulated_vz.max:.2f} m/s")
print(f"Real data: {actual_vz_filtered.max:.2f} m/s")
velocity_error = simulated_vz.max - actual_vz_filtered.max
print(f"Absolute error: {velocity_error:.2f} m/s")
relative_error = abs(velocity_error) / actual_vz_filtered.max * 100
print(f"Relative error: {relative_error:.2f} %")

Max Velocity
RocketPy:  171.85 m/s
Real data: 164.21 m/s
Absolute error: 7.63 m/s
Relative error: 4.65 %

print("Max Acceleration (before parachute deployment)")

# For some reason, the acceleration data gets a shock at the deployment of a parachute
# We will investigate the acceleration data before the parachute deployment
# For pragmatical reasons, we will consider the parachute deployment to be at the half of the flight

# Compute the maximum acceleration for the first half of the flight
simulated_half_length = len(simulated_az) // 2
max_simulated_az = np.max(simulated_az[:, 1][:simulated_half_length])

actual_half_length = len(actual_az_filtered) // 2
max_actual_az_filtered = np.max(actual_az_filtered[:, 1][:actual_half_length])

# Print the results
print(f"RocketPy: {max_simulated_az:.2f} m/s²")
print(f"Real data: {actual_az_filtered.max:.2f} m/s²")

# Compute and print the errors
acceleration_error = max_simulated_az - actual_az_filtered.max
print(f"Absolute error: {acceleration_error:.2f} m/s^2")
relative_error = abs(acceleration_error) / actual_az_filtered.max * 100
print(f"Relative error: {relative_error:.2f} %")

Max Acceleration (before parachute deployment)
RocketPy: 67.08 m/s²
Real data: 72.45 m/s²
Absolute error: -5.36 m/s^2
Relative error: 7.40 %

# 1. Measured Stability Margin 2.875 cal
# 2. Official Target Altitude 4,444 ft
# 3. Measured Altitude 4,320 ft or 1316.736 m
# 4. Drift: 2275 ft