In [1]:
%matplotlib inline
import math
from types import SimpleNamespace
import numpy as np
import sympy
from sympy import Symbol, solve
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from copy import deepcopy as DC
sympy.init_printing()
Origin: On the ground below plane at $t_0$.
X axis [0]: 90 degrees clockwise of Y axis.
Y axis [1]: Direction of helicopter ground track.
Z axis [2]: Pointed up vertically.
Absolute velocity (for determining position) is with respect to the ground.
Relative velocity (for determining drag) is with respect to the wind.
East: 0 degrees, i.e. direction of positive X axis.
North: 90 degrees, i.e. direction of positive Y axis.
Note, these conventions are somewhat arbitrary and selected for convenience. In the Monte Carlo simulations, we utilize a different set of conventions that corresponds to the perspective of an aviation pilot.
| Variable | Symbol | Meaning |
|---|---|---|
| g | $g$ | standard gravity (m/s^2) |
| rho | $\rho$ | density of air (kg/m^3) |
| epsilon | $\epsilon$ | an arbitrarily small number |
| m | $m$ | mass of system (kg) |
| drop_height | --- | height system is dropped from (m) |
| sl | $sl$ | length of static line (m) |
| d_deploy_t | --- | drogue deployment timer (s) |
| r_d | $r_d$ | radius of drogue chute (m) |
| dl | $dl$ | length of drogue line (m) |
| vt_d | $vt_d$ | terminal velocity of drogue (m/s) |
| c_deploy_t | --- | main chute deployment timer (s) |
| r_c | $r_c$ | radius of main chute (m) |
| C_d_c | $C_{d,c}$ | drag coefficient of main chute |
| C_d_dud | $C_{d,dud}$ | drag coefficient of dud system |
| A_dud | $A_{dud}$ | frontal area of dud system (m^2) |
You can modify assumptions on the fly in this document. It runs a sim with default values and graphs it automatically. Keep in mind that you can use main() with its arguments to specify some of the variables of the model, but some variables are hardcoded as parameters. Also remember that this tool has been broken into many functions for a reason. We will import it as a library in another document and run many different sims to compare, so it has been coded with that in mind.
In [2]:
# g = acceleration due to gravity [kg*m/s^2]
g = -9.80665
# density of air [kg/m^3]
rho = 1.2 # rough estimate
# simulation timestep [s]
epsilon = 10**(-10) # just in case we want a small number
time_step = 0.25
120in. Rocketman High Performance CD 2.2 Parachute. This is our main parachute.
https://the-rocketman.com/high-performance-parachutes/
We use this site to eyeball the terminal velocity of the drogue parachute.
https://web.archive.org/web/20181110101551/http://the-rocketman.com/drogue-decent-rate.html
In [3]:
# mass of payload [kg]
m = 30. # needs a sanity check but this determines almost everything to follow
drop_height = 609.6 # [m] (2000 ft)
# static line length [m]
sl = 2.
# drogue deployment timer [s]
d_deploy_t = 4.
# radius of drogue [m]
r_d = 0.762
# drogue line length [m]
dl = 50. * 0.3048 # feet, so convert to meters
# terminal velocity of drogue [m/s]
vt_d = -18.5 # extremely rough eyeball estimate from rocketman
# main chute deployment timer [s]
c_deploy_t = 10.
# radius of main chute [m]
r_c = 3.048 / 2 # [m] (120 in diameter)
# drag coefficients of main chute [m/s]
C_d_c = 2.2 # from rocketman, or your money back
C_d_dud = 0.3 # this is entirely an arbitrary guess. kinda close to a car, kinda close to lv4
A_dud = np.pi* 0.127**2 # [m] (5 in radius) i pulled this number out of a hat but we could probably just measure it
Source [http://www.pcprg.com/rocketre.htm]
In [4]:
weight = abs(m * g) # drag force at terminal velocity [N]
A_d = np.pi * r_d**2 # frontal area of drogue [m^2]
C_d_d = (2 * weight) / (rho * A_d * vt_d**2) # drag coeff. of drogue []
A_c = np.pi * r_c**2 # frontal area of main chute [m^2]
def term_vel(C_d, A):
return -np.sqrt(2* weight/(C_d * rho *A))
vt_c = term_vel(C_d_c, A_c)
In [5]:
def m_to_ft(x): # convert meters to feet
return 3.28084 * x
def ms_to_knots(v): # convert m/s to knots
return 1.943844 * v
# passively rotates a coordinate frame with a fixed vector. to go from relative to absolute
def rotate(x, theta):
R = np.array([[np.cos(theta), - np.sin(theta)],
[np.sin(theta), np.cos(theta)]])
return R.dot(x)
# https://en.wikipedia.org/wiki/Wind_profile_power_law
# assuming stable atmosphere over land without too rough
# of terrain and wind measurements at 10 m
def wind_profile(z, u_0, z_0=10, alpha=0.143):
if z < 0: z = 0
return u_0 * (z / z_0)**alpha
# calculates drag in all directions
# note that i am assuming that drag splits up in all dimensions equally
# idk if i should be assuming something about laminar flow in horizontal plane
# because i just learned what that even was
def drag_forces(v, C_d, A):
norm = np.linalg.norm(v)
return (-A/2 * C_d * rho * norm) * v
# solve a quartic to find the time at which the drogue chute becomes taut
# i'll make the names more iconic later, for now its fine
# at some point i'll do some actual math and make sure this can only return one positive real instead of assuming
def quartic(v_now):
Ps, vds, gs, Lvds, Lhds, tds, vcs, side = sympy.symbols('Ps vds gs Lvds Lhds tds vcs side')
Dparms = {side: v_now[0], Ps: v_now[1], vcs: v_now[2], gs: g}
tdEqn = (side*tds)**2 + (Ps*tds)**2 + (vcs*tds + 0.5*gs*tds**2)**2 - dl**2
tdSolns = sympy.solve(tdEqn.subs(Dparms))
for soln in [complex(x) for x in tdSolns]:
if (soln.imag > epsilon) or (soln.real <= 0):
#if (soln.imag != 0) or (soln.real <= 0):
pass
else:
td = soln.real
return td
# func for integrating. dis our ODE
# i know it's kinda weird to put logic inside an ode but it made sense at the time
def ode_f(x, v, t, phs, st8, v_w):
wind = np.array([wind_profile(x[2], v_wi) for v_wi in v_w])
if phs['live']: # drogue or chute active
if phs['steadystate']: # aerodynamic equilibrium
rel_v = np.zeros(3)
abs_v = wind
rel_v[2], abs_v[2] = (vt_c, vt_c) if phs['mainchute'] else (vt_d, vt_d)
else:
rel_v = v - wind
abs_v = v
# remember that drag depends on velocity of object relative to medium
D = drag_forces(rel_v, st8[0], st8[1]) / m
elif not phs['live']:
rel_v = v - wind
D = drag_forces(rel_v, C_d_dud, A_dud) / m # took a stab at guessing dud parameters, plz sanitize eventually
abs_v = v
a = D + np.array([0, 0, g])
return np.array([abs_v, a])
# runge kutta integrator
def runge_kutta(sim, next_step):
dt = sim.dt
k1 = ode_f(sim.x[-1],
sim.v[-1],
sim.t[-1],
sim.phase,
sim.state,
sim.v_w)
k2 = ode_f(sim.x[-1] + k1[0] * dt/2,
sim.v[-1] + k1[1] * dt/2,
sim.t[-1] + dt/2,
sim.phase,
sim.state,
sim.v_w)
k3 = ode_f(sim.x[-1] + k2[0] * dt/2,
sim.v[-1] + k2[1] * dt/2,
sim.t[-1] + dt/2,
sim.phase,
sim.state,
sim.v_w)
k4 = ode_f(sim.x[-1] + k3[0] * dt,
sim.v[-1] + k3[1] * dt,
sim.t[-1] + dt,
sim.phase,
sim.state,
sim.v_w)
next_step.x = sim.x[-1] + dt * (k1[0] + 2*k2[0] + 2*k3[0] + k4[0]) / 6
next_step.v = sim.v[-1] + dt * (k1[1] + 2*k2[1] + 2*k3[1] + k4[1]) / 6
next_step.t = sim.t[-1] + dt
next_step.a = (k1[1] + 2*k2[1] + 2*k3[1] + k4[1]) / 6
# keep eyebleed hidden
def update(sim, next_step):
sim.x = np.append(sim.x, [next_step.x], axis=0)
sim.v = np.append(sim.v, [next_step.v], axis=0)
sim.t = np.append(sim.t, next_step.t)
sim.a = np.append(sim.a, [next_step.a], axis=0)
# this is gonna get called during each phase for some amount of time
def iterate(sim, next_step, time):
start = sim.t[-1] # offset
# pretttttyy logic
while (not sim.phase['mainchute'] and sim.t[-1] < start+time) or (sim.phase['mainchute'] and sim.x[-1][2] > 0):
runge_kutta(sim, next_step)
update(sim, next_step)
if sim.x[-1][2] <= 0: break # kind of a kludge in case drogue fails and we impact before main deploys
if sim.assume_limit: # in case we want to check whether we're close to terminal velocity
terminal = vt_c if sim.phase['mainchute'] else vt_d
if abs(next_step.v[2] - terminal) < epsilon: sim.phase['steadystate'] = True # jump to steady state
# no interesting code in this func, just graphing stuff
def graph_sol(sim, start=0, end=None):
# break out the solutions into convenient names:
time = sim.t[start:end]
soln_d_x = [s[0] for s in sim.x][start:end]
soln_d_y = [s[1] for s in sim.x][start:end]
soln_d_z = [s[2] for s in sim.x][start:end]
soln_d_xdot = [s[0] for s in sim.v][start:end]
soln_d_ydot = [s[1] for s in sim.v][start:end]
soln_d_zdot = [s[2] for s in sim.v][start:end]
soln_d_xddot = np.diff(soln_d_xdot) # x acceleration
soln_d_yddot = np.diff(soln_d_ydot) # y acceleration
soln_d_zddot = np.diff(soln_d_zdot) # z acceleration
# plot da shiz:
plt.figure(1)
plt.plot(soln_d_y, soln_d_z)
plt.axis('equal')
plt.grid()
plt.title('down-range trajectory')
plt.xlabel('horizontal range (m)')
plt.ylabel('vertical range (m)')
plt.figure(2)
plt.plot(time, soln_d_ydot)
plt.xlabel('time (s)')
plt.ylabel('horizontal velocity (m/s)')
plt.figure(3)
plt.plot(time, soln_d_zdot)
plt.xlabel('time (s)')
plt.ylabel('vertical velocity (m/s)')
plt.figure(4)
plt.plot(time[:-1], soln_d_yddot)
plt.xlabel('time (s)')
plt.ylabel('horizontal acceleration (m/s^2)')
plt.figure(5)
plt.plot(time[:-1], soln_d_zddot)
plt.xlabel('time (s)')
plt.ylabel('vertical acceleration (m/s^2)')
plt.figure(6)
plt.plot(sim.x[0][0], sim.x[0][1],'r^')
plt.plot(sim.x[-1][0], sim.x[-1][1],'bs')
plt.xlabel('First and Last X pos')
plt.ylabel('First and Last Y pos')
plt.grid()
plt.axis('equal')
plt.figure(7)
plt.plot(soln_d_y, soln_d_ydot)
plt.xlabel('horizontal range (m)')
plt.ylabel('horizontal velocity (m/s)')
plt.figure(8)
plt.plot(soln_d_ydot[:-1], soln_d_yddot)
plt.xlabel('horizontal velocity (m/s)')
plt.ylabel('horizontal acceleration (m/s^2)')
plt.figure(9)
plt.plot(sim.t, sim.a)
plt.xlabel('time (s)')
plt.ylabel('acceleration (m/s^2)')
plt.legend(['x', 'y', 'z'])
In [6]:
phase = {'live': False, 'steadystate': False, 'mainchute': False}
fail_mode = {'drogue': False, 'mainchute': False, 'early': False}
# we're using this to set conditions in the model based on an abstract code (i.e. A,B,C,D)
def set_failure_mode(sim, mode):
sim.fail_mode = DC(fail_mode)
if mode == 'A':
sim.fail_mode['mainchute'] = True
if sim.prnt: print('Main Chute Failure!\n')
elif mode == 'B':
sim.fail_mode['drogue'] = True
if sim.prnt: print('Drogue Failure!\n')
elif mode == 'C':
sim.fail_mode['drogue'] = True
sim.fail_mode['early'] = True
if sim.prnt: print('Drogue Failure and Early Main Chute Deployment!\n')
elif mode == 'D': # d is for dud, it's good enough for me
sim.fail_mode['mainchute'] = True
sim.fail_mode['drogue'] = True
if sim.prnt: print('Recovery system is a dud! Now Im freee. Free fallin!\n')
else:
if sim.prnt: print('No Failures detected!\n')
# in case we want to track the simulation as it occurs
def test_print(next_step):
print('t: ',next_step.t, '. x: ', next_step.x, '. v: ', next_step.v)
# initialize a trajectory simulation with all relevant parameters
def init(v_plane=18, v_side=0.1, deg_plane=90+45, v_wind=3.1, deg_wind=360-45, x_0_offset=[0, 0, 0], t_0=0,
mode='', dt=time_step, prnt=False):
# v_plane = speed of plane (m/s)
# deg_plane = bearing of plane, east to north (degrees)
# v_side = speed system is thrown right out of airlock at (m/s)
# v_wind = average wind speed (m/s)
# deg_wind = bearing of wind, east to north (degrees)
# bearing of system referencing plane (degrees)
deg_side = deg_plane - 90
# bearing of plane in radians
theta_plane = np.radians(deg_plane)
# bearing of system in radians
theta_side = np.radians(deg_side)
# bearing of wind in radians
theta_wind = np.radians(deg_wind)
# angle between wind and plane (rad)
theta = theta_wind - theta_plane
# [cross-wind left to right (from pilot POV), wind in flight direction, vertical wind]
dir_w = np.array([np.sin(-theta), np.cos(theta), 0])
next_step = SimpleNamespace()
next_step.x, next_step.v, next_step.t, next_step.a = np.zeros(3), np.zeros(3), np.zeros(1), np.zeros(3)
sim = SimpleNamespace()
sim.v_w = dir_w * v_wind
# here is our state space!
sim.x = np.array([np.array([0, 0, drop_height]) + np.array(x_0_offset)])
sim.v = np.array([np.array([v_side, v_plane, 0])])
sim.t = np.array([t_0])
sim.a = np.array([np.array([0,0,g])])
sim.dt = dt
sim.prnt = prnt # silent by default
sim.state = [C_d_d, A_d] # drogue armed
sim.phase = DC(phase) # keep track of boolean states
set_failure_mode(sim, mode) # handle the logic of failure modes so we have high level simple controls
sim.assume_limit = False # True if we want to assume natural (not forced) convergence to terminal velocity
return sim, next_step
# half deprecated
def main_chute_steady(sim, next_step):
if sim.prnt: print('Main Chute Actively Dragging!\n')
# safety distance above the ground that the main chute should deploy at [m]
safety_dist= 304.8 # = 1000 ft
# time until ground impact from safety distance at main chute's terminal velocity
chute_t = abs(safety_dist / vt_c)
duration = chute_t
if not sim.fail_mode['mainchute']: sim.state = [C_d_c, A_c] # use main chute drag if it isn't failing
# if main chute fails, drag is only live if the drogue didn't fail. if drogue failed, we free fall
sim.phase['live'] = True if not sim.fail_mode['mainchute'] else not sim.fail_mode['drogue']
sim.phase['mainchute'] = True # either way we fall until impact
sim.phase['steadystate'] = True # jump straight to steady state
# kludge to account for jumping from one steady state to another
sim.v[-1][0] = sim.v_w[0]
sim.v[-1][1] = sim.v_w[1]
sim.v[-1][2] = vt_c
iterate(sim, next_step, duration)
### these are the important steps, cleanly presented
def static_line(sim, next_step):
duration = np.sqrt(abs(2 * sl / g))
iterate(sim, next_step, duration)
if sim.prnt: print('Static Line Taut!\n')
def drogue_deploy(sim, next_step): # and eNSR ring separation
duration = d_deploy_t
iterate(sim, next_step, duration)
if sim.prnt: print('Drogue Deployed!\n')
def drogue_line(sim, next_step):
duration = quartic(sim.v[-1])
iterate(sim, next_step, duration)
if sim.prnt: print('Drogue Line Taut!\n')
def drogue_action(sim, next_step):
if sim.prnt: print('Drogue Actively Dragging!\n')
duration = c_deploy_t
sim.phase['live'] = not sim.fail_mode['drogue'] # not a dud, unless it is
iterate(sim, next_step, duration)
def main_chute(sim, next_step):
if sim.prnt: print('Main Chute Actively Dragging!\n')
if not sim.fail_mode['mainchute']: sim.state = [C_d_c, A_c] # use main chute drag if it isn't failing
# if main chute fails, drag is only live if the drogue didn't fail. if drogue failed, we free fall
sim.phase['live'] = True if not sim.fail_mode['mainchute'] else not sim.fail_mode['drogue']
sim.phase['mainchute'] = True # either way we fall until impact
iterate(sim, next_step, 0) # time condition doesn't matter since determined by height
# run this sim! so pretty!
def trajectory(v_plane=18, v_side=0.0, deg_plane=90+45, v_wind=3.1, deg_wind=360-45, x_0_offset=[0, 0, 0], t_0=0,
mode='', dt=time_step, prnt=False):
sim, next_step = init(v_plane, v_side, deg_plane, v_wind, deg_wind, x_0_offset, t_0,
mode, dt, prnt)
static_line(sim, next_step)
drogue_deploy(sim, next_step)
if not sim.fail_mode['early']: # skip drogue and go straight to main chute if failure mode is early deployment
drogue_line(sim, next_step)
drogue_action(sim, next_step)
main_chute(sim, next_step)
return (np.array([sim.x[-1][0], sim.x[-1][1]]), sim)
# this block is so that our notebook won't automatically run a sim when it is imported ;)
if __name__ == '__main__' and not '__file__' in globals():
coords, results = trajectory(v_plane=15, v_wind=5, dt=0.005, prnt=True)
print('TOTAL HORIZONTAL DISTANCE FROM ORIGIN:', round(np.linalg.norm(results.x[-1][0:2]), 2),
'm', '=', round(np.linalg.norm(results.x[-1][0:2]) * 3.28084, 2), 'ft')
print('TOTAL VERTICAL DISTANCE DESCENDED:', round(drop_height - results.x[-1][2], 2),
'm', '=', round((drop_height - results.x[-1][2]) * 3.28084, 2), 'ft')
print('TOTAL TIME FOR DESCENT:', round(results.t[-1], 2), 's =', round(results.t[-1] / 60, 2), 'min')
print('Final coordinates: (x, y) = (', coords[0], ",", coords[1],')')
graph_sol(results, start=0, end=None)
Recall that this will be used by a simulator upstream. It will be run many times, so that simulator will stipulate a new time step that is much rougher. To verify that this doesn't significantly harm our model, we will run the above simulation which was run at dt=0.005 a second time at the Monte Carlo time step of dt=0.25 so we can compare the results and get a sense of how consistent the two numerical integrations are.
Then we will run a third time using the manual methods which should save some time without effecting the results.
In [7]:
if __name__ == '__main__' and not '__file__' in globals():
coords, results = trajectory(v_plane=15, v_wind=5, dt=0.05, prnt=True)
print('TOTAL HORIZONTAL DISTANCE FROM ORIGIN:', round(np.linalg.norm(results.x[-1][0:2]), 2),
'm', '=', round(np.linalg.norm(results.x[-1][0:2]) * 3.28084, 2), 'ft')
print('TOTAL VERTICAL DISTANCE DESCENDED:', round(drop_height - results.x[-1][2], 2),
'm', '=', round((drop_height - results.x[-1][2]) * 3.28084, 2), 'ft')
print('TOTAL TIME FOR DESCENT:', round(results.t[-1], 2), 's =', round(results.t[-1] / 60, 2), 'min')
print('Final coordinates: (x, y) = (', coords[0], ",", coords[1],')')
graph_sol(results, start=0, end=None)
Now we will investigate our failure modes. The cases are A (no main chute), B (no drogue), C (no drogue and early main deployment), and D (total dud). Note that a blank case will represent mission success. Also note that we are still ignoring drag from the body of the recovery system. I would prefer to account for it, but I have yet to estimate its coefficients. As such, there is no difference between case B and case D, since the main chute deploys after impact without a drogue to slow descent. Therefore, we will only check the success case, case A, case B, and case C in our Monte Carlo simulation.
I've checked the default wind/plane scenario and found that the variance between smooth and rough trajectories seems negligible for our failure modes as well. You are free to double check yourself, it only requires changing which trajectory function you run. I ran the smooth trajectory here because the graphs from trajectory_fast() are uninsightful.
In [8]:
# this block is so that our notebook won't automatically run a sim when it is imported ;)
if __name__ == '__main__' and not '__file__' in globals():
coords, results = trajectory(v_plane=15, v_wind=5, mode='D', dt=0.005, prnt=True)
print('TOTAL HORIZONTAL DISTANCE FROM ORIGIN:', round(np.linalg.norm(results.x[-1][0:2]), 2),
'm', '=', round(np.linalg.norm(results.x[-1][0:2]) * 3.28084, 2), 'ft')
print('TOTAL VERTICAL DISTANCE DESCENDED:', round(drop_height - results.x[-1][2], 2),
'm', '=', round((drop_height - results.x[-1][2]) * 3.28084, 2), 'ft')
print('TOTAL TIME FOR DESCENT:', round(results.t[-1], 2), 's =', round(results.t[-1] / 60, 2), 'min')
print('Final coordinates: (x, y) = (', coords[0], ",", coords[1],')')
graph_sol(results, start=0, end=None)
In [9]:
# this block is so that our notebook won't automatically run a sim when it is imported ;)
if __name__ == '__main__' and not '__file__' in globals():
coords, results = trajectory(v_plane=15, v_wind=5, mode='B', dt=0.005, prnt=True)
print('TOTAL HORIZONTAL DISTANCE FROM ORIGIN:', round(np.linalg.norm(results.x[-1][0:2]), 2),
'm', '=', round(np.linalg.norm(results.x[-1][0:2]) * 3.28084, 2), 'ft')
print('TOTAL VERTICAL DISTANCE DESCENDED:', round(drop_height - results.x[-1][2], 2),
'm', '=', round((drop_height - results.x[-1][2]) * 3.28084, 2), 'ft')
print('TOTAL TIME FOR DESCENT:', round(results.t[-1], 2), 's =', round(results.t[-1] / 60, 2), 'min')
print('Final coordinates: (x, y) = (', coords[0], ",", coords[1],')')
graph_sol(results, start=0, end=None)
In [10]:
# this block is so that our notebook won't automatically run a sim when it is imported ;)
if __name__ == '__main__' and not '__file__' in globals():
coords, results = trajectory(v_plane=15, v_wind=5, mode='A', dt=0.005, prnt=True)
print('TOTAL HORIZONTAL DISTANCE FROM ORIGIN:', round(np.linalg.norm(results.x[-1][0:2]), 2),
'm', '=', round(np.linalg.norm(results.x[-1][0:2]) * 3.28084, 2), 'ft')
print('TOTAL VERTICAL DISTANCE DESCENDED:', round(drop_height - results.x[-1][2], 2),
'm', '=', round((drop_height - results.x[-1][2]) * 3.28084, 2), 'ft')
print('TOTAL TIME FOR DESCENT:', round(results.t[-1], 2), 's =', round(results.t[-1] / 60, 2), 'min')
print('Final coordinates: (x, y) = (', coords[0], ",", coords[1],')')
graph_sol(results, start=0, end=None)
In [11]:
# this block is so that our notebook won't automatically run a sim when it is imported ;)
if __name__ == '__main__' and not '__file__' in globals():
coords, results = trajectory(v_plane=15, v_wind=5, mode='C', dt=0.005, prnt=True)
print('TOTAL HORIZONTAL DISTANCE FROM ORIGIN:', round(np.linalg.norm(results.x[-1][0:2]), 2),
'm', '=', round(np.linalg.norm(results.x[-1][0:2]) * 3.28084, 2), 'ft')
print('TOTAL VERTICAL DISTANCE DESCENDED:', round(drop_height - results.x[-1][2], 2),
'm', '=', round((drop_height - results.x[-1][2]) * 3.28084, 2), 'ft')
print('TOTAL TIME FOR DESCENT:', round(results.t[-1], 2), 's =', round(results.t[-1] / 60, 2), 'min')
print('Final coordinates: (x, y) = (', coords[0], ",", coords[1],')')
graph_sol(results, start=0, end=None)
In [ ]: