This is adapted from a script created by natronics, 7deeptide, and jalouke to create engine files for OpenRocket. Its original purpose was to determine some crucial properties of the engine system based on initial assumptions and specifications a few years before the multidisciplinary design optimization for LV4 began.
It has since been reformed to take its parameters from the optimization in order to utilize OpenRocket as a downstream verifier of potential designs, since the optimizer's internal model of aerodynamics is particularly low-fidelity. In addition, since this script is imported by our trajectory model and optimizer, we may define certain functions, variables, and parameters in one location. This eliminates redundancy and makes our lives easier if we would like to change any of our parameters based on improved assumptions or information.
In [1]:
%run ./inputs.ipynb # load all our inputs, parameters, and constants
rkt_prefix = "rocket_farm/" # the rockets live on a cute little farm upstate where they frolic in fields
### PARAMETERS END HERE ###
## Utility Functions
# unpack rocket template temporarily
def unzip():
with ZipFile('psas_rocket.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
I will explain the proportion function algebraically.
'A' is the total amount, 'O' is the proportion that is LOX, 'F' is the proportion that is IPA/H2O. $$A = O + F$$ $$\frac{O + F}{F + O} = \frac{A}{F + O} = 1$$ $$\frac{O}{F + O} + \frac{F}{F + O} = 1 \space\space\space(*)$$ $$\frac{F}{F}\frac{\frac{O}{F}}{1 + \frac{O}{F}} + \frac{F}{F}\frac{1}{1 + \frac{O}{F}} = 1$$ $$\frac{\frac{O}{F}}{1 + \frac{O}{F}} + \frac{1}{1 + \frac{O}{F}} = 1$$ $$O = A\frac{\frac{O}{F}}{1 + \frac{O}{F}}$$ $$F = A\frac{1}{1 + \frac{O}{F}}$$ These last two identities should be clear to see from the starred equation, and the preceding lines are to demystify all the fractions.
The upstream optimization doesn't account for the distinction between the radius of the airframe and of the propellant tank, so we do that here.
Because the radius of a tank changes while the mass and density of its contents remain constant, we must derive a new length so that its volume is unchanged. $$Area = \pi r^2$$ $$Volume = \frac{m}{\rho}$$ $$Length = \frac{Volume}{Area}$$
The upstream optimization does not account for there being two separate tanks with different requirements, so we do that here. After proportioning the propellant masses and determining the tank radii, we determine the length of each tank based on the mass and density of its contents. For a layer of realism, we add a percentage to the tank lengths to account for our inability to perfectly pack a tank without wasted space.
To be clear, tank_thickness is engineering black magic and beyond my understanding. I have changed it to account for a different tank pressure than the original MDO assumed. Nowhere do we account for the difference in volume that the inner radius entails, but this is hopefully a small change that is absorbed by the assumptions regarding ullage and factors of safety. It would be relatively trivial to account for that difference explicitly in split_tanks if we did some elementary arithmetic. If someone asks politely, I'll do it. tank_mass is more engineering black magic but straight-forward; its just a uniformly dense thicc cylinder with lids, and the estimate is doubled for good pessimism.
Georges demanded I account for the mass of bulkheads in the gaps between subsystems, so we made up a rough estimate of how heavy they might be based off our guesses about their geometry.
Be aware that split_tanks, system_length, and system_mass are used heavily by the upstream optimization and trajectory simulation where the fine detail is necessary. The airframe diameter and original total tank length are pretty much just used to determine the volume and thus total propellant mass.
In [2]:
# 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):
stuff_o = amount * OF/(1 + OF)
stuff_f = amount * 1/(1 + OF)
return stuff_o, stuff_f
# outer Radius of the tanks offset from airframe inner, converts in to m
def tank_r(total_dia):
return (total_dia/2 * 0.0254) - airframe_offset - insulation_layers * insulation_thickness
# inner Radius of airframe, m. Inverse of above function. depends on outer radius of tank
def body_r(eng_r):
return eng_r + airframe_offset + insulation_layers * insulation_thickness
# length of tank based on its inner radius, and mass and density of contents
def tank_length(m, rho, r):
return m / (rho * pi * r**2)
# tank thickness, depends on outer radius
# minimum tank thickness calculations defunct, since manufacturing constraints call for a higher minimum
# note i'm returning a value from inputs so that none of the code depending on this changes
def tank_thickness(tank, r):
#P_i = 689475.7 # Tank pressure in Pa (~100 PSI), assuming pressurized by N2 and ramped up later by EFS
#design_stress = tank['Sy']/factor_of_safety
#radius_i = sqrt(design_stress * (r**2) / ((2*P_i) + design_stress)) # inner radius
#return r - radius_i
return min_tank_thickness
# turn total propellant mass and inner airframe diameter into two tanks of propellants
def split_tanks(prop_mass, total_dia):
m_o, m_f = proportion(prop_mass)
r = tank_r(total_dia) # tank outer radius
thick_o = tank_thickness(Al, r)
thick_f = tank_thickness(Al, r)
l_o = tank_length(m_o, rho_lox, r - thick_o)
l_f = tank_length(m_f, rho_ipa, r - thick_f)
l_o *= ullage # add 10% for ullage
l_f *= ullage # add 10% for ullage
return r, l_o, l_f
# Tank Mass, depends on outer radius
def tank_mass(l, tank, r, insulated=False):
thick = tank_thickness(tank, r)
s_area = 2*pi*r*(l + r) # outer surface area of tank
extra = 0 if not insulated else insulation_layers * insulation_thickness * rho_cryo
return mass_fudger * s_area * (isogrid_reduction * thick * tank['rho'] + extra)
# bulkhead mass, depends on airframe inner radius
def bulkhead(r):
thicc = 0.00635 # m (0.25")
perforation = 0.4 # percentage of holiness
return perforation * Al['rho'] * thicc * pi * r**2
## System level functions
# Returns mass of each tank, depends on outer radius
def tank_builder(r, l_o, l_f):
return tank_mass(l_o, Al, r, insulated=False), tank_mass(l_f, Al, r)
# Total length of engine system
def system_length(l_o, l_f):
return sum([l_o, l_f, l_feed, l_ems, l_engine, 4*gaps])
# total dry mass of engine system, depends on outer radius
def system_mass(r, l_o, l_f):
bulkheads = bulkhead(body_r(r))
t1, t2 = tank_builder(r, l_o, l_f)
return sum([m_engine, m_ems, m_feed, t1, t2, 8*bulkheads])
This is a nice picture from some previous work that gives an idea of what the next code block accounts for.
Basically, we will eventually hand OpenRocket a black box for an engine system that will account for its own center of mass as well as changes in that. I don't believe these equations require explanation, they're just ugly routine computations, but if asked I will provide further detail. There may be a prettier way of writing this code that involves less eye-bleeding, but I'm fairly sure it's correct.
In [1]:
from IPython.display import SVG
SVG(filename='archive/motor-system.svg')
Out[1]:
In [ ]:
# For openrocket thrust Curve, since the mass and cm change over time.
# dry center of mass of engine system, note these are positions not lengths
def dry_c_of_m(r, l_o, l_f):
m_tank_o, m_tank_f = tank_builder(r, l_o, l_f)
bulkheads = bulkhead(body_r(r))
# fuel tank cm is in the center of the tank
cm_tank_f = l_f / 2.0
# including gaps since they have bulkheads now
cm_gap1 = l_f + gaps/2.
# next tank down (lox) has a gap.
cm_tank_o = l_f + gaps + dist_after_f + (l_o/2.0)
# next gap
cm_gap2 = l_f + gaps + dist_after_f + l_o + gaps/2.
# now feedsystem
cm_feed = l_f + gaps + dist_after_f + l_o + gaps + dist_after_o + (l_feed/2.0)
#next gap
cm_gap3 = l_f + gaps + dist_after_f + l_o + gaps + dist_after_o + l_feed + gaps/2.
#ems
cm_ems = l_f + gaps + dist_after_f + l_o + gaps + dist_after_o + l_feed + gaps + (l_ems/2.)
#last gap
cm_gap4 = l_f + gaps + dist_after_f + l_o + gaps + dist_after_o + l_feed + gaps + l_ems + gaps/2.
# finally the engine
cm_engine = l_f + gaps + dist_after_f + l_o + gaps + dist_after_o + \
l_feed + gaps + l_ems + gaps + (l_engine/2.0)
# sum cm
dry_cm = sum([cm_tank_f*m_tank_f, cm_tank_o*m_tank_o, cm_feed*m_feed, cm_engine*m_engine, \
cm_gap1*2*bulkheads, cm_gap2*2*bulkheads, cm_gap3*2*bulkheads, cm_gap4*2*bulkheads])
dry_cm /= system_mass(r, l_o, l_f)
return dry_cm
# Total center of mass (including propellants) at a specific time
def c_of_m(r, l_o, l_f, prop_mass):
thick_o = tank_thickness(Al, r)
thick_f = tank_thickness(Al, r)
# propellant masses
m_o, m_f = proportion(prop_mass)
# dry stuff
dry_cm = dry_c_of_m(r, l_o, l_f)
dry_mass = system_mass(r, l_o, l_f)
#accounts for gravity pushing each propellent to bottom of its tank
cm_f_l = l_f - tank_length(m_f, rho_ipa, r - thick_f)/2.0
cm_o_l = l_f + gaps + dist_after_f + l_o - tank_length(m_o, rho_lox, r - thick_o)/2.0
weighted_cm_prop = m_f*cm_f_l + m_o*cm_o_l
return (dry_cm*dry_mass + weighted_cm_prop) / (dry_mass + m_f + m_o)
Our optimization will spit out a list of relevant details, which is passed to this function so we can save everything potentially useful in one nice way.
Our rocket template is an xml file set to use one node to determine the diameter of each subsystem so we just have to climb that tree to change out the airframe radius and the two tank lengths and viola! Note, this effects just the rocket file, it is the next function that creates the engine file.
There is probably a nice way to use the imported XML library instead of the code from the previous MDO, but I don't quite care enough since it's perfectly functional. This takes its parameters straight from the trajectory simulation of an optimized rocket. It doesn't do anything particularly exciting, just creates a thrust curve and calculates some odds and ends that OpenRocket likes. It's unfortunate that thrust in OpenRocket doesn't account for altitude, but the discrepancy in final apogee is usually no more than 2 or 3 km (our trajectory simulation seems to be conservative).
Be aware that I have not yet tested this in Windows or MACOSX, or outside of Ubuntu 18 for that matter.
In [ ]:
# Takes rocket properties from optimized trajectory and creates a list of all relevant properties
def print_characteristics(mdot, prop_mass, r, l_o, l_f, index, res_text):
res_text.append("\n")
res_text.append("\nENGINE SYSTEM DETAILS")
res_text.append("\n-----------------------------")
# Mass flow for each propllent
mdot_o, mdot_f = proportion(mdot)
res_text.append("\nOx flow: . . . . . . . . . . %7.3f kg/s" % mdot_o)
res_text.append("\nFuel flow: %7.3f kg/s" % mdot_f)
# Propellent Mass for each propllent
mprop_o, mprop_f = proportion(prop_mass)
res_text.append("\nOx mass: . . . . . . . . . . . %5.3f kg" % mprop_o)
res_text.append("\nFuel mass: %5.3f kg" % mprop_f)
# dimensions of each tank
res_text.append("\nTank outer diameters: . . . . . . . %7.3f m" % (2*r))
res_text.append("\nOx tank length + ullage: %7.3f m" % l_o)
res_text.append("\nFuel tank length + ullage: %7.3f m" % l_f)
# Tank thickness for each tank (mm)
thickness_o = tank_thickness(Al, r)
thickness_f = tank_thickness(Al, r)
res_text.append("\nOx tank thickness: %5.3f mm" % (thickness_o*1000))
res_text.append("\nFuel tank thickness: %5.3f mm" % (thickness_f*1000))
# Mass of each tank
m_tank_o = tank_mass(l_o, Al, r, insulated=True)
m_tank_f = tank_mass(l_f, Al, r)
res_text.append("\nOx tank mass: . . . . . . . . %5.3f kg" % m_tank_o)
res_text.append("\nFuel tank mass: %5.3f kg" % m_tank_f)
# create a file with all this info in it
with open(rkt_prefix+'psas_rocket_'+index+'_traj.txt', 'w') as traj:
for line in res_text:
traj.write(line)
# create a rocket file for our engine's dimensions and characteristics
def update_body(index, eng_r, l_o, l_f):
unzip() # unpack template
with open('rocket.ork', 'rb') as xml_file:
tree = ET.parse(xml_file)
root = tree.getroot()
for child in root.iter():
#set radius, this node propagates itself in openrocket
for kid in child.iterfind('aftradius'):
kid.text = str(body_r(eng_r) + airframe_thickness)
for kid in child.iterfind("*[name='Fuel Tank']"):
kid.find('length').text = str(gaps + l_f) # set fuel tank length
for kid in child.iterfind("*[name='LOX Tank']"):
kid.find('length').text = str(gaps + l_o) # set lox tank length
tree.write('rocket.ork')
zipit(index) # repack template
print("New rocket generated from template!")
# when a father engine and a mother engine love each other very much...
# they make an openrocket engine
# that is approximately equivalent to the trajectory profile from the optimimization
def make_engine(mdot, prop_mass, total_dia, Thrust, throttle, Burn_time, Isp, res_text):
index = str(get_index())
r, l_o, l_f = split_tanks(prop_mass[0], total_dia)
print_characteristics(mdot, prop_mass[0], r, l_o, l_f, index, res_text)
length = system_length(l_o, l_f) + dist_after_f + dist_after_o # in case we want to shift subsystems around
dry_mass = system_mass(r, l_o, l_f)
n = len(Thrust)
peak = max(Thrust)
average = float(sum(Thrust) / n)
file_head = """<engine-database>
<engine-list>
<engine mfg="PSAS" code="{code}" Type="Liquid" dia="{diameter}" len="{length}"
initWt="{total_mass}" propWt="{M_prop}" delays="0" auto-calc-mass="0" auto-calc-cg="0"
avgThrust="{a_thrust}" peakThrust="{p_thrust}" throatDia="0." exitDia="0." Itot="{impulse}"
burn-time="{burn_time}" massFrac="{m_frac}" Isp="{Isp}" tDiv="10" tStep="-1." tFix="1"
FDiv="10" FStep="-1." FFix="1" mDiv="10" mStep="-1." mFix="1" cgDiv="10"
cgStep="-1." cgFix="1">
<comments>Optimized engine</comments>
<data>
""".format(**{'code': 'PSAS'+index,
'diameter': r*2*1000,
'length': length*1000,
'total_mass': (prop_mass[0] + dry_mass)*1000,
'M_prop': prop_mass[0]*1000,
'a_thrust': average,
'p_thrust': peak,
'burn_time': Burn_time,
'm_frac': dry_mass/(dry_mass + prop_mass[0]),
'impulse': average*Burn_time,
'Isp': Isp,
})
data = [] # this is going to be our thrust curve!
# may be slightly inaccurate if/when altitudes vary between trajectory.py and openrocket
resolution = float(Burn_time/(n-1)) # sec per step
for i in range(n):
data.append(' <eng-data t="{t}" f="{thrust}" m="{mass}" cg="{cg}"/>\n'.format(**{
't': i * resolution, # sec
'thrust': Thrust[i],
'mass': (dry_mass + prop_mass[i]) * 1000,
'cg': c_of_m(r, l_o, l_f, prop_mass[i]) * 1000,
}))
file_tail = """
</data>
</engine>
</engine-list>
</engine-database>"""
# we're gonna put the engine file in the default location
prefix = "./"
if 'linux' in _platform:
home = os.path.expanduser("~")
prefix = os.path.join(home, '.openrocket/ThrustCurves/')
elif "darwin" in _platform:
home = os.path.expanduser("~")
prefix = os.path.join(home, 'Library/Application Support/OpenRocket/ThrustCurves/')
elif "win" in _platform:
home = os.getenv("APPDATA")
prefix = os.path.join(home, "OpenRocket/ThrustCurves/")
if not os.path.isdir(prefix):
os.mkdirs(prefix) # in case openrocket not installed to prevent crash
# now write the file, great job!
with open(os.path.join(prefix, 'psas_motor_'+index+'.rse'), 'w') as eng:
eng.write(file_head)
for d in data:
eng.write(d)
eng.write(file_tail)
update_body(index, r, l_o, l_f) # make a rocket to correspond with our new engine
print("Reopen OpenRocket to run simulation.")
In [ ]: