WIP: N3hybrid

environment.yml

@@ -12,7 +12,7 @@ dependencies:
   - pip
   - pip:
     - pytest
-    - openmdao
+    - openmdao<=3.2.1
     - six
     - sphinx_rtd_theme
-    - redbaron
+    - redbaron

examples/B738_Dymos.py

@@ -0,0 +1,156 @@
+from __future__ import division
+import sys
+import os
+import numpy as np
+
+sys.path.insert(0, os.getcwd())
+import openmdao.api as om
+import openconcept.api as oc
+# imports for the airplane model itself
+from openconcept.analysis.aerodynamics import PolarDrag
+from examples.aircraft_data.B738 import data as acdata
+from openconcept.analysis.performance.mission_profiles import MissionWithReserve
+from openconcept.components.cfm56 import CFM56
+
+class B738AirplaneModel(oc.IntegratorGroup):
+    """
+    A custom model specific to the Boeing 737-800 airplane.
+    This class will be passed in to the mission analysis code.
+
+    """
+    def initialize(self):
+        self.options.declare('num_nodes', default=1)
+        self.options.declare('flight_phase', default=None)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        flight_phase = self.options['flight_phase']
+
+
+        # a propulsion system needs to be defined in order to provide thrust
+        # information for the mission analysis code
+        # propulsion_promotes_outputs = ['fuel_flow', 'thrust']
+        propulsion_promotes_inputs = ["fltcond|*", "throttle"]
+
+        self.add_subsystem('propmodel', CFM56(num_nodes=nn, plot=False),
+                           promotes_inputs=propulsion_promotes_inputs)
+
+        doubler = om.ExecComp(['thrust=2*thrust_in', 'fuel_flow=2*fuel_flow_in'], 
+                  thrust_in={'value': 1.0*np.ones((nn,)),
+                     'units': 'kN'},
+                  thrust={'value': 1.0*np.ones((nn,)),
+                       'units': 'kN'},
+                  fuel_flow={'value': 1.0*np.ones((nn,)),
+                     'units': 'kg/s',
+                     'tags': ['integrate', 'state_name:fuel_used', 'state_units:kg', 'state_val:1.0', 'state_promotes:True']},
+                  fuel_flow_in={'value': 1.0*np.ones((nn,)),
+                       'units': 'kg/s'})
+
+        self.add_subsystem('doubler', doubler, promotes_outputs=['*'])
+        self.connect('propmodel.thrust', 'doubler.thrust_in')
+        self.connect('propmodel.fuel_flow', 'doubler.fuel_flow_in')
+
+        # use a different drag coefficient for takeoff versus cruise
+        if flight_phase not in ['v0v1', 'v1v0', 'v1vr', 'rotate']:
+            cd0_source = 'ac|aero|polar|CD0_cruise'
+        else:
+            cd0_source = 'ac|aero|polar|CD0_TO'
+        self.add_subsystem('drag', PolarDrag(num_nodes=nn),
+                           promotes_inputs=['fltcond|CL', 'ac|geom|*', ('CD0', cd0_source),
+                                            'fltcond|q', ('e', 'ac|aero|polar|e')],
+                           promotes_outputs=['drag'])
+
+        # generally the weights module will be custom to each airplane
+        passthru = om.ExecComp('OEW=x',
+                  x={'value': 1.0,
+                     'units': 'kg'},
+                  OEW={'value': 1.0,
+                       'units': 'kg'})
+        self.add_subsystem('OEW', passthru,
+                           promotes_inputs=[('x', 'ac|weights|OEW')],
+                           promotes_outputs=['OEW'])
+
+        self.add_subsystem('weight', oc.AddSubtractComp(output_name='weight',
+                                                     input_names=['ac|weights|MTOW', 'fuel_used'],
+                                                     units='kg', vec_size=[1, nn],
+                                                     scaling_factors=[1, -1]),
+                           promotes_inputs=['*'],
+                           promotes_outputs=['weight'])
+
+class B738AnalysisGroup(om.Group):
+    def setup(self):
+        # Define number of analysis points to run pers mission segment
+        nn = 11
+
+        # Define a bunch of design varaiables and airplane-specific parameters
+        dv_comp = self.add_subsystem('dv_comp',  oc.DictIndepVarComp(acdata),
+                                     promotes_outputs=["*"])
+        dv_comp.add_output_from_dict('ac|aero|CLmax_TO')
+        dv_comp.add_output_from_dict('ac|aero|polar|e')
+        dv_comp.add_output_from_dict('ac|aero|polar|CD0_TO')
+        dv_comp.add_output_from_dict('ac|aero|polar|CD0_cruise')
+
+        dv_comp.add_output_from_dict('ac|geom|wing|S_ref')
+        dv_comp.add_output_from_dict('ac|geom|wing|AR')
+        dv_comp.add_output_from_dict('ac|geom|wing|c4sweep')
+        dv_comp.add_output_from_dict('ac|geom|wing|taper')
+        dv_comp.add_output_from_dict('ac|geom|wing|toverc')
+        dv_comp.add_output_from_dict('ac|geom|hstab|S_ref')
+        dv_comp.add_output_from_dict('ac|geom|hstab|c4_to_wing_c4')
+        dv_comp.add_output_from_dict('ac|geom|vstab|S_ref')
+
+        dv_comp.add_output_from_dict('ac|geom|nosegear|length')
+        dv_comp.add_output_from_dict('ac|geom|maingear|length')
+
+        dv_comp.add_output_from_dict('ac|weights|MTOW')
+        dv_comp.add_output_from_dict('ac|weights|W_fuel_max')
+        dv_comp.add_output_from_dict('ac|weights|MLW')
+        dv_comp.add_output_from_dict('ac|weights|OEW')
+
+        dv_comp.add_output_from_dict('ac|propulsion|engine|rating')
+
+        dv_comp.add_output_from_dict('ac|num_passengers_max')
+        dv_comp.add_output_from_dict('ac|q_cruise')
+
+        # Run a full mission analysis including takeoff, reserve_, cruise,reserve_ and descereserve_nt
+        analysis = self.add_subsystem('analysis',
+                                      MissionWithReserve(num_nodes=nn,
+                                                          aircraft_model=B738AirplaneModel),
+                                      promotes_inputs=['*'], promotes_outputs=['*'])
+
+# def configure_problem():
+#     prob = om.Problem()
+#     prob.model = B738AnalysisGroup()
+#     prob.model.nonlinear_solver = om.NewtonSolver(iprint=2,solve_subsystems=True)
+#     prob.model.linear_solver = om.DirectSolver()
+#     prob.model.nonlinear_solver.options['maxiter'] = 20
+#     prob.model.nonlinear_solver.options['atol'] = 1e-6
+#     prob.model.nonlinear_solver.options['rtol'] = 1e-6
+#     prob.model.nonlinear_solver.linesearch = om.BoundsEnforceLS(bound_enforcement='scalar', print_bound_enforce=False)
+#     return prob
+
+# def set_values(prob, num_nodes):
+#     # set some (required) mission parameters. Each pahse needs a vertical and air-speed
+#     # the entire mission needs a cruise altitude and range
+#     prob.set_val('climb.fltcond|vs', np.linspace(2300.,  600.,num_nodes), units='ft/min')
+#     prob.set_val('climb.fltcond|Ueas', np.linspace(230, 220,num_nodes), units='kn')
+#     prob.set_val('cruise.fltcond|vs', np.ones((num_nodes,)) * 4., units='ft/min')
+#     prob.set_val('cruise.fltcond|Ueas', np.linspace(265, 258, num_nodes), units='kn')
+#     prob.set_val('descent.fltcond|vs', np.linspace(-1000, -150, num_nodes), units='ft/min')
+#     prob.set_val('descent.fltcond|Ueas', np.ones((num_nodes,)) * 250, units='kn')
+#     prob.set_val('reserve_climb.fltcond|vs', np.linspace(3000.,  2300.,num_nodes), units='ft/min')
+#     prob.set_val('reserve_climb.fltcond|Ueas', np.linspace(230, 230,num_nodes), units='kn')
+#     prob.set_val('reserve_cruise.fltcond|vs', np.ones((num_nodes,)) * 4., units='ft/min')
+#     prob.set_val('reserve_cruise.fltcond|Ueas', np.linspace(250, 250, num_nodes), units='kn')
+#     prob.set_val('reserve_descent.fltcond|vs', np.linspace(-800, -800, num_nodes), units='ft/min')
+#     prob.set_val('reserve_descent.fltcond|Ueas', np.ones((num_nodes,)) * 250, units='kn')
+#     prob.set_val('loiter.fltcond|vs', np.linspace(0.0, 0.0, num_nodes), units='ft/min')
+#     prob.set_val('loiter.fltcond|Ueas', np.ones((num_nodes,)) * 200, units='kn')
+#     prob.set_val('cruise|h0',33000.,units='ft')
+#     prob.set_val('reserve|h0',15000.,units='ft')
+#     prob.set_val('mission_range',2050,units='NM')
+
+
+
+if __name__ == "__main__":
+    run_738_analysis(plots=True)    

examples/tests/test_example_aircraft.py

@@ -115,6 +115,8 @@ def setUp(self):
     def test_values_B738(self):
         prob = self.prob
         # block fuel
-        assert_near_equal(prob.get_val('descent.fuel_used_final', units='lbm'), 28688.32933661591, tolerance=2e-5)
+        assert_near_equal(prob.get_val('descent.fuel_used_final', units='lbm'), 28549.432517, tolerance=3e-4)
+        # changelog: 9/2020 - previously 28688.329, updated CFM surrogate model to reject spurious high Mach, low altitude points
         # total fuel
-        assert_near_equal(prob.get_val('loiter.fuel_used_final', units='lbm'), 34555.31347454542, tolerance=2e-5)
+        assert_near_equal(prob.get_val('loiter.fuel_used_final', units='lbm'), 34424.68533072, tolerance=3e-4)
+        # changelog: 9/2020 - previously 34555.313, updated CFM surrogate model to reject spurious high Mach, low altitude points

openconcept/analysis/atmospherics/compute_atmos_props.py

@@ -20,6 +20,8 @@ class ComputeAtmosphericProperties(Group):
         Altitude (vector, km)
     fltcond|Ueas : float
         Equivalent airspeed (vector, m/s)
+    fltcond|TempIncrement : float
+        Temperature increment for non-standard day (vector, degC)
 
     Outputs
     -------
@@ -50,115 +52,13 @@ def initialize(self):
     def setup(self):
         nn = self.options['num_nodes']
         tas_in = self.options['true_airspeed_in']
-        self.add_subsystem('inputconv', InputConverter(num_nodes=nn),promotes_inputs=['*'])
-        self.add_subsystem('temp', TemperatureComp(num_nodes=nn))
-        self.add_subsystem('pressure',PressureComp(num_nodes=nn))
-        self.add_subsystem('density',DensityComp(num_nodes=nn))
-        self.add_subsystem('speedofsound',SpeedOfSoundComp(num_nodes=nn))
-        self.add_subsystem('outputconv',OutputConverter(num_nodes=nn),promotes_outputs=['*'])
+        self.add_subsystem('temp', TemperatureComp(num_nodes=nn), promotes_inputs=['fltcond|h','fltcond|TempIncrement'], promotes_outputs=['fltcond|T'])
+        self.add_subsystem('pressure',PressureComp(num_nodes=nn), promotes_inputs=['fltcond|h'], promotes_outputs=['fltcond|p'])
+        self.add_subsystem('density',DensityComp(num_nodes=nn), promotes_inputs=['fltcond|p', 'fltcond|T'], promotes_outputs=['fltcond|rho'])
+        self.add_subsystem('speedofsound',SpeedOfSoundComp(num_nodes=nn), promotes_inputs=['fltcond|T'], promotes_outputs=['fltcond|a'])
         if tas_in:
             self.add_subsystem('equivair',EquivalentAirspeedComp(num_nodes=nn),promotes_inputs=['*'],promotes_outputs=['*'])
         else:
             self.add_subsystem('trueair',TrueAirspeedComp(num_nodes=nn),promotes_inputs=['*'],promotes_outputs=['*'])
         self.add_subsystem('dynamicpressure',DynamicPressureComp(num_nodes=nn),promotes_inputs=["*"],promotes_outputs=["*"])
-        self.add_subsystem('machnumber',MachNumberComp(num_nodes=nn),promotes_inputs=["*"],promotes_outputs=["*"])
-
-        self.connect('inputconv.h_km','temp.h_km')
-        self.connect('inputconv.h_km','pressure.h_km')
-        self.connect('pressure.p_MPa','density.p_MPa')
-        self.connect('temp.T_1e2_K',['density.T_1e2_K','speedofsound.T_1e2_K'])
-        self.connect('pressure.p_MPa','outputconv.p_MPa')
-        self.connect('temp.T_1e2_K','outputconv.T_1e2_K')
-        self.connect('speedofsound.a_1e2_ms','outputconv.a_1e2_ms')
-        self.connect('density.rho_kg_m3','outputconv.rho_kg_m3')
-
-
-class InputConverter(ExplicitComponent):
-    """
-    This component adds a unitized interface to the Hwang and Jasa model.
-
-    Inputs
-    ------
-    fltcond|h : float
-        Altitude (vector, km)
-
-    Outputs
-    -------
-    h_km : float
-        Altitude in km to pass to the standard atmosphere modules (vector, unitless)
-
-    Options
-    -------
-    num_nodes : int
-        Number of analysis points to run (sets vec length) (default 1)
-    """
-    def initialize(self):
-        self.options.declare('num_nodes', default=1, desc='Number of flight/control conditions')
-    def setup(self):
-        nn = self.options['num_nodes']
-        self.add_input('fltcond|h', units='km', desc='Flight condition altitude',shape=(nn,))
-        #outputs and partials
-        self.add_output('h_km', desc='Height in kilometers with no units',shape=(nn,))
-        self.declare_partials('h_km','fltcond|h', rows=range(nn), cols=range(nn))
-
-    def compute(self, inputs, outputs):
-        outputs['h_km'] = inputs['fltcond|h']
-
-    def compute_partials(self, inputs, J):
-        nn = self.options['num_nodes']
-        J['h_km','fltcond|h'] = np.ones(nn)
-
-
-class OutputConverter(ExplicitComponent):
-    """
-    This component adds a unitized interface to the Hwang and Jasa model.
-
-    Inputs
-    ------
-    p_MPa : float
-        Pressure in megapascals from the standard atm model (vector, unitless)
-    T_1e2_K : float
-        Tempreature in 100K units from the std atm model (vector, unitless)
-    rho_kg_m3 : float
-        Density in kg / m3 from the std atm model (vector, unitless)
-
-    Outputs
-    -------
-    fltcond|p : float
-        Pressure with units (vector, Pa)
-    fltcond|rho : float
-        Density with units (vector, kg/m3)
-    fltcond|T : float
-        Temperature with units (vector, K)
-
-    Options
-    -------
-    num_nodes : int
-        Number of analysis points to run (sets vec length) (default 1)
-    """
-    def initialize(self):
-        self.options.declare('num_nodes', default=1, desc='Number of flight/control conditions')
-    def setup(self):
-        nn = self.options['num_nodes']
-        self.add_input('p_MPa', desc='Flight condition pressures',shape=(nn,))
-        self.add_input('T_1e2_K', desc='Flight condition temp',shape=(nn,))
-        self.add_input('rho_kg_m3', desc='Flight condition density',shape=(nn,))
-        self.add_input('a_1e2_ms', desc='Flight condition speed of sound',shape=(nn,))
-
-        #outputs and partials
-        self.add_output('fltcond|p', units='Pa', desc='Flight condition pressure with units',shape=(nn,))
-        self.add_output('fltcond|rho', units='kg * m**-3', desc='Flight condition density with units',shape=(nn,))
-        self.add_output('fltcond|T', units='K', desc='Flight condition temp with units',shape=(nn,))
-        self.add_output('fltcond|a', units='m * s**-1', desc='Flight condition speed of sound with units',shape=(nn,))
-
-
-        self.declare_partials(['fltcond|p'], ['p_MPa'], rows=range(nn), cols=range(nn), val=1e6*np.ones(nn))
-        self.declare_partials(['fltcond|rho'], ['rho_kg_m3'], rows=range(nn), cols=range(nn), val=np.ones(nn))
-        self.declare_partials(['fltcond|T'], ['T_1e2_K'], rows=range(nn), cols=range(nn), val=100*np.ones(nn))
-        self.declare_partials(['fltcond|a'], ['a_1e2_ms'], rows=range(nn), cols=range(nn), val=100*np.ones(nn))
-
-    def compute(self, inputs, outputs):
-        outputs['fltcond|p'] = inputs['p_MPa'] * 1e6
-        outputs['fltcond|rho'] = inputs['rho_kg_m3']
-        outputs['fltcond|T'] = inputs['T_1e2_K'] * 100
-        outputs['fltcond|a'] = inputs['a_1e2_ms'] * 100
+        self.add_subsystem('machnumber',MachNumberComp(num_nodes=nn),promotes_inputs=["*"],promotes_outputs=["*"])

openconcept/analysis/atmospherics/density_comp.py

@@ -21,26 +21,26 @@ def initialize(self):
     def setup(self):
         num_points = self.options['num_nodes']
 
-        self.add_input('p_MPa', shape=num_points)
-        self.add_input('T_1e2_K', shape=num_points)
-        self.add_output('rho_kg_m3', shape=num_points)
+        self.add_input('fltcond|p', shape=num_points, units='Pa')
+        self.add_input('fltcond|T', shape=num_points, units='K')
+        self.add_output('fltcond|rho', shape=num_points, units='kg / m**3')
 
         arange = np.arange(num_points)
-        self.declare_partials('rho_kg_m3', 'p_MPa', rows=arange, cols=arange)
-        self.declare_partials('rho_kg_m3', 'T_1e2_K', rows=arange, cols=arange)
+        self.declare_partials('fltcond|rho', 'fltcond|p', rows=arange, cols=arange)
+        self.declare_partials('fltcond|rho', 'fltcond|T', rows=arange, cols=arange)
 
     def compute(self, inputs, outputs):
-        p_Pa = inputs['p_MPa'] * 1e6
-        T_K = inputs['T_1e2_K'] * 1e2
+        p_Pa = inputs['fltcond|p']
+        T_K = inputs['fltcond|T']
 
-        outputs['rho_kg_m3'] = p_Pa / R / T_K
+        outputs['fltcond|rho'] = p_Pa / R / T_K
 
     def compute_partials(self, inputs, partials):
-        p_Pa = inputs['p_MPa'] * 1e6
-        T_K = inputs['T_1e2_K'] * 1e2
+        p_Pa = inputs['fltcond|p']
+        T_K = inputs['fltcond|T']
 
         data = 1.0 / R / T_K
-        partials['rho_kg_m3', 'p_MPa'] = data * 1e6
+        partials['fltcond|rho', 'fltcond|p'] = data
 
         data = -p_Pa / R / T_K ** 2
-        partials['rho_kg_m3', 'T_1e2_K'] = data * 1e2
+        partials['fltcond|rho', 'fltcond|T'] = data

openconcept/analysis/atmospherics/pressure_comp.py

@@ -22,26 +22,26 @@ def initialize(self):
     def setup(self):
         num_points = self.options['num_nodes']
 
-        self.add_input('h_km', shape=num_points)
-        self.add_output('p_MPa', shape=num_points, lower=0.)
+        self.add_input('fltcond|h', shape=num_points, units='m')
+        self.add_output('fltcond|p', shape=num_points, lower=0., units='Pa')
 
         arange = np.arange(num_points)
-        self.declare_partials('p_MPa', 'h_km', rows=arange, cols=arange)
+        self.declare_partials('fltcond|p', 'fltcond|h', rows=arange, cols=arange)
 
     def compute(self, inputs, outputs):
         num_points = self.options['num_nodes']
 
-        h_m = inputs['h_km'] * 1e3
+        h_m = inputs['fltcond|h']
         self.tropos_mask, self.strato_mask, self.smooth_mask = get_mask_arrays(h_m)
         p_Pa = compute_pressures(h_m, self.tropos_mask, self.strato_mask, self.smooth_mask)
 
-        outputs['p_MPa'] = p_Pa / 1e6
+        outputs['fltcond|p'] = p_Pa
 
     def compute_partials(self, inputs, partials):
         num_points = self.options['num_nodes']
 
-        h_m = inputs['h_km'] * 1e3
+        h_m = inputs['fltcond|h']
 
         derivs = compute_pressure_derivs(h_m, self.tropos_mask, self.strato_mask, self.smooth_mask)
 
-        partials['p_MPa', 'h_km'] = derivs * 1e3 / 1e6
+        partials['fltcond|p', 'fltcond|h'] = derivs