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
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.
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:
derivs_method
use_jit
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'))
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.
The following variables have been overridden:
'aircraft:design:touchdown_mass 152800 lbm
'aircraft:engine:mass [7400.] lbm
'aircraft:fins:mass 0 lbm
'aircraft:fuel:auxiliary_fuel_capacity 0 lbm
'aircraft:fuel:fuselage_fuel_capacity 0 lbm
'aircraft:fuel:total_capacity 45694 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:landing_gear:main_gear_oleo_length 102 inch
'aircraft:landing_gear:nose_gear_oleo_length 67 inch
'aircraft:vertical_tail:wetted_area 581.13 ft**2
'aircraft:wing:aspect_ratio 11.22091 unitless
'aircraft:wing:control_surface_area 137 ft**2
'aircraft:wing:wetted_area 2396.56 ft**2
The following variables have been overridden by an external subsystem:
'aircraft:wing:mass
/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_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:2323: 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:2323: 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:2323: 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:2323: 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:2323: 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 "
--- 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]