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')
User has specified Design.NUM_* passenger values but CrewPyaload.NUM_* has been left blank or set to zero.
Assuming they are equal to maintain backwards compatibility with GASP and FLOPS output files.
If you intended to have no passengers on this flight, please set Aircraft.CrewPayload.TOTAL_PAYLOAD_MASS to zero in aviary_values.
/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.
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
/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 "
--- 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.solver_sub
====================================
+
+ ======================================================
+ traj.phases.climb.rhs_all.solver_sub.core_aerodynamics
+ ======================================================
+ NL: Newton 0 ; 2.58915687 1
+ NL: Newton 1 ; 0.000434175902 0.000167690072
+ NL: Newton 2 ; 9.80418557e-11 3.78663251e-11
+ NL: Newton Converged
NL: Newton 0 ; 0.0742164013 1
+
+ ======================================================
+ traj.phases.climb.rhs_all.solver_sub.core_aerodynamics
+ ======================================================
+ NL: Newton 0 ; 3.43770326e-15 1
+ NL: Newton Converged
NL: Newton 1 ; 1.4421144e-05 0.00019431209
+
+ ======================================================
+ traj.phases.climb.rhs_all.solver_sub.core_aerodynamics
+ ======================================================
+ NL: Newton 0 ; 3.22409496e-15 1
+ NL: Newton Converged
NL: Newton 2 ; 3.22427199e-15 4.34441974e-14
NL: Newton Converged
===============================
traj.phases.cruise.indep_states
===============================
NL: Newton Converged in 0 iterations
=====================================
traj.phases.cruise.rhs_all.solver_sub
=====================================
+
+ =======================================================
+ traj.phases.cruise.rhs_all.solver_sub.core_aerodynamics
+ =======================================================
+ NL: Newton 0 ; 2.18415672 1
+ NL: Newton 1 ; 0.000370762304 0.000169750778
+ NL: Newton 2 ; 1.27843307e-11 5.85321126e-12
+ NL: Newton Converged
NL: Newton 0 ; 0.0114651758 1
+
+ =======================================================
+ traj.phases.cruise.rhs_all.solver_sub.core_aerodynamics
+ =======================================================
+ NL: Newton 0 ; 9.59985337e-16 1
+ NL: Newton Converged
NL: Newton 1 ; 5.58856097e-07 4.87437879e-05
+
+ =======================================================
+ traj.phases.cruise.rhs_all.solver_sub.core_aerodynamics
+ =======================================================
+ NL: Newton 0 ; 1.06696405e-15 1
+ NL: Newton Converged
NL: Newton 2 ; 1.06702762e-15 9.30668345e-14
NL: Newton Converged
================================
traj.phases.descent.indep_states
================================
NL: Newton Converged in 0 iterations
======================================
traj.phases.descent.rhs_all.solver_sub
======================================
+
+ ========================================================
+ traj.phases.descent.rhs_all.solver_sub.core_aerodynamics
+ ========================================================
+ NL: Newton 0 ; 2.24938732 1
+ NL: Newton 1 ; 0.00021530306 9.57163129e-05
+ NL: Newton 2 ; 1.942e-11 8.63346202e-12
+ NL: Newton Converged
NL: Newton 0 ; 0.123934084 1
+
+ ========================================================
+ traj.phases.descent.rhs_all.solver_sub.core_aerodynamics
+ ========================================================
+ NL: Newton 0 ; 1.20982351e-15 1
+ NL: Newton Converged
NL: Newton 1 ; 0.00368262399 0.0297143761
+
+ ========================================================
+ traj.phases.descent.rhs_all.solver_sub.core_aerodynamics
+ ========================================================
+ NL: Newton 0 ; 1.10441271e-15 1
+ NL: Newton Converged
NL: Newton 2 ; 1.77757596e-05 0.000143429144
+
+ ========================================================
+ traj.phases.descent.rhs_all.solver_sub.core_aerodynamics
+ ========================================================
+ NL: Newton 0 ; 4.65661287e-16 1
+ NL: Newton Converged
NL: Newton 3 ; 4.65941674e-16 3.75959266e-15
NL: Newton Converged