Externally Computed Polars

Externally Computed Polars#

This example shows how to build, using the level-2 interface, an aviary model that includes an external susbsystem that computes a lift and drag polar and passes them into the mission aerodynamics for a 3-phase mission (climb, cruise, descent). During the mission, Aviary will interpolate on the computed polars to compute actual lift and drag for a given flight condition.

We start with the assumption that we have an external component called ExternalAero that can compute the lift and drag at any given altitude, mach number, and angle of attack. The details of such a component may be highly complicated and not important for the purposes of this example. We will be using a structured grid, which assumes the data table is regularly spaced in all dimensions. We want to compute lift and drag over a grid of altitudes (in ‘ft’), mach numbers, and angles of attack given by:

Altitude (ft) [    0.  3000.  6000.  9000. 12000. 15000. 18000. 21000. 24000. 27000.
 30000. 33000. 36000. 38000. 42000.]
Mach [0.   0.2  0.4  0.5  0.6  0.7  0.75 0.8  0.85 0.9 ]
Angle of Attack (deg) [-2.  0.  2.  4.  6.  8. 10.]
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/pyoptsparse/pyOpt_MPI.py:68: UserWarning: mpi4py could not be imported. mpi4py is required to use the parallel gradient analysis and parallel objective analysis for non-gradient based optimizers. Continuing using a dummy MPI module from pyOptSparse.
  warnings.warn(warn)

In a structured grid, interpolation data must be present for every combination of inputs. In other words, our ExternalAero component must run a full factorial of points spanning those 3 variables. The Aviary variable hierarchy includes two variables for the polars: Aircraft.Design.LIFT_POLAR, and Aircraft.Design.DRAG_POLAR. The data in each of these polars should be a n x m x k numpy array, where n is the number of altitudes, m is the number of mach numbers, and k is the number of angles of attack. The ExternalAero will need to compute these values and place them into an array of this shape.

If use of a structured grid is not desirable, then the data does not need to meet these formatting requirements. In that case, the data table does not have to be regularly spaced, and each variable (Altitude, Mach, angle_of_attack, LIFT_POLAR, and DRAG_POLAR) must be 1-dimensional numpy arrays of equal length.

Using the level-2 interface, we create a builder for our external ExternalAero subsystem. In this example, the component produces outputs drag_table and lift_table, but we can specify an alias to Aircraft.Design.DRAG_POLAR and Aircraft.Design.LIFT_POLAR respectively. It is important that we inherit from the AerodynamicsBuilderBase to let Aviary know this is builder produces aerodynamics components. Some mission analysis methods require special handling of aerodynamics components that will not occur if we skip this step.

class ExternalAeroBuilder(av.AerodynamicsBuilderBase):
    """
    An example subsystem builder that adds an external aerodynamics component

    Parameters
    ----------
    aero_data : NamedValues
        Altitude, Mach number, and angle of attack data, all in ascending order.
    """

    def __init__(self, name='aero', altitude=None, mach=None,
                 angle_of_attack=None):
        super().__init__(name)
        self.altitude = altitude
        self.mach = mach
        self.angle_of_attack = angle_of_attack
        
    def build_pre_mission(self, aviary_inputs):
        """
        Build an OpenMDAO system for the pre-mission computations of the subsystem.

        Returns
        -------
        pre_mission_sys : openmdao.core.Group
            An OpenMDAO group containing all computations that need to happen in
            the pre-mission part of the Aviary problem. This includes sizing, design,
            and other non-mission parameters.
        """
        aero_group = om.Group()
        aero = ExternalAero(altitude=self.altitude, mach=self.mach, angle_of_attack=self.angle_of_attack)
        aero_group.add_subsystem(
            'premission_aero', 
            aero,
            promotes_inputs=['*'],
            promotes_outputs=[
                ('drag_table', Aircraft.Design.DRAG_POLAR),
                ('lift_table', Aircraft.Design.LIFT_POLAR)
            ]
        )
        return aero_group

Notice that we have passed the altitude, Mach, and angle of attack arrays into the builder so that the ExternalAero component can use them as instantiation arguments.

Next, we add the builder to our phase_info as usual. We are using a single-aisle commercial transport aircraft and mission.

phase_info = av.default_height_energy_phase_info.copy()

external_aero = ExternalAeroBuilder(name='external_aero',
                                    altitude=altitude, mach=mach, angle_of_attack=angle_of_attack)

phase_info['pre_mission']['external_subsystems'] = [external_aero]

Next, the existing mission phases need to be given the information they need to set up our aerodynamics analysis using phase_info. We use the solved_alpha method of Aviary’s included aerodynamics for this, which can accept the input passed from our external subsystem. Since we are using Aviary’s built-in aerodynamics methods, we use the default name “core_aerodynamics”. Don’t forget to update the subsystem_options for each phase. We must specify the method, the aero_data that contains our altitude, Mach, and angle of attack data, as well as the connect_training_data flag to denote we are passing our drag polars via openMDAO connections.

subsystem_options = {'method': 'solved_alpha',
                     'aero_data': aero_data,
                     'connect_training_data': True}

phase_info['climb']['subsystem_options'] = {'core_aerodynamics': subsystem_options}
phase_info['cruise']['subsystem_options'] = {'core_aerodynamics': subsystem_options}
phase_info['descent']['subsystem_options'] = {'core_aerodynamics': subsystem_options}

Finally, we can instantiate the AviaryProblem like normal. However, we need to tell Aviary the size of our lift and drag polars so that it can allocate the right shape for the connection.

from aviary.utils.functions import get_aviary_resource_path

input_file = get_aviary_resource_path("models/test_aircraft/aircraft_for_bench_GwFm.csv")

prob = av.AviaryProblem()
prob.load_inputs(input_file, phase_info)

# Preprocess inputs
prob.check_and_preprocess_inputs()

# Add correctly-sized polar to aviary_inputs so that input components are sized correctly.
shape = (altitude.size, mach.size, angle_of_attack.size)

prob.aviary_inputs.set_val(Aircraft.Design.LIFT_POLAR, np.zeros(shape), units='unitless')
prob.aviary_inputs.set_val(Aircraft.Design.DRAG_POLAR, np.zeros(shape), units='unitless')
/home/runner/work/Aviary/Aviary/aviary/utils/process_input_decks.py:175: UserWarning: Variable 'aircraft:engine:engine_mass_specific' is not in meta_data nor in 'guess_names'. It will be ignored.
  warnings.warn(
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:897: OMDeprecationWarning:None: The method `add_polynomial_control` is deprecated and will be removed in Dymos 2.1. Please use `add_control` with the appropriate options to define a polynomial control.
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:2323: RuntimeWarning: Invalid options for non-optimal control 'mach' in phase 'climb': lower, upper, ref
  warnings.warn(f"Invalid options for non-optimal control '{name}' in phase "
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:2323: RuntimeWarning: Invalid options for non-optimal control 'altitude' in phase 'climb': lower, upper, ref
  warnings.warn(f"Invalid options for non-optimal control '{name}' in phase "
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:2323: RuntimeWarning: Invalid options for non-optimal control 'mach' in phase 'cruise': lower, upper, ref
  warnings.warn(f"Invalid options for non-optimal control '{name}' in phase "
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:2323: RuntimeWarning: Invalid options for non-optimal control 'altitude' in phase 'cruise': lower, upper, ref
  warnings.warn(f"Invalid options for non-optimal control '{name}' in phase "
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:2323: RuntimeWarning: Invalid options for non-optimal control 'mach' in phase 'descent': lower, upper, ref
  warnings.warn(f"Invalid options for non-optimal control '{name}' in phase "
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:2323: RuntimeWarning: Invalid options for non-optimal control 'altitude' in phase 'descent': lower, upper, ref
  warnings.warn(f"Invalid options for non-optimal control '{name}' in phase "
The following variables have been overridden:
  'aircraft:design:lift_curve_slope  7.1765  1/rad
  'aircraft:fuel:auxiliary_fuel_capacity  0  lbm
  'aircraft:fuel:total_capacity  45694  lbm
  'aircraft:fuel:wing_volume_geometric_max  1114  ft**3
  'aircraft:fuselage:avg_diameter  12.75  ft
  'aircraft:fuselage:length  128  ft
  'aircraft:fuselage:mass  18357.13345514  lbm
  'aircraft:fuselage:planform_area  1578.24  ft**2
  'aircraft:fuselage:wetted_area  4158.62  ft**2
  'aircraft:horizontal_tail:wetted_area  592.65  ft**2
  'aircraft:nacelle:avg_diameter  [7.94]  ft
  'aircraft:nacelle:avg_length  [12.3]  ft
  'aircraft:nacelle:fineness  [2.]  unitless
  'aircraft:nacelle:surface_area  [329.615]  ft**2
  'aircraft:vertical_tail:wetted_area  581.13  ft**2
  'aircraft:wing:area  1370  ft**2
  'aircraft:wing:aspect_ratio  10.13  unitless
  'aircraft:wing:average_chord  12.615  ft
  'aircraft:wing:span  117.83  ft
  'aircraft:wing:thickness_to_chord_unweighted  0.1397  unitless
  'aircraft:wing:ultimate_load_factor  3.91650835  unitless
  'aircraft:wing:wetted_area  2396.56  ft**2

--- Constraint Report [traj] ---
    --- climb ---
        [path]    0.0000e+00 <= throttle <= 1.0000e+00  [unitless]
    --- cruise ---
        [initial] 0.0000e+00 <= throttle <= 1.0000e+00  [unitless]
        [final]   0.0000e+00 <= throttle <= 1.0000e+00  [unitless]
    --- descent ---
        [path]    0.0000e+00 <= throttle <= 1.0000e+00  [unitless]
=============================================================
pre_mission.core_subsystems.core_mass.wing_mass.isolated_mass
=============================================================
NL: Newton 0 ; 1709.38381 1
NL: Newton 1 ; 0.0712464535 4.16796117e-05
NL: Newton 2 ; 1.20053301e-10 7.02319162e-14
NL: Newton Converged

===============================================
pre_mission.core_subsystems.core_mass.fuel_mass
===============================================
Warning: req_fuel_mass > max_wingfuel_mass, adding a body tank
Warning: req_fuel_mass > max_wingfuel_mass, adding a body tank
NL: Newton 0 ; 45744.3838 1
Warning: req_fuel_mass > max_wingfuel_mass, adding a body tank
Warning: req_fuel_mass > max_wingfuel_mass, adding a body tank
NL: Newton 1 ; 10.8955866 0.000238184139
Warning: req_fuel_mass > max_wingfuel_mass, adding a body tank
Warning: req_fuel_mass > max_wingfuel_mass, adding a body tank
NL: Newton 2 ; 0 0
NL: Newton Converged

=========================
traj.phases.climb.rhs_all
=========================
+  
+  ===========================================
+  traj.phases.climb.rhs_all.core_aerodynamics
+  ===========================================
+  NL: Newton 0 ; 2.59001986 1
+  NL: Newton 1 ; 0.000434550464 0.000167778816
+  NL: Newton 2 ; 9.81087013e-11 3.7879517e-11
+  NL: Newton Converged
NL: Newton 0 ; 0.0742154392 1
+  
+  ===========================================
+  traj.phases.climb.rhs_all.core_aerodynamics
+  ===========================================
+  NL: Newton 0 ; 3.46518961e-15 1
+  NL: Newton Converged
NL: Newton 1 ; 1.44208449e-05 0.000194310578
+  
+  ===========================================
+  traj.phases.climb.rhs_all.core_aerodynamics
+  ===========================================
+  NL: Newton 0 ; 4.12741161e-15 1
+  NL: Newton Converged
NL: Newton 2 ; 4.1275495e-15 5.56157794e-14
NL: Newton Converged

===============================
traj.phases.cruise.indep_states
===============================
NL: Newton Converged in 0 iterations

==========================
traj.phases.cruise.rhs_all
==========================
+  
+  ============================================
+  traj.phases.cruise.rhs_all.core_aerodynamics
+  ============================================
+  NL: Newton 0 ; 2.18733565 1
+  NL: Newton 1 ; 0.000368185841 0.000168326174
+  NL: Newton 2 ; 1.26385647e-11 5.7780637e-12
+  NL: Newton Converged
NL: Newton 0 ; 0.0114704127 1
+  
+  ============================================
+  traj.phases.cruise.rhs_all.core_aerodynamics
+  ============================================
+  NL: Newton 0 ; 1.0215403e-15 1
+  NL: Newton Converged
NL: Newton 1 ; 5.59006288e-07 4.87346274e-05
+  
+  ============================================
+  traj.phases.cruise.rhs_all.core_aerodynamics
+  ============================================
+  NL: Newton 0 ; 7.63386316e-16 1
+  NL: Newton Converged
NL: Newton 2 ; 7.63522833e-16 6.65645478e-14
NL: Newton Converged

================================
traj.phases.descent.indep_states
================================
NL: Newton Converged in 0 iterations

===========================
traj.phases.descent.rhs_all
===========================
+  
+  =============================================
+  traj.phases.descent.rhs_all.core_aerodynamics
+  =============================================
+  NL: Newton 0 ; 2.25049651 1
+  NL: Newton 1 ; 0.000216424917 9.61676306e-05
+  NL: Newton 2 ; 1.97108295e-11 8.75843595e-12
+  NL: Newton Converged
NL: Newton 0 ; 0.123935965 1
+  
+  =============================================
+  traj.phases.descent.rhs_all.core_aerodynamics
+  =============================================
+  NL: Newton 0 ; 5.4603626e-16 1
+  NL: Newton Converged
NL: Newton 1 ; 0.00368270888 0.0297146101
+  
+  =============================================
+  traj.phases.descent.rhs_all.core_aerodynamics
+  =============================================
+  NL: Newton 0 ; 5.70316273e-16 1
+  NL: Newton Converged
NL: Newton 2 ; 1.80675877e-05 0.000145781635
+  
+  =============================================
+  traj.phases.descent.rhs_all.core_aerodynamics
+  =============================================
+  NL: Newton 0 ; 3.0800599e-16 1
+  NL: Newton Converged
NL: Newton 3 ; 3.08104732e-16 2.48599937e-15
NL: Newton Converged