Bella Lui - EPFL - 2020

Bella Lui Kaltbrunn Mission from ERT (EPFL Rocket Team) Permission to use flight data given by Antoine Scardigli, 2020

%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.227 - 0.001, 0.010),  # propellant mass = 1.373
    # propulsion details
    "impulse": (2157, 0.03 * 2157),
    "burn_time": (2.43, 0.1),
    "nozzle_radius": (44.45 / 1000, 0.001),
    "throat_radius": (21.4376 / 1000, 0.001),
    "grain_separation": (3 / 1000, 1 / 1000),
    "grain_density": (782.4, 30),
    "grain_outer_radius": (85.598 / 2000, 0.001),
    "grain_initial_inner_radius": (33.147 / 1000, 0.002),
    "grain_initial_height": (152.4 / 1000, 0.001),
    # Aerodynamic Details
    "inertia_i": (0.78267, 0.03 * 0.78267),
    "inertia_z": (0.064244, 0.03 * 0.064244),
    "radius": (156 / 2000, 0.001),
    "distance_rocket_nozzle": (-1.1356, 0.100),
    "distance_rocket_propellant": (-1, 0.100),
    "power_off_drag": (1, 0.05),
    "power_on_drag": (1, 0.05),
    "nose_length": (0.242, 0.001),
    "nose_distance_to_cm": (1.3, 0.100),
    "fin_span": (0.200, 0.001),
    "fin_root_chord": (0.280, 0.001),
    "fin_tip_chord": (0.125, 0.001),
    "fin_distance_to_cm": (-0.75, 0.100),
    "tail_top_radius": (156 / 2000, 0.001),
    "tail_bottom_radius": (135 / 2000, 0.001),
    "tail_length": (0.050, 0.001),
    "tail_distance_to_cm": (-1.0856, 0.001),
    # Launch and Environment Details
    "wind_direction": (0, 5),
    "wind_speed": (1, 0.05),
    "inclination": (89, 1),
    "heading": (45, 5),
    "rail_length": (4.2, 0.001),
    # Parachute Details
    "CdS_drogue": (np.pi / 4, 0.20 * np.pi / 4),
    "lag_rec": (1, 0.020),
}

Environment

Define the Environment object

# Environment conditions
env = Environment(
    gravity=9.81,
    latitude=47.213476,
    longitude=9.003336,
    date=(2020, 2, 22, 13),
    elevation=407,
)

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

Visualize the Environment object

env.info()

Gravity Details

Acceleration of Gravity at Lauch Site: 9.81 m/s²


Launch Site Details

Launch Date: 2020-02-22 13:00:00 UTC
Launch Site Latitude: 47.21348°
Launch Site Longitude: 9.00334°
Reference Datum: SIRGAS2000
Launch Site UTM coordinates: 500252.61 E    5228887.37 N
Launch Site UTM zone: 32T
Launch Site Surface Elevation: 407.0 m


Atmospheric Model Details

Atmospheric Model Type: Reanalysis
Reanalysis Maximum Height: 4.378 km
Reanalysis Time Period: From  2020-02-22 00:00:00  to  2020-02-22 18:00:00  UTC
Reanalysis Hour Interval: 4  hrs
Reanalysis Latitude Range: From  48.0 ° To  46.0 °
Reanalysis Longitude Range: From  8.0 ° To  10.0 °


Surface Atmospheric Conditions

Surface Wind Speed: 1.26 m/s
Surface Wind Direction: 213.21°
Surface Wind Heading: 33.21°
Surface Pressure: 980.43 hPa
Surface Temperature: 286.63 K
Surface Air Density: 1.192 kg/m³
Surface Speed of Sound: 339.39 m/s


Atmospheric Model Plots

Motor

Define the SolidMotor object

k828fj = SolidMotor(
    thrust_source="../../tests/fixtures/acceptance/EPFL_Bella_Lui/bella_lui_motor_AeroTech_K828FJ.eng",
    burn_time=parameters.get("burn_time")[0],
    dry_mass=0.001,
    dry_inertia=(0, 0, 0),
    center_of_dry_mass_position=0.3,
    grains_center_of_mass_position=0.3,
    grain_number=3,
    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=0,
    coordinate_system_orientation="nozzle_to_combustion_chamber",
)

k828fj.plots.draw()

k828fj.info()

Nozzle Details
Nozzle Radius: 0.04445 m
Nozzle Throat Radius: 0.0214376 m

Grain Details
Number of Grains: 3
Grain Spacing: 0.003 m
Grain Density: 782.4 kg/m3
Grain Outer Radius: 0.042799 m
Grain Inner Radius: 0.033146999999999996 m
Grain Height: 0.1524 m
Grain Volume: 0.000 m3
Grain Mass: 0.275 kg

