This notebook will encompass all calculations regarding the LV3 Recovery/eNSR Drop Test.
[http://www.usma.edu/math/Military%20Math%20Modeling/C5.pdf]
[http://www.the-rocketman.com/drogue-decent-rate.html]
[http://wind.willyweather.com/or/crook-county/prineville-reservoir.html]
In [1]:
import math
import sympy
from sympy import Symbol, solve
from scipy.integrate import odeint
from types import SimpleNamespace
import numpy as np
import matplotlib.pyplot as plt
sympy.init_printing()
%matplotlib inline
In [2]:
# P = speed of plane [m/s]
P = 38
# wind speed [m/s]
w = 2.59283
# wind bearing, measured from east to north [degrees]
theta_w_deg= 360-45
# aircraft bearing, measured from east to north [degrees]
theta_a_deg= 90+45
# safety distance above the ground that the main chute should deploy at [m]
mainSafetyDist= 304.8 # 1000 ft = 304.8 m
In [3]:
# worst case wind speed [m/s]
# w_worst= 12.86
# mass of payload [kg]
m = 28
# g = acceleration due to gravity [kg*m/s^2]
g = 9.81
# density of air [kg/m^3]
rho = 1.2
# terminal velocity of drogue [m/s]
vt_d= 18.5 # according to rocketman
# radius of drogue [m]
R_d= 0.762
# static line length [m]
sl = 2
# drogue line length [m]
dl= 50
# terminal velocity of main chute [m/s]
vt_m= 4.83108 # according to rocketman
# radius of main chute [m]
R_m= 5.4864
In [4]:
# wind speed in the direction of flight
theta_a= theta_a_deg*2*np.pi/360
theta_w= theta_w_deg*2*np.pi/360
wx= w*np.cos(theta_w-theta_a)
# cross-wind speed from left to right (pilot's perspective)
wz= -w*np.sin(theta_w-theta_a)
In [5]:
va = 0
## vertical distance gained
## since the static line is 2m, assuming it falls only in the vertical direction:
Lva = sl
# horizontal distance gained
## speed of plane times time to drop 2m static line
## 1/2*g*ta**2 = sl
ta = math.sqrt(2*sl/g)
Lha = P*ta
print('step a (from drop to static line disconnect):')
print('time to free fall fall 2 m = ', round(ta,4), ' s')
print('vertical length gained = ', round(Lva,4), ' m')
print('horizontal length gained = ', round(Lha,4), ' m')
deployment timer is a 2 sec timer
In [6]:
# vertical velocity at end of static line, beginning of timer
vb = va + (g*ta)
# since the deployment is controlled by a 2 sec timer:
tb = 2
# vertical length gained
Lvb = (vb*tb) + (0.5*g*(tb**2))
# horizontal length gained
Lhb = P*tb
print('step b (from static line disconnect to timer runout):')
print('vertical velocity at beginning of step b = ', round(vb,4), ' m/s')
print('vertical length gained = ', round(Lvb,4), ' m')
print('horizontal length gained = ', round(Lhb,4), ' m')
In [7]:
# velocity at time of ring separation, end of timer
vc = vb + g*tb
Lhc = 0
Lvc = 0
print('vertical velocity at ring separation = ', round(vc,4), ' m/s')
vertical length gained from step d
$Lvd = vc*td + 0.5*g*(td^2)$
horizontal length gained from step d
$Lhd = P*td$
calculate td by replacing x and y in the above equation
$dl^2 = (P*td)^2 + (vc*td + g*td^2)^2$
In [8]:
Ps, vds, gs, Lvds, Lhds, tds, vcs = sympy.symbols('Ps vds gs Lvds Lhds tds vcs')
Dparms= {Ps: P, gs: g, vcs: vc}
tdEqn= (Ps*tds)**2 + (vcs*tds + 0.5*gs*tds**2)**2 - dl**2
tdSolns= sympy.solve(tdEqn.subs(Dparms))
print('possible solutions:', tdSolns)
for soln in [complex(x) for x in tdSolns]:
if (soln.imag != 0) or (soln.real <= 0):
pass
else:
print(soln, 'seems fine')
td= soln.real
# now go back and calculate x and y
Lhd = P*td
Lvd = vc*td + g*(td**2)
# vertical velocity gained after the 50' drop
vd = vc + g*td
print()
print('time to pull out drogue:', round(td,4), 's')
print('horizontal distance gained = ', round(Lhd,4), 'm')
print('vertical distance gained = ', round(Lvd,4), 'm')
print('vertical velocity at instant line becomes taught = ', round(vd,4), 'm/s')
print('horizontal velocity: ', P, 'm/s')
Just start with Newton's 2nd law. The $-1/2\rho$ stuff is the drag force. It's negative because it opposes the motion. The biz with the $|\dot{\vec r}|\dot{\vec r}$ is to get a vector that has the magnitude of $r^2$ and the direction of $\vec r$.
$
m*\ddot{\vec r} = -1/2*\rho*A_d*Cd*|\dot{\vec r}|\dot{\vec r} +m*\vec g\\
$
Break it out into components. (This is where we see that it's an ugly coupled diffeq.)
$
m* \ddot r_x = -1/2*\rho*A_d*Cd*\sqrt{\dot r_x^2+\dot r_y^2}*\dot r_x \\
m*\ddot r_y = -1/2*\rho*A_d*Cd*\sqrt{\dot r_x**2+\dot r_y**2}*\dot r_y -m*g
$
In [9]:
# make a function that translates our equtions into odeint() format
def dragFunc(y, t0, p):
# map the positions and velocities to convenient names:
r_x= y[0]
r_y= y[1]
r_z= y[2]
rdot_x= y[3]
rdot_y= y[4]
rdot_z= y[5]
# calculate the accelerations:
rddot_x= 1/p.m*(-1/2*p.rho*p.A*p.Cd*np.sqrt((rdot_x-p.wx)**2+rdot_y**2+(rdot_z-p.wz)**2)*(rdot_x-p.wx))
rddot_y= 1/p.m*(-1/2*p.rho*p.A*p.Cd*np.sqrt((rdot_x-p.wx)**2+rdot_y**2+(rdot_z-p.wz)**2)*rdot_y -p.m*p.g)
rddot_z= 1/p.m*(-1/2*p.rho*p.A*p.Cd*np.sqrt((rdot_x-p.wx)**2+rdot_y**2+(rdot_z-p.wz)**2)*(rdot_z-p.wz))
# return the velocities and accelerations:
return([rdot_x, rdot_y, rdot_z, rddot_x, rddot_y, rddot_z])
In [10]:
D_d = m*g # drag force on drogue at terminal velocity [N]
A_d = math.pi*(R_d**2) # frontal area of drogue [m^2]
cd_d = (2*D_d)/(rho*A_d*vt_d**2) # drag coeff. of drogue []
# bundle up the parameters needed by dragFunc():
pd = SimpleNamespace()
pd.rho = rho
pd.A = A_d
pd.Cd = cd_d
pd.m = m
pd.g = g
pd.wx = wx
pd.wz = wz
# set the boundary conditions for the solver:
y0 = [0,0,0, P, -vd, 0]
t_step = 0.001
t_start = 0
t_final = 10
times_d = np.linspace(t_start, t_final, (t_final-t_start)/t_step)
# run the simulation:
soln_d = odeint(func= dragFunc, y0= y0, t= times_d, args= (pd,))
# find the time when it's okay to deploy the main chute:
# for i in range(0, len(soln)):
# if (soln_d_xddot[i] < 0.01*soln_d_xddot[0]) and (soln_d_yddot[i] < 0.01*soln_d_yddot[0]):
# print('At time', round(times_d[i],4), 'x and y acceleration are below 1% their original values.')
# tcr_d= times_d[i]
# break
# chop of the stuff after the critical time:
#soln= soln[range(0,i)]
#times= times[range(0,i)]
In [11]:
# break out the solutions into convenient names:
soln_d_x= [s[0] for s in soln_d]
soln_d_y= [s[1] for s in soln_d]
soln_d_z= [s[2] for s in soln_d]
soln_d_xdot= [s[3] for s in soln_d]
soln_d_ydot= [s[4] for s in soln_d]
soln_d_zdot= [s[5] for s in soln_d]
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_x, soln_d_y)
plt.axis('equal')
plt.xlabel('horizontal range (m)')
plt.ylabel('vertical range (m)')
plt.figure(2)
plt.plot(times_d, soln_d_xdot)
plt.xlabel('time (s)')
plt.ylabel('horizontal velocity (m/s)')
plt.figure(3)
plt.plot(times_d, soln_d_ydot)
plt.xlabel('time (s)')
plt.ylabel('vertical velocity (m/s)')
plt.figure(4)
plt.plot(times_d[range(0, len(soln_d_xddot))], soln_d_xddot)
plt.xlabel('time (s)')
plt.ylabel('horizontal acceleration (m/s^2)')
plt.figure(5)
plt.plot(times_d[range(0, len(soln_d_yddot))], soln_d_yddot)
plt.xlabel('time (s)')
plt.ylabel('vertical acceleration (m/s^2)')
Out[11]:
In [12]:
Lhe = soln_d_x[-1]
Lve = -soln_d_y[-1]
Lle = soln_d_z[-1]
te = times_d[-1]
print('horizontal distance travelled in step e:', Lhe)
print('vertical distance travelled in step e:', Lve)
print('lateral distance travelled in step e:', Lle)
print('time taken in step e:', te)
In [13]:
# x-direction calculations
##########################
# from usma:
# mx" + cd*x' = cd*w
####### need python help here #######
# ugh, I have to go learn how to use scipy... 1 sec -- Joe
# mx" + cd*x
## homogeneous equation mx" + rho*x' = 0
## characteristic equation for the homogeneous differential equation is:
## mr^2 + rho*r = 0
## where the roots are:
## r1 = 0, r2 = -(rho/m)
## complementary solution:
## xc = C1*e^0 + C2* e^(-(rho*t/m))
## non-homogeneous equation mx" + rho*x' = rho*w
## complete solution x = C1 + C2*e^(-(rho*t/m)) + wt
## solving for C1 and C2 using results from step d as initial conditions
## except time = 0 since we are making calculations just for this step
## i.e. x(0) = x_curr_tot and vx(0) = P
## therefore C1 = and C2 =
# x_0 = Lha + Lhb + Lhc + Lhd
# t = 0
# vx_0 = P
# C1 = Symbol('C1')
# C2 = Symbol('C2')
# C_1 = solve(C1 + C2*math.exp(-(rho*t/m)) + w*t - x_0, C1)
# C_1
# print(C_1)
# C_2 = solve(C2*(-(rho/m)) + w - vx_0, C2)
# print(C_2)
# ## NEEEED HELLLPPP should be using piecewise to solve this
# ## copying C_1 output from just above with the C_2 value
# calc_C1 = 147.560492558936 + 586.6
# print(calc_C1)
#
# ## therefore the complete solution is:
# ## x = 734.1605 - 586.6*exp(-(rho*t/m)) + w*t
#
# ## if the drogue falls for 3 seconds, then
# t = 3
# Lhe = 734.1605 - 586.6*math.exp(-(rho*t/m)) + w*t
#
# print('horizontal distance gained = ', round(Lhe,4), 'm')
# print(' ')
#
# # y-direction calculations
# ##########################
#
# ## from usma
# ## characteristic equation:
# ## m*r^2 + rho*r = 0
# ## where the roots are r = 0 and r = (-b/m)
#
# ## complete solution:
# ## y = C1 + C2*exp(-(rho*t)/m)
# ## solving for C1 and C2 using results from step d as initial conditions
# ## except time = 0 since we are making calculations just for this step
#
# y_0 = Lva + Lvb + Lvc + Lvd
# print('y_0 = ', y_0)
# vy_0 = vd
# print('vy_0 = ',vy_0)
# t_0 = 0
# C1 = Symbol('C1')
# C2 = Symbol('C2')
# ## NEEEED HELLLPPP should be using piecewise to solve this
# # C1 equation
# C_1 = solve(C1 + C2*math.exp(-(rho*t_0/m)) - y_0, C1)
# print('C1 equation: ', C_1)
# # C2 equation/value
# C_2 = solve(C2*(-(rho/m)*math.exp(-(rho*t_0/m))) - vy_0, C2)
# print('C2 = ', C_2)
# ## copying C_1 output from just above with the C_2 value
# calc_C1 = 793.253769802079 + 62.2619406518579 #62.2619406518579 + (0.879350749407306*793.253769802079)
# print('C1 = ', calc_C1)
#
# # NEED HELP: need to make C_2 a number (int, float)
# ## if the drogue falls for 3 seconds, then
# t = 3
# Lve = calc_C1 + (-793.253769802079*math.exp(-(rho/m)*t))
#
# print('vertical distance gained = ', Lve, 'm')
#
# ## Maayybbbeee
#
# vert_length = v*t
# print(vert_length)
If you want to justify to yourself that the main chute hits terminal velocity [almost] instantly, you can mess with the inputs for the numerical solution in step e.
In [14]:
Lvf= mainSafetyDist
# step f time = vertical distance / main terminal velocity
tf= Lvf/vt_m
# horizontal distance= wind speed * step f time
Lhf= wx*tf
Llf= wz*tf
In [15]:
print('horizontal distance travelled in step f:', Lhf, 'm')
print('vertical distance travelled in step f:', Lvf, 'm')
print('time taken in step f:', tf, 's')
In [16]:
# TOTAL HORIZONTAL DISTANCE TRAVELED
X_TOT = Lha + Lhb + Lhc + Lhd + Lhe + Lhf
X_TOT_ft = X_TOT*3.28084
print('TOTAL HORIZONTAL DISTANCE TRAVELED = ', round(X_TOT,2), 'm ', ' = ', round(X_TOT_ft,2), 'ft')
# TOTAL VERTICAL DISTANCE DESCENDED
Y_TOT = Lva + Lvb + Lvc + Lvd + Lve + Lvf
Y_TOT_ft = Y_TOT*3.28084
print('TOTAL VERTICAL DISTANCE DESCENDED = ', round(Y_TOT,2), 'm ', ' = ', round(Y_TOT_ft,2), 'ft')
# TOTAL TIME FOR DESCENT
T_TOT = ta + tb + td + te + tf
# in minutes
t_tot_min = T_TOT/60
print('TOTAL TIME FOR DESCENT', round(T_TOT,2), 's = ', round(t_tot_min,2), 'min')
In [17]:
delta_xs= np.array([0, Lha, Lhb, Lhc, Lhd, Lhe, Lhf])
delta_ys= -np.array([0, Lva, Lvb, Lvc, Lvd, Lve, Lvf])
delta_zs= np.array([0, 0, 0, 0, 0, Lle, Llf])
xs= np.cumsum(delta_xs)
ys= np.cumsum(delta_ys)
zs= np.cumsum(delta_zs)
plt.close('all')
plt.figure(1)
plt.plot(xs,ys)
_= plt.axis('equal')
plt.grid()
plt.title('down-range trajectory')
plt.xlabel('down-range distance from drop (m)')
plt.ylabel('altitude relative to drop (m)')
plt.figure(2)
plt.plot(zs, ys)
_= plt.axis('equal')
plt.grid()
plt.title('lateral trajectory')
plt.xlabel('lateral (left to right) distance from drop (m)')
plt.ylabel('altitude relative to drop (m)')
print('xs:', xs)
print('ys:', ys)
print('zs:', zs)
print('note that Y is up and Z is to the right of the aircraft... because I don\'t want to change my code.')
In [18]:
Es= xs*np.cos(theta_a) +zs*np.sin(theta_a)
Ns= xs*np.sin(theta_a) -zs*np.cos(theta_a)
plt.figure(2)
plt.plot(Es,ys)
_= plt.axis('equal')
plt.grid()
plt.title('east trajectory')
plt.xlabel('eastern distance from drop (m)')
plt.ylabel('altitude relative to drop (m)')
plt.figure(3)
plt.plot(Ns, ys)
_= plt.axis('equal')
plt.grid()
plt.title('north trajectory')
plt.xlabel('northern distance from drop (m)')
plt.ylabel('altitude relative to drop (m)')
print('Es:', Es)
print('ys:', ys)
print('Ns:', Ns)
In [ ]:
In [ ]: