Getting Started¶
Here we go through a simplified rocket trajectory simulation to get you started. Let’s start by importing the rocketpy module.
[1]:
from rocketpy import Environment, SolidMotor, Rocket, Flight
If you are using Jupyter Notebooks, it is recommended to run the following line to make matplotlib plots which will be shown later interactive and higher quality.
[2]:
%matplotlib widget
Setting Up a Simulation¶
Creating an Environment for Spaceport America¶
[3]:
Env = Environment(
railLength=5.2, latitude=32.990254, longitude=-106.974998, elevation=1400
)
To get weather data from the GFS forecast, available online, we run the following lines.
First, we set tomorrow’s date.
[4]:
import datetime
tomorrow = datetime.date.today() + datetime.timedelta(days=1)
Env.setDate((tomorrow.year, tomorrow.month, tomorrow.day, 12)) # Hour given in UTC time
Then, we tell Env to use a GFS forecast to get the atmospheric conditions for flight.
Don’t mind the warning, it just means that not all variables, such as wind speed or atmospheric temperature, are available at all altitudes given by the forecast.
[5]:
Env.setAtmosphericModel(type="Forecast", file="GFS")
We can see what the weather will look like by calling the info method!
[6]:
Env.info()
Launch Site Details
Launch Rail Length: 5.2 m
Launch Date: 2022-11-13 12:00:00 UTC
Launch Site Latitude: 32.99025°
Launch Site Longitude: -106.97500°
Reference Datum: SIRGAS2000
Launch Site UTM coordinates: 315468.64 W 3651938.65 N
Launch Site UTM zone: 13S
Launch Site Surface Elevation: 1471.5 m
Atmospheric Model Details
Atmospheric Model Type: Forecast
Forecast Maximum Height: 78.547 km
Forecast Time Period: From 2022-11-12 18:00:00 to 2022-11-28 18:00:00 UTC
Forecast Hour Interval: 3 hrs
Forecast Latitude Range: From -90.0 ° To 90.0 °
Forecast Longitude Range: From 0.0 ° To 359.75 °
Surface Atmospheric Conditions
Surface Wind Speed: 3.68 m/s
Surface Wind Direction: 149.63°
Surface Wind Heading: 329.63°
Surface Pressure: 851.87 hPa
Surface Temperature: 280.20 K
Surface Air Density: 1.059 kg/m³
Surface Speed of Sound: 335.57 m/s
Atmospheric Model Plots
Creating a Motor¶
A solid rocket motor is used in this case. To create a motor, the SolidMotor class is used and the required arguments are given.
The SolidMotor class requires the user to have a thrust curve ready. This can come either from a .eng file for a commercial motor, such as below, or a .csv file from a static test measurement.
Besides the thrust curve, other parameters such as grain properties and nozzle dimensions must also be given.
[7]:
Pro75M1670 = SolidMotor(
thrustSource="../../data/motors/Cesaroni_M1670.eng",
burnOut=3.9,
grainNumber=5,
grainSeparation=5 / 1000,
grainDensity=1815,
grainOuterRadius=33 / 1000,
grainInitialInnerRadius=15 / 1000,
grainInitialHeight=120 / 1000,
nozzleRadius=33 / 1000,
throatRadius=11 / 1000,
interpolationMethod="linear",
)
To see what our thrust curve looks like, along with other import properties, we invoke the info method yet again. You may try the allInfo method if you want more information all at once!
[8]:
Pro75M1670.info()
Motor Details
Total Burning Time: 3.9 s
Total Propellant Mass: 2.956 kg
Propellant Exhaust Velocity: 2038.745 m/s
Average Thrust: 1545.218 N
Maximum Thrust: 2200.0 N at 0.15 s after ignition.
Total Impulse: 6026.350 Ns
Plots
Creating a Rocket¶
A rocket is composed of several components. Namely, we must have a motor (good thing we have the Pro75M1670 ready), a couple of aerodynamic surfaces (nose cone, fins and tail) and parachutes (if we are not launching a missile).
Let’s start by initializing our rocket, named Calisto, supplying it with the Pro75M1670 engine, entering its inertia properties, some dimensions and also its drag curves.
[9]:
Calisto = Rocket(
motor=Pro75M1670,
radius=127 / 2000,
mass=19.197 - 2.956,
inertiaI=6.60,
inertiaZ=0.0351,
distanceRocketNozzle=-1.255,
distanceRocketPropellant=-0.85704,
powerOffDrag="../../data/calisto/powerOffDragCurve.csv",
powerOnDrag="../../data/calisto/powerOnDragCurve.csv",
)
Calisto.setRailButtons([0.2, -0.5])
Adding Aerodynamic Surfaces¶
Now we define the aerodynamic surfaces. They are really straight forward.
[10]:
NoseCone = Calisto.addNose(length=0.55829, kind="vonKarman", distanceToCM=0.71971)
FinSet = Calisto.addTrapezoidalFins(
n=4,
rootChord=0.120,
tipChord=0.040,
span=0.100,
distanceToCM=-1.04956,
cantAngle=0,
radius=None,
airfoil=None,
)
Tail = Calisto.addTail(
topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656
)
[11]:
Calisto.allInfo()
Inertia Details
Rocket Mass: 16.241 kg (No Propellant)
Rocket Mass: 19.197 kg (With Propellant)
Rocket Inertia I: 6.600 kg*m2
Rocket Inertia Z: 0.035 kg*m2
Geometrical Parameters
Rocket Maximum Radius: 0.0635 m
Rocket Frontal Area: 0.012668 m2
Rocket Distances
Rocket Center of Mass - Nozzle Exit Distance: -1.255 m
Rocket Center of Mass - Motor reference point: -0.85704 m
Rocket Center of Mass - Rocket Loaded Center of Mass: -0.132 m
Aerodynamic Components Parameters
Currently not implemented.
Aerodynamics Lift Coefficient Derivatives
Nose Cone Lift Coefficient Derivative: 2.000/rad
Fins Lift Coefficient Derivative: 5.603/rad
Tail Lift Coefficient Derivative: -1.061/rad
Aerodynamics Center of Pressure
Nose Cone Center of Pressure to CM: 0.999 m
Fins Center of Pressure to CM: -1.084 m
Tail Center of Pressure to CM: -1.223 m
Distance - Center of Pressure to CM: -0.425 m
Initial Static Margin: 2.306 c
Final Static Margin: 3.345 c
Mass Plots
Aerodynamics Plots
Adding Parachutes¶
Finally, we have parachutes! Calisto will have two parachutes, Drogue and Main.
Both parachutes are activated by some special algorithm, which is usually really complex and a trade secret. Most algorithms are based on pressure sampling only, while some also use acceleration info.
RocketPy allows you to define a trigger function which will decide when to activate the ejection event for each parachute. This trigger function is supplied with pressure measurement at a predefined sampling rate. This pressure signal is usually noisy, so artificial noise parameters can be given. Call help(Rocket.addParachute) for more details. Furthermore, the trigger function also receives the complete state vector of the rocket, allowing us to use velocity, acceleration or even attitude to decide when the parachute event should be triggered.
Here, we define our trigger functions rather simply using Python. However, you can call the exact code which will fly inside your rocket as well.
[12]:
def drogueTrigger(p, y):
# p = pressure
# y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]
# activate drogue when vz < 0 m/s.
return True if y[5] < 0 else False
def mainTrigger(p, y):
# p = pressure
# y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]
# activate main when vz < 0 m/s and z < 800 + 1400 m (+1400 due to surface elevation).
return True if y[5] < 0 and y[2] < 800 + 1400 else False
Main = Calisto.addParachute(
"Main",
CdS=10.0,
trigger=mainTrigger,
samplingRate=105,
lag=1.5,
noise=(0, 8.3, 0.5),
)
Drogue = Calisto.addParachute(
"Drogue",
CdS=1.0,
trigger=drogueTrigger,
samplingRate=105,
lag=1.5,
noise=(0, 8.3, 0.5),
)
Just be careful if you run this last cell multiple times! If you do so, your rocket will end up with lots of parachutes which activate together, which may cause problems during the flight simulation. We advise you to re-run all cells which define our rocket before running this, preventing unwanted old parachutes. Alternatively, you can run the following lines to remove parachutes.
Calisto.parachutes.remove(Drogue)
Calisto.parachutes.remove(Main)
Simulating a Flight¶
Simulating a flight trajectory is as simple as initializing a Flight class object givin the rocket and environnement set up above as inputs. The launch rail inclination and heading are also given here.
[13]:
TestFlight = Flight(rocket=Calisto, environment=Env, inclination=85, heading=0)
Analyzing the Results¶
RocketPy gives you many plots, thats for sure! They are divided into sections to keep them organized. Alternatively, see the Flight class documentation to see how to get plots for specific variables only, instead of all of them at once.
[14]:
TestFlight.allInfo()
Initial Conditions
Position - x: 0.00 m | y: 0.00 m | z: 1471.47 m
Velocity - Vx: 0.00 m/s | Vy: 0.00 m/s | Vz: 0.00 m/s
Attitude - e0: 0.999 | e1: -0.044 | e2: -0.000 | e3: 0.000
Euler Angles - Spin φ : 0.00° | Nutation θ: -5.00° | Precession ψ: 0.00°
Angular Velocity - ω1: 0.00 rad/s | ω2: 0.00 rad/s| ω3: 0.00 rad/s
Launch Rail Orientation
Launch Rail Inclination: 85.00°
Launch Rail Heading: 0.00°
Surface Wind Conditions
Frontal Surface Wind Speed: 3.17 m/s
Lateral Surface Wind Speed: 1.86 m/s
Rail Departure State
Rail Departure Time: 0.363 s
Rail Departure Velocity: 25.800 m/s
Rail Departure Static Margin: 2.389 c
Rail Departure Angle of Attack: 8.179°
Rail Departure Thrust-Weight Ratio: 10.143
Rail Departure Reynolds Number: 1.980e+05
BurnOut State
BurnOut time: 3.900 s
Altitude at burnOut: 659.413 m (AGL)
Rocket velocity at burnOut: 279.843 m/s
Freestream velocity at burnOut: 279.497 m/s
Mach Number at burnOut: 0.835
Kinetic energy at burnOut: 6.359e+05 J
Apogee
Apogee Altitude: 4819.613 m (ASL) | 3348.147 m (AGL)
Apogee Time: 26.043 s
Apogee Freestream Speed: 10.755 m/s
Events
Drogue Ejection Triggered at: 26.048 s
Drogue Parachute Inflated at: 27.548 s
Drogue Parachute Inflated with Freestream Speed of: 17.781 m/s
Drogue Parachute Inflated at Height of: 3337.285 m (AGL)
Main Ejection Triggered at: 162.790 s
Main Parachute Inflated at: 164.290 s
Main Parachute Inflated with Freestream Speed of: 18.079 m/s
Main Parachute Inflated at Height of: 701.402 m (AGL)
Impact
X Impact: 1661.987 m
Y Impact: 1716.832 m
Time of Impact: 288.863 s
Velocity at Impact: -5.483 m/s
Maximum Values
Maximum Speed: 285.979 m/s at 3.39 s
Maximum Mach Number: 0.852 Mach at 3.39 s
Maximum Reynolds Number: 2.073e+06 at 3.30 s
Maximum Dynamic Pressure: 4.073e+04 Pa at 3.34 s
Maximum Acceleration: 105.100 m/s² at 0.15 s
Maximum Gs: 10.717 g at 0.15 s
Maximum Upper Rail Button Normal Force: 0.452 N
Maximum Upper Rail Button Shear Force: 1.738 N
Maximum Lower Rail Button Normal Force: 0.257 N
Maximum Lower Rail Button Shear Force: 0.863 N
Numerical Integration Information
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: 702
Number of Derivative Functions Evaluation: 2050
Average Function Evaluations per Time Step: 2.920228
Trajectory 3d Plot
Trajectory Kinematic Plots
Angular Position Plots
Path, Attitude and Lateral Attitude Angle plots
Trajectory Angular Velocity and Acceleration Plots
Trajectory Force Plots
Trajectory Energy Plots
Trajectory Fluid Mechanics Plots
Trajectory Stability and Control Plots
Export Flight Trajectory to a .kml file so it can be opened on Google Earth
[15]:
TestFlight.exportKML(
fileName="trajectory.kml",
extrude=True,
altitudeMode="relativetoground",
)
File trajectory.kml saved with success!
Using Simulation for Design¶
Here, we go through a couple of examples which make use of RocketPy in cool ways to help us design our rocket.
Dynamic Stability Analysis¶
Ever wondered how static stability translates into dynamic stability? Different static margins result in different dynamic behavior, which also depends on the rocket’s rotational inertial.
Let’s make use of RocketPy’s helper class called Function to explore how the dynamic stability of Calisto varies if we change the fins span by a certain factor.
[16]:
# Helper class
from rocketpy import Function
# Prepare Rocket Class
Calisto = Rocket(
motor=Pro75M1670,
radius=127 / 2000,
mass=19.197 - 2.956,
inertiaI=6.60,
inertiaZ=0.0351,
distanceRocketNozzle=-1.255,
distanceRocketPropellant=-0.85704,
powerOffDrag="../../data/calisto/powerOffDragCurve.csv",
powerOnDrag="../../data/calisto/powerOnDragCurve.csv",
)
Calisto.setRailButtons([0.2, -0.5])
Nose = Calisto.addNose(length=0.55829, kind="vonKarman", distanceToCM=0.71971)
FinSet = Calisto.addTrapezoidalFins(
4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956
)
Tail = Calisto.addTail(
topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656
)
# Prepare Environment Class
Env = Environment(5.2, 9.8)
Env.setAtmosphericModel(type="CustomAtmosphere", wind_v=-5)
# Simulate Different Static Margins by Varying Fin Position
simulation_results = []
for factor in [0.5, 0.7, 0.9, 1.1, 1.3]:
# Modify rocket fin set by removing previous one and adding new one
Calisto.aerodynamicSurfaces.remove(FinSet)
FinSet = Calisto.addTrapezoidalFins(
4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956*factor
)
# Simulate
print(
"Simulating Rocket with Static Margin of {:1.3f}->{:1.3f} c".format(
Calisto.staticMargin(0), Calisto.staticMargin(Calisto.motor.burnOutTime)
)
)
TestFlight = Flight(
rocket=Calisto,
environment=Env,
inclination=90,
heading=0,
maxTimeStep=0.01,
maxTime=5,
terminateOnApogee=True,
verbose=True,
)
# Post process flight data
TestFlight.postProcess()
# Store Results
staticMarginAtIgnition = Calisto.staticMargin(0)
staticMarginAtOutOfRail = Calisto.staticMargin(TestFlight.outOfRailTime)
staticMarginAtSteadyState = Calisto.staticMargin(TestFlight.tFinal)
simulation_results += [
(
TestFlight.attitudeAngle,
"{:1.2f} c | {:1.2f} c | {:1.2f} c".format(
staticMarginAtIgnition,
staticMarginAtOutOfRail,
staticMarginAtSteadyState,
),
)
]
Function.comparePlots(
simulation_results,
lower=0,
upper=1.5,
xlabel="Time (s)",
ylabel="Attitude Angle (deg)",
)
Simulating Rocket with Static Margin of -1.444->-0.405 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of -0.046->0.993 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of 1.352->2.391 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of 2.750->3.789 c
Simulation Completed at Time: 5.0000 s
Simulating Rocket with Static Margin of 4.147->5.186 c
Simulation Completed at Time: 5.0000 s
Characteristic Frequency Calculation¶
Here we analyze the characteristic frequency of oscillation of our rocket just as it leaves the launch rail. Note that when we ran TestFlight.allInfo(), one of the plots already showed us the frequency spectrum of our flight. Here, however, we have more control of what we are plotting.
[17]:
import numpy as np
import matplotlib.pyplot as plt
Env = Environment(
railLength=5.2, latitude=32.990254, longitude=-106.974998, elevation=1400
)
Env.setAtmosphericModel(type="CustomAtmosphere", wind_v=-5)
# Prepare Motor
Pro75M1670 = SolidMotor(
thrustSource="../../data/motors/Cesaroni_M1670.eng",
burnOut=3.9,
grainNumber=5,
grainSeparation=5 / 1000,
grainDensity=1815,
grainOuterRadius=33 / 1000,
grainInitialInnerRadius=15 / 1000,
grainInitialHeight=120 / 1000,
nozzleRadius=33 / 1000,
throatRadius=11 / 1000,
interpolationMethod="linear",
)
# Prepare Rocket
Calisto = Rocket(
motor=Pro75M1670,
radius=127 / 2000,
mass=19.197 - 2.956,
inertiaI=6.60,
inertiaZ=0.0351,
distanceRocketNozzle=-1.255,
distanceRocketPropellant=-0.85704,
powerOffDrag="../../data/calisto/powerOffDragCurve.csv",
powerOnDrag="../../data/calisto/powerOnDragCurve.csv",
)
Calisto.setRailButtons([0.2, -0.5])
Nose = Calisto.addNose(length=0.55829, kind="vonKarman", distanceToCM=0.71971)
FinSet = Calisto.addTrapezoidalFins(
4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956
)
Tail = Calisto.addTail(
topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656
)
# Simulate first 5 seconds of Flight
TestFlight = Flight(
rocket=Calisto,
environment=Env,
inclination=90,
heading=0,
maxTimeStep=0.01,
maxTime=5,
)
TestFlight.postProcess()
# Perform a Fourier Analysis
Fs = 100.0
# sampling rate
Ts = 1.0 / Fs
# sampling interval
t = np.arange(1, 400, Ts) # time vector
ff = 5
# frequency of the signal
y = TestFlight.attitudeAngle(t) - np.mean(TestFlight.attitudeAngle(t))
n = len(y) # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T # two sides frequency range
frq = frq[range(n // 2)] # one side frequency range
Y = np.fft.fft(y) / n # fft computing and normalization
Y = Y[range(n // 2)]
fig, ax = plt.subplots(2, 1)
ax[0].plot(t, y)
ax[0].set_xlabel("Time")
ax[0].set_ylabel("Signal")
ax[0].set_xlim((0, 5))
ax[1].plot(frq, abs(Y), "r") # plotting the spectrum
ax[1].set_xlabel("Freq (Hz)")
ax[1].set_ylabel("|Y(freq)|")
ax[1].set_xlim((0, 5))
plt.subplots_adjust(hspace=0.5)
plt.show()
Apogee as a Function of Mass¶
This one is a classic one! We always need to know how much our rocket’s apogee will change when our payload gets heavier.
[18]:
def apogee(mass):
# Prepare Environment
Env = Environment(
railLength=5.2,
latitude=32.990254,
longitude=-106.974998,
elevation=1400,
date=(2018, 6, 20, 18),
)
Env.setAtmosphericModel(type="CustomAtmosphere", wind_v=-5)
# Prepare Motor
Pro75M1670 = SolidMotor(
thrustSource="../../data/motors/Cesaroni_M1670.eng",
burnOut=3.9,
grainNumber=5,
grainSeparation=5 / 1000,
grainDensity=1815,
grainOuterRadius=33 / 1000,
grainInitialInnerRadius=15 / 1000,
grainInitialHeight=120 / 1000,
nozzleRadius=33 / 1000,
throatRadius=11 / 1000,
interpolationMethod="linear",
)
# Prepare Rocket
Calisto = Rocket(
motor=Pro75M1670,
radius=127 / 2000,
mass=mass,
inertiaI=6.60,
inertiaZ=0.0351,
distanceRocketNozzle=-1.255,
distanceRocketPropellant=-0.85704,
powerOffDrag="../../data/calisto/powerOffDragCurve.csv",
powerOnDrag="../../data/calisto/powerOnDragCurve.csv",
)
Calisto.setRailButtons([0.2, -0.5])
Nose = Calisto.addNose(length=0.55829, kind="vonKarman", distanceToCM=0.71971)
FinSet = Calisto.addTrapezoidalFins(
4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956
)
Tail = Calisto.addTail(
topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656
)
# Simulate Flight until Apogee
TestFlight = Flight(
rocket=Calisto,
environment=Env,
inclination=85,
heading=0,
terminateOnApogee=True,
)
return TestFlight.apogee
apogeebymass = Function(apogee, inputs="Mass (kg)", outputs="Estimated Apogee (m)")
apogeebymass.plot(8, 20, 20)
c:\users\ghceo\repos\rocketpy\rocketpy\Function.py:1043: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
fig = plt.figure()
Out of Rail Speed as a Function of Mass¶
To finish off, lets make a really important plot. Out of rail speed is the speed our rocket has when it is leaving the launch rail. This is crucial to make sure it can fly safely after leaving the rail. A common rule of thumb is that our rocket’s out of rail speed should be 4 times the wind speed so that it does not stall and become unstable.
[19]:
def speed(mass):
# Prepare Environment
Env = Environment(
railLength=5.2,
latitude=32.990254,
longitude=-106.974998,
elevation=1400,
date=(2018, 6, 20, 18),
)
Env.setAtmosphericModel(type="CustomAtmosphere", wind_v=-5)
# Prepare Motor
Pro75M1670 = SolidMotor(
thrustSource="../../data/motors/Cesaroni_M1670.eng",
burnOut=3.9,
grainNumber=5,
grainSeparation=5 / 1000,
grainDensity=1815,
grainOuterRadius=33 / 1000,
grainInitialInnerRadius=15 / 1000,
grainInitialHeight=120 / 1000,
nozzleRadius=33 / 1000,
throatRadius=11 / 1000,
interpolationMethod="linear",
)
# Prepare Rocket
Calisto = Rocket(
motor=Pro75M1670,
radius=127 / 2000,
mass=mass,
inertiaI=6.60,
inertiaZ=0.0351,
distanceRocketNozzle=-1.255,
distanceRocketPropellant=-0.85704,
powerOffDrag="../../data/calisto/powerOffDragCurve.csv",
powerOnDrag="../../data/calisto/powerOnDragCurve.csv",
)
Calisto.setRailButtons([0.2, -0.5])
Nose = Calisto.addNose(length=0.55829, kind="vonKarman", distanceToCM=0.71971)
FinSet = Calisto.addTrapezoidalFins(
4, span=0.100, rootChord=0.120, tipChord=0.040, distanceToCM=-1.04956
)
Tail = Calisto.addTail(
topRadius=0.0635, bottomRadius=0.0435, length=0.060, distanceToCM=-1.194656
)
# Simulate Flight until Apogee
TestFlight = Flight(
rocket=Calisto,
environment=Env,
inclination=85,
heading=0,
terminateOnApogee=True,
)
return TestFlight.outOfRailVelocity
speedbymass = Function(speed, inputs="Mass (kg)", outputs="Out of Rail Speed (m/s)")
speedbymass.plot(8, 20, 20)