Using Aviary and OpenAeroStruct Together#

This is an example of an external subsystem using the OpenAeroStruct (OAS) structural analysis system to perform aerostructural optimization of a typical large single aisle transport aircraft wing. The subsystem is based on the OpenAeroStruct aerostructural optimization with wingbox example problem.

This example performs a sub-optimization for minimum wing weight that is then used by Aviary. Another use case would be to perform a structural analysis only. Structural design variables would be passed to the subsystem from Aviary and wing weight and a constraint or constraints representing the structural analysis would be passed back to Aviary for use by the top level optimization.

Motivation#

There may be a need for a higher fidelity tool to compute wing weight instead of relying on the empirical methods in core Aviary. A structural analysis external tool is usually used because of an unusual aircraft configuration that may not be predicted my Aviary empirical weight estimation methods, but in this example case it is simply a demonstration of an external capability to compute wing weight.

External Dependencies#

The user must install OpenAeroStruct into their Python environment using the command ‘pip install openaerostruct’. The user must also install the ambiance package using the command ‘pip install ambiance’.

Subsystem Details#

There are two parts to building an external subsystem – the analysis tool and the Aviary external subsystem builder for that tool. The analysis tool takes inputs and parameters from Aviary and return outputs that Aviary can use to override existing variables. The subsystem builder uses the Aviary external subsystem builder template to connect the analysis tool to Aviary as either a pre-mission, mission or post-mission subsystem.

For this case, the analysis tool will compute a wing weight in the pre-mission portion of the Aviary analysis and return its value to Aviary to override the empirical wing weight value. Fuel weight is passed in from Aviary as the only input currently, but other inputs may also be passed in through the subsystem builder, OAS_wing_weight_builder, by the promotes_inputs parameter. Other Aviary variables can also be added as additional inputs based on user needs.

Note

Some modifications of the OAS_wing_weight_analysis code will be necessary to add new inputs not already defined.

Here is the builder object for the OAS wing weight analysis example:

# %load ../../examples/external_subsystems/OAS_weight/OAS_wing_weight_builder.py
"""
Builder for an OpenAeroStruct component that computes a new wing mass.

"""
import openmdao.api as om
import aviary.api as av
from aviary.examples.external_subsystems.OAS_weight.OAS_wing_weight_analysis import OAStructures


class OASWingWeightBuilder(av.SubsystemBuilderBase):
    def __init__(self, name='wing_weight'):
        super().__init__(name)

    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.System
            An OpenMDAO system 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.
        '''

        wing_group = om.Group()
        wing_group.add_subsystem("aerostructures",
                                 OAStructures(
                                     symmetry=True,
                                     wing_weight_ratio=1.0,
                                     S_ref_type='projected',
                                     n_point_masses=1,
                                     num_twist_cp=4,
                                     num_box_cp=51),
                                 promotes_inputs=[
                                     ('fuel', av.Mission.Design.FUEL_MASS),
                                 ],
                                 promotes_outputs=[('wing_weight', av.Aircraft.Wing.MASS)])

        return wing_group
/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)

Analysis Model Details#

This analysis is based on the Aviary benchmark aircraft_for_bench_FwFm input data representing a typical large single aisle class transport aircraft. The analysis code OAS_wing_weight_analysis contains the OAStructures class which performs a structural analysis of the wing. The image below shows a simplified XDSM diagram of the pre-mission data flow in this example.

OAS XDSM

We’ll now discuss this code in more detail.

First, we create a mesh that defines the shape of the wing based on the span, the location of the wing break typical of a transport aircraft, the dihedral and the wing root, wing break, and wing tip chord lengths. The twist of the wing is defined along the span at a set of control points and must be present as it is used in the optimization problem.

We then use this mesh to define a simplified finite element model of the wing structural box. We also define the airfoil shapes as an input to OpenAeroStruct for a given wing thickness to chord ratio (t/c) to represent the wing box thickness. We then set initial values for the wing skin thickness and spar thickness are set, along with material properties and stress allowables for a metal material, typically aluminum. OpenAeroStruct will then calculate aeroelastic loads for a 2.5g maneuver condition and apply those loads to the finite element model of the wing structure.

Results of the structural optimization determine the optimum wing skin thickness, spar cap thickness, wing twist, wing t/c and maneuver angle of attack that satisfies strength constraints while minimizing wing weight. The ‘OAStructures’ class returns the optimized wing mass and the fuel mass burned but currently only the wing mass is used to override the Aviary variable ‘Aircraft.Wing.MASS’.

The OAS_wing_weight_analysis code may be executed in Python to test the OpenAeroStruct analysis outside of the Aviary subsystem interface. Default values for each of the inputs and options are included at the bottom of the analysis code file. This can be a useful test to demonstrate that the OpenAeroStruct analysis model has been properly defined and the model returns reasonable results.

Once the user is satisfied that the results are acceptable, the analysis tool can then be added as an external subsystem and tested in a mission.

Subsystem variables#

A variety of parameters may be defined for an OpenAeroStruct model. These allow the user to control how the aerodynamic and finite element meshes are subdivided, give details about the aerodynamic solution and provide structural material properties and structural scaling factors. The input variables passed in from Aviary may include the fuel mass, reserve fuel mass, airfoil description, engine mass and its location, lift and drag coefficients and the cruise conditions of Mach, altitude, thrust-specific fuel consumption (SFC) and range.

This is a list of the available options defined for the structural analysis:

  • symmetry

  • chord_cos_spacing

  • span_cos_spacing

  • num_box_cp

  • num_twist_cp

  • S_ref_type

  • fem_model_type

  • with_viscous

  • with_wave

  • k_lam

  • c_max_t

  • E

  • G

  • yield

  • mrho

  • strength_factor_for_upper_skin

  • wing_weight_ratio

  • exact_failure_constraint

  • struct_weight_relief

  • distributed_fuel_weight

  • n_point_masses

  • fuel_density

This is a list of the inputs defined for the structural analysis:

  • box_upper_x

  • box_lower_x

  • box_upper_y

  • box_lower_y

  • twist_cp

  • spar_thickness_cp

  • skin_thickness_cp

  • t_over_c_cp

  • airfoil_t_over_c

  • fuel

  • fuel_reserve

  • CL0

  • CD0

  • cruise_Mach

  • cruise_altitude

  • cruise_range

  • cruise_SFC

  • engine_mass

  • engine_location

The 2 outputs from the analysis tool are:

  • wing_weight

  • fuel_burn

See OAS_wing_weight_analysis and the OpenAeroStruct aerostructural optimization with wingbox documentation for descriptions of these variables.

Test Case#

A simple Aviary mission is defined to test the inclusion of the OpenAeroStruct wing weight subsystem during the pre-mission phase. The test mission is defined in run_simple_OAS_mission and is a mission based on input data read from the benchmark data file aircraft_for_bench_FwFm.

The OpenAeroStruct subsystem is used to compute an optimum wing mass which will override the Aviary computed wing mass value. The value of the Aviary variable Aircraft.Wing.MASS is printed at the conclusion of the mission to verify that the wing weight from the subsystem is overriding the Aviary computed wing weight.

A variable in the run_simple_OAS_mission file named use_OAS may be set by the user to True or False to run the simple mission with or without the OpenAeroStruct subsystem included. This will allow the mission to be flown either using the Aviary empirical wing weight estimation (use_OAS=False) or with the OpenAeroStruct subsystem (use_OAS=True).

Wing weight optimization of this type usually does not have knowledge of non-optimum wing mass values such as leading and training edge structure, actuators, stiffeners, etc. The optimum wing mass computed by the OAStructures class can be scaled using the option wing_weight_ratio to better match either the Aviary empirical wing weight value or a known fly-away weight estimate for your wing model. One method to determine the wing_weight_ratio would be to run the mission to calculate the Aviary empirical wing weight and then run OAS_wing_weight_analysis by itself using its default input values and compare wing weights. The wing_weight_ratio value may then be set to calibrate the OpenAeroStruct wing weight to the expected fly-away weight.

This calibration step has already been performed for this model, so the user may run the simple mission with or without the OpenAeroStruct subsystem active and compare the results.

Example Run Script#

Here is the full run script used to run the simple mission with the OpenAeroStruct subsystem active. This run script is also available in the run_simple_OAS_mission file.

Note

We do not actually perform the optimization below. Instead, we define and set up the model and the call to run the optimization is commented you. You can uncomment this and run the code block to perform a full optimization.

# %load ../../examples/external_subsystems/OAS_weight/run_simple_OAS_mission.py
'''
This is a simple test mission to demonstrate the inclusion of a
pre-mission user defined external subsystem.  The simple mission
is based on input data read from the benchmark data file bench_4.csv,
which represents a single-aisle commercial transport aircraft.  The
OpenAeroStruct (OAS) subsystem is used to compute an optimum wing
mass which will override the Aviary computed wing mass value.

The flag 'use_OAS' is set to 'True' to include the OAS subsystem in
the mission, or set to 'False' to run the mission without the
subsystem so that wing mass values between the 2 methods may be
compared.

'''

import numpy as np
import openmdao.api as om
import aviary.api as av
from aviary.examples.external_subsystems.OAS_weight.OAS_wing_weight_builder import OASWingWeightBuilder

# flag to turn on/off OpenAeroStruct subsystem for comparison testing
use_OAS = True

wing_weight_builder = OASWingWeightBuilder()

# Load the phase_info and other common setup tasks
phase_info = {
    'climb_1': {
        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},
        'user_options': {
            'optimize_mach': False,
            'optimize_altitude': False,
            'polynomial_control_order': 1,
            'num_segments': 5,
            'order': 3,
            'solve_for_distance': False,
            'initial_mach': (0.2, 'unitless'),
            'final_mach': (0.72, 'unitless'),
            'mach_bounds': ((0.18, 0.74), 'unitless'),
            'initial_altitude': (0.0, 'ft'),
            'final_altitude': (32000.0, 'ft'),
            'altitude_bounds': ((0.0, 34000.0), 'ft'),
            'throttle_enforcement': 'path_constraint',
            'fix_initial': True,
            'constrain_final': False,
            'fix_duration': False,
            'initial_bounds': ((0.0, 0.0), 'min'),
            'duration_bounds': ((64.0, 192.0), 'min'),
        },
        'initial_guesses': {'time': ([0, 128], 'min')},
    },
    'climb_2': {
        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},
        'user_options': {
            'optimize_mach': False,
            'optimize_altitude': False,
            'polynomial_control_order': 1,
            'num_segments': 5,
            'order': 3,
            'solve_for_distance': False,
            'initial_mach': (0.72, 'unitless'),
            'final_mach': (0.72, 'unitless'),
            'mach_bounds': ((0.7, 0.74), 'unitless'),
            'initial_altitude': (32000.0, 'ft'),
            'final_altitude': (34000.0, 'ft'),
            'altitude_bounds': ((23000.0, 38000.0), 'ft'),
            'throttle_enforcement': 'boundary_constraint',
            'fix_initial': False,
            'constrain_final': False,
            'fix_duration': False,
            'initial_bounds': ((64.0, 192.0), 'min'),
            'duration_bounds': ((56.5, 169.5), 'min'),
        },
        'initial_guesses': {'time': ([128, 113], 'min')},
    },
    'descent_1': {
        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},
        'user_options': {
            'optimize_mach': False,
            'optimize_altitude': False,
            'polynomial_control_order': 1,
            'num_segments': 5,
            'order': 3,
            'solve_for_distance': False,
            'initial_mach': (0.72, 'unitless'),
            'final_mach': (0.36, 'unitless'),
            'mach_bounds': ((0.34, 0.74), 'unitless'),
            'initial_altitude': (34000.0, 'ft'),
            'final_altitude': (500.0, 'ft'),
            'altitude_bounds': ((0.0, 38000.0), 'ft'),
            'throttle_enforcement': 'path_constraint',
            'fix_initial': False,
            'constrain_final': True,
            'fix_duration': False,
            'initial_bounds': ((120.5, 361.5), 'min'),
            'duration_bounds': ((29.0, 87.0), 'min'),
        },
        'initial_guesses': {'time': ([241, 58], 'min')},
    },
    'post_mission': {
        'include_landing': False,
        'constrain_range': True,
        'target_range': (1800., 'nmi'),
    },
}

phase_info['pre_mission'] = {'include_takeoff': False, 'optimize_mass': True}
if use_OAS:
    phase_info['pre_mission']['external_subsystems'] = [wing_weight_builder]

aircraft_definition_file = 'models/test_aircraft/aircraft_for_bench_FwFm.csv'
mission_method = 'simple'
mass_method = 'FLOPS'
make_plots = False
max_iter = 100
optimizer = 'IPOPT'

prob = av.AviaryProblem()

prob.load_inputs(aircraft_definition_file, phase_info)
prob.check_and_preprocess_inputs()
prob.add_pre_mission_systems()
prob.add_phases()
prob.add_post_mission_systems()
prob.link_phases()
prob.add_driver(optimizer=optimizer, max_iter=max_iter)
prob.add_design_variables()
prob.add_objective()
prob.setup()

if use_OAS:
    OAS_sys = 'pre_mission.wing_weight.aerostructures.'
    prob.set_val(OAS_sys + 'box_upper_x', np.array([0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32,
                 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6]), units='unitless')
    prob.set_val(OAS_sys + 'box_lower_x', np.array([0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32,
                 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6]), units='unitless')
    prob.set_val(OAS_sys + 'box_upper_y', np.array([0.0447,  0.046,  0.0472,  0.0484,  0.0495,  0.0505,  0.0514,  0.0523,  0.0531,  0.0538, 0.0545,  0.0551,  0.0557, 0.0563,  0.0568, 0.0573,  0.0577,  0.0581,  0.0585,  0.0588,  0.0591,  0.0593,  0.0595,  0.0597,
                 0.0599,  0.06,    0.0601,  0.0602,  0.0602,  0.0602,  0.0602,  0.0602,  0.0601,  0.06,    0.0599,  0.0598,  0.0596,  0.0594,  0.0592,  0.0589,  0.0586,  0.0583,  0.058,   0.0576,  0.0572,  0.0568,  0.0563,  0.0558,  0.0553,  0.0547,  0.0541]), units='unitless')
    prob.set_val(OAS_sys + 'box_lower_y', np.array([-0.0447, -0.046, -0.0473, -0.0485, -0.0496, -0.0506, -0.0515, -0.0524, -0.0532, -0.054, -0.0547, -0.0554, -0.056, -0.0565, -0.057, -0.0575, -0.0579, -0.0583, -0.0586, -0.0589, -0.0592, -0.0594, -0.0595, -0.0596, -
                 0.0597, -0.0598, -0.0598, -0.0598, -0.0598, -0.0597, -0.0596, -0.0594, -0.0592, -0.0589, -0.0586, -0.0582, -0.0578, -0.0573, -0.0567, -0.0561, -0.0554, -0.0546, -0.0538, -0.0529, -0.0519, -0.0509, -0.0497, -0.0485, -0.0472, -0.0458, -0.0444]), units='unitless')
    prob.set_val(OAS_sys + 'twist_cp', np.array([-6., -6., -4., 0.]), units='deg')
    prob.set_val(OAS_sys + 'spar_thickness_cp',
                 np.array([0.004, 0.005, 0.008, 0.01]), units='m')
    prob.set_val(OAS_sys + 'skin_thickness_cp',
                 np.array([0.005, 0.01, 0.015, 0.025]), units='m')
    prob.set_val(OAS_sys + 't_over_c_cp',
                 np.array([0.08, 0.08, 0.10, 0.08]), units='unitless')
    prob.set_val(OAS_sys + 'airfoil_t_over_c', 0.12, units='unitless')
    prob.set_val(OAS_sys + 'fuel', 40044.0, units='lbm')
    prob.set_val(OAS_sys + 'fuel_reserve', 3000.0, units='lbm')
    prob.set_val(OAS_sys + 'CD0', 0.0078, units='unitless')
    prob.set_val(OAS_sys + 'cruise_Mach', 0.785, units='unitless')
    prob.set_val(OAS_sys + 'cruise_altitude', 11303.682962301647, units='m')
    prob.set_val(OAS_sys + 'cruise_range', 3500, units='nmi')
    prob.set_val(OAS_sys + 'cruise_SFC', 0.53 / 3600, units='1/s')
    prob.set_val(OAS_sys + 'engine_mass', 7400, units='lbm')
    prob.set_val(OAS_sys + 'engine_location', np.array([25, -10.0, 0.0]), units='m')

prob.set_initial_guesses()
# prob.run_aviary_problem('dymos_solution.db', make_plots=False)
# print('wing mass = ', prob.model.get_val(av.Aircraft.Wing.MASS, units='lbm'))
/usr/share/miniconda/envs/test/lib/python3.12/site-packages/dymos/phase/phase.py:902: 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:2328: RuntimeWarning: Invalid options for non-optimal control 'mach' in phase 'climb_1': 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:2328: RuntimeWarning: Invalid options for non-optimal control 'altitude' in phase 'climb_1': 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:2328: RuntimeWarning: Invalid options for non-optimal control 'mach' in phase 'climb_2': 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:2328: RuntimeWarning: Invalid options for non-optimal control 'altitude' in phase 'climb_2': 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:2328: RuntimeWarning: Invalid options for non-optimal control 'mach' in phase 'descent_1': 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:2328: RuntimeWarning: Invalid options for non-optimal control 'altitude' in phase 'descent_1': lower, upper, ref
  warnings.warn(f"Invalid options for non-optimal control '{name}' in phase "
The following variables have been overridden:
  'aircraft:design:touchdown_mass
  'aircraft:engine:mass
  'aircraft:fins:mass
  'aircraft:fuel:auxiliary_fuel_capacity
  'aircraft:fuel:fuselage_fuel_capacity
  'aircraft:fuel:total_capacity
  'aircraft:fuselage:planform_area
  'aircraft:fuselage:wetted_area
  'aircraft:horizontal_tail:wetted_area
  'aircraft:landing_gear:main_gear_oleo_length
  'aircraft:landing_gear:nose_gear_oleo_length
  'aircraft:vertical_tail:wetted_area
  'aircraft:wing:aspect_ratio
  'aircraft:wing:control_surface_area
  'aircraft:wing:wetted_area

The following variables have been overridden by an external subsystem:
  'aircraft:wing:mass

--- Constraint Report [traj] ---
    --- climb_1 ---
        [path]    0.0000e+00 <= throttle <= 1.0000e+00  [unitless]
    --- climb_2 ---
        [initial] 0.0000e+00 <= throttle <= 1.0000e+00  [unitless]
        [final]   0.0000e+00 <= throttle <= 1.0000e+00  [unitless]
    --- descent_1 ---
        [path]    0.0000e+00 <= throttle <= 1.0000e+00  [unitless]