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 System_Definition.ipynb
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]:
# 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)
# Total length of engine system
def system_length(l_o, l_f):
return sum([l_o, l_f, l_feed, l_ems, l_engine, 0.152]) # plus pass through
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 [3]:
from IPython.display import SVG
SVG(filename='../archive/motor-system.svg')
Out[3]:
In [4]:
'''# 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(l_o, l_f, m_tank_o, m_tank_f, system_mass):
# fuel tank cm is in the center of the tank
cm_tank_f = l_f / 2.0
cm_passthru= l_f + ers_l/2
# next tank down (lox)
cm_tank_o = l_f + ers_l + (l_o/2.0)
# ems
cm_ems = l_f + ers_l + l_o + (l_ems/2.0)
# now feedsystem
cm_feed = l_f + ers_l + l_o + l_ems + (l_feed/2.)
# finally the engine
cm_engine = l_f + ers_l + l_o + l_feed + l_ems + (l_engine/2.0)
# get center of mass
dry_cm = sum([cm_tank_f*m_tank_f, cm_tank_o*m_tank_o, cm_feed*m_feed, cm_ems*m_ems, cm_engine*m_engine])
dry_cm /= system_mass
return dry_cm
# Total center of mass (including propellants) at a specific time
def c_of_m(dry_cm, dry_mass, r, l_o, l_f, prop_mass):
# propellant masses
m_o, m_f = proportion(prop_mass, OF)
#accounts for gravity pushing each propellent to bottom of its tank
cm_f_l = l_f - tank_length(m_f, IPA['rho'], r)/2.0
cm_o_l = l_f + ers_l + l_o - tank_length(m_o, LOX['rho'], r)/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 [5]:
# create a rocket file for our engine's dimensions and characteristics
def update_body(index, out_r, l_o, l_f, ballast, fin):
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(out_r)
for kid in child.iterfind("*[name='Tip Ballast']"):
kid.find('mass').text = str(ballast) # set ballast mass
for kid in child.iterfind("*[name='Fuel Tank']"):
kid.find('length').text = str(l_f) # set fuel tank length
for kid in child.iterfind("*[name='LOX Tank']"):
kid.find('length').text = str(l_o) # set lox tank length
for kid in child.iterfind("*[name='Trapezoidal fin set']"):
kid.find('rootchord').text = str(fin.root) # set fin geometry
kid.find('tipchord').text = str(fin.tip)
kid.find('sweeplength').text = str(fin.sweep_length)
kid.find('height').text = str(fin.height)
kid.find('thickness').text = str(fin.thickness)
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, Thrust,
inr_radius, thickness,
l_o, l_f, m_tank_o, m_tank_f,
Burn_time, Isp, sys_mass, length, CoM,
ballast, fin):
index = str(get_index())
radius = inr_radius + thickness
#length = system_length(l_o, l_f)
eng_sys_dry_mass = sys_mass
#eng_sys_dry_mass = sum([m_tank_o, m_tank_f, m_feed, m_ems, m_engine])
#dry_cm = dry_c_of_m(l_o, l_f, m_tank_o, m_tank_f, eng_sys_dry_mass)
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': radius*2*1000,
'length': length*1000,
'total_mass': (prop_mass[0] + eng_sys_dry_mass)*1000,
'M_prop': prop_mass[0]*1000,
'a_thrust': average,
'p_thrust': peak,
'burn_time': Burn_time,
'm_frac': eng_sys_dry_mass/(eng_sys_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) # 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': (eng_sys_dry_mass + prop_mass[i]) * 1000,
'cg': (length - CoM[i]) * 1000
#'cg': c_of_m(dry_cm, eng_sys_dry_mass, inr_radius, 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, radius, l_o, l_f, ballast, fin) # make a rocket to correspond with our new engine
print("Reopen OpenRocket to run simulation.")
In [ ]: