Trajectory Simulation

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.

  1. There are model parameters at the top (and at the top of openrocket_interface) where they are easily accessible in case some basic assumptions change.
  2. The engine is itself a subsystem at the bottom of the hierarchy, so it is the most general class.
    • It only initializes engine-specific variables.
      • mdot, p_e, p_ch, T_ch, k_e, R_e
    • It contains the thrust function (even though thrust technically depends on ambient pressure at the altitude).
  3. The rocket is the main system whose class inherits an engine.
    • This initializes relevant variables.
      • t, alt, v, a, m_prop, m_dry, length, airframe frontal area
    • Rocket is concerned only with location, geometry, and mass.
  4. The physical environment is a supersystem in which the rocket is embedded (i.e. the environment inherits a rocket).
    • This determines air characteristics, drag profile (based off .csv either from OpenRocket or an actual Aerobee), and initial thrust (since that depends on initial location).
    • This keeps track of the atmosphere, drag, and gravitational force.
  5. Then there is the trajectory itself and its initialization, as well as the logical and mathematical details of each discrete step of the numerical integration.
  6. Finally, there is the overall simulation from start to finish, as well as the calculation of variables that are relevant to the optimization. trajectory returns the whole simulation as an object rather than a list of variables from it.

Assumptions

Model

  • The assumptions in openrocket_interface are consistent.
  • The launch tower is 60' tall.
  • Actual thrust is 100% of ideal thrust.
    • At a later date, we may wish to model inefficiencies.
  • The nose cone is 1.5 m long.
  • The ERS and RCS subsystems are 6" long.
  • The Avionics and N2 subsystems are 18" long.
  • The gaps between those subsystems is negligible. Note that the gaps between subsystems on the lower half of the rocket are included, so perhaps this is a poor assumption.

Engine

  • The standard equation for thrust isn't garbage.
  • Accounting for the thrust cosine loss is a good idea.
  • Thrust goes straight up z axis.
  • Exit pressure remains constant through-out flight.
  • We may throttle our mass flow rate up to 60%.

Rocket

  • The airframe can be considered a cylinder to calculate its mass if we pretend it's a little thinner than it is (to account for the nose cone tapering). This assumption could be replaced by some geometry, but it is close enough to OpenRocket to hand-wave.
  • There is a host of assumptions about the masses of each subsystem. Nothing too crazy though.

Environment

  • The drag coefficients with respect to Mach are equivalent to the sample used and provided by OpenRocket for a rocket with the same profile but different length. If you don't like this assumption you are free to assume the drag of an empirical Aerobee-150 is equivalent, but until we get a wind tunnel or good CFD sims, this is fine...ish.
  • There are some arcane but reasonable assumptions about the thermodynamic properties of dry air and the US 1976 Standard Atmosphere model.
  • The equation for drag is sane.
  • The equation for gravity at the latitude/altitude of the launch site is a fine approximation.

Trajectory

  • The average velocity around departing the launch tower is a good enough approximation of the actual velocity at that moment.
  • The loss of mass as well as the force from expelling N2 for the RCS system is negligible.
  • Trajectory is 1 dimensional, loss from other 5 degrees of freedom is irrelevant.
  • Internal steps of Runge-Kutta only need to vary the time, altitude, and velocity. In the python file, there are other functions that vary more of the state variables, but that both takes longer and is more optimistic, so this simple integration is better for our purposes.
  • Our discretization of the actual physics involved is reasonable.

Parameters and Helper Functions


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

Engine Subsystem

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

Rocket System

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

Environment Supersystem

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)')

Trajectory Simulation Logic and Math

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)

Top-Level Trajectory Simulator

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