Motor Details
Total Burning Time: 2.43 s
Total Propellant Mass: 0.824 kg
Average Propellant Exhaust Velocity: 2514.035 m/s
Average Thrust: 852.260 N
Maximum Thrust: 1303.79 N at 0.04 s after ignition.
Total Impulse: 2070.992 Ns

Rocket

Create the Rocket object

bella_lui = 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=0.43,
    power_on_drag=0.43,
    center_of_mass_without_motor=0,
)
bella_lui.set_rail_buttons(0.1, -0.5)
bella_lui.add_motor(motor=k828fj, position=parameters.get("distance_rocket_nozzle")[0])

Adding aerodynamic surfaces

nose_cone = bella_lui.add_nose(
    length=parameters.get("nose_length")[0],
    kind="tangent",
    position=parameters.get("nose_distance_to_cm")[0]
    + parameters.get("nose_length")[0],
)
fin_set = bella_lui.add_trapezoidal_fins(
    3,
    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_distance_to_cm")[0],
)
tail = bella_lui.add_tail(
    top_radius=parameters.get("tail_top_radius")[0],
    bottom_radius=parameters.get("tail_bottom_radius")[0],
    length=parameters.get("tail_length")[0],
    position=parameters.get("tail_distance_to_cm")[0],
)

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

Adding Parachute

Drogue = bella_lui.add_parachute(
    "Drogue",
    cd_s=parameters.get("CdS_drogue")[0],
    trigger="apogee",
    sampling_rate=105,
    lag=parameters.get("lag_rec")[0],
    noise=(0, 8.3, 0.5),
)

Modify the Drag Coefficient curve

# Define aerodynamic drag coefficients
bella_lui.power_off_drag = Function(
    [
        (0.01, 0.51),
        (0.02, 0.46),
        (0.04, 0.43),
        (0.28, 0.43),
        (0.29, 0.44),
        (0.45, 0.44),
        (0.49, 0.46),
    ],
    "Mach Number",
    "Drag Coefficient with Power Off",
    "linear",
    "constant",
)
bella_lui.power_on_drag = Function(
    [
        (0.01, 0.51),
        (0.02, 0.46),
        (0.04, 0.43),
        (0.28, 0.43),
        (0.29, 0.44),
        (0.45, 0.44),
        (0.49, 0.46),
    ],
    "Mach Number",
    "Drag Coefficient with Power On",
    "linear",
    "constant",
)
bella_lui.power_off_drag *= parameters.get("power_off_drag")[0]
bella_lui.power_on_drag *= parameters.get("power_on_drag")[0]

bella_lui.info()

Inertia Details

Rocket Mass: 18.226 kg
Rocket Dry Mass: 18.227 kg (With Motor)
Rocket Mass: 19.051 kg (With Propellant)
Rocket Inertia (with motor, but without propellant) 11: 0.785 kg*m2
Rocket Inertia (with motor, but without propellant) 22: 0.785 kg*m2
Rocket Inertia (with motor, but without propellant) 33: 0.064 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.078 m
Rocket Frontal Area: 0.019113 m2

Rocket Distances
Rocket Center of Dry Mass - Center of Mass withour Motor: 0.000 m
Rocket Center of Dry Mass - Nozzle Exit Distance: 1.136 m
Rocket Center of Dry Mass - Center of Propellant Mass: 1.436 m
Rocket Center of Mass - Rocket Loaded Center of Mass: 0.062 m


Aerodynamics Lift Coefficient Derivatives

Nose Cone Lift Coefficient Derivative: 2.000/rad
Fins Lift Coefficient Derivative: 10.281/rad
Tail Lift Coefficient Derivative: -0.502/rad

Aerodynamics Center of Pressure

Nose Cone Center of Pressure to CM: 1.433 m
Fins Center of Pressure to CM: -0.871 m
Tail Center of Pressure to CM: -1.110 m
Distance - Center of Pressure to Center of Dry Mass: 0.407 m
Initial Static Margin: 2.610 c
Final Static Margin: 3.008 c


Parachute Details

Parachute Drogue with a cd_s of 0.7854 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

bella_lui.draw()

Flight

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

test_flight.info()
test_flight.plots.trajectory_3d()

Initial Conditions

Position - x: 0.00 m | y: 0.00 m | z: 407.00 m
Velocity - Vx: 0.00 m/s | Vy: 0.00 m/s | Vz: 0.00 m/s
Attitude - e0: 0.924 | e1: -0.008 | e2: 0.003 | e3: -0.383
Euler Angles - Spin φ : 0.00° | Nutation θ: -1.00° | Precession ψ: -45.00°
Angular Velocity - ω1: 0.00 rad/s | ω2: 0.00 rad/s| ω3: 0.00 rad/s


Surface Wind Conditions

