Add aerostructural support for multi-section surfaces (#461)

* Add aerostructural support for multi-section

* Add support for struct design variables and fixed uni_t_over_c and uni_mesh output naming

* Added new test for aerostructural multi-section

* Update multisection aerostructural tests

* Fix code for latest refactor and clean up tests

* Add docs and tutorial

* Added support for rotational VLM to the aerostructural groups

* Update docs and tests

* Fix flake8

* Fix flake 8 again

Author: sabakhshi

openaerostruct/docs/advanced_features/multi_section_surfaces.rst

@@ -1,14 +1,12 @@
 .. _Custom_Mesh:
 
-Multi-section Surfaces
+Multi-section surfaces for aerodynamic problems
 ==========================
 
 OpenAeroStruct features the ability to specify surfaces as a series of sequentially connected sections. 
 Rather than controling the geometry of the surface as a whole, the optimizer can control the geometric parameters of each section individually using any of the geometric transformations available in OpenAeroStruct. 
 This feature was developed for the modular wing morphing applications but can be useful in situations where the user wishes to optimize a particular wing section while leaving the others fixed. 
 
-*Please note that aerostructural optimization with multi-section wings is currently not supported*
-
 This example script demonstrate the usage of the multi-section wing geometry features in OpenAeroStruct.
 We first start with the induced drag minimization of a simple two-section symmetrical wing.
 
@@ -125,4 +123,50 @@ The steps required to setup the multi-section geometry group by constraint featu
    :end-before: checkpoint 4
 
 
+Multi-section surfaces for aerostructural problems
+==========================
+
+Both multi-section geometery parameterization approaches in OpenAeroStruct are fully compatible with the aerostructural simulation capability through the multi-section aerostructural geometery class.
+Both the tube and wingbox models can be used however note that the structural model itself is not multi-section.
+The structural member geometery is not defined for each section but rather the entire surface as it normally is.
+The finite-element beam model is therefore coupled with the unified surface mesh which is where the aerodynamic simulation is done.
+
+
+We start by making the necessary imports for an aerostructural optimization in OpenAeroStruct.
+Note that we also import a special geometry group for multi-section aerostructural problems.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_AS.py
+   :start-after: checkpoint 0
+   :end-before: checkpoint 1
+
+Next, we define the surface dictionary as we did for the construction-based multi-section geometry parameterization but also include the structural parameters.
+
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_AS.py
+   :start-after: checkpoint 1
+   :end-before: checkpoint 2
+
+The independent variables and unified twist spline are setup in the same way as in the construction-based multi-section geometry parameterization.
 
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_AS.py
+   :start-after: checkpoint 2
+   :end-before: checkpoint 3
+
+Next, we add the multi-section aerostructural geometery and aerostructural analysis point groups and make the necessary connections.
+Note that this step is very similar to the single section aerostructural problem setup with only the geometery group being replaced.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_AS.py
+   :start-after: checkpoint 3
+   :end-before: checkpoint 4
+
+We can now setup our optimization problem, configure the optimizer and run the optimization.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_AS.py
+   :start-after: checkpoint 4
+   :end-before: checkpoint 5
+
+Lastly, we plot and print the results.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_AS.py
+   :start-after: checkpoint 5
+   :end-before: checkpoint 6

openaerostruct/docs/advanced_features/scripts/basic_2_sec_AS.py

@@ -0,0 +1,184 @@
+"""Optimizes the section twist distribution of a two section symmetrical wing using the construction-based approach for section
+joining and the aerostructural tube model. This example is referenced as part of the multi-section tutorial."""
+
+# docs checkpoint 0
+import numpy as np
+import openmdao.api as om
+from openaerostruct.integration.aerostruct_groups import MultiSecAerostructGeometry, AerostructPoint
+from openaerostruct.utils.constants import grav_constant
+from openaerostruct.geometry.utils import build_section_dicts, unify_mesh, build_multi_spline, connect_multi_spline
+
+# docs checkpoint 1
+
+# The geometry parameterization used here is identical to the one in the two section contruction based example. However,
+# instead of chord we apply the principle to the twist B-spline.
+
+# Set-up B-splines for each section. Done here since this information will be needed multiple times.
+sec_twist_cp = [np.zeros(2), np.zeros(2)]
+
+# Note the additional of structural variables to the surface dictionary
+surface = {
+    # Wing definition
+    # Basic surface parameters
+    "name": "surface",
+    "is_multi_section": True,
+    "num_sections": 2,  # The number of sections in the multi-section surface
+    "sec_name": ["sec0", "sec1"],  # names of the individual sections
+    "symmetry": True,  # if true, model one half of wing. reflected across the midspan of the root section
+    "S_ref_type": "wetted",  # how we compute the wing area, can be 'wetted' or 'projected'
+    "root_section": 1,
+    # Geometry Parameters
+    "taper": [1.0, 1.0],  # Wing taper for each section
+    "span": [10.0, 10.0],  # Wing span for each section
+    "sweep": [0.0, 0.0],  # Wing sweep for each section
+    "twist_cp": sec_twist_cp,
+    "t_over_c_cp": [np.array([0.15]), np.array([0.15])],  # thickness over chord ratio (NACA0015)
+    "root_chord": 5.0,  # Wing root chord for each section
+    # Mesh Parameters
+    "meshes": "gen-meshes",  # Supply a mesh for each section or "gen-meshes" for automatic mesh generation
+    "nx": 2,  # Number of chordwise points. Same for all sections
+    "ny": [3, 3],  # Number of spanwise points for each section
+    # Aerodynamic Parameters
+    "CL0": 0.0,  # CL of the surface at alpha=0
+    "CD0": 0.015,  # 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
+    "c_max_t": 0.303,  # chordwise location of maximum (NACA0015)
+    # thickness
+    "with_viscous": True,  # if true, compute viscous drag
+    "with_wave": False,  # if true, compute wave drag
+    "groundplane": False,
+    # Structural
+    "fem_model_type": "tube",
+    "thickness_cp": 0.1 * np.ones((2)),
+    "E": 70.0e9,  # [Pa] Young's modulus of the spar
+    "G": 30.0e9,  # [Pa] shear modulus of the spar
+    "yield": 500.0e6 / 2.5,  # [Pa] yield stress divided by 2.5 for limiting case
+    "mrho": 3.0e3,  # [kg/m^3] material density
+    "fem_origin": 0.35,  # normalized chordwise location of the spar
+    "wing_weight_ratio": 2.0,
+    "struct_weight_relief": False,  # True to add the weight of the structure to the loads on the structure
+    "distributed_fuel_weight": False,
+    # Constraints
+    "exact_failure_constraint": False,  # if false, use KS function
+}
+
+# docs checkpoint 2
+
+# Create the problem and assign the model group
+prob = om.Problem(reports=False)
+
+# Add problem information as an independent variables component
+indep_var_comp = om.IndepVarComp()
+indep_var_comp.add_output("v", val=248.136, units="m/s")
+indep_var_comp.add_output("alpha", val=9.0, units="deg")
+indep_var_comp.add_output("Mach_number", val=0.84)
+indep_var_comp.add_output("re", val=1.0e6, units="1/m")
+indep_var_comp.add_output("rho", val=0.38, units="kg/m**3")
+indep_var_comp.add_output("CT", val=grav_constant * 17.0e-6, units="1/s")
+indep_var_comp.add_output("R", val=11.165e6, units="m")
+indep_var_comp.add_output("W0", val=0.4 * 3e5, units="kg")
+indep_var_comp.add_output("speed_of_sound", val=295.4, units="m/s")
+indep_var_comp.add_output("load_factor", val=1.0)
+indep_var_comp.add_output("empty_cg", val=np.zeros((3)), units="m")
+
+prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"])
+
+
+# Generate the sections and unified mesh here in addition to adding the components.
+# This has to also be done here since AeroPoint has to know the unified mesh size.
+section_surfaces = build_section_dicts(surface)
+uniMesh = unify_mesh(section_surfaces)
+surface["mesh"] = uniMesh
+
+
+# Build a component with B-spline control points that joins the sections by construction
+twist_comp = build_multi_spline("twist_cp", len(section_surfaces), sec_twist_cp)
+prob.model.add_subsystem("twist_bspline", twist_comp)
+
+# Connect the B-spline component to the section B-splines
+connect_multi_spline(prob, section_surfaces, sec_twist_cp, "twist_cp", "twist_bspline", surface["name"])
+
+# docs checkpoint 3
+
+# Create and add a group that handles the geometry for the
+# aerostructual multi-section lifting surface
+multi_geom_group = MultiSecAerostructGeometry(surface=surface)
+prob.model.add_subsystem(surface["name"], multi_geom_group)
+
+name = surface["name"]
+
+point_name = "AS_point_0"
+
+# Create the aero point group and add it to the model
+AS_point = AerostructPoint(surfaces=[surface])
+
+prob.model.add_subsystem(
+    point_name,
+    AS_point,
+    promotes_inputs=[
+        "v",
+        "alpha",
+        "Mach_number",
+        "re",
+        "rho",
+        "CT",
+        "R",
+        "W0",
+        "speed_of_sound",
+        "empty_cg",
+        "load_factor",
+    ],
+)
+
+com_name = point_name + "." + name + "_perf"
+prob.model.connect(name + ".local_stiff_transformed", point_name + ".coupled." + name + ".local_stiff_transformed")
+prob.model.connect(name + ".nodes", point_name + ".coupled." + name + ".nodes")
+
+# Connect aerodyamic mesh to coupled group mesh
+prob.model.connect(name + ".mesh", point_name + ".coupled." + name + ".mesh")
+
+# Connect performance calculation variables
+prob.model.connect(name + ".radius", com_name + ".radius")
+prob.model.connect(name + ".thickness", com_name + ".thickness")
+prob.model.connect(name + ".nodes", com_name + ".nodes")
+prob.model.connect(name + ".cg_location", point_name + "." + "total_perf." + name + "_cg_location")
+prob.model.connect(name + ".structural_mass", point_name + "." + "total_perf." + name + "_structural_mass")
+prob.model.connect(name + ".t_over_c", com_name + ".t_over_c")
+
+
+# docs checkpoint 4
+
+prob.driver = om.ScipyOptimizeDriver()
+prob.driver.options["optimizer"] = "SLSQP"
+prob.driver.options["tol"] = 1e-6
+prob.driver.options["disp"] = True
+prob.driver.options["maxiter"] = 1000
+
+
+# Setup problem and add design variables, constraint, and objective
+prob.model.add_design_var("twist_bspline.twist_cp_spline", lower=-10.0, upper=15.0)
+prob.model.add_design_var("surface.thickness_cp", lower=0.01, upper=0.5, scaler=1e2)
+prob.model.add_constraint("AS_point_0.surface_perf.failure", upper=0.0)
+prob.model.add_constraint("AS_point_0.surface_perf.thickness_intersects", upper=0.0)
+
+# Add design variables, constraints, and objective on the problem
+prob.model.add_design_var("alpha", lower=-10.0, upper=10.0)
+prob.model.add_constraint("AS_point_0.L_equals_W", equals=0.0)
+prob.model.add_objective("AS_point_0.fuelburn", scaler=1e-5)
+
+
+# Set up the problem
+prob.setup(check=True)
+
+# Run the optimization
+optResult = prob.run_driver()
+# om.n2(prob, show_browser=False)
+
+
+# docs checkpoint 5
+
+# Print fuelburn
+print(prob.get_val("AS_point_0.fuelburn"))
+# docs checkpoint 6

openaerostruct/docs/advanced_features/scripts/basic_2_sec_visc.py

@@ -167,6 +167,7 @@
 meshUni = prob.get_val(name + "." + unification_name + "." + name + "_uni_mesh")
 
 
+# Plot the results
 def plot_meshes(meshes):
     """this function plots to plot the mesh"""
     plt.figure(figsize=(8, 4))