At first blush, this will appear to be much more complicated than the trajectory simulation of the original MDO. That is because it is. However, that is just the trade-off we pay for an increase in modularity and flexibility, as well as (I hope) conceptual clarity.
Because we are doing systems engineering, I have imposed a hierarchical organization by making use of sequencing and Python classes.
In [1]:
# A numerical integration for rocket trajectories
%run "./openrocket_interface.ipynb"
# helper function for unit conversions
def unit_conv(p_ch, dia, p_e):
p_ch = p_ch*6894.76 # Chamber pressure, convert psi to Pa
# LV4 design variables
dia = dia*0.0254 # Convert in. to m
p_e = p_e*1000 # Exit pressure, convert kPa to Pa
return p_ch, dia, p_e
Copying (with minor changes) from the previous MDO,
"The rocket thrust force, given a huge number of simplifying assumptions, is a function of the mass flow rate design variable, the chamber pressure constant, and the ambient pressure state variable. It is given by $$F_t = \dot{m}V_e + A_e(p_e - p_a) $$
where $V_e$ is the rocket exhaust velocity, $A_e$ is the rocket nozzle exit area. The exhaust velocity, again given a large number of simplifying assumptions, is given by
$$V_e = \sqrt{\frac{2kRT_{ch}}{k - 1}\left(1 - \left(\frac{p_e}{p_{ch}}\right)^{((k - 1)/k)}\right)}$$where $T_{ch}$ is the chamber temperature, which along with the combustion chamber gas ratio of specific heats $k$, and gas constant $R$ are functions of the mixture ratio of the propellants, and the chamber pressure $p_{ch}$ determined using the NASA Chemical Equilibrium Analysis (CEA) tool and are design constants. It should be noted that the $p_e$ is a design variable tied to the expansion ratio of the rocket engine, but that $p_a$ is constantly changing with altitude. The thrust function is maximized when $p_e = p_a$, however this occurs at one, and only one, altitude. Thus the design expansion ratio selection directly benefits from trajectory optimization."
To make our function "clearer", I'll typeset each equation from below. $$p_t = p_{ch} \times (1 + \frac{k - 1}{2})^{\frac{-k}{k - 1}}$$ $$T_t = T_{ch} \times \frac{1}{1 + \frac{k - 1}{2}}$$ $$A_t = \frac{\dot{m}}{p_t} \times \sqrt{R \times \frac{T_t}{k}}$$ $$A_e = \frac{A_t \times (\frac{2}{k + 1})^{\frac{1}{k - 1}} \times (\frac{p_{ch}}{p_e})^{\frac{1}{k}}}{\sqrt{\frac{k + 1}{k - 1} \times (1 - (\frac{p_e}{p_{ch}})^{\frac{k - 1}{k}})}}$$ To be quite frank, I have not closely examined the consistency of this set of equations, but it apparently worked for the previous MDO, so I've been saving the hand-calc derivations for a rainy day.
In [2]:
# this class is for conceptual clarity, to keep rocket internals organized
class EngineSys:
# ***
# let the model keep track of its own properties
def __init__(self, mdot, p_e, p_ch, T_ch, ke, Re, lower, upper, minimum):
self.mdot = mdot
self.p_e = p_e
self.p_ch = p_ch
self.T_ch = T_ch
self.ke = ke
self.Re = Re
self.throttle = [1.]
self.minimum = minimum
self.lower = lower
self.upper = upper
self.slope = (minimum - 1.) / (upper - lower)
self.bias = 1 - self.slope * lower
# Throat pressure [Pa]
self.p_t = p_ch * (1 + (ke - 1)/2)**(-ke /(ke - 1))
# Throat temperature [K]
self.T_t = T_ch*(1 / (1 + (ke - 1)/2))
# Throat area [m^2]
self.A_t = (mdot / self.p_t) * sqrt(Re * self.T_t / ke)
# Expansion ratio
self.ex = (2 / (ke + 1))**(1 / (ke - 1)) * (p_ch / p_e)**(1/ke) / sqrt(
(ke + 1) / (ke - 1) * (1 - (p_e / p_ch)**(((ke - 1) / ke))))
# Exit area [m^2]
self.A_e = self.ex * self.A_t
# Lookup table of divergence angles, assuming 80% bell length
alpha_t = [14, 11, 10, 9]
# Lookup table of expansion ratios from alpha_t
ex_t = [5, 10, 15, 20]
alpha= np.interp(self.ex, ex_t, alpha_t)
# Thrust cosine loss correction,
lam = 0.5 * (1 + cos(alpha * pi/180))
# even in extreme cases this is definitely not an O(1) effect (that comment is from previous MDO)
# Exhaust velocity [m/s]
self.Ve = lam * sqrt(2*ke / (ke - 1) * Re * T_ch *
(1 - (p_e/p_ch)**((ke - 1)/ke))) * loss_factor
# linear throttle for simplicity
def throttle_engine(self, D):
if can_throttle == False:
return 1.
elif D <= self.lower:
return 1.
elif D >= self.upper:
return self.minimum
else:
return D * self.slope + self.bias
# calculates engine thrust at a height
def thrust(self, p_a, thrtl):
# Thrust force [N]
return thrtl*self.mdot*self.Ve + (self.p_e - p_a)*self.A_e
Here we can see that when a rocket is initialized, the design variable L only determines the total propellant masses (based on the volume given by the airframe diameter). This mass is then used to determine the actual lengths and masses of each tank. The algebra involved is pretty simple so I left out its derivation.
It may be a cleaner idea to simply use the total propellant mass as the design variable, but it shouldn't make a difference. Trajectory is sensitive to mass, so be wary of the assumptions made in dry_mass and update them as better information becomes available.
In [3]:
# this class is for conceptual clarity to keep rocket externals and internals organized
class Rocket(EngineSys):
# let the model keep track of its own properties
def __init__(self, m_prop, dia, mdot, p_e, p_ch, T_ch, ke, Re, lower, upper, minimum, x_init):
# rocket's got an engine don't it
EngineSys.__init__(self, mdot, p_e, p_ch, T_ch, ke, Re, lower, upper, minimum)
# rocket trajectory variables
self.t = [0.]
self.alt = [x_init]
self.v = [0.]
self.a = [0.]
## relevant measurements
self.dia = dia
# Airframe frontal area projected onto a circle of diameter given by dia
self.A = pi * (dia/2 + airframe_thickness)**2
# propellant mass
self.m_prop = [m_prop]
# tank geometry to properly contain propellant, split_tanks() expects inches, not meters
self.r, self.l_o, self.l_f = split_tanks(m_prop, dia/0.0254)
# we're doing this to make sure we have at least enough N2 to keep tanks pressurized
# one day we will add more N2 for RCS system. This all follows from ideal gas law
tube_area = pi*(self.r - min_tank_thickness)**2
volume_o = self.l_o * tube_area
volume_f = self.l_f * tube_area
self.N2_mass = rcs_req_mass + prop_tank_pressure * (volume_o + volume_f) * N2_M / (N2_temp * R_univ) # kg
shirt_l = sum([nose_l, ers_l, rcs_l, av_l, n2_l, 4*gaps]) # m
pants_l = system_length(self.l_o, self.l_f)
# length of entire rocket
self.total_length = shirt_l + pants_l
# dry mass of whole engine system
m_eng_sys = system_mass(self.r, self.l_o, self.l_f)
# Total dry mass of rocket
self.m_dry = self.dry_mass(m_eng_sys, dia)
# GLOW
self.m = [self.m_dry + m_prop]
# calculates length-diameter ratio
self.ld_ratio = self.total_length / (dia + airframe_thickness*2)
# Subsystem dry masses, needs sanity check from dirty ME's
# currently not accounting for different material density of avionics module
def dry_mass(self, engine_sys_mass, dia):
bulkheads = bulkhead(body_r(self.r)) # bulkheads in rocket's torso
# carbon fiber composite Airframe mass [kg] of cylinder
m_airframe_cf = (cf_thickness * CF['rho'] + nomex_thickness *
rho_nomex) *pi*dia * (n2_l + l_feed + l_ems + l_engine)
# fiberglass module for avionics system
m_av_frame = airframe_thickness * rho_fiber * pi * dia * av_l
m_ers_rcs_fuelgap_frames = (ers_l + rcs_l + rcs_l) * airframe_thickness * Al['rho'] * pi * dia
m_airframe = m_av_frame + m_nose_frame + m_ers_rcs_fuelgap_frames
return (m_ringsclamps + m_ballast + m_recovery +
m_rcs + self.N2_mass + m_avionics + m_n2 + 7*bulkheads +
engine_sys_mass + m_airframe + m_fins) * airframe_fudger
The previous MDO used empirical drag coefficient data from an Aerobee-150, now I am using drag coefficient data taken straight from OpenRocket for a rocket with nearly equivalent geometry.
We are pretending there is no wind. We can expect about 3 m/s winds at the launch site in reality.
We're accounting for the change in gravity because it is noticable at high altitudes, and this brings our simulation closer to what we will see in an OpenRocket simulation.
From the previous MDO,
"We usually express drag as $$F_d = \frac{1}{2} \rho V^2 C_d A$$
where $\rho$ is the local air density, $A$ is the frontal area, $C_d$ is the drag coefficient which for a given airframe geometry is a function of angle-of-attack (which for a 1-DoF system is identically zero), and Mach number. The Mach number is a function given by
$$ \frac{V}{c} = \frac{V}{\sqrt{kRT}} $$where $c$ is the local speed of sound, $k$ is the gas ratio of specific heats, $R$ is the specific gas constant (assumed constant in altitude), and $T$ is the gas temperature. These numbers, as well as $\rho$, and $p_a$ the ambient pressure, are supplied by the 1976 U.S. Standard Atmosphere model."
In [4]:
# this class is for conceptual clarity to keep our rocket externals organized
class Environment(Rocket):
# let the model keep track of its own properties
def __init__(self, m_prop, dia, mdot, p_e, p_ch, T_ch, ke, Re, lower, upper, minimum, x_init):
# an environment has a rocket in it, ideally
Rocket.__init__(self, m_prop, dia, mdot, p_e, p_ch, T_ch, ke, Re, lower, upper, minimum, x_init)
# Drag coefficient look up
self.C_d_t = []
self.Ma_t = []
csv_file = open('lv4drag.csv') # use estimated drag data from openrocket for one of our lv4 designs
drag_data = csv.reader(csv_file, delimiter=',')
for row in drag_data:
self.Ma_t.append(float(row[0]))
self.C_d_t.append(float(row[1]))
# air
self.g = [self.gravity(self.alt[0])]
self.ka = 1.4 # Ratio of specific heats, air
self.Ra = 287.1 # Avg. specific gas constant (dry air)
p_a, rho, T_a = self.std_at(self.alt[0])
self.p_a = [p_a]
self.rho = [rho]
self.T_a = [T_a]
# drag and thrust
self.D = [0.]
self.q = [0.]
self.Ma = [0.]
self.F = [self.thrust(self.p_a[0], self.throttle[0])+60]
# all your atmospheric needs are here, probably should sanity check, matches openrocket at least
def std_at(self, h): # U.S. 1976 Standard Atmosphere
if h < 11000:
T = 15.04 - 0.00649*h
p = 101.29*((T + 273.1)/288.08)**5.256
elif 11000 <= h and h <25000:
T = -56.46
p = 22.65*exp(1.73 - 0.000157*h)
else:
T = -131.21 + 0.00299*h
p = 2.488 * ((T + 273.1)/216.6)**(-11.388)
p_a = p*1000 # Ambient air pressure [Pa]
rho = p/(0.2869*(T + 273.1)) # Ambient air density [kg/m^3]
T_a = T + 273.1 # Ambient air temperature [K]
return p_a, rho, T_a
# calculates drag force
def drag(self, x, v):
# Check Knudsen number and switch drag models (e.g. rarified gas dyn vs. quadratic drag)
# ^ comment from previous MDO; I'm not sure what the above comment means for us, so I ignored it
p_a, rho, T_a = self.std_at(x)
Ma = v / sqrt(self.ka * self.Ra * T_a)
C_d = np.interp(Ma, self.Ma_t, self.C_d_t) # Drag coefficient function
q = 0.5 * rho * v**2 # Dynamic pressure [Pa]
D = q * C_d * self.A # Drag force [N]
return D, q, Ma
# https://www.sensorsone.com/local-gravity-calculator/
# International Gravity Formula (IGF) 1980
# Geodetic Reference System 1980 (GRS80)
# Free Air Correction (FAC)
def gravity(self, h):
lat = radians(32) # degrees north, launch site
igf = 9.780327 * (1 + 0.0053024 * sin(lat)**2 - 0.0000058 * sin(2*lat)**2)
fac = -3.086 * 10**(-6) * h
return igf + fac
def save_results(self):
arr = np.array([self.t, self.alt, self.v, self.a, self.F, self.m, self.Ma, self.rho]).T
np.savetxt(rkt_prefix +'psas_rocket_'+str(get_index()-1)+'_traj.csv', arr, delimiter=',',
header='time (s),altitude (m),velocity (m/s),acceleration (m/s^2),thrust (N),mass (kg),Mach,air density (kg/m^3)')
I have endeavored to write self-explanatory code in this section. Where I have failed, please inform me and I'll go into detail. If you want more information about Runge-Kutta methods, see https://lpsa.swarthmore.edu/NumInt/NumIntFourth.html
From the previous MDO,
"The trajectory of the rocket, in a single degree of freedom, can be best described by altitude $h(t)$, velocity $V(t)$, and total mass $m(t)$ state variables which are functions of time. The initial values of the state variables are given by $h_0$, $V_0$, and $m_0$ respectively. The governing equation of the rocket trajectory is Newton's 2nd law of motion
$$F = \frac{d(mV)}{dt}.$$This can be expanded in terms of the sum of forces acting on the rocket during free flight going up to apogee
$$Thrust - Drag - m g = \frac{d}{dt} \left(m \frac{dh}{dt}\right)."$$
In [5]:
class Simulation(Environment):
# initialize trajectory sim
def __init__(self, m_prop, mdot, dia, p_e, lower, upper, minimum,
x_init, p_ch, T_ch, ke, Re):
p_ch, dia, p_e = unit_conv(p_ch, dia, p_e) # gently change units
Environment.__init__(self, m_prop, dia, mdot, p_e, p_ch, T_ch, ke, Re, lower, upper, minimum, x_init)
# launch site
self.off_rail = False
self.rail_height = x_init + (launch_tower * 0.3048) # 60 ft, convert to m, scale for site altitude
self.launch_speed = 0. # in case rocket fails to launch
# engine cut-off
self.motor_burnout = False
# by default, use the runge kutta integrator. make this choosable later
self.method = self.rungekutta
# abstracting some of the logic of which regime the rocket is in.
def phase_driver(self, next_step):
# Check if the propellant tanks are currently non-empty.
# If they are, we have a fueled phase
if self.m_prop[-1] > 0:
# check if we have just left launch rail, in order to capture current speed
# this is to ensure aerodynamic stability so nothing goes sideways
if (self.off_rail == False) and (self.alt[-1] >= self.rail_height):
# take average velocity
self.launch_speed = (self.v[-1] + self.v[-2]) / 2
self.off_rail = True # don't check again
# throttle percent of mass flow rate depending on drag force
# to reduce max q and save a couple kilos
# this might make the MEs sad but it might make them happier later
next_step.throttle = self.throttle_engine(self.D[-1])
# burn fuel
next_step.m_prop = self.m_prop[-1] - self.throttle[-1]*self.mdot*self.dt
# calculate thrust at current height
next_step.F = self.thrust(self.p_a[-1], self.throttle[-1])
# If the tanks are empty, not much happens
elif self.m_prop[-1] <= 0:
# there's no propellant now, did we have it last moment (if there was a last moment)?
if (self.motor_burnout == False) and (len(self.m_prop) >= 2) and (self.m_prop[-2] > 0):
# need in order to calculate burn time, not including this unfueled step
self.F_index = len(self.m_prop) - 1
self.motor_burnout = True
# no throttle,
next_step.throttle = 0
# no fuel,
next_step.m_prop = 0
# no thrust.
next_step.F = 0
# keep track of our mass
next_step.m = self.m_dry + self.m_prop[-1]
return next_step
# just a reminder. breathe deep. you're almost there!
def air_exists(self, next_step):
# update air properties
next_step.p_a, next_step.rho, next_step.T_a = self.std_at(self.alt[-1])
# calculate current drag force
next_step.D, next_step.q, next_step.Ma = self.drag(self.alt[-1], self.v[-1])
return next_step
# overkill integrator, in case you want to use the sledgehammer.
# this may be more optimistic than openrocket, so I wouldn't, pending further investigation.
# this is also currently broken until updated to new model
def rungekutta(self, next_step):
# helper func for runge kutta, this is our ODE :)
# v = dx/dt, a = dv/dt
def ode_f(t, x, v, F, D, m, g):
v = v
a = (F - D)/m - g
return np.array([v, a])
k1 = ode_f(self.t[-1],
self.alt[-1],
self.v[-1],
self.F[-1],
self.D[-1],
self.m[-1],
self.g[-1])
alt = self.alt[-1] + k1[0]*self.dt/2
v = self.v[-1] + k1[1]*self.dt/2
drag = self.drag(alt, v)[0]
throttle = self.throttle_engine(drag)
k2 = ode_f(self.t[-1] + self.dt/2,
alt,
v,
0 if self.F[-1]==0 else self.thrust(self.std_at(alt)[0],
throttle),
drag,
self.m_dry + self.m_prop[-1] - (0 if self.F[-1]==0 else throttle*self.mdot*self.dt/2),
self.gravity(alt))
alt = self.alt[-1] + k2[0]*self.dt/2
v = self.v[-1] + k2[1]*self.dt/2
drag = self.drag(alt, v)[0]
throttle = self.throttle_engine(drag)
k3 = ode_f(self.t[-1] + self.dt/2,
alt,
v,
0 if self.F[-1]==0 else self.thrust(self.std_at(alt)[0],
throttle),
drag,
self.m_dry + self.m_prop[-1] - (0 if self.F[-1]==0 else throttle*self.mdot*self.dt/2),
self.gravity(alt))
alt = self.alt[-1] + k3[0]*self.dt
v = self.v[-1] + k3[1]*self.dt
drag = self.drag(alt, v)[0]
throttle = self.throttle_engine(drag)
k4 = ode_f(self.t[-1] + self.dt,
alt,
v,
0 if self.F[-1]==0 else self.thrust(self.std_at(alt)[0],
throttle),
drag,
self.m_dry + self.m_prop[-1] - (0 if self.F[-1]==0 else throttle*self.mdot*self.dt),
self.gravity(alt))
next_step.a = (k1[1] + 2*k2[1] + 2*k3[1] + k4[1]) / 6
next_step.alt = self.alt[-1] + self.dt * (k1[0] + 2*k2[0] + 2*k3[0] + k4[0]) / 6
next_step.v = self.v[-1] + self.dt * next_step.a
return next_step
# this is ugly-looking but efficient
# truly a proper integrator for us engineering folk
# essentially linearizing acceleration and doing a slightly higher-order approximation of velocity
def trapezoidal_int(self, next_step):
k1 = self.v[-1]
k2 = self.v[-1] + self.a[-1] * self.dt/2
k3 = self.v[-1] + self.a[-1] * self.dt
next_step.a = (self.F[-1] - self.D[-1]) / self.m[-1] - self.g[-1]
next_step.alt = self.alt[-1] + self.dt * (k1 + 4*k2 + k3) / 6
next_step.v = self.v[-1] + next_step.a * self.dt
return next_step
# deprecated function for basic monkey integration. please don't call it.
def forwardeuler(self, next_step):
next_step.a = (self.F[-1] - self.D[-1]) / self.m[-1] - self.g[-1]
next_step.alt = self.alt[-1] + self.dt * self.v[-1]
next_step.v = self.v[-1] + self.dt * self.a[-1]
return next_step
# calculate trajectory!!
def do_physics(self, next_step):
next_step.t = self.t[-1] + self.dt
next_step.g = self.gravity(self.alt[-1])
next_step = self.method(next_step)
return next_step
# convenience function, hide the eye-bleed!
def step_forward(self, next_step):
self.t.append(next_step.t)
self.alt.append(next_step.alt)
self.v.append(next_step.v)
self.a.append(next_step.a)
self.g.append(next_step.g)
self.F.append(next_step.F)
self.D.append(next_step.D)
self.m.append(next_step.m)
self.throttle.append(next_step.throttle)
self.m_prop.append(next_step.m_prop)
self.q.append(next_step.q)
self.Ma.append(next_step.Ma)
self.p_a.append(next_step.p_a)
self.rho.append(next_step.rho)
self.T_a.append(next_step.T_a)
Here, all of our abstraction has paid off. As we can see, the integration loop has been reduced to four lines, three of which handle related state variables, and the fourth of which simply moves onto the next step.
There are several variables to calculate for our optimization, which also expects some numpy arrays, so we do that at the end. Note that all our relevant variables are assigned to the simulation object itself, which is returned to the optimization.
In [6]:
# calculates trajectory of a design
# dt is time step [s]
def trajectory(m_prop, mdot, dia, p_e, bounds=(500, 1500), x_init=0, dt=.1,
p_ch=p_ch, T_ch=T_ch, ke=ke, Re=Re):
sim = Simulation(m_prop, mdot, dia, p_e, bounds[0], bounds[1], min_throttle, x_init,
p_ch, T_ch, ke, Re) # get this puppy rolling
sim.dt = dt
next_step = SimpleNamespace() # this is to prevent interference between steps
# now step our rocket through time until we reach apogee
while (sim.v[-1] > 0) or (sim.off_rail == False):
next_step = sim.phase_driver(next_step)
next_step = sim.air_exists(next_step)
next_step = sim.do_physics(next_step)
sim.step_forward(next_step)
# DONE!
# numpy arrays are nicer for some reasons
sim.alt = np.array(sim.alt)
sim.v = np.array(sim.v)
sim.a = np.array(sim.a)
sim.F = np.array(sim.F)
sim.D = np.array(sim.D)
sim.q = np.array(sim.q)
sim.p_a = np.array(sim.p_a)
sim.rho = np.array(sim.rho)
sim.t = np.array(sim.t)
# this will run after we reach apogee to do last-minute calculations and conversions
sim.TWR = sim.a[1] / g_n # Thrust-to-weight ratio constraint
sim.S_crit = sim.p_e / sim.p_a[0] # Sommerfield criterion constraint
sim.massratio = (sim.m_prop[0] + sim.m_dry) / sim.m_dry # Mass ratio
sim.dV1 = sim.Ve * log(sim.massratio) / 1000 # Tsiolkovsky's bane (delta-V)
sim.max_g_force = max(abs(sim.a)) / g_n # calculates top g-force
sim.maxq = max(sim.q)
# ends our suffering and returns all the stats of the launch
return sim