Merge pull request #162 from timryanb/tacs_oas_example

Adding aerostructural example w/ tacs oas coupling

examples/aerostructural/oas_tacs_wing/README.md

@@ -0,0 +1,7 @@
+This example performs an aerostructural optimization using a TACS structural analysis coupled to an OpenAerostruct aerodynamic analysis.
+The optimization model is based off a Boeing 777 class aircraft geometry (publicly available [here](https://hangar.openvsp.org/vspfiles/375)).
+The goal of the optimization is to minimize the wing structural weight under a 2.5G maneuver flight condition subject to an aggregated failure constraint.
+
+The optimization can be run from the `run_oas_tacs_wing.py` script.
+In order to run this script the user is required to [OpenVSP Python API](https://openvsp.org/api_docs/latest/) installed, 
+in addition to the standard MPhys libraries (TACS, OpenAeroStruct, FunToFEM, etc.).

examples/aerostructural/oas_tacs_wing/run_oas_tacs_wing.py

@@ -0,0 +1,124 @@
+import os
+
+import numpy as np
+import openmdao.api as om
+
+from mphys import Multipoint
+from mphys.scenario_aerostructural import ScenarioAeroStructural
+from openaerostruct.geometry.utils import generate_vsp_surfaces
+from openaerostruct.mphys import AeroBuilder
+from funtofem.mphys import MeldBuilder
+from tacs.mphys import TacsBuilder
+
+import tacs_setup
+
+class Top(Multipoint):
+    def setup(self):
+        # Input files
+        vsp_file = os.path.join(os.path.dirname(__file__), "twin_aisle_aircraft.vsp3")
+        bdf_file = os.path.join(os.path.dirname(__file__), "wingbox.bdf")
+        # Read the geometry.
+        # VSP components to include in VLM mesh
+        vsp_comps = ["Wing"]
+        # Generate half-wing mesh of wing
+        surfaces = generate_vsp_surfaces(vsp_file, symmetry=True, include=vsp_comps)
+
+        # Default surface dictionary values
+        default_dict = {
+            # Wing definition
+            "type": "aero",
+            # reflected across the plane y = 0
+            "S_ref_type": "wetted",  # how we compute the wing area,
+            # can be 'wetted' or 'projected'
+            # "twist_cp": np.zeros(3),  # Define twist using 3 B-spline cp's
+            "CL0": 0.0,  # CL of the surface at alpha=0
+            "CD0": 0.0,  # CD of the surface at alpha=0
+            # Airfoil properties for viscous drag calculation
+            "k_lam": 0.05,  # percentage of chord with laminar
+            # flow, used for viscous drag
+            "t_over_c": 0.12,  # thickness over chord ratio (NACA0015)
+            "c_max_t": 0.303,  # chordwise location of maximum (NACA0015)
+            # thickness
+            "with_viscous": False,  # if true, compute viscous drag,
+            "with_wave": False}
+
+        #  Update our surface dict with defaults
+        for surface in surfaces:
+            surface.update(default_dict)
+
+        mach = 0.6415
+        aoa = 1.0
+        rho = 1.23
+        yaw = 0.0
+        vel = 220.
+        re = 1e6
+
+        dvs = self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
+        dvs.add_output('aoa', val=aoa, units='deg')
+        dvs.add_output('yaw', val=yaw, units='deg')
+        dvs.add_output('rho', val=rho, units='kg/m**3')
+        dvs.add_output('mach', mach)
+        dvs.add_output('v', vel, units='m/s')
+        dvs.add_output('reynolds', re, units="1/m")
+
+        # OpenAeroStruct
+        aero_builder = AeroBuilder(surfaces)
+        aero_builder.initialize(self.comm)
+
+        self.add_subsystem('mesh_aero', aero_builder.get_mesh_coordinate_subsystem())
+
+        # TACS setup
+        struct_builder = TacsBuilder(mesh_file=bdf_file, element_callback=tacs_setup.element_callback,
+                                     problem_setup=tacs_setup.problem_setup,
+                                     constraint_setup=tacs_setup.constraint_setup, coupled=True)
+
+        struct_builder.initialize(self.comm)
+        ndv_struct = struct_builder.get_ndv()
+
+        self.add_subsystem('mesh_struct', struct_builder.get_mesh_coordinate_subsystem())
+        dvs.add_output('dv_struct', np.array(ndv_struct*[0.02]))
+
+        # MELD setup
+        ldxfer_builder = MeldBuilder(aero_builder, struct_builder, isym=1)
+        ldxfer_builder.initialize(self.comm)
+
+        # Scenario
+        self.mphys_add_scenario('maneuver', ScenarioAeroStructural(aero_builder=aero_builder,
+                                                                 struct_builder=struct_builder,
+                                                                 ldxfer_builder=ldxfer_builder))
+
+        for discipline in ['aero', 'struct']:
+            self.mphys_connect_scenario_coordinate_source(
+                'mesh_%s' % discipline, 'maneuver', discipline)
+
+        for dv in ['aoa', 'yaw', 'rho', 'mach', 'v', 'reynolds']:
+            self.connect(dv, f'maneuver.{dv}')
+        self.connect('dv_struct', 'maneuver.dv_struct')
+
+prob = om.Problem()
+prob.model = Top()
+
+model = prob.model
+
+# Add wingbox panel thicknesses and angle of attack as design variables
+model.add_design_var('dv_struct', lower=0.002, upper=0.2, scaler=1000.0)
+model.add_design_var('aoa', lower=-10, upper=10.0, scaler=0.1)
+# Structural mass of half wing
+model.add_objective('maneuver.mass', scaler=1.0/1000.0)
+# Max stress constraint
+model.add_constraint('maneuver.ks_vmfailure', upper=1.0, scaler=1.0)
+# Lift constraint
+model.add_constraint('maneuver.CL', equals=0.625, scaler=1.0)
+# Wingbox panel adjacency constraints
+model.add_constraint('maneuver.adjacency.LE_SPAR', lower=-2.5e-3, upper=2.5e-3, scaler=1e3, linear=True)
+model.add_constraint('maneuver.adjacency.TE_SPAR', lower=-2.5e-3, upper=2.5e-3, scaler=1e3, linear=True)
+model.add_constraint('maneuver.adjacency.UPPER_SKIN', lower=-2.5e-3, upper=2.5e-3, scaler=1e3, linear=True)
+model.add_constraint('maneuver.adjacency.LOWER_SKIN', lower=-2.5e-3, upper=2.5e-3, scaler=1e3, linear=True)
+
+prob.driver = om.ScipyOptimizeDriver(debug_print=['objs', 'nl_cons'], maxiter=200)
+prob.driver.options['optimizer'] = 'SLSQP'
+
+prob.setup()
+om.n2(prob, show_browser=False, outfile='tacs_oas_aerostruct.html')
+
+prob.run_driver()

examples/aerostructural/oas_tacs_wing/tacs_setup.py

@@ -0,0 +1,65 @@
+import numpy as np
+from tacs import elements, constitutive, functions
+
+# Callback function used to setup TACS element objects and DVs
+def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs):
+    rho = 2780.0  # density, kg/m^3
+    E = 73.1e9  # elastic modulus, Pa
+    nu = 0.33  # poisson's ratio
+    ys = 324.0e6  # yield stress, Pa
+    thickness = 0.003
+    min_thickness = 0.002
+    max_thickness = 0.05
+
+    # Setup (isotropic) property and constitutive objects
+    prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys)
+    # Set one thickness dv for every component
+    con = constitutive.IsoShellConstitutive(prop, t=thickness, tNum=dvNum, tlb=min_thickness, tub=max_thickness)
+
+    # For each element type in this component,
+    # pass back the appropriate tacs element object
+    transform = None
+
+    elemList = []
+    for elemDescript in elemDescripts:
+        if elemDescript in ["CQUAD4", "CQUADR"]:
+            elem = elements.Quad4Shell(transform, con)
+        elif elemDescript in ["CTRIA3", "CTRIAR"]:
+            elem = elements.Tri3Shell(transform, con)
+        else:
+            raise ValueError("Uh oh, '%s' not recognized" % (elemDescript))
+        elemList.append(elem)
+
+    return elemList
+
+def problem_setup(scenario_name, fea_assembler, problem):
+    """
+    Helper function to add fixed forces and eval functions
+    to structural problems used in tacs builder
+    """
+    # Add TACS Functions
+    # Only include mass from elements that belong to pytacs components (i.e. skip concentrated masses)
+    problem.addFunction('mass', functions.StructuralMass)
+    problem.addFunction('ks_vmfailure', functions.KSFailure, safetyFactor=1.0, ksWeight=50.0)
+
+    # Add 2.5G gravity load
+    g = np.array([0.0, 0.0, -9.81])  # m/s^2
+    problem.addInertialLoad(2.5*g)
+
+def constraint_setup(scenario_name, fea_assembler, constraint_list):
+    """
+    Helper function to setup tacs constraint classes
+    """
+    if scenario_name == "maneuver":
+        # Setup adjacency constraints for skin and spar panel thicknesses
+        constr = fea_assembler.createAdjacencyConstraint("adjacency")
+        compIDs = fea_assembler.selectCompIDs(include="UPPER_SKIN")
+        constr.addConstraint("UPPER_SKIN", compIDs=compIDs)
+        compIDs = fea_assembler.selectCompIDs(include="LOWER_SKIN")
+        constr.addConstraint("LOWER_SKIN", compIDs=compIDs)
+        compIDs = fea_assembler.selectCompIDs(include="LE_SPAR")
+        constr.addConstraint("LE_SPAR", compIDs=compIDs)
+        compIDs = fea_assembler.selectCompIDs(include="TE_SPAR")
+        constr.addConstraint("TE_SPAR", compIDs=compIDs)
+        constraint_list.append(constr)
+