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)
Measured Stability Margin 2.875 cal
Official Target Altitude 4,444 ft
Measured Altitude 4,320 ft or 1316.736 m
Drift: 2275 ft
[275]:
%load_ext autoreload
%autoreload 2
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
[276]:
# 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
[277]:
plt.style.use("seaborn-v0_8-dark-palette")
RocketPy Simulation#
Define a dictionary with the inputs for the simulation
[278]:
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
[279]:
# 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
[280]:
# env.info()
Motor#
Define the SolidMotor
object
[281]:
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"
)
[282]:
motor_l1395.plots.draw()
[283]:
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
[284]:
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
[285]:
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
[286]:
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),
)
[287]:
ndrt2020.draw()
[289]:
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#
[290]:
# Flight
flight = Flight(
rocket=ndrt2020,
environment=env,
rail_length=parameters.get("rail_length")[0],
inclination=parameters.get("inclination")[0],
heading=parameters.get("heading")[0],
)
[291]:
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)
[292]:
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
[293]:
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
[294]:
# 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",
)
[295]:
Function.compare_plots(
plot_list=[(actual_vz, "Actual"), (actual_vz_filtered, "Filtered")],
ylabel="Vertical velocity (m/s)",
title="Vertical velocity",
xlabel="Time (s)",
)
[296]:
# 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",
)
[297]:
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
[298]:
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#
[299]:
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()
[300]:
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()
[301]:
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#
[302]:
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 %
[303]:
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 %
[304]:
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 %
[305]:
# 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