Frontal Surface Wind Speed: 1.23 m/s
Lateral Surface Wind Speed: 0.26 m/s


Launch Rail

Launch Rail Length: 4.2  m
Launch Rail Inclination: 89.00°
Launch Rail Heading: 45.00°


Rail Departure State

Rail Departure Time: 0.359 s
Rail Departure Velocity: 16.185 m/s
Rail Departure Static Margin: 2.680 c
Rail Departure Angle of Attack: 4.457°
Rail Departure Thrust-Weight Ratio: 5.416
Rail Departure Reynolds Number: 1.691e+05


Burn out State

Burn out time: 2.430 s
Altitude at burn out: 122.980 m (AGL)
Rocket velocity at burn out: 85.150 m/s
Freestream velocity at burn out: 85.192 m/s
Mach Number at burn out: 0.251
Kinetic energy at burn out: 6.608e+04 J


Apogee State

Apogee Altitude: 867.503 m (ASL) | 460.503 m (AGL)
Apogee Time: 10.607 s
Apogee Freestream Speed: 3.568 m/s


Parachute Events

Drogue Ejection Triggered at: 10.610 s
Drogue Parachute Inflated at: 11.610 s
Drogue Parachute Inflated with Freestream Speed of: 10.104 m/s
Drogue Parachute Inflated at Height of: 455.668 m (AGL)


Impact Conditions

X Impact: 2.017 m
Y Impact: -5.890 m
Time of Impact: 35.819 s
Velocity at Impact: -19.573 m/s


Maximum Values

Maximum Speed: 86.201 m/s at 2.25 s
Maximum Mach Number: 0.254 Mach at 2.25 s
Maximum Reynolds Number: 8.915e+05 at 2.24 s
Maximum Dynamic Pressure: 4.384e+03 Pa at 2.25 s
Maximum Acceleration During Motor Burn: 58.470 m/s² at 0.04 s
Maximum Gs During Motor Burn: 5.962 g at 0.04 s
Maximum Acceleration After Motor Burn: -0.000 m/s² at 0.00 s
Maximum Gs After Motor Burn: -0.000 g at 0.00 s
Maximum Upper Rail Button Normal Force: 0.307 N
Maximum Upper Rail Button Shear Force: 0.468 N
Maximum Lower Rail Button Normal Force: 1.195 N
Maximum Lower Rail Button Shear Force: 1.826 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:  336
Number of Derivative Functions Evaluation:  1188
Average Function Evaluations per Time Step: 3.535714

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/EPFL_Bella_Lui/bella_lui_flight_data_filtered.csv",
    skiprows=1,
    delimiter=",",
    usecols=(2, 3, 4),
)

Convert to Function objects

actual_z = Function(
    source=np.column_stack((flight_data[:573, 0], flight_data[:573, 1])),
    inputs="Time (s)",
    outputs="Altitude above ground level (m)",
    interpolation="linear",
    extrapolation="zero",
)

actual_vz = Function(
    source=np.column_stack((flight_data[:573, 0], flight_data[:573, 2])),
    inputs="Time (s)",
    outputs="Vertical velocity (m/s)",
    interpolation="linear",
    extrapolation="zero",
)

# the actual acceleration will be calculated as the derivative of the actual velocity
actual_az = actual_vz.derivative_function()

We may need to filter the acceleration data to reduce the noise

az_filtered = savgol_filter(x=actual_az.source[:, 1], window_length=51, polyorder=3)

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

Function.compare_plots(
    [(actual_az, "actual"), (actual_az_filtered, "filtered")],
    xlabel="Time (s)",
    ylabel="Vertical acceleration (m/s^2)",
)

Get the simulated results

simulated_z = test_flight.z - env.elevation
simulated_vz = test_flight.vz
simulated_az = test_flight.az
simulated_t_final = test_flight.t_final
simulated_apogee = test_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, round(simulated_apogee, -2))
plt.xlim(0, round(simulated_t_final, -1))
plt.legend()
plt.grid()
plt.show()

plt.plot(actual_vz[:, 0], actual_vz[:, 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.ylim()
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.ylim()
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: 460.50 m
Real data: 458.97 m
Absolute error: 1.53 m
Relative error: 0.33 %

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

Max Velocity
RocketPy:  86.17 m/s
Real data: 90.00 m/s
Absolute error: -3.83 m/s
Relative error: 4.26 %

print("Max Acceleration")
print(f"RocketPy:  {simulated_az.max:.2f} m/s²")
print(f"Real data (derivative): {actual_az_filtered.max:.2f} m/s²")
acceleration_error = simulated_az.max - 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
RocketPy:  58.46 m/s²
Real data (derivative): 58.73 m/s²
Absolute error: -0.27 m/s^2
Relative error: 0.46 %