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