diff --git a/openaerostruct/__init__.py b/openaerostruct/__init__.py index 892994aa6..43ce13db0 100644 --- a/openaerostruct/__init__.py +++ b/openaerostruct/__init__.py @@ -1 +1 @@ -__version__ = "2.8.0" +__version__ = "2.9.0" diff --git a/openaerostruct/aerodynamics/aero_groups.py b/openaerostruct/aerodynamics/aero_groups.py index 961e32fa2..0b5981d24 100644 --- a/openaerostruct/aerodynamics/aero_groups.py +++ b/openaerostruct/aerodynamics/aero_groups.py @@ -30,6 +30,37 @@ def setup(self): surfaces = self.options["surfaces"] rotational = self.options["rotational"] + # Check for multi-section surfaces and create suitable surface dictionaries for them + for i, surface in enumerate(surfaces): + # If multisection mesh then build a single surface with the unified mesh data + if "is_multi_section" in surface.keys(): + import copy + + target_keys = [ + # Essential Info + "name", + "symmetry", + "S_ref_type", + "ref_axis_pos", + "mesh", + # aerodynamics + "CL0", + "CD0", + "with_viscous", + "with_wave", + "groundplane", + "k_lam", + "t_over_c_cp", + "c_max_t", + ] + + # Constructs a surface dictionary and adds the specified supported keys and values from the mult-section surface dictionary. + aeroSurface = {} + for k in set(surface).intersection(target_keys): + aeroSurface[k] = surface[k] + # print(aeroSurface["name"]) + surfaces[i] = copy.deepcopy(aeroSurface) + # Loop through each surface and connect relevant parameters for surface in surfaces: name = surface["name"] diff --git a/openaerostruct/docs/advanced_features.rst b/openaerostruct/docs/advanced_features.rst index 2beee7262..1cda14b93 100644 --- a/openaerostruct/docs/advanced_features.rst +++ b/openaerostruct/docs/advanced_features.rst @@ -12,6 +12,7 @@ These examples show some advanced features in OpenAeroStruct. advanced_features/geometry_manipulation.rst advanced_features/custom_mesh_example.rst advanced_features/openvsp_mesh_example.rst + advanced_features/multi_section_surfaces.rst advanced_features/multiple_surfaces.rst advanced_features/multipoint.rst advanced_features/multipoint_parallel.rst diff --git a/openaerostruct/docs/advanced_features/figs/multi_section_2_sym.png b/openaerostruct/docs/advanced_features/figs/multi_section_2_sym.png new file mode 100644 index 000000000..1fb988451 Binary files /dev/null and b/openaerostruct/docs/advanced_features/figs/multi_section_2_sym.png differ diff --git a/openaerostruct/docs/advanced_features/multi_section_surfaces.rst b/openaerostruct/docs/advanced_features/multi_section_surfaces.rst new file mode 100644 index 000000000..f21d88d90 --- /dev/null +++ b/openaerostruct/docs/advanced_features/multi_section_surfaces.rst @@ -0,0 +1,107 @@ +.. _Custom_Mesh: + +Multi-section Surfaces +========================== + +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. + +Let's start with the necessary imports. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 0 + :end-before: checkpoint 1 + + +Then we are ready to discuss the multi-section parameterization and multi-section surface dictionary. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 1 + :end-before: checkpoint 2 + +Next, we setup the flow variables. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 2 + :end-before: checkpoint 3 + + +Giving the optimizer control of each wing section without any measure to constrain its motion can lead to undesirable geometries, such as wing sections separating from each other or large kinks appearing along the span, that cause numerical deficiencies in the optimization. +This concern is addressed by enforcing C0 continuity along the span-wise junctions between the sections. +C0 continuity can be enforced for any geometric parameter OAS controls for a given section. +For example, the tip chord or twist of a given section should match the root chord or twist of the subsequent section. +There are two approaches for enforcing C0 continuity between sections: a constraint-based approach and a construction-based approach. + +The constaint-based approach involves explicitly setting a position constraint that joins the surface points at section junctions to a specified tolerance. +This constraint is useful when varying section parameters like span or dihedral, where entire sections could collide or separate if C0 continuity is not strictly enforced. +A fully differentiated implementation that facilitates setting this constraint is incorporated into OAS. This approach is robust but introduces at least two linear constraints per segment junction. +The number of linear constraints can quickly grow for problems with many sections that this study anticipates for morphing applications. +In cases where geometric parameters can be specified as a function of span, however, it is possible to eliminate these additional constraints and maintain C0 continuity using the construction-based approach. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 3 + :end-before: checkpoint 4 + + +The construction-based approach forgoes the constraint entirely. +Continuity by construction is enforced by assigning the B-spline control point located at section edges to the same independent variable controlled by the optimizer. +Enforcing C0 continuity by construction solely applies to geometric parameters that employ B-spline parametrization in OAS. +These include chord distribution, twist distribution, shear distribution in all three directions, and thickness and radius distributions for structural spars. +Geometric transformations parameterized by a single value, such as sweep, taper, and dihedral, cannot be used with C0 continuity enforcement by construction. + + +.. literalinclude:: /advanced_features/scripts/basic_2_sec_construction.py + :start-after: checkpoint 0 + :end-before: checkpoint 1 + +.. literalinclude:: /advanced_features/scripts/basic_2_sec_construction.py + :start-after: checkpoint 2 + :end-before: checkpoint 3 + +We can now create the aerodynamic analysis group. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 4 + :end-before: checkpoint 5 + +Connecting the geometry and analysis groups requires care when using the multi-section parameterization. +While the multi-section geometry group is effectively a drop-in replacement for the standard geometery group for multi-section wings, it does require that the unified mesh component be connected to the AeroPoint analysis. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 5 + :end-before: checkpoint 6 + +*The following section applies to user considering cases with thickess to chord B-splines and viscous effects.* + +When using a thickness to chord ratio B-spline to account for viscous effects the user should be careful to connect the unified thickess to chord ratio B-spline automatically generated by the multi-section geometery group. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec_visc.py + :start-after: checkpoint 0 + :end-before: checkpoint 1 + +We can now setup our optimization problem and run it. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 6 + :end-before: checkpoint 7 + +We then finish by plotting the result. + +.. literalinclude:: /advanced_features/scripts/basic_2_sec.py + :start-after: checkpoint 7 + :end-before: checkpoint 8 + + + + + +The following shows the resulting optimized mesh. +The result is identical regardless of if the constraint-based or construction-based joining approahces are used. + +.. image:: /advanced_features/figs/multi_section_2_sym.png diff --git a/openaerostruct/docs/advanced_features/scripts/basic_2_sec.py b/openaerostruct/docs/advanced_features/scripts/basic_2_sec.py new file mode 100644 index 000000000..cc1a063ea --- /dev/null +++ b/openaerostruct/docs/advanced_features/scripts/basic_2_sec.py @@ -0,0 +1,241 @@ +"""Optimizes the section chord distribution of a two section symmetrical wing using the constraint-based approach for section +joining. 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.geometry.geometry_group import MultiSecGeometry +from openaerostruct.aerodynamics.aero_groups import AeroPoint +from openaerostruct.geometry.geometry_group import build_sections +from openaerostruct.geometry.geometry_unification import unify_mesh +import matplotlib.pyplot as plt + + +# docs checkpoint 1 + +# The multi-section geometry parameterization number section from left to right starting with section #0. A two-section symmetric wing parameterization appears as follows. +# For a symmetrical wing the last section in the sequence will always be marked as the "root section" as it's adjacent to the geometric centerline of the wing. +# Geometeric parameters must be specified for each section using lists with values corresponding in order of the surface numbering. Section section supports all the +# standard OpenAeroStruct geometery transformations including B-splines. + + +""" + +----------------------------------------------- ^ +| | | | +| | | | +| sec 0 | sec 1 | | root symmetrical BC +| | "root section" | | chord +|______________________|_______________________| | + _ + y = 0 ------------------> + y + +""" + + +# A multi-section surface dictionary is very similar to the standard one. However, it features some additional options and requires that the user specify +# parameters for each desired section. The multi-section geometery group also features an built in mesh generator so the wing mesh parameters can be specified right +# in the surface dictionary. Let's create a dictionary with info and options for a two-section aerodynamic lifting surface + +surface = { + # Wing definition + # Basic surface parameters + "name": "surface", + "is_multi_section": True, # This key must be present for the AeroPoint to correctly interpret this surface as multi-section + "num_sections": 2, # The number of sections in the multi-section surface + "sec_name": [ + "sec0", + "sec1", + ], # names of the individual sections. Each section must be named and the list length must match the specified number of 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' + # Geometry Parameters + "taper": [1.0, 1.0], # Wing taper for each section. The list length must match the specified number of sections. + "span": [2.0, 2.0], # Wing span for each section. The list length must match the specified number of sections. + "sweep": [0.0, 0.0], # Wing sweep for each section. The list length must match the specified number of sections. + "chord_cp": [ + np.array([1, 1]), + np.array([1, 1]), + ], # The chord B-spline parameterization for EACH SECTION. The list length must match the specified number of sections. + "twist_cp": [ + np.zeros(2), + np.zeros(2), + ], # The twist B-spline parameterization for EACH SECTION. The list length must match the specified number of sections. + "root_chord": 1.0, # Root chord length of the section indicated as "root section"(required if using the built-in mesh generator) + # Mesh Parameters + "meshes": "gen-meshes", # Supply a list of meshes for each section or "gen-meshes" for automatic mesh generation + "nx": 2, # Number of chordwise points. Same for all sections.(required if using the built-in mesh generator) + "ny": [ + 21, + 21, + ], # Number of spanwise points for each section. The list length must match the specified number of sections. (required if using the built-in mesh generator) + # 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": False, # if true, compute viscous drag + "with_wave": False, # if true, compute wave drag + "groundplane": False, +} + +# docs checkpoint 2 + +# Create the OpenMDAO problem +prob = om.Problem() + +# Create an independent variable component that will supply the flow +# conditions to the problem. +indep_var_comp = om.IndepVarComp() +indep_var_comp.add_output("v", val=1.0, units="m/s") +indep_var_comp.add_output("alpha", val=10.0, units="deg") +indep_var_comp.add_output("Mach_number", val=0.3) +indep_var_comp.add_output("re", val=1.0e5, units="1/m") +indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") +indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + +# Add this IndepVarComp to the problem model +prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + +# docs checkpoint 3 + +# Instead of creating a standard geometery group, here we will create a multi-section group that will accept our multi-section surface +# dictionary and allow us to specify any C0 continuity constraints between the sections. In this example we will constrain the sections +# into a C0 continuous surface using a component that the optimizer can use as a constraint. The joining constraint component returns the +# distance between the leading edge and trailing edge points at section interections. Any combination of the x,y, and z distances can be returned +# to constrain the surface in a particular direction. + + +""" + LE1 LE2 cLE = [LE2x-LE1x,LE2y-LE1y,LE2z-LE1z] +------------------------- <-------------> ------------------------- +| | | | +| | | | +| sec 0 | | sec 1 | +| | | "root section" | +|______________________ | <-------------> |_______________________| + TE1 TE2 cTE = [TE2x-TE1x,TE2y-TE1y,TE2z-TE1z] + + + +""" + +# We pass in the multi-section surface dictionary to the MultiSecGeometry geometery group. We also enabled joining_comp and pass two array to dim_contr +# These two arrays should only consists of 1 and 0 and tell the joining component which of the x,y, and z distance constraints we wish to enforce at the LE and TE +# In this example, we only wish to constraint the x-distance between the sections at both the leading and trailing edge. + +multi_geom_group = MultiSecGeometry( + surface=surface, joining_comp=True, dim_constr=[np.array([1, 0, 0]), np.array([1, 0, 0])] +) +prob.model.add_subsystem(surface["name"], multi_geom_group) + +# docs checkpoint 4 + +# In this next part, we will setup the aerodynamics group. First we use a utility function called build_sections which takes our multi-section surface dictionary and outputs a +# surface dictionary for each individual section. We then inputs these dictionaries into the mesh unification function unify_mesh to produce a single mesh array for the the entire surface. +# We then add this mesh to the multi-section surface dictionary +section_surfaces = build_sections(surface) +uniMesh = unify_mesh(section_surfaces) +surface["mesh"] = uniMesh + +# Create the aero point group, which contains the actual aerodynamic +# analyses. This step is exactly as it's normally done except the surface dictionary we pass in is the multi-surface one +aero_group = AeroPoint(surfaces=[surface]) +point_name = "aero_point_0" +prob.model.add_subsystem(point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"]) + +# docs checkpoint 5 + +# The following steps are similar to a normal OAS surface script but note the differences in surface naming. Note that +# unified surface created by the multi-section geometry group needs to be connected to AeroPoint(be careful with the naming) + +# Get name of surface and construct the name of the unified surface mesh +name = surface["name"] +unification_name = "{}_unification".format(surface["name"]) + +# Connect the mesh from the mesh unification component to the analysis point. +prob.model.connect(name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh") + +# Perform the connections with the modified names within the +# 'aero_states' group. +prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + ".aero_states." + "surface" + "_def_mesh" +) + +# docs checkpoint 6 + +# Next, we add the DVs to the OpenMDAO problem. Note that each surface's geometeric parameters are under the given section names specified in the multi-surface dictionary earlier. +# Here we use the chord B-spline that we specified earlier for each section and the angle-of-attack as DVs. +prob.model.add_design_var("surface.sec0.chord_cp", lower=0.1, upper=10.0, units=None) +prob.model.add_design_var("surface.sec1.chord_cp", lower=0.1, upper=10.0, units=None) +prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + +# Next, we add the C0 continuity constraint for this problem by constraining the x-distance between sections to 0. +# NOTE: SLSQP optimizer does not handle the joining equality constraint properly so the constraint needs to be specified as an inequality constraint +# All other optimizers like SNOPT can handle the equality constraint as is. + +prob.model.add_constraint("surface.surface_joining.section_separation", upper=0, lower=0) # FOR SLSQP +# prob.model.add_constraint('surface.surface_joining.section_separation',equals=0.0,scaler=1e-4) #FOR OTHER OPTIMIZERS + +# Add CL constraint +prob.model.add_constraint(point_name + ".CL", equals=0.3) + +# Add Wing total area constraint +prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + +# Add objective +prob.model.add_objective(point_name + ".CD", scaler=1e4) + + +prob.driver = om.ScipyOptimizeDriver() +prob.driver.options["optimizer"] = "SLSQP" +prob.driver.options["tol"] = 1e-3 +prob.driver.options["disp"] = True +prob.driver.options["maxiter"] = 1000 +prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + +# Set up and run the optimization problem +prob.setup() + +# prob.run_model() +prob.run_driver() +# om.n2(prob) + +# docs checkpoint 7 + +# Get each section mesh +mesh1 = prob.get_val("surface.sec0.mesh", units="m") +mesh2 = prob.get_val("surface.sec1.mesh", units="m") + +# Get the unified mesh +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)) + for i, mesh in enumerate(meshes): + mesh_x = mesh[:, :, 0] + mesh_y = mesh[:, :, 1] + color = "k" + for i in range(mesh_x.shape[0]): + plt.plot(mesh_y[i, :], 1 - mesh_x[i, :], color, lw=1) + plt.plot(-mesh_y[i, :], 1 - mesh_x[i, :], color, lw=1) # plots the other side of symmetric wing + for j in range(mesh_x.shape[1]): + plt.plot(mesh_y[:, j], 1 - mesh_x[:, j], color, lw=1) + plt.plot(-mesh_y[:, j], 1 - mesh_x[:, j], color, lw=1) # plots the other side of symmetric wing + plt.axis("equal") + plt.xlabel("y (m)") + plt.ylabel("x (m)") + plt.savefig("opt_planform_constraint.png") + + +# plot_meshes([mesh1,mesh2]) +plot_meshes([meshUni]) +# plt.show() +# docs checkpoint 8 diff --git a/openaerostruct/docs/advanced_features/scripts/basic_2_sec_construction.py b/openaerostruct/docs/advanced_features/scripts/basic_2_sec_construction.py new file mode 100644 index 000000000..b7f42a220 --- /dev/null +++ b/openaerostruct/docs/advanced_features/scripts/basic_2_sec_construction.py @@ -0,0 +1,215 @@ +"""Optimizes the section chord distribution of a two section symmetrical wing using the construction-based approach for section +joining. This example is referenced as part of the multi-section tutorial.""" + +import numpy as np + +import openmdao.api as om + +from openaerostruct.geometry.geometry_group import MultiSecGeometry +from openaerostruct.aerodynamics.aero_groups import AeroPoint +from openaerostruct.geometry.geometry_group import build_sections +from openaerostruct.geometry.geometry_unification import unify_mesh +import matplotlib.pyplot as plt + +# docs checkpoint 0 + +# To use the construction based approach an additional import is required. +from openaerostruct.geometry.multi_unified_bspline_utils import build_multi_spline, connect_multi_spline + +# docs checkpoint 1 + +# Set-up B-splines for each section. Done here since this information will be needed multiple times. +sec_chord_cp = [np.array([1, 1]), np.array([1.0, 0.2])] + + +# Create a dictionary with info and options about the multi-section aerodynamic +# lifting surface +surface = { + # Wing definition + # Basic surface parameters + "name": "surface", + "is_multi_section": True, # This key must be present for the AeroPoint to correctly interpret this surface as multi-section + "num_sections": 2, # The number of sections in the multi-section surface + "sec_name": [ + "sec0", + "sec1", + ], # names of the individual sections. Each section must be named and the list length must match the specified number of 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' + # Geometry Parameters + "taper": [1.0, 1.0], # Wing taper for each section. The list length must match the specified number of sections. + "span": [2.0, 2.0], # Wing span for each section. The list length must match the specified number of sections. + "sweep": [0.0, 0.0], # Wing sweep for each section. The list length must match the specified number of sections. + "chord_cp": [ + np.array([1, 1]), + np.array([1, 1]), + ], # The chord B-spline parameterization for EACH SECTION. The list length must match the specified number of sections. + "twist_cp": [ + np.zeros(2), + np.zeros(2), + ], # The twist B-spline parameterization for EACH SECTION. The list length must match the specified number of sections. + "root_chord": 1.0, # Root chord length of the section indicated as "root section"(required if using the built-in mesh generator) + # Mesh Parameters + "meshes": "gen-meshes", # Supply a list of meshes for each section or "gen-meshes" for automatic mesh generation + "nx": 2, # Number of chordwise points. Same for all sections.(required if using the built-in mesh generator) + "ny": [ + 21, + 21, + ], # Number of spanwise points for each section. The list length must match the specified number of sections. (required if using the built-in mesh generator) + # 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)# docs checkpoint 1 + # thickness + "with_viscous": False, # if true, compute viscous drag + "with_wave": False, # if true, compute wave drag + "groundplane": False, +} + +# Create the OpenMDAO problem +prob = om.Problem() + +# Create an independent variable component that will supply the flow +# conditions to the problem. +indep_var_comp = om.IndepVarComp() +indep_var_comp.add_output("v", val=1.0, units="m/s") +indep_var_comp.add_output("alpha", val=10.0, units="deg") +indep_var_comp.add_output("Mach_number", val=0.3) +indep_var_comp.add_output("re", val=1.0e5, units="1/m") +indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") +indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + +# Add this IndepVarComp to the problem model +prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + +# Generate the sections and unified mesh here. It's needed to join the sections by construction. +section_surfaces = build_sections(surface) +uniMesh = unify_mesh(section_surfaces) +surface["mesh"] = uniMesh + +# docs checkpoint 2 + +"""This functions builds an OpenMDAO B-spline component for the surface with the correct number of control points +corresponding to each section junction on the surface. Refer to the functions documentions for input details. After +the compnent has been generated it needs to be added to the model.""" +chord_comp = build_multi_spline("chord_cp", len(section_surfaces), sec_chord_cp) +prob.model.add_subsystem("chord_bspline", chord_comp) + +"""In order to properly transform the surface geometry the surface B-spline's control points need to be connected +to the corresponding control points of the local B-spline component on each section. This function automates this +process as it can be tedious. + +The figure below explains how the surface B-spline's control points are connected to the control point of the local +section B-spline. In this example, each section features a two point B-spline with control points at the section tips +however the principle is the same for B-splines with more points. + + + surface B-spline +0;;;;;;;;;;;;;;;;;;;;;;1;;;;;;;;;;;;;;;;;;;;;;;2 +^ ^ ^ +| | | +| | | +| | sec 1 B-spline | + sec 0 B-spline c:::::::::::::::::::::::d +a::::::::::::::::::::::b +----------------------------------------------- ^ +| | | | +| | | | +| sec 0 | sec 1 | | root +| | | | chord +|______________________|_______________________| | + _ + y = 0 ------------------> + y + + +An edge case in this process is when a section features a B-spline with a single control point. The same control point +cannot be assigned to two different control points on the surface B-spline. In these situations a constraint will need +to be used to maintain C0 continuity. See the connect_multi_spline documentation for details. +""" +connect_multi_spline(prob, section_surfaces, sec_chord_cp, "chord_cp", "chord_bspline") + + +""" With the surface B-spline connected we can add the multi-section geometry group. Note that in this case the joining +component does not need to be specified as we are not joining the sections by constraint.""" +multi_geom_group = MultiSecGeometry(surface=surface) +prob.model.add_subsystem(surface["name"], multi_geom_group) + +# docs checkpoint 3 + +# Create the aero point group, which contains the actual aerodynamic +# analyses +aero_group = AeroPoint(surfaces=[surface]) +point_name = "aero_point_0" +prob.model.add_subsystem(point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"]) + +# Get name of surface and construct unified mesh name +name = surface["name"] +unification_name = "{}_unification".format(surface["name"]) + +# Connect the mesh from the mesh unification component to the analysis point +prob.model.connect(name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh") + +# Perform the connections with the modified names within the +# 'aero_states' group. +prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + ".aero_states." + "surface" + "_def_mesh" +) + +# Add DVs +prob.model.add_design_var("chord_bspline.chord_cp_spline", lower=0.1, upper=10.0, units=None) +prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + +# Add CL constraint +prob.model.add_constraint(point_name + ".CL", equals=0.3) + +# Add Area constraint +prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + +# Add objective +prob.model.add_objective(point_name + ".CD", scaler=1e4) + +prob.driver = om.ScipyOptimizeDriver() +prob.driver.options["optimizer"] = "SLSQP" +prob.driver.options["tol"] = 1e-7 +prob.driver.options["disp"] = True +prob.driver.options["maxiter"] = 1000 +# prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + +# Set up and run the optimization problem +prob.setup() +prob.run_driver() +# om.n2(prob) + + +mesh1 = prob.get_val("surface.sec0.mesh", units="m") +mesh2 = prob.get_val("surface.sec1.mesh", units="m") + +meshUni = prob.get_val(name + "." + unification_name + "." + name + "_uni_mesh") + + +def plot_meshes(meshes): + """this function plots to plot the mesh""" + plt.figure(figsize=(8, 4)) + for i, mesh in enumerate(meshes): + mesh_x = mesh[:, :, 0] + mesh_y = mesh[:, :, 1] + color = "w" + for i in range(mesh_x.shape[0]): + plt.plot(mesh_y[i, :], 1 - mesh_x[i, :], color, lw=1) + plt.plot(-mesh_y[i, :], 1 - mesh_x[i, :], color, lw=1) # plots the other side of symmetric wing + for j in range(mesh_x.shape[1]): + plt.plot(mesh_y[:, j], 1 - mesh_x[:, j], color, lw=1) + plt.plot(-mesh_y[:, j], 1 - mesh_x[:, j], color, lw=1) # plots the other side of symmetric wing + plt.axis("equal") + plt.xlabel("y (m)") + plt.ylabel("x (m)") + plt.savefig("opt_planform_construction.pdf") + + +plot_meshes([meshUni]) +# plt.show() diff --git a/openaerostruct/docs/advanced_features/scripts/basic_2_sec_visc.py b/openaerostruct/docs/advanced_features/scripts/basic_2_sec_visc.py new file mode 100644 index 000000000..119ee8838 --- /dev/null +++ b/openaerostruct/docs/advanced_features/scripts/basic_2_sec_visc.py @@ -0,0 +1,190 @@ +"""Optimizes the section chord distribution of a two section symmetrical wing with thickenss and viscous effects accounted +for. This examples is similar to the inviscid case but explains how to connect the unified t/c B-spline to the OAS performance +component for correct viscous drag computation. This example is referenced as part of the multi-section tutorial.""" + +import numpy as np + +import openmdao.api as om + +from openaerostruct.geometry.geometry_group import MultiSecGeometry +from openaerostruct.aerodynamics.aero_groups import AeroPoint +from openaerostruct.geometry.geometry_group import build_sections +from openaerostruct.geometry.geometry_unification import unify_mesh +from openaerostruct.geometry.multi_unified_bspline_utils import build_multi_spline, connect_multi_spline +import matplotlib.pyplot as plt + + +# Set-up B-splines for each section. Done here since this information will be needed multiple times. +sec_chord_cp = [np.array([1.0, 1.0]), np.array([1.0, 1.0])] + + +# Create a dictionary with info and options about the multi-section aerodynamic +# lifting surface +surface = { + # Wing definition + # Basic surface parameters + "name": "surface", + "is_multi_section": True, # This key must be present for the AeroPoint to correctly interpret this surface as multi-section + "num_sections": 2, # The number of sections in the multi-section surface + "sec_name": [ + "sec0", + "sec1", + ], # names of the individual sections. Each section must be named and the list length must match the specified number of 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' + # Geometry Parameters + "taper": [1.0, 1.0], # Wing taper for each section. The list length must match the specified number of sections. + "span": [2.0, 2.0], # Wing span for each section. The list length must match the specified number of sections. + "sweep": [0.0, 0.0], # Wing sweep for each section. The list length must match the specified number of sections. + "chord_cp": [ + np.array([1, 1]), + np.array([1, 1]), + ], # The chord B-spline parameterization for EACH SECTION. The list length must match the specified number of sections. + "twist_cp": [ + np.zeros(2), + np.zeros(2), + ], # The twist B-spline parameterization for EACH SECTION. The list length must match the specified number of sections. + "t_over_c_cp": [ + np.array([0.15]), + np.array([0.15]), + ], # The thickenss to chord ratio B-spline parameterization for EACH SECTION. The list length must match the specified number of sections. + "root_chord": 1.0, # Root chord length of the section indicated as "root section"(required if using the built-in mesh generator) + # Mesh Parameters + "meshes": "gen-meshes", # Supply a list of meshes for each section or "gen-meshes" for automatic mesh generation + "nx": 2, # Number of chordwise points. Same for all sections.(required if using the built-in mesh generator) + "ny": [ + 21, + 21, + ], # Number of spanwise points for each section. The list length must match the specified number of sections. (required if using the built-in mesh generator) + # 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, +} + +# Create the OpenMDAO problem +prob = om.Problem() + +# Create an independent variable component that will supply the flow +# conditions to the problem. +indep_var_comp = om.IndepVarComp() +indep_var_comp.add_output("v", val=1.0, units="m/s") +indep_var_comp.add_output("alpha", val=10.0, units="deg") +indep_var_comp.add_output("Mach_number", val=0.3) +indep_var_comp.add_output("re", val=1.0e5, units="1/m") +indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") +indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + +# Add this IndepVarComp to the problem model +prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + +# Generate the sections and unified mesh here. It's needed to join the sections by construction. +section_surfaces = build_sections(surface) +uniMesh = unify_mesh(section_surfaces) +surface["mesh"] = uniMesh + +# Build a component with B-spline control points that joins the sections by construction +chord_comp = build_multi_spline("chord_cp", len(section_surfaces), sec_chord_cp) +prob.model.add_subsystem("chord_bspline", chord_comp) + +# Connect the B-spline component to the section B-splines +connect_multi_spline(prob, section_surfaces, sec_chord_cp, "chord_cp", "chord_bspline") + +# Create and add a group that handles the geometry for the +# aerodynamic lifting surface +multi_geom_group = MultiSecGeometry(surface=surface) +prob.model.add_subsystem(surface["name"], multi_geom_group) + + +# Create the aero point group, which contains the actual aerodynamic +# analyses +aero_group = AeroPoint(surfaces=[surface]) +point_name = "aero_point_0" +prob.model.add_subsystem(point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"]) + +# Get name of surface and construct unified mesh name +name = surface["name"] +unification_name = "{}_unification".format(surface["name"]) + +# Connect the mesh from the mesh unification component to the analysis point +prob.model.connect(name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh") + +# Perform the connections with the modified names within the +# 'aero_states' group. +prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + ".aero_states." + "surface" + "_def_mesh" +) + + +# docs checkpoint 0 + +"""If the surface features a thickness to chord ratio B-spline, then the Muli-section geometry component will automatically +combine them together while enforcing C0 continuity. This process will occurs regardless of whether the constraint or construction-based +section joining method are used as otherwise OAS will not be able to calculate the viscous drag correctly. Pay close attention to how the +unified t_over_c B-spline(given a unique name by the multi-section geometry group) needs to be connect to the AeroPoint performance component.""" +prob.model.connect( + name + "." + unification_name + "." + name + "_uni_t_over_c", point_name + "." + name + "_perf." + "t_over_c" +) + +# docs checkpoint 1 + +# Add DVs +prob.model.add_design_var("chord_bspline.chord_cp_spline", lower=0.1, upper=10.0, units=None) +prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + +# Add CL constraint +prob.model.add_constraint(point_name + ".CL", equals=0.3) + +# Add Area constraint +prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + +# Add objective +prob.model.add_objective(point_name + ".CD", scaler=1e4) + +prob.driver = om.ScipyOptimizeDriver() +prob.driver.options["optimizer"] = "SLSQP" +prob.driver.options["tol"] = 1e-7 +prob.driver.options["disp"] = True +prob.driver.options["maxiter"] = 1000 +# prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + +# Set up and run the optimization problem +prob.setup() +prob.run_driver() +# om.n2(prob) + + +mesh1 = prob.get_val("surface.sec0.mesh", units="m") +mesh2 = prob.get_val("surface.sec1.mesh", units="m") + +meshUni = prob.get_val(name + "." + unification_name + "." + name + "_uni_mesh") + + +def plot_meshes(meshes): + """this function plots to plot the mesh""" + plt.figure(figsize=(8, 4)) + for i, mesh in enumerate(meshes): + mesh_x = mesh[:, :, 0] + mesh_y = mesh[:, :, 1] + color = "k" + for i in range(mesh_x.shape[0]): + plt.plot(mesh_y[i, :], mesh_x[i, :], color, lw=1) + plt.plot(-mesh_y[i, :], mesh_x[i, :], color, lw=1) # plots the other side of symmetric wing + for j in range(mesh_x.shape[1]): + plt.plot(mesh_y[:, j], mesh_x[:, j], color, lw=1) + plt.plot(-mesh_y[:, j], mesh_x[:, j], color, lw=1) # plots the other side of symmetric wing + plt.axis("equal") + plt.xlabel("y (m)") + plt.ylabel("x (m)") + # plt.savefig('opt_planform_visc_construction.png') + + +plot_meshes([meshUni]) +# plt.show() diff --git a/openaerostruct/docs/user_reference/mesh_surface_dict.rst b/openaerostruct/docs/user_reference/mesh_surface_dict.rst index 9f53ccbcd..a3cce0750 100644 --- a/openaerostruct/docs/user_reference/mesh_surface_dict.rst +++ b/openaerostruct/docs/user_reference/mesh_surface_dict.rst @@ -128,6 +128,59 @@ The surface dict will be provided to Groups, including ``Geometry``, ``AeroPoint - - Position of reference axis along the chord about which to apply twist, chord, taper, and span geometry transformations. 1 is the trailing edge, 0 is the leading edge. +.. list-table:: Multi-section Surface definition + :widths: 20 20 5 55 + :header-rows: 1 + + * - Key + - Example value + - Units + - Description + * - is_multi_section + - True + - + - This key must be present and set to True for the AeroPoint to correctly interpret this surface as multi-section. + * - num_sections + - 2 + - + - The number of sections in the multi-section surface. + * - sec_name + - ["sec0","sec1","sec2"] + - + - Names of the individual sections. Each section must be named and the list length must match the specified number of sections. + * - meshes + - "gen-meshes" or [mesh1,mesh2,...] + - + - Supply a list of meshes for each section or "gen-meshes" for automatic mesh generation + * - root_chord + - 1.0 + - m + - Root chord length of the section indicated as "root section"(required if using the built-in mesh generator) + * - span + - [10.0,10.0] + - m + - Wing span for each section. The list length must match the specified number of sections. + * - ny + - [21,21] + - + - Number of spanwise points for each section. The list length must match the specified number of sections. (required if using the built-in multi-section mesh generator) + * - nx + - 10 + - + - Number of chordwise points. Same for all sections.(required if using the built-in multi-section mesh generator) + * - bpanels + - [10,10] + - + - Number of spanwise panels for each section. The list length must match the specified number of sections. An alternative to specifying nx. + * - cpanels + - [10,10] + - + - Number of chordwise panels for each section. An alternative to specifying ny. + * - root_section + - 1 + - + - Root chord length of the section indicated as "root section"(required if using the built-in mesh generator) + .. list-table:: Aerodynamics definitions :widths: 20 20 5 55 :header-rows: 1 diff --git a/openaerostruct/geometry/geometry_group.py b/openaerostruct/geometry/geometry_group.py index b7563e006..8a633e16a 100644 --- a/openaerostruct/geometry/geometry_group.py +++ b/openaerostruct/geometry/geometry_group.py @@ -2,6 +2,9 @@ import openmdao.api as om from openaerostruct.utils.check_surface_dict import check_surface_dict_keys +import openaerostruct.geometry.geometry_mesh_gen as meshGen +from openaerostruct.geometry.geometry_unification import GeomMultiUnification +from openaerostruct.geometry.geometry_multi_join import GeomMultiJoin from openaerostruct.utils.interpolation import get_normalized_span_coords @@ -190,3 +193,155 @@ def setup(self): self.add_subsystem( "mesh", GeometryMesh(surface=surface), promotes_inputs=bsp_inputs, promotes_outputs=["mesh"] ) + + +# Function that constructs the individual section surface data dictionaries +def build_sections(surface): + """This function returns an OpenMDAO Independent Variable Component with an output vector appropriately + named and sized to function as an unified B-spline that joins multiple sections by construction. + + Parameters + ---------- + surface: dict + OpenAeroStruct multi-section surface dictionary + + Returns + ------- + section_surfaces : list + List of OpenAeroStruct surface dictionaries for each individual surface + + """ + # Get number of sections + num_sections = surface["num_sections"] + + if surface["meshes"] == "gen-meshes": + # Verify that all required inputs for automatic mesh generation are provided for each section + if len(surface["ny"]) != num_sections: + raise ValueError("Number of spanwise points needs to be provided for each section") + if len(surface["taper"]) != num_sections: + raise ValueError("Taper needs to be provided for each section") + if len(surface["span"]) != num_sections: + raise ValueError("Span needs to be provided for each section") + if len(surface["sweep"]) != num_sections: + raise ValueError("Sweep needs to be provided for each section") + + # Generate unified and individual section meshes + mesh, sec_meshes = meshGen.generate_mesh(surface) + else: + # Allow user to provide mesh for each section + if len(surface["meshes"]) != num_sections: + raise ValueError("A mesh needs to be provided for each section.") + sec_meshes = surface["meshes"] + + if len(surface["sec_name"]) != num_sections: + raise ValueError("A name needs to be provided for each section.") + + # List of support keys for multi-section wings + # NOTE: make sure this is consistent to the documentation's surface dict page + target_keys = [ + # Essential Info + "num_section", + "symmetry", + "S_ref_type", + "ref_axis_pos", + # wing definition + "span", + "taper", + "sweep", + "dihedral", + "twist_cp", + "chord_cp", + "xshear_cp", + "yshear_cp", + "zshear_cp", + # aerodynamics + "CL0", + "CD0", + "with_viscous", + "with_wave", + "groundplane", + "k_lam", + "t_over_c_cp", + "c_max_t", + ] + + # Constructs a list of section dictionaries and adds the specified supported keys and values from the mult-section surface dictionary. + surface_sections = [] + num_sections = surface["num_sections"] + + for i in range(num_sections): + section = {} + for k in set(surface).intersection(target_keys): + if type(surface[k]) is list: + section[k] = surface[k][i] + else: + section[k] = surface[k] + section["mesh"] = sec_meshes[i] + section["name"] = surface["sec_name"][i] + surface_sections.append(section) + return surface_sections + + +class MultiSecGeometry(om.Group): + """ + Group that contains the section geometery groups for the multi-section surface + + + This group handles the creation of each section geometry group based on parameters + supplied in the multi-section surface dictionary. Meshes for each section can be + provided by the user or automatically generated based on parameters supplied in the + surface dictionary. The group also adds a mesh unification component that combines the + individual section for each mesh into a singular unified mesh for use in aero components. + Optionally, the joining component can be added that computes the edge distances between sections. + This information can be used to set a distance constraint along the specified axes if needed. + """ + + def initialize(self): + self.options.declare("surface", types=dict) # Multi-section surface dictionary + self.options.declare( + "joining_comp", types=bool, default=False + ) # Specify if a distance computation component should be added + self.options.declare( + "dim_constr", types=list, default=[] + ) # List of arrays corresponding to each shared edge between section along the surface. Each array inidicates along which axes the distance constarint is applied([x y z]) + + def setup(self): + surface = self.options["surface"] + joining_comp = self.options["joining_comp"] + dc = self.options["dim_constr"] + + # key validation of the surface dict + check_surface_dict_keys(surface) + + sec_dicts = build_sections(surface) + + section_names = [] + for sec in sec_dicts: + geom_group = Geometry(surface=sec) + self.add_subsystem(sec["name"], geom_group) + section_names.append(sec["name"]) + + # Add the mesh unification component + unification_name = "{}_unification".format(surface["name"]) + + uni_mesh = GeomMultiUnification(sections=sec_dicts, surface_name=surface["name"]) + self.add_subsystem(unification_name, uni_mesh) + + # Connect each section mesh to mesh unification component inputs + for sec_name in section_names: + self.connect("{}.mesh".format(sec_name), "{}.{}_def_mesh".format(unification_name, sec_name)) + + # Connect each section t over c B-spline to t over c unification component if needed + if "t_over_c_cp" in surface.keys(): + for sec_name in section_names: + self.connect("{}.t_over_c".format(sec_name), "{}.{}_t_over_c".format(unification_name, sec_name)) + + if joining_comp: + # Add section joining component to output edge distances + joining_name = "{}_joining".format(surface["name"]) + + join = GeomMultiJoin(sections=sec_dicts, dim_constr=dc) + self.add_subsystem(joining_name, join) + + for sec_name in section_names: + self.connect("{}.mesh".format(sec_name), "{}.{}_join_mesh".format(joining_name, sec_name)) diff --git a/openaerostruct/geometry/geometry_mesh_gen.py b/openaerostruct/geometry/geometry_mesh_gen.py new file mode 100644 index 000000000..70bba2a32 --- /dev/null +++ b/openaerostruct/geometry/geometry_mesh_gen.py @@ -0,0 +1,320 @@ +""" Utility for quickly generating a multi-section user specificed OAS mesh +""" + +import numpy as np +import matplotlib.pyplot as plt + + +def generate_mesh(surface): + """ + Automatically generates meshes for single and multi-section wings using the OpenAeroStruct surface dictionary. + + Parameters + ---------- + surface : dict + OpenAeroStruct surface or multi-section surface dictionary. + + Returns + ------- + mesh[nx, ny, 3] : numpy array + Nodal mesh defining the aerodynamic surface. + + sec_mesh : list + List of nodal meshes corresponding to each section in a multi-section surface. + """ + + # Get Data + num_sections = surface["num_sections"] + symmetry = surface["symmetry"] + + if symmetry or num_sections == 1: + root_section = num_sections - 1 + else: + if "root_section" not in surface.keys(): + raise Exception("The root section of an asymmetrical mesh needs to be identified") + else: + root_section = surface["root_section"] + + # Geometry data dictionary + section_data = { + "taper": surface["taper"], + "sweep": surface["sweep"], + "span": surface["span"], + "root_chord": surface["root_chord"], + } + + if "bpanels" in surface.keys(): + ny = surface["bpanels"] + 1 + else: + ny = surface["ny"] + + if "cpanels" in surface.keys(): + nx = surface["cpanels"] + 1 + else: + nx = surface["nx"] + + # Generate the mesh data + panel_gx, panel_gy = generate_section_geometry(num_sections, symmetry, section_data, ny, nx, root_section) + + # Sitch the mesh into a unified mesh and out in OAS format + panel_geom_x, panel_geom_y = stitch_section_geometry(num_sections, panel_gy, panel_gx) + mesh = output_oas_mesh(panel_geom_x, panel_geom_y) + + # Produce meshes for each section in OAS format + sec_meshes = [] + for section in range(num_sections): + sec_mesh = output_oas_mesh(panel_gx[section], panel_gy[section]) + sec_meshes.append(sec_mesh) + + return mesh, sec_meshes + + +def generate_section_geometry(sections, symmetry, section_data, ny, nx, root_section): + """ + Constructs the multi-section wing geometry specified by the user and generates a mesh for each section. + + Parameters + ---------- + sections : int + Integer for number of wing sections specified + symmetry : bool + Bool inidicating if the funciton should only generate the left span of the wing + section_data : dict + Dictionary with arrays corresponding the taper, span, sweep, and root chord of each section + ny : numpy array + Array with ints correponding to the number of spanwise points per section + nx : int + Number of chordwise points + root_section: + The section number that should be treated as the root section(y=0 origin) + + + Returns + ------- + panel_gx : List + List containing the mesh x-coordinates for each section + panel_gy : List + List containing the mesh y-coordinates for each section + """ + + # Preallocate the lists + panel_gy = [None] * sections + panel_gx = [None] * sections + + # Jump to root section and build left wing + for sec in np.arange(root_section, -1, -1): + taper = section_data["taper"][sec] + b = section_data["span"][sec] + le_lambda = section_data["sweep"][sec] + + if sec == root_section: + root_c = section_data["root_chord"] + root_te = 0 + root_y = 0 + else: + root_c = np.abs(panel_gx[sec + 1][0, 0] - panel_gx[sec + 1][nx - 1, 0]) + root_te = panel_gx[sec + 1][nx - 1, 0] + root_y = panel_gy[sec + 1][0] + tip_c = root_c * taper + + root_le = root_c + root_te + tip_le = root_le - b * np.tan(le_lambda) + tip_te = tip_le - tip_c + + root_x = np.linspace(root_le, root_te, nx) + + if tip_le == tip_te: + tip_x = tip_le * np.ones(len(root_x)) + else: + tip_x = np.linspace(tip_le, tip_te, len(root_x)) + + panel_geom_y = np.zeros(ny[sec]) + panel_geom_x = np.zeros([nx, ny[sec]]) + + panel_geom_y = np.linspace(root_y - b, root_y, ny[sec]) + + for i in range(len(root_x)): + panel_geom_x[i, :] = root_x[i] - ((tip_x[i] - root_x[i]) / b) * (panel_geom_y - root_y) + + panel_gy[sec] = panel_geom_y + panel_gx[sec] = panel_geom_x + + # Build the right wing if asymmetrical + if not symmetry: + for sec in np.arange(root_section + 1, sections): + taper = section_data["taper"][sec] + b = section_data["span"][sec] + le_lambda = section_data["sweep"][sec] + + root_c = np.abs(panel_gx[sec - 1][0, -1] - panel_gx[sec - 1][nx - 1, -1]) + tip_c = root_c * taper + root_y = panel_gy[sec - 1][-1] + + root_te = panel_gx[sec - 1][nx - 1, 0] + root_le = root_c + root_te + tip_le = root_le + b * np.tan(le_lambda) + tip_te = tip_le - tip_c + + root_x = np.linspace(root_le, root_te, nx) + + if tip_le == tip_te: + tip_x = tip_le * np.ones(len(root_x)) + else: + tip_x = np.linspace(tip_le, tip_te, len(root_x)) + + panel_geom_y = np.zeros(ny[sec]) + panel_geom_x = np.zeros([nx, ny[sec]]) + + panel_geom_y[0 : ny[sec]] = np.linspace(root_y, root_y + b, ny[sec]) + + for i in range(len(root_x)): + panel_geom_x[i, :] = root_x[i] + ((tip_x[i] - root_x[i]) / (b / 2)) * (panel_geom_y - root_y) + + panel_gy[sec] = panel_geom_y + panel_gx[sec] = panel_geom_x + + return panel_gx, panel_gy + + +def stitch_section_geometry(sections, panel_gy, panel_gx): + """ + Combines the split section array into singular unified mesh + + Parameters + ---------- + sections : int + Integer for number of wing sections specified + panel_gx : List + List containing the mesh x-coordinates for each section + panel_gy : List + List containing the mesh y-coordinates for each section + + + Returns + ------- + panel_geom_x : numpy array + Array of the mesh x-coordinates + panel_geom_y : numpy array + Array of the mesh y-coordinates + """ + # Stitch the results into a singular mesh + + if sections > 1: + panel_geom_y = panel_gy[0][:-1] + panel_geom_x = panel_gx[0][:, :-1] + for i in np.arange(1, sections - 1): + panel_geom_y = np.concatenate((panel_geom_y, panel_gy[i][:-1])) + panel_geom_x = np.concatenate((panel_geom_x, panel_gx[i][:, :-1]), axis=1) + panel_geom_y = np.concatenate((panel_geom_y, panel_gy[sections - 1])) + panel_geom_x = np.concatenate((panel_geom_x, panel_gx[sections - 1]), axis=1) + else: + panel_geom_y = panel_gy[0] + panel_geom_x = panel_gx[0] + return panel_geom_x, panel_geom_y + + +def reflect_symmetric(panel_geom_x, panel_geom_y): + """ + Reflects the mesh over y=0 + + Parameters + ---------- + panel_geom_x : numpy array + Array of the mesh x-coordinates + panel_geom_y : numpy array + Array of the mesh y-coordinates + + Returns + ------- + panel_geom_x : numpy array + Array of the mesh x-coordinates + panel_geom_y : numpy array + Array of the mesh y-coordinates + """ + panel_geom_x = np.hstack((panel_geom_x, panel_geom_x)) + panel_geom_y = np.hstack((panel_geom_y, -panel_geom_y)) + return panel_geom_x, panel_geom_y + + +def output_oas_mesh(panel_geom_x, panel_geom_y): + """ + Outputs the mesh in OAS format + + Parameters + ---------- + panel_geom_x : numpy array + 2D array of the mesh x-coordinates + panel_geom_y : numpy array + 1D array of the mesh y-coordinates + + Returns + ------- + mesh : numpy array + 3-D array with the OAS format mesh + """ + panel_geom_y = np.broadcast_to(panel_geom_y, (panel_geom_x.shape[0], len(panel_geom_y))) + mesh = np.zeros((panel_geom_x.shape[0], panel_geom_y.shape[1], 3)) + mesh[:, :, 0] = panel_geom_x[::-1] + mesh[:, :, 1] = panel_geom_y + return mesh + + +if __name__ == "__main__": + """Runs the mesh generator to be run independently of OAS. Example 2 section mesh provided.""" + + surface = { + # Wing definition + # Basic surface parameters + "name": "surface", + "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": np.array([1.0, 1.0]), # Wing taper for each section + "span": np.array([1.0, 1.0]), # Wing span for each section + "sweep": np.array([0.0, 0.0]), # Wing sweep for each section + "chord_cp": [np.array([1, 1]), np.array([1.0, 0.2])], + # "sec_chord_cp": [np.ones(1),2*np.ones(1),3*np.ones(1)], #Chord B-spline control points for each section + "root_chord": 1.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": np.array([2, 2]), # 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 + "sec_t_over_c_cp": [np.array([0.15]), np.array([0.15])], # thickness over chord ratio (NACA0015) + "sec_c_max_t": 0.303, # chordwise location of maximum (NACA0015) + # thickness + "with_viscous": False, # if true, compute viscous drag + "with_wave": False, # if true, compute wave drag + "groundplane": False, + } + + meshT, sec_meshes = generate_mesh(surface) + + def plot_meshes(meshes): + """this function plots to plot the mesh""" + plt.figure(figsize=(8, 4)) + for i, mesh in enumerate(meshes): + mesh_x = mesh[:, :, 0] + mesh_y = mesh[:, :, 1] + color = "k" + for i in range(mesh_x.shape[0]): + plt.plot(mesh_y[i, :], mesh_x[i, :], color, lw=1) + # plt.plot(-mesh_y[i, :], mesh_x[i, :], color, lw=1) # plots the other side of symmetric wing + for j in range(mesh_x.shape[1]): + plt.plot(mesh_y[:, j], mesh_x[:, j], color, lw=1) + # plt.plot(-mesh_y[:, j], mesh_x[:, j], color, lw=1) # plots the other side of symmetric wing + plt.axis("equal") + plt.xlabel("y (m)") + plt.ylabel("x (m)") + + plot_meshes(sec_meshes) + plt.show() diff --git a/openaerostruct/geometry/geometry_multi_join.py b/openaerostruct/geometry/geometry_multi_join.py new file mode 100644 index 000000000..cafe3b56c --- /dev/null +++ b/openaerostruct/geometry/geometry_multi_join.py @@ -0,0 +1,186 @@ +import numpy as np +import openmdao.api as om + +default_vec = np.ones(3) + + +def get_section_edge_left(mesh, v=default_vec, edge_cur=0, edges_all_constraints=default_vec): + """ + Function that gets the coordinates of the leading and trailing edge points of the left edge of a section. The output can be masked to only retreive the x,y, or z coordinate. + The function also returns the row and column vectors of non zero entries in the edge coordinate jacobian. + + Parameters + ---------- + mesh : numpy array + OAS mesh of a given section + v : numpy array[3] + Numpy array of length three populuted with ones and zeros used to mask the output so that only the specified of the x, y, and z coordinates are returned. + edge_cur : int + Integer indicating which unique intersection of edges this particular edge is associated with. Required to return the correct jacobian sparsity pattern. + edges_all_constraints: list + See dim_constr in the GeomMultiJoin component. This array needs to passed into this function from the component in order to return the correct jacobian sparsity pattern. + + Returns + ------- + edge points : numpy array + Array of points corresponding to the leading and trailing edges of the left section edge. Masked according to input. + rows : numpy array + Array of the rows of the non-zero jacobian entries + cols : numpy array + Array of the columns of the non-zero jacobian entries + + """ + nx = mesh.shape[0] + le_index = 0 + te_index = np.ravel_multi_index((nx - 1, 0, 0), mesh.shape) + mask = np.array(v, dtype="bool") + + rows = np.arange(0, 2 * np.sum(v)) + 2 * int(np.sum(edges_all_constraints[:edge_cur])) + cols = np.concatenate([np.arange(le_index, le_index + 3)[mask], np.arange(te_index, te_index + 3)[mask]]) + + return mesh[[0, -1], 0][:, np.arange(0, 3)[mask]], rows, cols + + +def get_section_edge_right(mesh, v=default_vec, edge_cur=0, edges_all_constraints=default_vec): + """ + Function that gets the coordinates of the leading and trailing edge points of the right edge of a section. The output can be masked to only retreive the x,y, or z coordinate. + The function also returns the row and column vectors of non zero entries in the edge coordinate jacobian. + + Parameters + ---------- + mesh : numpy array + OAS mesh of a given section + v : numpy array[3] + Numpy array of length three populuted with ones and zeros used to mask the output so that only the specified of the x, y, and z coordinates are returned. + edge_cur : int + Integer indicating which unique intersection of edges this particular edge is associated with. Required to return the correct jacobian sparsity pattern. + edges_all_constraints: list + See dim_constr in the GeomMultiJoin component. This array needs to passed into this function from the component in order to return the correct jacobian sparsity pattern. + + Returns + ------- + edge points : numpy array + Array of points corresponding to the leading and trailing edges of the right section edge. Masked according to input. + rows : numpy array + Array of the rows of the non-zero jacobian entries + cols : numpy array + Array of the columns of the non-zero jacobian entries + """ + + nx = mesh.shape[0] + ny = mesh.shape[1] + le_index = np.ravel_multi_index((0, ny - 1, 0), mesh.shape) + te_index = np.ravel_multi_index((nx - 1, ny - 1, 0), mesh.shape) + mask = np.array(v, dtype="bool") + + rows = np.arange(0, 2 * np.sum(v)) + 2 * int(np.sum(edges_all_constraints[:edge_cur])) + cols = np.concatenate([np.arange(le_index, le_index + 3)[mask], np.arange(te_index, te_index + 3)[mask]]) + + return mesh[[0, -1], -1][:, np.arange(0, 3)[mask]], rows, cols + + +class GeomMultiJoin(om.ExplicitComponent): + """ + OpenMDAO component that outputs the distance between the leading and trailing edge corners of each section corresponding to each shared edge. + Users can specify along which axes the distance should be computed([x y z]). + + Parameters + ---------- + sections : list + List of section OpenAeroStruct surface dictionaries + dim_constr : list + A list of vectors of length three corresponding to each edge. Entries correspond to the dimension([x,y,z]) the user wishes to constraint should be set to 1. Remaining entries should be zero. + + Returns + ------- + section_separation[2*count_nonzero(concatenate(dim_constr)] : numpy array + Numpy array of distances along each axis between corresponding leading and trailing edge points between subsequent sections. + """ + + def initialize(self): + """ + Declare options. + """ + self.options.declare("sections", types=list, desc="A list of section surface dictionaries to be joined.") + self.options.declare( + "dim_constr", + types=list, + default=[np.ones(3)], + desc="A list of vectors of length three corresponding to each edge. Entries corresponding the dimension([x,y,z]) the user wishes to constraint should be set to 1. Remaining entries should be zero.", + ) + + def setup(self): + sections = self.options["sections"] + self.num_sections = len(sections) + self.dim_constr = self.options["dim_constr"] + + # Compute total number of unique intersecting edges between each section + edge_total = self.num_sections - 1 + + # Defaults to distane along x-axis for each intersecting edge + if len(self.dim_constr) != (edge_total): + self.dim_constr = [np.array([1, 0, 0]) for i in range(edge_total)] + + # Compute size of and add output + constr_size = 2 * np.count_nonzero(np.concatenate(self.dim_constr)) + self.add_output("section_separation", val=np.zeros(constr_size)) + + """Generate the Jacobian of the edge seperation distance with respect to the section mesh. + Jacobian is just ones, zeros, and negative ones so we can declare here. However the sparsity pattern is complicated and requires two helper functions and the loop below.""" + + # Counter used to track the current unique edge interection between section being processed. + edge_cur = 0 + + for i_sec, section in enumerate(sections): + mesh = section["mesh"] + nx = mesh.shape[0] + ny = mesh.shape[1] + name = section["name"] + + # Add the input + mesh_name = "{}_join_mesh".format(name) + self.add_input(mesh_name, shape=(nx, ny, 3), units="m") + + # Get the sparsity patterns for each section. First and last sections only have one edge intersection. + if i_sec == 0: + rows, cols = get_section_edge_right(mesh, self.dim_constr[i_sec], edge_cur, self.dim_constr)[1:] + vals = -1 * np.ones_like(rows) + elif i_sec < len(sections) - 1: + rows1, cols1 = get_section_edge_left(mesh, self.dim_constr[i_sec - 1], edge_cur, self.dim_constr)[1:] + vals1 = np.ones_like(rows1) + + edge_cur += 1 + rows2, cols2 = get_section_edge_right(mesh, self.dim_constr[i_sec], edge_cur, self.dim_constr)[1:] + vals2 = -1 * np.ones_like(rows2) + + rows = np.concatenate([rows1, rows2]) + cols = np.concatenate([cols1, cols2]) + vals = np.concatenate([vals1, vals2]) + else: + rows, cols = get_section_edge_left(mesh, self.dim_constr[i_sec - 1], edge_cur, self.dim_constr)[1:] + vals = np.ones_like(rows) + + # Declare partials for the current section + self.declare_partials("section_separation", mesh_name, rows=rows, cols=cols, val=vals) + + def compute(self, inputs, outputs): + # Compute the distances between the corresponding leading and trailing edges along the edge interection between each section + sections = self.options["sections"] + edges = [] + edge_constraints = [] + for i_sec, section in enumerate(sections): + name = section["name"] + mesh_name = "{}_join_mesh".format(name) + + if i_sec == 0: + edges.append(get_section_edge_right(inputs[mesh_name], self.dim_constr[i_sec])[0]) + elif i_sec < len(sections) - 1: + edges.append(get_section_edge_left(inputs[mesh_name], self.dim_constr[i_sec - 1])[0]) + edges.append(get_section_edge_right(inputs[mesh_name], self.dim_constr[i_sec])[0]) + else: + edges.append(get_section_edge_left(inputs[mesh_name], self.dim_constr[i_sec - 1])[0]) + + for i in range(self.num_sections - 1): + edge_constraints.append((edges[2 * i + 1] - edges[2 * i]).flatten()) + + outputs["section_separation"] = np.array(edge_constraints).flatten() diff --git a/openaerostruct/geometry/geometry_unification.py b/openaerostruct/geometry/geometry_unification.py new file mode 100644 index 000000000..e2c1d8b1f --- /dev/null +++ b/openaerostruct/geometry/geometry_unification.py @@ -0,0 +1,231 @@ +import numpy as np +import openmdao.api as om + + +def compute_uni_mesh_dims(sections): + """ + Function computes the dimensions of the unified mesh as well as an index array + + Parameters + ---------- + sections : list + List of section OpenAeroStruct surface dictionaries + + Returns + ------- + uni_mesh_indices : numpy array + Array of indicies matching the size of the mesh array + uni_nx : int + Number of chordwise points in the unified mesh + uni_ny : int + Number of spanwise points in the unified mesh + """ + uni_nx = sections[0]["mesh"].shape[0] + uni_ny = 0 + for i_sec in range(len(sections)): + if i_sec == len(sections) - 1: + uni_ny += sections[i_sec]["mesh"].shape[1] + else: + uni_ny += sections[i_sec]["mesh"].shape[1] - 1 + uni_mesh_indices = np.arange(uni_nx * uni_ny * 3).reshape((uni_nx, uni_ny, 3)) + + return uni_mesh_indices, uni_nx, uni_ny + + +def compute_uni_mesh_index_blocks(sections, uni_mesh_indices): + """ + Function that computes the index block that corresponds to each individual wing section + + Parameters + ---------- + sections : list + List of section OpenAeroStruct surface dictionaries + uni_mesh_indices : numpy array + Array of indicies matching the size of the mesh array + + Returns + ------- + blocks : list + List of numpy arrays containing the incidies corresponding to each section. Basically breaks the uni_mesh_indices + array up into pieces that correspond to each section. + """ + blocks = [] + + # cursor to track the y position of each section along the unified mesh + y_curr = 0 + + for i_sec in range(len(sections)): + mesh = sections[i_sec]["mesh"] + ny = mesh.shape[1] + + if i_sec == len(sections) - 1: + block = uni_mesh_indices[:, y_curr:, :] + y_curr += ny + else: + block = uni_mesh_indices[:, y_curr : y_curr + (ny - 1), :] + y_curr += ny - 1 + + blocks.append(block.flatten()) + + return blocks + + +def unify_mesh(sections): + """ + Function that produces a unified mesh from all the individual wing section meshes. + + Parameters + ---------- + sections : list + List of section OpenAeroStruct surface dictionaries + + Returns + ------- + uni_mesh : numpy array + Unfied surface mesh in OAS format + """ + for i_sec in np.arange(0, len(sections) - 1): + mesh = sections[i_sec]["mesh"] + + import copy + + if i_sec == 0: + uni_mesh = copy.deepcopy(mesh[:, :-1, :]) + else: + uni_mesh = np.concatenate([uni_mesh, mesh[:, :-1, :]], axis=1) + # Stitch the results into a singular mesh + mesh = sections[len(sections) - 1]["mesh"] + if len(sections) == 1: + import copy + + uni_mesh = copy.deepcopy(mesh) + else: + uni_mesh = np.concatenate([uni_mesh, mesh], axis=1) + + return uni_mesh + + +class GeomMultiUnification(om.ExplicitComponent): + """ + OpenMDAO component that combines the meshes associated with each individual section + of a multi-section wing into a unified mesh. Section 0 as defined by the user(typically + the most inboard section) will be added to the unified mesh unchanged. The remaining sections + are then sequencially added outboard to section 0 along the y-direction. The mesh unification + component operates under the assumption that it is unifying section that are C0 continous along + the span. The most inboard column of points of section i+1 will be coincident with the + most outbord column of points of section i. As result, the inboard column of points for each section + after section 0 are removed prior to the section being added to the unified mesh. The component also + produces the appropriate sparse Jacobian matrix for the column removal so that the partials can be computed + appropriately. + + Parameters + ---------- + sections : list + A list of section surface dictionaries in OAS format. + + Returns + ------- + mesh[nx, ny, 3] : numpy array + Nodal mesh defining the unified surface + """ + + def initialize(self): + """ + Declare options. + """ + self.options.declare("sections", types=list, desc="A list of section surface dictionaries to be unified.") + self.options.declare("surface_name", types=str, desc="The name of the multi-section surface") + + def setup(self): + sections = self.options["sections"] + name = self.options["surface_name"] + + # Get the unified mesh size, index array, and block of indicies for each section + [uni_mesh_indices, uni_nx, uni_ny] = compute_uni_mesh_dims(sections) + uni_mesh_blocks = compute_uni_mesh_index_blocks(sections, uni_mesh_indices) + uni_mesh_name = "{}_uni_mesh".format(name) + + # Loop through each section to build the sparsity pattern for the Jacobian. This Jacobian is fixed based on the mesh size so can be declared here + for i_sec, section in enumerate(sections): + mesh = section["mesh"] + nx = mesh.shape[0] + ny = mesh.shape[1] + name = section["name"] + + mesh_name = "{}_def_mesh".format(name) + + self.add_input(mesh_name, shape=(nx, ny, 3), units="m", tags=["mphys_coupling"]) + + # Generate index array + mesh_indices = np.arange(nx * ny * 3).reshape((nx, ny, 3)) + + if i_sec == len(sections) - 1: + cols = mesh_indices.flatten() + else: + # If not section N then the most inboard column is disregarded + cols = mesh_indices[:, :-1, :].flatten() + + # Get data from section block in unified mesh + rows = uni_mesh_blocks[i_sec] + + # Fill non zero Jacobian entries with ones + data = np.ones_like(rows) + + self.declare_partials(uni_mesh_name, mesh_name, val=data, rows=rows, cols=cols) + + self.add_output(uni_mesh_name, shape=(uni_nx, uni_ny, 3), units="m") + + # Unify the t/c output of each section if that has been specified + if "t_over_c_cp" in sections[0].keys(): + uni_tc_name = "{}_uni_t_over_c".format(self.options["surface_name"]) + + n_acc = 0 + for section in sections: + name = section["name"] + t_over_c_name = "{}_t_over_c".format(name) + n = int(ny - 1) + self.add_input(t_over_c_name, shape=(n), tags=["mphys_coupling"]) + + self.declare_partials( + uni_tc_name, t_over_c_name, rows=np.arange(0, n) + n_acc, cols=np.arange(0, n), val=np.ones(n) + ) + + n_acc += n + + self.add_output(uni_tc_name, shape=(uni_ny - 1)) + + def compute(self, inputs, outputs): + sections = self.options["sections"] + surface_name = self.options["surface_name"] + uni_mesh_name = "{}_uni_mesh".format(surface_name) + + # Loop through all sections to unify the mesh + for i_sec in np.arange(0, len(sections) - 1): + name = sections[i_sec]["name"] + mesh_name = "{}_def_mesh".format(name) + + if i_sec == 0: + uni_mesh = inputs[mesh_name][:, :-1, :] + else: + uni_mesh = np.concatenate([uni_mesh, inputs[mesh_name][:, :-1, :]], axis=1) + + mesh_name = "{}_def_mesh".format(sections[len(sections) - 1]["name"]) + if len(sections) == 1: + uni_mesh = inputs[mesh_name] + else: + uni_mesh = np.concatenate([uni_mesh, inputs[mesh_name]], axis=1) + + if "t_over_c_cp" in sections[0].keys(): + uni_tc_name = "{}_uni_t_over_c".format(self.options["surface_name"]) + for i_sec, section in enumerate(sections): + name = section["name"] + t_over_c_name = "{}_t_over_c".format(name) + t_over_c = inputs[t_over_c_name] + + if i_sec == 0: + uni_t_over_c = inputs[t_over_c_name] + else: + uni_t_over_c = np.concatenate([uni_t_over_c, t_over_c]) + outputs[uni_tc_name] = uni_t_over_c + + outputs[uni_mesh_name] = uni_mesh diff --git a/openaerostruct/geometry/multi_unified_bspline_utils.py b/openaerostruct/geometry/multi_unified_bspline_utils.py new file mode 100644 index 000000000..deca793bf --- /dev/null +++ b/openaerostruct/geometry/multi_unified_bspline_utils.py @@ -0,0 +1,85 @@ +import numpy as np +import openmdao.api as om + + +def build_multi_spline(out_name, num_sections, control_points): + """This function returns an OpenMDAO Independent Variable Component with an output vector appropriately + named and sized to function as an unified B-spline that joins multiple sections by construction. + + Parameters + ---------- + out_name: string + Name of the output to assign to the B-spline + num_sections : int + Number of sections + control_points: list + List of B-spline control point arrays corresponding to each section + + Returns + ------- + spline_control : OpenMDAO component object + The unified B-spline indpendent variable component + + """ + if len(control_points) != num_sections: + raise Exception("Target sections need to match with control points!") + + single_sections = len([cp for cp in control_points if len(cp) == 1]) + + control_poin_vec = np.ones(len(np.concatenate(control_points)) - (num_sections - 1 - single_sections)) + + spline_control = om.IndepVarComp() + spline_control.add_output("{}_spline".format(out_name), val=control_poin_vec) + + return spline_control + + +def connect_multi_spline(prob, section_surfaces, sec_cp, out_name, comp_name, return_bind_inds=False): + """This function connects the the unified B-spline component with the individual B-splines + of each section. There is a point of overlap at each section so that each edge control point control the edge + controls points of each section's B-spline. This is how section joining by consturction is acheived. + An issue occurs however when a B-spline in a particular section only has one control point. In this case the one + section control point is bound to the left edge B-spline component control point. As result, there is nothing to + maintain C0 continuity with the next section. As result a constraint will need to be manually set. To facilitate this, + the array bind_inds will contain a list of the B-spline control point indicies that will need to be manually constrained to + their previous sections. + + + Parameters + ---------- + prob : OpenMDAO problem object + The OpenAeroStruct problem object with the unified B-spline component added. + section_surfaces : list + List of the surface dictionaries for each section. + sec_cp : list + List of B-spline control point arrays for each section. + out_name: string + Name of the unified B-spline component output to connect from + comp_name: string + Name of the unified B-spline component added to the problem object + return_bind_inds: bool + Return list of unjoined unified B-spline inidices. Default is False. + + Returns + ------- + bind_inds : list + List of unified B-spline control point indicies not connected due to the presence of a single control point section.(Only if return bind_inds specified) + + """ + acc = 0 + bind_inds = [] + for i, section in enumerate(section_surfaces): + point_count = len(sec_cp[i]) + src_inds = np.arange(acc, acc + point_count) + acc += point_count - 1 + if point_count == 1: + acc += 1 + bind_inds.append(acc) + prob.model.connect( + "{}.{}".format(comp_name, out_name) + "_spline", + "surface." + section["name"] + ".{}".format(out_name), + src_indices=src_inds, + ) + + if return_bind_inds: + return bind_inds diff --git a/openaerostruct/geometry/tests/multisection/test_joining_comp.py b/openaerostruct/geometry/tests/multisection/test_joining_comp.py new file mode 100644 index 000000000..0c4857bfa --- /dev/null +++ b/openaerostruct/geometry/tests/multisection/test_joining_comp.py @@ -0,0 +1,20 @@ +import numpy as np +import unittest + +from openaerostruct.geometry.geometry_multi_join import GeomMultiJoin +from openaerostruct.geometry.geometry_group import build_sections +from openaerostruct.utils.testing import run_test, get_three_section_surface + + +class Test(unittest.TestCase): + def test(self): + (surface, chord_bspline) = get_three_section_surface() + sec_dicts = build_sections(surface) + + comp = GeomMultiJoin(sections=sec_dicts, dim_constr=[np.ones(3), np.ones(3), np.ones(3)]) + + run_test(self, comp, complex_flag=True, method="cs") + + +if __name__ == "__main__": + unittest.main() diff --git a/openaerostruct/geometry/tests/multisection/test_unification_comp.py b/openaerostruct/geometry/tests/multisection/test_unification_comp.py new file mode 100644 index 000000000..fbbac154e --- /dev/null +++ b/openaerostruct/geometry/tests/multisection/test_unification_comp.py @@ -0,0 +1,19 @@ +import unittest + +from openaerostruct.geometry.geometry_unification import GeomMultiUnification +from openaerostruct.geometry.geometry_group import build_sections +from openaerostruct.utils.testing import run_test, get_three_section_surface + + +class Test(unittest.TestCase): + def test(self): + (surface, chord_bspline) = get_three_section_surface() + sec_dicts = build_sections(surface) + + comp = GeomMultiUnification(sections=sec_dicts, surface_name=surface["name"]) + + run_test(self, comp, complex_flag=True, method="cs") + + +if __name__ == "__main__": + unittest.main() diff --git a/openaerostruct/utils/check_surface_dict.py b/openaerostruct/utils/check_surface_dict.py index 17c3a10e2..c3a331165 100755 --- a/openaerostruct/utils/check_surface_dict.py +++ b/openaerostruct/utils/check_surface_dict.py @@ -67,6 +67,18 @@ def check_surface_dict_keys(surface): # FFD "mx", "my", + # Multisection + "is_multi_section", + "num_sections", + "sec_name", + "meshes", + "root_chord", + "span", + "ny", + "nx", + "bpanels", + "cpanels", + "root_section", ] for key in surface.keys(): diff --git a/openaerostruct/utils/testing.py b/openaerostruct/utils/testing.py index cccd2a9d4..61cf7417d 100644 --- a/openaerostruct/utils/testing.py +++ b/openaerostruct/utils/testing.py @@ -217,6 +217,155 @@ def get_ground_effect_surfaces(): return surfaces +def get_three_section_surface(sym=True, visc=False): + # Outputs a three section wing surface + # Set-up B-splines for each section. Done here since this information will be needed multiple times. + sec_chord_cp = [np.ones(2), np.ones(2), np.ones(2)] + + surface_dict = { + # Wing definition + # Basic surface parameters + "name": "surface", + "is_multi_section": True, + "num_sections": 3, # The number of sections in the multi-section surface + "sec_name": ["sec0", "sec1", "sec2"], # names of the individual sections + "symmetry": sym, # 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' + # Geometry Parameters + "taper": [1.0, 1.0, 1.0], # Wing taper for each section + "span": [1.0, 1.0, 1.0], # Wing span for each section + "sweep": [0.0, 0, 0.0], # Wing sweep for each section + "twist_cp": [np.zeros(2), np.zeros(2), np.zeros(2)], + "chord_cp": sec_chord_cp, + "root_chord": 1.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": [11, 11, 11], # 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 + "t_over_c_cp": [ + np.array([0.15]), + np.array([0.15]), + np.array([0.15]), + ], # thickness over chord ratio (NACA0015) + "c_max_t": 0.303, # chordwise location of maximum (NACA0015) + # thickness + "with_viscous": visc, # if true, compute viscous drag + "with_wave": False, # if true, compute wave drag + "groundplane": False, + } + + if sym is False: + surface_dict["root_section"] = 1 + + if visc is True: + surface_dict["t_over_c_cp"] = [np.array([0.15]), np.array([0.15])] + surface_dict["nx"] = 3 + + return surface_dict, sec_chord_cp + + +def get_two_section_surface(sym=True, visc=False): + # Outputs a symmetric two section wing surface + # Set-up B-splines for each section. Done here since this information will be needed multiple times. + sec_chord_cp = [np.array([1.0, 1.0]), np.array([1.0, 1.0])] + + # Create a dictionary with info and options about the multi-section aerodynamic + # lifting surface + surface_dict = { + # 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": sym, # 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": [1.0, 1.0], # Wing span for each section + "sweep": [0.0, 0.0], # Wing sweep for each section + "chord_cp": sec_chord_cp, + "twist_cp": [np.zeros(2), np.zeros(2)], + # "chord_cp": [np.ones(1),2*np.ones(1),3*np.ones(1)], #Chord B-spline control points for each section + "root_chord": 1.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": [21, 21], # 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 + # "t_over_c_cp": [np.array([0.15]),np.array([0.15])], # thickness over chord ratio (NACA0015) + "c_max_t": 0.303, # chordwise location of maximum (NACA0015) + # thickness + "with_viscous": visc, # if true, compute viscous drag + "with_wave": False, # if true, compute wave drag + "groundplane": False, + } + + if sym is False: + surface_dict["root_section"] = 1 + + if visc is True: + surface_dict["t_over_c_cp"] = [np.array([0.15]), np.array([0.15])] + surface_dict["nx"] = 3 + + return surface_dict, sec_chord_cp + + +def get_single_section_surface(): + """Create a dictionary with info and options about the aerodynamic + single section lifting surface""" + + # Create a dictionary to store options about the mesh + mesh_dict = { + "num_y": 81, + "num_x": 2, + "wing_type": "rect", + "span": 2.0, + "root_chord": 1.0, + "symmetry": True, + "span_cos_spacing": 0, + "chord_cos_spacing": 0, + } + + # Generate the aerodynamic mesh based on the previous dictionary + mesh = generate_mesh(mesh_dict) + surface_dict = { + # Wing definition + "name": "surface", # name of the surface + "symmetry": True, # if true, model one half of wing + # 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(2), + "mesh": mesh, + "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": False, # if true, compute viscous drag + "with_wave": False, # if true, compute wave drag + "groundplane": False, + } + + return surface_dict + + def assert_check_totals(totals, atol=1e-6, rtol=1e-6): for of_wrt_pair in totals.keys(): of, wrt = of_wrt_pair diff --git a/tests/integration_tests/test_multi_2_sec.py b/tests/integration_tests/test_multi_2_sec.py new file mode 100644 index 000000000..6740121c9 --- /dev/null +++ b/tests/integration_tests/test_multi_2_sec.py @@ -0,0 +1,410 @@ +from openmdao.utils.assert_utils import assert_near_equal +import unittest + + +class Test(unittest.TestCase): + def assert_check_results_sym(self, prob): + """We put the assert_near_equal for the test_constraint and test construction cases here since they should both return the same results.""" + assert_near_equal(prob["aero_point_0.surface_perf.CD"][0], 0.02920058, 1e-3) + assert_near_equal(prob["aero_point_0.surface_perf.CL"][0], 0.2999996, 1e-3) + assert_near_equal(prob["aero_point_0.CM"][1], -0.07335945, 1e-3) + + def test_constraint(self): + import numpy as np + + import openmdao.api as om + + from openaerostruct.geometry.geometry_group import MultiSecGeometry + from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openaerostruct.geometry.geometry_group import build_sections + from openaerostruct.geometry.geometry_unification import unify_mesh + from openaerostruct.utils.testing import get_two_section_surface + + surface, sec_chord_cp = get_two_section_surface() + + # Create the OpenMDAO problem for the constrained version + prob = om.Problem() + + # Create an independent variable component that will supply the flow + # conditions to the problem. + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=1.0, units="m/s") + indep_var_comp.add_output("alpha", val=10.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.3) + indep_var_comp.add_output("re", val=1.0e5, units="1/m") + indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") + indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + + # Add this IndepVarComp to the problem model + 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_sections(surface) + uni_mesh = unify_mesh(section_surfaces) + surface["mesh"] = uni_mesh + + # Create and add a group that handles the geometry for the + # aerodynamic lifting surface + multi_geom_group = MultiSecGeometry( + surface=surface, joining_comp=True, dim_constr=[np.array([1, 0, 0]), np.array([1, 0, 0])] + ) + prob.model.add_subsystem(surface["name"], multi_geom_group) + + # Create the aero point group, which contains the actual aerodynamic + # analyses + aero_group = AeroPoint(surfaces=[surface]) + point_name = "aero_point_0" + prob.model.add_subsystem( + point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"] + ) + + # Get name of surface and construct unified mesh name + name = surface["name"] + unification_name = "{}_unification".format(surface["name"]) + + # Connect the mesh from the mesh unification component to the analysis point + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh" + ) + + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", + point_name + ".aero_states." + "surface" + "_def_mesh", + ) + + # Add DVs + prob.model.add_design_var("surface.sec0.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("surface.sec1.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + # Add joined mesh constraint + prob.model.add_constraint("surface.surface_joining.section_separation", upper=0, lower=0) + + # Add CL constraint + prob.model.add_constraint(point_name + ".CL", equals=0.3) + + # Add Area constraint + prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + + # Add objective + prob.model.add_objective(point_name + ".CD", scaler=1e4) + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["tol"] = 1e-7 + prob.driver.options["disp"] = True + prob.driver.options["maxiter"] = 1000 + # prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up and run the optimization problem + prob.setup() + prob.run_driver() + + self.assert_check_results_sym(prob) + + def test_construction(self): + import numpy as np + + import openmdao.api as om + + from openaerostruct.geometry.geometry_group import MultiSecGeometry + from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openaerostruct.geometry.geometry_group import build_sections + from openaerostruct.geometry.geometry_unification import unify_mesh + from openaerostruct.geometry.multi_unified_bspline_utils import build_multi_spline, connect_multi_spline + from openaerostruct.utils.testing import get_two_section_surface + + surface, sec_chord_cp = get_two_section_surface() + + # Create the OpenMDAO problem + prob = om.Problem() + + # Create an independent variable component that will supply the flow + # conditions to the problem. + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=1.0, units="m/s") + indep_var_comp.add_output("alpha", val=10.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.3) + indep_var_comp.add_output("re", val=1.0e5, units="1/m") + indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") + indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + + # Add this IndepVarComp to the problem model + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Generate the sections and unified mesh here. It's needed to join the sections by construction. + section_surfaces = build_sections(surface) + uni_mesh = unify_mesh(section_surfaces) + surface["mesh"] = uni_mesh + + # Build a component with B-spline control points that joins the sections by construction + chord_comp = build_multi_spline("chord_cp", len(section_surfaces), sec_chord_cp) + prob.model.add_subsystem("chord_bspline", chord_comp) + + # Connect the B-spline component to the section B-splines + connect_multi_spline(prob, section_surfaces, sec_chord_cp, "chord_cp", "chord_bspline") + + # Create and add a group that handles the geometry for the + # aerodynamic lifting surface + multi_geom_group = MultiSecGeometry(surface=surface) + prob.model.add_subsystem(surface["name"], multi_geom_group) + + # Create the aero point group, which contains the actual aerodynamic + # analyses + aero_group = AeroPoint(surfaces=[surface]) + point_name = "aero_point_0" + prob.model.add_subsystem( + point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"] + ) + + # Get name of surface and construct unified mesh name + name = surface["name"] + unification_name = "{}_unification".format(surface["name"]) + + # Connect the mesh from the mesh unification component to the analysis point + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh" + ) + + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", + point_name + ".aero_states." + "surface" + "_def_mesh", + ) + + # Add DVs + prob.model.add_design_var("chord_bspline.chord_cp_spline", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + # Add CL constraint + prob.model.add_constraint(point_name + ".CL", equals=0.3) + + # Add Area constraint + prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + + # Add objective + prob.model.add_objective(point_name + ".CD", scaler=1e4) + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["tol"] = 1e-7 + prob.driver.options["disp"] = True + prob.driver.options["maxiter"] = 1000 + # prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up and run the optimization problem + prob.setup() + prob.run_driver() + + self.assert_check_results_sym(prob) + + def test_asymmetrical(self): + import numpy as np + + import openmdao.api as om + + from openaerostruct.geometry.geometry_group import MultiSecGeometry + from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openaerostruct.geometry.geometry_group import build_sections + from openaerostruct.geometry.geometry_unification import unify_mesh + from openaerostruct.utils.testing import get_two_section_surface + + surface, sec_chord_cp = get_two_section_surface(sym=False, visc=False) + + # Create the OpenMDAO problem + prob = om.Problem() + + # Create an independent variable component that will supply the flow + # conditions to the problem. + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=1.0, units="m/s") + indep_var_comp.add_output("alpha", val=10.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.3) + indep_var_comp.add_output("re", val=1.0e5, units="1/m") + indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") + indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + + # Add this IndepVarComp to the problem model + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Create and add a group that handles the geometry for the + # aerodynamic lifting surface + multi_geom_group = MultiSecGeometry( + surface=surface, joining_comp=True, dim_constr=[np.array([1, 0, 0]), np.array([1, 0, 0])] + ) + prob.model.add_subsystem(surface["name"], multi_geom_group) + + # 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_sections(surface) + uni_mesh = unify_mesh(section_surfaces) + surface["mesh"] = uni_mesh + + # Create the aero point group, which contains the actual aerodynamic + # analyses + aero_group = AeroPoint(surfaces=[surface]) + point_name = "aero_point_0" + prob.model.add_subsystem( + point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"] + ) + + # Get name of surface and construct unified mesh name + name = surface["name"] + unification_name = "{}_unification".format(surface["name"]) + + # Connect the mesh from the mesh unification component to the analysis point + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh" + ) + + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", + point_name + ".aero_states." + "surface" + "_def_mesh", + ) + + prob.model.add_design_var("surface.sec0.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("surface.sec1.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + # Add joined mesh constraint + prob.model.add_constraint("surface.surface_joining.section_separation", upper=0, lower=0) + # prob.model.add_constraint('surface.surface_joining.section_separation',equals=0.0,scaler=1e-4) + + # Add CL constraint + prob.model.add_constraint(point_name + ".CL", equals=0.3) + + # Add Area constraint + prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + + # Add objective + prob.model.add_objective(point_name + ".CD", scaler=1e4) + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["tol"] = 1e-3 + prob.driver.options["disp"] = True + prob.driver.options["maxiter"] = 1000 + # prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up and run the optimization problem + prob.setup() + prob.run_driver() + + assert_near_equal(prob["aero_point_0.surface_perf.CD"][0], 0.02902683, 1e-6) + assert_near_equal(prob["aero_point_0.surface_perf.CL"][0], 0.29999999, 1e-6) + assert_near_equal(prob["aero_point_0.CM"][1], -0.07375071, 1e-6) + + def test_visc(self): + import numpy as np + + import openmdao.api as om + + from openaerostruct.geometry.geometry_group import MultiSecGeometry + from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openaerostruct.geometry.geometry_group import build_sections + from openaerostruct.geometry.geometry_unification import unify_mesh + from openaerostruct.geometry.multi_unified_bspline_utils import build_multi_spline, connect_multi_spline + from openaerostruct.utils.testing import get_two_section_surface + + surface, sec_chord_cp = get_two_section_surface(sym=True, visc=True) + + # Create the OpenMDAO problem + prob = om.Problem() + + # Create an independent variable component that will supply the flow + # conditions to the problem. + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=1.0, units="m/s") + indep_var_comp.add_output("alpha", val=10.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.3) + indep_var_comp.add_output("re", val=1.0e5, units="1/m") + indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") + indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + + # Add this IndepVarComp to the problem model + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Generate the sections and unified mesh here. It's needed to join the sections by construction. + section_surfaces = build_sections(surface) + uni_mesh = unify_mesh(section_surfaces) + surface["mesh"] = uni_mesh + + # Build a component with B-spline control points that joins the sections by construction + chord_comp = build_multi_spline("chord_cp", len(section_surfaces), sec_chord_cp) + prob.model.add_subsystem("chord_bspline", chord_comp) + + # Connect the B-spline component to the section B-splines + connect_multi_spline(prob, section_surfaces, sec_chord_cp, "chord_cp", "chord_bspline") + + # Create and add a group that handles the geometry for the + # aerodynamic lifting surface + multi_geom_group = MultiSecGeometry(surface=surface) + prob.model.add_subsystem(surface["name"], multi_geom_group) + + # Create the aero point group, which contains the actual aerodynamic + # analyses + aero_group = AeroPoint(surfaces=[surface]) + point_name = "aero_point_0" + prob.model.add_subsystem( + point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"] + ) + + # Get name of surface and construct unified mesh name + name = surface["name"] + unification_name = "{}_unification".format(surface["name"]) + + # Connect the mesh from the mesh unification component to the analysis point + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh" + ) + + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", + point_name + ".aero_states." + "surface" + "_def_mesh", + ) + + # Connect t over c B-spline for viscous drag + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_t_over_c", + point_name + "." + name + "_perf." + "t_over_c", + ) + + # Add DVs + prob.model.add_design_var("chord_bspline.chord_cp_spline", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + # Add CL constraint + prob.model.add_constraint(point_name + ".CL", equals=0.3) + + # Add Area constraint + prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + + # Add objective + prob.model.add_objective(point_name + ".CD", scaler=1e4) + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["tol"] = 1e-7 + prob.driver.options["disp"] = True + prob.driver.options["maxiter"] = 1000 + # prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up and run the optimization problem + prob.setup() + prob.run_driver() + + assert_near_equal(prob["aero_point_0.surface_perf.CD"][0], 0.04946752, 1e-6) + assert_near_equal(prob["aero_point_0.surface_perf.CL"][0], 0.3, 1e-6) + assert_near_equal(prob["aero_point_0.CM"][1], -0.05045758, 1e-6) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/integration_tests/test_multi_3_sec.py b/tests/integration_tests/test_multi_3_sec.py new file mode 100644 index 000000000..fd1660c3d --- /dev/null +++ b/tests/integration_tests/test_multi_3_sec.py @@ -0,0 +1,213 @@ +from openmdao.utils.assert_utils import assert_near_equal +import unittest + + +class Test(unittest.TestCase): + def test_symmetric(self): + import numpy as np + + import openmdao.api as om + + from openaerostruct.geometry.geometry_group import MultiSecGeometry + from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openaerostruct.geometry.geometry_group import build_sections + from openaerostruct.geometry.geometry_unification import unify_mesh + from openaerostruct.utils.testing import get_three_section_surface + + # Create a dictionary with info and options about the multi-section aerodynamic + # lifting surface + surface, sec_chord_cp = get_three_section_surface(sym=True, visc=False) + + # Create the OpenMDAO problem + prob = om.Problem() + + # Create an independent variable component that will supply the flow + # conditions to the problem. + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=1.0, units="m/s") + indep_var_comp.add_output("alpha", val=10.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.3) + indep_var_comp.add_output("re", val=1.0e5, units="1/m") + indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") + indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + + # Add this IndepVarComp to the problem model + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Create and add a group that handles the geometry for the + # aerodynamic lifting surface + multi_geom_group = MultiSecGeometry( + surface=surface, + joining_comp=True, + dim_constr=[np.array([1, 0, 0]), np.array([1, 0, 0]), np.array([1, 0, 0])], + ) + prob.model.add_subsystem(surface["name"], multi_geom_group) + + # 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_sections(surface) + uni_mesh = unify_mesh(section_surfaces) + surface["mesh"] = uni_mesh + + # Create the aero point group, which contains the actual aerodynamic + # analyses + aero_group = AeroPoint(surfaces=[surface]) + point_name = "aero_point_0" + prob.model.add_subsystem( + point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"] + ) + + # Get name of surface and construct unified mesh name + name = surface["name"] + unification_name = "{}_unification".format(surface["name"]) + + # Connect the mesh from the mesh unification component to the analysis point + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh" + ) + + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", + point_name + ".aero_states." + "surface" + "_def_mesh", + ) + + # Add DVs + prob.model.add_design_var("surface.sec0.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("surface.sec1.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("surface.sec2.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + # Add joined mesh constraint + prob.model.add_constraint("surface.surface_joining.section_separation", upper=0, lower=0) + # prob.model.add_constraint('surface.surface_joining.section_separation',equals=0.0) + + # Add CL constraint + prob.model.add_constraint(point_name + ".CL", equals=0.3) + + # Add Area constraint + prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0) + + # Add objective + prob.model.add_objective(point_name + ".CD", scaler=1e4) + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["tol"] = 1e-7 + prob.driver.options["disp"] = True + prob.driver.options["maxiter"] = 1000 + # prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up and run the optimization problem + prob.setup() + prob.run_driver() + + assert_near_equal(prob["aero_point_0.surface_perf.CD"][0], 0.02126492, 1e-6) + assert_near_equal(prob["aero_point_0.surface_perf.CL"][0], 0.3, 1e-6) + assert_near_equal(prob["aero_point_0.CM"][1], -0.10778177, 1e-6) + + def test_asymmetric(self): + import numpy as np + + import openmdao.api as om + + from openaerostruct.geometry.geometry_group import MultiSecGeometry + from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openaerostruct.geometry.geometry_group import build_sections + from openaerostruct.geometry.geometry_unification import unify_mesh + from openaerostruct.utils.testing import get_three_section_surface + + surface, sec_chord_cp = get_three_section_surface(sym=False, visc=False) + + # Create the OpenMDAO problem + prob = om.Problem() + + # Create an independent variable component that will supply the flow + # conditions to the problem. + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=1.0, units="m/s") + indep_var_comp.add_output("alpha", val=10.0, units="deg") + indep_var_comp.add_output("Mach_number", val=0.3) + indep_var_comp.add_output("re", val=1.0e5, units="1/m") + indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") + indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + + # Add this IndepVarComp to the problem model + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Create and add a group that handles the geometry for the + # aerodynamic lifting surface + multi_geom_group = MultiSecGeometry( + surface=surface, + joining_comp=True, + dim_constr=[np.array([1, 0, 0]), np.array([1, 0, 0]), np.array([1, 0, 0])], + ) + prob.model.add_subsystem(surface["name"], multi_geom_group) + + # 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_sections(surface) + uni_mesh = unify_mesh(section_surfaces) + surface["mesh"] = uni_mesh + + # Create the aero point group, which contains the actual aerodynamic + # analyses + aero_group = AeroPoint(surfaces=[surface]) + point_name = "aero_point_0" + prob.model.add_subsystem( + point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"] + ) + + # Get name of surface and construct unified mesh name + name = surface["name"] + unification_name = "{}_unification".format(surface["name"]) + + # Connect the mesh from the mesh unification component to the analysis point + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh" + ) + + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", + point_name + ".aero_states." + "surface" + "_def_mesh", + ) + + prob.model.add_design_var("surface.sec0.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("surface.sec1.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("surface.sec2.chord_cp", lower=0.1, upper=10.0, units=None) + prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + # Add joined mesh constraint + prob.model.add_constraint("surface.surface_joining.section_separation", upper=0, lower=0) + # prob.model.add_constraint('surface.surface_joining.section_separation',equals=0.0,scaler=1e-4) + + # Add CL constraint + prob.model.add_constraint(point_name + ".CL", equals=0.3) + + # Add Area constraint + prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.5) + + # Add objective + prob.model.add_objective(point_name + ".CD", scaler=1e4) + + prob.driver = om.ScipyOptimizeDriver() + prob.driver.options["optimizer"] = "SLSQP" + prob.driver.options["tol"] = 1e-7 + prob.driver.options["disp"] = True + prob.driver.options["maxiter"] = 1000 + # prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up and run the optimization problem + prob.setup() + prob.run_driver() + + assert_near_equal(prob["aero_point_0.surface_perf.CD"][0], 0.02270839, 1e-6) + assert_near_equal(prob["aero_point_0.surface_perf.CL"][0], 0.3, 1e-6) + assert_near_equal(prob["aero_point_0.CM"][1], -0.08717908, 1e-6) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/integration_tests/test_multi_single.py b/tests/integration_tests/test_multi_single.py new file mode 100644 index 000000000..85b337c43 --- /dev/null +++ b/tests/integration_tests/test_multi_single.py @@ -0,0 +1,141 @@ +from openmdao.utils.assert_utils import assert_near_equal +import unittest + + +class Test(unittest.TestCase): + def test(self): + import numpy as np + + import openmdao.api as om + + from openaerostruct.geometry.geometry_group import Geometry, MultiSecGeometry + from openaerostruct.aerodynamics.aero_groups import AeroPoint + from openaerostruct.geometry.geometry_group import build_sections + from openaerostruct.geometry.geometry_unification import unify_mesh + from openaerostruct.utils.testing import get_two_section_surface, get_single_section_surface + + """Create a dictionary with info and options about the aerodynamic + single section lifting surface""" + surfaceSingle = get_single_section_surface() + + """Create a dictionary with info and options about the multi-section aerodynamic + lifting surface""" + surfaceMulti, sec_chord_cp = get_two_section_surface() + surfaceMulti["name"] = "surfaceMulti" + + # Create the OpenMDAO problem + prob = om.Problem() + + # Create an independent variable component that will supply the flow + # conditions to the problem. + indep_var_comp = om.IndepVarComp() + indep_var_comp.add_output("v", val=1.0, units="m/s") + indep_var_comp.add_output("alpha", val=10.0 * np.ones(2), units="deg") + indep_var_comp.add_output("Mach_number", val=0.3) + indep_var_comp.add_output("re", val=1.0e5, units="1/m") + indep_var_comp.add_output("rho", val=0.38, units="kg/m**3") + indep_var_comp.add_output("cg", val=np.zeros((3)), units="m") + + # Add this IndepVarComp to the problem model + prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"]) + + # Create and add a group that handles the geometry for the + # aerodynamic single lifting surface + geom_group = Geometry(surface=surfaceSingle) + prob.model.add_subsystem(surfaceSingle["name"], geom_group) + + # Create the aero point group, which contains the actual aerodynamic + # analyses for the single section surface + aero_group = AeroPoint(surfaces=[surfaceSingle]) + point_name0 = "aero_point_0" + prob.model.add_subsystem(point_name0, aero_group) + name = surfaceSingle["name"] + + # Connect Flight conditions + prob.model.connect("v", point_name0 + ".v") + prob.model.connect("alpha", point_name0 + ".alpha", src_indices=[0]) + prob.model.connect("Mach_number", point_name0 + ".Mach_number") + prob.model.connect("re", point_name0 + ".re") + prob.model.connect("rho", point_name0 + ".rho") + prob.model.connect("cg", point_name0 + ".cg") + + # Connect the mesh from the geometry component to the analysis point + prob.model.connect(name + ".mesh", point_name0 + "." + name + ".def_mesh") + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect(name + ".mesh", point_name0 + ".aero_states." + name + "_def_mesh") + + # Create and add a group that handles the geometry for the + # multi-section aerodynamic lifting surface + multi_geom_group = MultiSecGeometry( + surface=surfaceMulti, joining_comp=True, dim_constr=[np.array([1, 0, 0]), np.array([1, 0, 0])] + ) + prob.model.add_subsystem(surfaceMulti["name"], multi_geom_group) + + # 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_sections(surfaceMulti) + uni_mesh = unify_mesh(section_surfaces) + surfaceMulti["mesh"] = uni_mesh + + # Create the aero point group, which contains the actual aerodynamic + # analyses + aero_group = AeroPoint(surfaces=[surfaceMulti]) + point_name1 = "aero_point_1" + prob.model.add_subsystem(point_name1, aero_group) + + # Connect Flight conditions + prob.model.connect("v", point_name1 + ".v") + prob.model.connect("alpha", point_name1 + ".alpha", src_indices=[1]) + prob.model.connect("Mach_number", point_name1 + ".Mach_number") + prob.model.connect("re", point_name1 + ".re") + prob.model.connect("rho", point_name1 + ".rho") + prob.model.connect("cg", point_name1 + ".cg") + + # Get name of surface and construct unified mesh name + name = surfaceMulti["name"] + unification_name = "{}_unification".format(surfaceMulti["name"]) + + # Connect the mesh from the mesh unification component to the analysis point + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", point_name1 + "." + name + ".def_mesh" + ) + + # Perform the connections with the modified names within the + # 'aero_states' group. + prob.model.connect( + name + "." + unification_name + "." + name + "_uni_mesh", + point_name1 + ".aero_states." + name + "_def_mesh", + ) + + # Add DVs + prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg") + + # Add CL constraint + prob.model.add_constraint(point_name0 + ".CL", equals=0.3) + prob.model.add_constraint(point_name1 + ".CL", equals=0.3) + + # Add objective + prob.model.add_objective(point_name0 + ".CD", scaler=1e4) + + 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 + # prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"] + + # Set up and run the optimization problem + prob.setup() + + prob.run_driver() + + CDDiff = np.abs(prob["aero_point_1.surfaceMulti_perf.CD"][0] - prob["aero_point_0.surface_perf.CD"][0]) + CLDiff = np.abs(prob["aero_point_1.surfaceMulti_perf.CL"][0] - prob["aero_point_0.surface_perf.CL"][0]) + + assert_near_equal(CDDiff, 0.0, 1e-6) + assert_near_equal(CLDiff, 0.0, 1e-6) + + +if __name__ == "__main__": + unittest.main()