In [1]:
#openrocket_interface imports
from math import pi, log, sqrt, exp, cos, sin, radians
import os
from sys import platform as _platform
import xml.etree.ElementTree as ET # xml library
from zipfile import ZipFile
#trajectory imports
from types import SimpleNamespace
import numpy as np
import csv
import copy
#optimizer imports
from scipy.optimize import minimize, differential_evolution
import pylab
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import rc
global OF, rho_ipa # this is technically evil
np.set_printoptions(precision=3) # this line may be deprecated
In [ ]:
# this is hamilton's quaternion product
def product(a, b):
v = b[0] * a[1:] + a[0] * b[1:] + np.cross(a[1:], b[1:])
return np.array([a[0] * b[0] - np.dot(a[1:], b[1:]), v[0], v[1], v[2]])
# this is the inverse of a unit quat
def conjugate(q):
return np.array([q[0], -q[1], -q[2], -q[3]])
# this rotates a vector with a fixed frame
def sandwich(q, v):
return product(q, product(np.array([0, v[0], v[1], v[2]]), conjugate(q)))[1:]
# this rotates a frame with a fixed vector
def frame_rotation(q, v):
return sandwich(conjugate(q), v)
# this constrains a quat to S3 or vector to S2
def normalize(q):
norm = np.linalg.norm(q)
norm = norm if norm !=0 else 1
return q / norm
Gathering all input variables, design parameters, and actual constants in one place.
In [2]:
# Physics, chemistry, and materials
g_n = 9.80665 # kg.m/s^2 Standard gravity
R_univ = 8314.46261815324 # universal gas constant, J/ (K * kmol)
# Nitrogen characteristics
# https://github.com/psas/reaction-control/blob/master/pubs/AIAA%20RCS%20Manuscript_FINAL2.pdf
N2_temp = 293 # K, holding this constant is sketchy but easy
N2_MM = 28.01 # nitrogen molecular mass [g/mol]
N2_ke = 1.4 # nitrogen specific heat ratio
# combustion gas properties ke, Re, T_ch, determined from CEArun
# with chamber pressure=350 psi, fuel temp=419.15 K,
# lox temp=90 K, OF=1.3 for fuel = 64.8% IPA (2propanol) / 35.2% H20
OF = 1.3 # O/F ratio, this is somewhat arbitrary but CEA says its good.
p_ch = 2413166 # chamber pressure, PSI
T_ch = 3097.82 # chamber temperature, K
ke = 1.1251 # specific heat ratio, propellant (aka gammas)
MM = 23.196 # molar mass
def Material(name, rho, Sy=None):
return { 'name': name,
'rho': rho, # kg/m^3 Density
'Sy': Sy} # Pa Yield strength
# Materials
# https://www.aircraftspruce.com/catalog/cmpages/anh4120honeycomb01-01574.php?clickkey=5444217
# note: i've calculated ~395 kg/m^3 for uniform airframe density based on our density/thickness estimates.
# however, measurement of LV3.1 module says module is around 130 kg/m^3. Not sure why these don't agree.
Nomex = Material('Nomex', 48.06)
Cryogel = Material('Cryogel', 160)
Fiberglass = Material('Fiberglass', 1850)
Alum = Material('Aluminum', 2800.0, 0.270e9)
CFiber = Material('Carbon Fiber', 1550.0, 0.450e9)
LOX = Material('LOX', 1141.0) # kg/m^3 Density of LOX
IPA = Material('IPA/H20', 849.28) # kg/m^3 Density of 64.8% IPA / 35.2% H20
# Consider that there are two tanks, and we will want to divide total mass flow rate and propellant mass
# oxygen first, fuel second
def proportion(amount, OF):
stuff_o = amount * OF/(1 + OF)
stuff_f = amount * 1/(1 + OF)
return stuff_o, stuff_f
In [3]:
# Launch constants
launch_site_alt = 1401 # m, altitude of launch site above sea level
launch_tower = 18.288 # launch rail height in m
In [4]:
# System Definition
# Rocket parameters
airfrm_in_rad = 0.1524 # rocket inner radius, m
throttle_window = (100., 500.) # lower and upper bounds of drag force (N) for throttling
min_throttle = 1. # the internet says between 60 - 70% is doable without ruining our lives.
ballast = 0.0
# Fin Geometry
root = 0.7 # Root length, m
tip = 0.45 # tip length, m
sweep = np.radians(60) # sweep angle, degrees
span = 0.42 # fin span/height, m
# RCS parameters
# low priority, find constraints for these and optimize
rcs_control = False # whether sims have RCS enabled
rcs_mdot = 0.03784/3.15 # kg/s RCS nozzle mass flow rate
rcs_p_e = 101300 # Pa, RCS nozzle exit pressure
rcs_p_ch = 1.72*10**7 # Pa, N2 tank chamber pressure
# upper subsystem module dimensions, if these change you will have to make changes in structure.ipynb
airframe_thickness = 0.008255 # m, (0.325 in)
nose_l = 1.50 # m
ers_l = 6 * 0.0254 # m (converted from in)
rcs_l = 6 * 0.0254 # m (converted from in)
av_l = 24 * 0.0254 # m (converted from in)
n2_l = 24 * 0.0254 # m (converted from in)
pas_l = 6 * 0.0254 # m (converted from in)
# Engine system dimensions. 'm_' = mass, 'l_' = length
# when you have time, these masses should transfer out of structures.ipynb and go into openrocket_interface.ipynb
ullage = 1.1 # percentage of length added to a tank to account for not filling
#m_feed = 10.0 # kg
l_feed = 0.4572 # m, this is 18"
#m_ems = 1 # kg
l_ems = 0.1016 # m, this is 4"
#m_engine = 6.0 # kg
l_engine = 0.300 # m
In [5]:
# SIMULATION AND OPTIMIZATION PARAMETERS
delta = 10**(-3) # a guess at a good margin for design "convergence"
mu_0 = 0.0025 # barrier parameter, this value lets altitudes approach lower bound pretty quickly
rho_0 = 5 # penalty parameter, i'm still playing with this value.
dt = 0.04 # change starting time-step for trajectory simulation
iterations = 3 # number of escalating iteration sequences
# INITIAL DESIGN GUESS
# be sure that you start with a feasible design, otherwise the problem will be ill-conditioned
m_prop = 155.127 # propellant mass (kg)
mdot = 3.355 # Propellant mass flow rate (kg/s)
p_e = 66786.346 # Exit Pressure (Pa)
x0 = np.array([m_prop, mdot, p_e, ballast, root, tip, sweep, span]) # initial design vector
# OPTIMIZATION CONSTRAINTS
cons_mass = 270. # GLOW constraint, kg, somewhat arbitrary
cons_ls = 22. # min launch speed from 60' tower constraint, m/s
cons_TWR = 2. # TWR constraint
cons_stbl = 2. # minimum stability margin calibers throughout flight
cons_S_crit = 0.35 # Critical pressure ratio constraint
cons_accel = 15. # Max acceleration constraint, g's
cons_LD = 20. # L/D ratio constraint, slightly arbitrary
cons_alt = 110000. + launch_site_alt # Min altitude constraint, m (adjusted to overshoot)
cons_thrust = 6. # max ground-level thrust, kN
cons_ceiling = 150000. + launch_site_alt # base-11 maximum apogee requirement, km
cons_stblty = 2.0 # minimum in flight stability margin caliber
In [ ]:
rkt_prefix = "../rocket_farm/" # the rockets live on a cute little farm upstate where they frolic in fields
## Utility Functions
# unpack rocket template temporarily
def unzip():
with ZipFile('../LV4_canonical/template.ork') as myzip:
myzip.extract('rocket.ork')
# package our new rocket and remove temporary template
def zipit(index):
with ZipFile(rkt_prefix+'psas_rocket_'+index+'.ork', 'w') as myzip:
myzip.write('rocket.ork')
if 'linux' in _platform:
os.system('rm rocket.ork')
elif "darwin" in _platform:
os.system('rm rocket.ork')
elif "win" in _platform:
os.system('del rocket.ork')
# pulls ALL file references from given directory
def all_files(directory):
for path, dirs, files in os.walk(directory):
for f in sorted(files):
yield os.path.join(path, f)
# counts how many rockets are in our directory and then increments by 1
def get_index():
ork_files = [f for f in all_files(rkt_prefix)
if f.endswith('.ork')]
return len(ork_files) + 1