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

# 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="../../data/weather/ndrt_2020_weather_data_ERA5.nc",
    dictionary="ECMWF",
)

Visualize the Environment object

# env.info()

Motor

Define the SolidMotor object

motor_l1395 = SolidMotor(
    thrust_source="../../data/motors/cesaroni/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],
)

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 (without motor)
Rocket Dry Mass: 20.846 kg (with unloaded motor)
Rocket Loaded Mass: 23.321 kg
Rocket Inertia (with unloaded motor) 11: 78.421 kg*m2
Rocket Inertia (with unloaded motor) 22: 78.421 kg*m2
Rocket Inertia (with unloaded motor) 33: 0.160 kg*m2
Rocket Inertia (with unloaded motor) 12: 0.000 kg*m2
Rocket Inertia (with unloaded motor) 13: 0.000 kg*m2
Rocket Inertia (with unloaded motor) 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 without Motor: 0.154 m
Rocket Center of Dry Mass - Nozzle Exit: 1.937 m
Rocket Center of Dry Mass - Center of Propellant Mass: 1.587 m
Rocket Center of Mass - Rocket Loaded Center of Mass: 0.168 m


Aerodynamics Lift Coefficient Derivatives

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

Center of Pressure

Nose Cone Center of Pressure position: 0.282 m
Tail Center of Pressure position: 1.261 m
Fins Center of Pressure position: 3.097 m

Stability

Center of Mass position (time=0): 1.623 m
Center of Pressure position (time=0): 2.438 m
Initial Static Margin (mach=0, time=0): 4.016 c
Final Static Margin (mach=0, time=burn_out): 4.846 c
Rocket Center of Mass (time=0) - Center of Pressure (mach=0): 0.815 m


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

Initial time: 0.000 s
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 (quaternions) - e0: -0.391 | e1: 0.000 | e2: 0.000 | e3: -0.921
Euler Angles - Spin φ : -293.00° | Nutation θ: 0.00° | Precession ψ: 67.00°
Angular Velocity - ω1: 0.00 rad/s | ω2: 0.00 rad/s | ω3: 0.00 rad/s
Initial Stability Margin: 4.017 c


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.253 s
Rail Departure Velocity: 13.044 m/s
Rail Departure Stability Margin: 4.072 c
Rail Departure Angle of Attack: 23.546°
Rail Departure Thrust-Weight Ratio: 6.657
Rail Departure Reynolds Number: 2.061e+05


Burn out State

Burn out time: 3.433 s
Altitude at burn out: 533.391 m (ASL) | 327.391 m (AGL)
Rocket speed at burn out: 170.838 m/s
Freestream velocity at burn out: 172.574 m/s
Mach Number at burn out: 0.518
Kinetic energy at burn out: 3.042e+05 J


Apogee State

Apogee Time: 16.650 s
Apogee Altitude: 1499.430 m (ASL) | 1293.430 m (AGL)
Apogee Freestream Speed: 23.981 m/s
Apogee X position: -88.586 m
Apogee Y position: -204.583 m
Apogee latitude: 41.7736065°
Apogee longitude: -86.5735356°


Parachute Events

Parachute: Drogue
        Ejection time: 16.657 s
        Inflation time: 17.657 s
        Freestream speed at inflation: 25.573 m/s
        Altitude at inflation: 1494.688 m (ASL) | 1293.429 m (AGL)
Parachute: Main
        Ejection time: 59.362 s
        Inflation time: 60.362 s
        Freestream speed at inflation: 27.777 m/s
        Altitude at inflation: 346.028 m (ASL) | 167.710 m (AGL)


Impact Conditions

Time of impact: 90.112 s
X impact: 416.230 m
Y impact: 134.939 m
Altitude impact: 206.000 m (ASL) | -0.000 m (AGL)
Latitude: 41.7766609°
Longitude: -86.5674457°
Vertical velocity at impact: -4.528 m/s
Number of parachutes triggered until impact: 2


Stability Margin

Initial Stability Margin: 4.017 c at 0.00 s
Out of Rail Stability Margin: 4.072 c at 0.25 s
Maximum Stability Margin: 5.010 c at 3.34 s
Minimum Stability Margin: 4.017 c at 0.00 s


Maximum Values

Maximum Speed: 173.142 m/s at 3.26 s
Maximum Mach Number: 0.525 Mach at 3.26 s
Maximum Reynolds Number: 2.469e+06 at 3.26 s
Maximum Dynamic Pressure: 1.840e+04 Pa at 3.26 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: 13.675 m/s² at 7.63 s
Maximum Gs After Motor Burn: 1.394 Gs at 7.63 s
Maximum Stability Margin: 5.010 c at 3.34 s
Maximum Upper Rail Button Normal Force: 1.278 N
Maximum Upper Rail Button Shear Force: 2.512 N
Maximum Lower Rail Button Normal Force: 5.994 N
Maximum Lower Rail Button Shear Force: 11.749 N


Numerical Integration Settings

Maximum Allowed Flight Time: 600.00 s
Maximum Allowed Time Step: inf s
Minimum Allowed Time Step: 0.00e+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: 487
Number of Derivative Functions Evaluation: 163
Average Function Evaluations per Time Step: 0.335

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(
    "../../data/rockets/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: 1293.43 m
Real data: 1320.29 m
Absolute error: 26.86 m
Relative error: 2.03 %

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.65 m/s
Real data: 164.21 m/s
Absolute error: 7.43 m/s
Relative error: 4.53 %

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